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 void mpmul_add(mpmul *b, mp *x)
  65 {
  66   size_t i = b->i;
  67
  68   /* --- Now do the reduction step --- */
  69
  70   x = MP_COPY(x);
  71
  72   while (i > 0) {
  73     if (MP_LEN(b->v[i - 1]) > MP_LEN(x))
  74       break;
  75     i--;
  76     x = mp_mul(x, x, b->v[i]);
  77     MP_DROP(b->v[i]);
  78   }
  79
  80   if (i > HWM) {
  81     while (i > LWM || (i > 0 && MP_LEN(b->v[i - 1]) <= MP_LEN(x))) {
  82       i--;
  83       x = mp_mul(x, x, b->v[i]);
  84       MP_DROP(b->v[i]);
  85     }
  86   }
  87
  88   b->v[i++] = x;
  89   b->i = i;
  90 }
  91
  92 /* --- @mpmul_done@ --- *
  93  *
  94  * Arguments:   @mpmul *b@ = pointer to big multiplication context
  95  *
  96  * Returns:     The product of all the numbers contributed.
  97  *
  98  * Use:         Returns a (large) product of numbers.  The context is
  99  *              deallocated.
 100  */
 101
 102 mp *mpmul_done(mpmul *b)
 103 {
 104   size_t i = b->i;
 105   mp *x;
 106
 107   if (!i)
 108     return (MP_ONE);
 109   i--;
 110   x = b->v[i];
 111   while (i > 0) {
 112     i--;
 113     x = mp_mul(x, x, b->v[i]);
 114     MP_DROP(b->v[i]);
 115   }
 116   return (x);
 117 }
 118
 119 /* --- @mp_factorial@ --- *
 120  *
 121  * Arguments:   @unsigned long i@ = number whose factorial should be
 122  *                      computed.
 123  *
 124  * Returns:     The requested factorial.
 125  */
 126
 127 mp *mp_factorial(unsigned long i)
 128 {
 129   unsigned long j;
 130   mp *x = MP_NEW;
 131   mpmul b = MPMUL_INIT;
 132
 133   for (j = 1; j <= i; j++) {
 134     x = mp_fromulong(x, j);
 135     mpmul_add(&b, x);
 136   }
 137   mp_drop(x);
 138   return (mpmul_done(&b));
 139 }
 140
 141 /*----- Test rig ----------------------------------------------------------*/
 142
 143 #ifdef TEST_RIG
 144
 145 #include <mLib/testrig.h>
 146
 147 static int vfact(dstr *v)
 148 {
 149   unsigned long x = *(unsigned long *)v[0].buf;
 150   mp *fx = *(mp **)v[1].buf;
 151   mp *y = mp_factorial(x);
 152   int ok = 1;
 153   if (!MP_EQ(fx, y)) {
 154     fprintf(stderr, "factorial failed\n");
 155     MP_FPRINTF(stderr, (stderr, "%lu! = ", x), fx);
 156     MP_EPRINT("result", y);
 157     ok = 0;
 158   }
 159   mp_drop(fx);
 160   mp_drop(y);
 161   assert(mparena_count(MPARENA_GLOBAL) == 0);
 162   return (ok);
 163 }
 164
 165 static test_chunk tests[] = {
 166   { "factorial", vfact, { &type_ulong, &type_mp, 0 } },
 167   { 0, 0, { 0 } }
 168 };
 169
 170 int main(int argc, char *argv[])
 171 {
 172   test_run(argc, argv, tests, SRCDIR "/tests/mp");
 173   return (0);
 174 }
 175
 176 #endif
 177
 178 /*----- That's all, folks -------------------------------------------------*/