chiark - git - mdw - catacomb/blob - mpmul.h

   1 /* -*-c-*-
   2  *
   3  * $Id: mpmul.h,v 1.1 2000/07/01 11:21:39 mdw Exp $
   4  *
   5  * Multiply many small numbers together
   6  *
   7  * (c) 2000 Straylight/Edgeware
   8  */
   9
  10 /*----- Licensing notice --------------------------------------------------*
  11  *
  12  * This file is part of Catacomb.
  13  *
  14  * Catacomb is free software; you can redistribute it and/or modify
  15  * it under the terms of the GNU Library General Public License as
  16  * published by the Free Software Foundation; either version 2 of the
  17  * License, or (at your option) any later version.
  18  *
  19  * Catacomb is distributed in the hope that it will be useful,
  20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22  * GNU Library General Public License for more details.
  23  *
  24  * You should have received a copy of the GNU Library General Public
  25  * License along with Catacomb; if not, write to the Free
  26  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  27  * MA 02111-1307, USA.
  28  */
  29
  30 /*----- Revision history --------------------------------------------------*
  31  *
  32  * $Log: mpmul.h,v $
  33  * Revision 1.1  2000/07/01 11:21:39  mdw
  34  * New interface for computing products of many (small) integers.
  35  *
  36  */
  37
  38 #ifndef CATACOMB_MPMUL_H
  39 #define CATACOMB_MPMUL_H
  40
  41 #ifdef __cplusplus
  42   extern "C" {
  43 #endif
  44
  45 /*----- Header files ------------------------------------------------------*/
  46
  47 #ifndef CATACOMB_MP_H
  48 #  include "mp.h"
  49 #endif
  50
  51 /*----- Magic numbers -----------------------------------------------------*/
  52
  53 /* --- How the algorithm works --- *
  54  *
  55  * Multiplication on large integers is least wasteful when the numbers
  56  * multiplied are approximately the same size.  When a new multiplier is
  57  * added to the system, we push it onto a stack.  Then we `reduce' the stack:
  58  * while the value on the top of the stack is not shorter than the value
  59  * below it, replace the top two elements by their product.
  60  *
  61  * Let %$b$% be the radix of our multiprecision integers, and let %$Z$% be
  62  * the maximum number of digits.  Then the largest integer we can represent
  63  * is %$M - 1 = b^Z - 1$%.  We could assume that all of the integers we're
  64  * given are about the same size.  This would give us the same upper bound as
  65  * that derived in `mptext.c'.
  66  *
  67  * However, we're in less control over our inputs.  In particular, if a
  68  * sequence of integers with strictly decreasing lengths is input then we're
  69  * sunk.  Suppose that the stack contains, from top to bottom, %$b^i$%,
  70  * %$b^{i+1}$%, ..., %$b^n$%.  The final product will therefore be
  71  * %$p = b^{(n+i)(n-i+1)/2}$%.  We must now find the maximum stack depth
  72  * %$d = n - i$% such that %$p > M$%.
  73  *
  74  * Taking logs of both sides gives that %$(d + 2 i)(d + 1) > 2 Z$%.  We can
  75  * maximize %$d$% by taking %$i = 0$%, which gives that %$d^2 + d > 2 Z$%, so
  76  * %$d$% must be approximately %$(\sqrt{8 Z + 1} - 1)/2$%, which is
  77  * uncomfortably large.
  78  *
  79  * We compromise by choosing double the `mptext' bound and imposing high- and
  80  * low-water marks for forced reduction.
  81  */
  82
  83 #define MPMUL_DEPTH (2 * (CHAR_BIT * sizeof(size_t) + 10))
  84
  85 #define HWM (MPMUL_DEPTH - 20)
  86 #define LWM (MPMUL_DEPTH / 2)
  87
  88 /*----- Data structures ---------------------------------------------------*/
  89
  90 typedef struct mpmul {
  91   size_t i;
  92   mp *v[MPMUL_DEPTH];
  93 } mpmul;
  94
  95 #define MPMUL_INIT { 0 }
  96
  97 /*----- Functions provided ------------------------------------------------*/
  98
  99 /* --- @mpmul_init@ --- *
 100  *
 101  * Arguments:   @mpmul *b@ = pointer to multiplier context to initialize
 102  *
 103  * Returns:     ---
 104  *
 105  * Use:         Initializes a big multiplier context for use.
 106  */
 107
 108 extern void mpmul_init(mpmul */*b*/);
 109
 110 /* --- @mpmul_add@ --- *
 111  *
 112  * Arguments:   @mpmul *b@ = pointer to multiplier context
 113  *              @mp *x@ = the next factor to multiply in
 114  *
 115  * Returns:     ---
 116  *
 117  * Use:         Contributes another factor to the mix.  It's important that
 118  *              the integer lasts at least as long as the multiplication
 119  *              context; this sort of rules out @mp_build@ integers.
 120  */
 121
 122 extern void mpmul_add(mpmul */*b*/, mp */*x*/);
 123
 124 /* --- @mpmul_done@ --- *
 125  *
 126  * Arguments:   @mpmul *b@ = pointer to big multiplication context
 127  *
 128  * Returns:     The product of all the numbers contributed.
 129  *
 130  * Use:         Returns a (large) product of numbers.  The context is
 131  *              deallocated.
 132  */
 133
 134 extern mp *mpmul_done(mpmul */*b*/);
 135
 136 /* --- @mp_factorial@ --- *
 137  *
 138  * Arguments:   @unsigned long i@ = number whose factorial should be
 139  *                      computed.
 140  *
 141  * Returns:     The requested factorial.
 142  */
 143
 144 extern mp *mp_factorial(unsigned long /*i*/);
 145
 146 /*----- That's all, folks -------------------------------------------------*/
 147
 148 #ifdef __cplusplus
 149   }
 150 #endif
 151
 152 #endif