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

   1 /* -*-c-*-
   2  *
   3  * $Id: mpmul.h,v 1.2 2004/04/08 01:36:15 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 #ifndef CATACOMB_MPMUL_H
  31 #define CATACOMB_MPMUL_H
  32
  33 #ifdef __cplusplus
  34   extern "C" {
  35 #endif
  36
  37 /*----- Header files ------------------------------------------------------*/
  38
  39 #ifndef CATACOMB_MP_H
  40 #  include "mp.h"
  41 #endif
  42
  43 /*----- Magic numbers -----------------------------------------------------*/
  44
  45 /* --- How the algorithm works --- *
  46  *
  47  * Multiplication on large integers is least wasteful when the numbers
  48  * multiplied are approximately the same size.  When a new multiplier is
  49  * added to the system, we push it onto a stack.  Then we `reduce' the stack:
  50  * while the value on the top of the stack is not shorter than the value
  51  * below it, replace the top two elements by their product.
  52  *
  53  * Let %$b$% be the radix of our multiprecision integers, and let %$Z$% be
  54  * the maximum number of digits.  Then the largest integer we can represent
  55  * is %$M - 1 = b^Z - 1$%.  We could assume that all of the integers we're
  56  * given are about the same size.  This would give us the same upper bound as
  57  * that derived in `mptext.c'.
  58  *
  59  * However, we're in less control over our inputs.  In particular, if a
  60  * sequence of integers with strictly decreasing lengths is input then we're
  61  * sunk.  Suppose that the stack contains, from top to bottom, %$b^i$%,
  62  * %$b^{i+1}$%, ..., %$b^n$%.  The final product will therefore be
  63  * %$p = b^{(n+i)(n-i+1)/2}$%.  We must now find the maximum stack depth
  64  * %$d = n - i$% such that %$p > M$%.
  65  *
  66  * Taking logs of both sides gives that %$(d + 2 i)(d + 1) > 2 Z$%.  We can
  67  * maximize %$d$% by taking %$i = 0$%, which gives that %$d^2 + d > 2 Z$%, so
  68  * %$d$% must be approximately %$(\sqrt{8 Z + 1} - 1)/2$%, which is
  69  * uncomfortably large.
  70  *
  71  * We compromise by choosing double the `mptext' bound and imposing high- and
  72  * low-water marks for forced reduction.
  73  */
  74
  75 #define MPMUL_DEPTH (2 * (CHAR_BIT * sizeof(size_t) + 10))
  76
  77 /*----- Data structures ---------------------------------------------------*/
  78
  79 typedef struct mpmul {
  80   size_t i;
  81   mp *v[MPMUL_DEPTH];
  82 } mpmul;
  83
  84 #define MPMUL_INIT { 0 }
  85
  86 /*----- Functions provided ------------------------------------------------*/
  87
  88 /* --- @mpmul_init@ --- *
  89  *
  90  * Arguments:   @mpmul *b@ = pointer to multiplier context to initialize
  91  *
  92  * Returns:     ---
  93  *
  94  * Use:         Initializes a big multiplier context for use.
  95  */
  96
  97 extern void mpmul_init(mpmul */*b*/);
  98
  99 /* --- @mpmul_add@ --- *
 100  *
 101  * Arguments:   @mpmul *b@ = pointer to multiplier context
 102  *              @mp *x@ = the next factor to multiply in
 103  *
 104  * Returns:     ---
 105  *
 106  * Use:         Contributes another factor to the mix.  It's important that
 107  *              the integer lasts at least as long as the multiplication
 108  *              context; this sort of rules out @mp_build@ integers.
 109  */
 110
 111 extern void mpmul_add(mpmul */*b*/, mp */*x*/);
 112
 113 /* --- @mpmul_done@ --- *
 114  *
 115  * Arguments:   @mpmul *b@ = pointer to big multiplication context
 116  *
 117  * Returns:     The product of all the numbers contributed.
 118  *
 119  * Use:         Returns a (large) product of numbers.  The context is
 120  *              deallocated.
 121  */
 122
 123 extern mp *mpmul_done(mpmul */*b*/);
 124
 125 /* --- @mp_factorial@ --- *
 126  *
 127  * Arguments:   @unsigned long i@ = number whose factorial should be
 128  *                      computed.
 129  *
 130  * Returns:     The requested factorial.
 131  */
 132
 133 extern mp *mp_factorial(unsigned long /*i*/);
 134
 135 /*----- That's all, folks -------------------------------------------------*/
 136
 137 #ifdef __cplusplus
 138   }
 139 #endif
 140
 141 #endif