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

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