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

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