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

   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 /*----- Header files ------------------------------------------------------*/
  29
  30 #include "mp.h"
  31 #include "mpint.h"
  32 #include "mpmul.h"
  33
  34 /*----- Main code ---------------------------------------------------------*/
  35
  36 /* --- @mpmul_init@ --- *
  37  *
  38  * Arguments:   @mpmul *b@ = pointer to multiplier context to initialize
  39  *
  40  * Returns:     ---
  41  *
  42  * Use:         Initializes a big multiplier context for use.
  43  */
  44
  45 void mpmul_init(mpmul *b)
  46 {
  47   b->i = 0;
  48 }
  49
  50 /* --- @mpmul_add@ --- *
  51  *
  52  * Arguments:   @mpmul *b@ = pointer to multiplier context
  53  *              @mp *x@ = the next factor to multiply in
  54  *
  55  * Returns:     ---
  56  *
  57  * Use:         Contributes another factor to the mix.  It's important that
  58  *              the integer lasts at least as long as the multiplication
  59  *              context; this sort of rules out @mp_build@ integers.
  60  */
  61
  62 #define HWM (MPMUL_DEPTH - 20)
  63 #define LWM (MPMUL_DEPTH / 2)
  64
  65 void mpmul_add(mpmul *b, mp *x)
  66 {
  67   size_t i = b->i;
  68
  69   /* --- Now do the reduction step --- */
  70
  71   x = MP_COPY(x);
  72
  73   while (i > 0) {
  74     if (MP_LEN(b->v[i - 1]) > MP_LEN(x))
  75       break;
  76     i--;
  77     x = mp_mul(x, x, b->v[i]);
  78     MP_DROP(b->v[i]);
  79   }
  80
  81   if (i > HWM) {
  82     while (i > LWM || (i > 0 && MP_LEN(b->v[i - 1]) <= MP_LEN(x))) {
  83       i--;
  84       x = mp_mul(x, x, b->v[i]);
  85       MP_DROP(b->v[i]);
  86     }
  87   }
  88
  89   b->v[i++] = x;
  90   b->i = i;
  91 }
  92
  93 /* --- @mpmul_done@ --- *
  94  *
  95  * Arguments:   @mpmul *b@ = pointer to big multiplication context
  96  *
  97  * Returns:     The product of all the numbers contributed.
  98  *
  99  * Use:         Returns a (large) product of numbers.  The context is
 100  *              deallocated.
 101  */
 102
 103 mp *mpmul_done(mpmul *b)
 104 {
 105   size_t i = b->i;
 106   mp *x;
 107
 108   if (!i)
 109     return (MP_ONE);
 110   i--;
 111   x = b->v[i];
 112   while (i > 0) {
 113     i--;
 114     x = mp_mul(x, x, b->v[i]);
 115     MP_DROP(b->v[i]);
 116   }
 117   return (x);
 118 }
 119
 120 /* --- @mp_factorial@ --- *
 121  *
 122  * Arguments:   @unsigned long i@ = number whose factorial should be
 123  *                      computed.
 124  *
 125  * Returns:     The requested factorial.
 126  */
 127
 128 mp *mp_factorial(unsigned long i)
 129 {
 130   unsigned long j;
 131   mp *x = MP_NEW;
 132   mpmul b = MPMUL_INIT;
 133
 134   for (j = 1; j <= i; j++) {
 135     x = mp_fromulong(x, j);
 136     mpmul_add(&b, x);
 137   }
 138   mp_drop(x);
 139   return (mpmul_done(&b));
 140 }
 141
 142 /*----- Test rig ----------------------------------------------------------*/
 143
 144 #ifdef TEST_RIG
 145
 146 #include <mLib/testrig.h>
 147
 148 static int vfact(dstr *v)
 149 {
 150   unsigned long x = *(unsigned long *)v[0].buf;
 151   mp *fx = *(mp **)v[1].buf;
 152   mp *y = mp_factorial(x);
 153   int ok = 1;
 154   if (!MP_EQ(fx, y)) {
 155     fprintf(stderr, "factorial failed\n");
 156     MP_FPRINTF(stderr, (stderr, "%lu! = ", x), fx);
 157     MP_EPRINT("result", y);
 158     ok = 0;
 159   }
 160   mp_drop(fx);
 161   mp_drop(y);
 162   assert(mparena_count(MPARENA_GLOBAL) == 0);
 163   return (ok);
 164 }
 165
 166 static test_chunk tests[] = {
 167   { "factorial", vfact, { &type_ulong, &type_mp, 0 } },
 168   { 0, 0, { 0 } }
 169 };
 170
 171 int main(int argc, char *argv[])
 172 {
 173   test_run(argc, argv, tests, SRCDIR "/t/mp");
 174   return (0);
 175 }
 176
 177 #endif
 178
 179 /*----- That's all, folks -------------------------------------------------*/