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

   1 /* -*-c-*-
   2  *
   3  * $Id: mpmul.c,v 1.1 2000/07/01 11:21:39 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 /*----- Revision history --------------------------------------------------*
  31  *
  32  * $Log: mpmul.c,v $
  33  * Revision 1.1  2000/07/01 11:21:39  mdw
  34  * New interface for computing products of many (small) integers.
  35  *
  36  */
  37
  38 /*----- Header files ------------------------------------------------------*/
  39
  40 #include "mp.h"
  41 #include "mpint.h"
  42 #include "mpmul.h"
  43
  44 /*----- Main code ---------------------------------------------------------*/
  45
  46 /* --- @mpmul_init@ --- *
  47  *
  48  * Arguments:   @mpmul *b@ = pointer to multiplier context to initialize
  49  *
  50  * Returns:     ---
  51  *
  52  * Use:         Initializes a big multiplier context for use.
  53  */
  54
  55 void mpmul_init(mpmul *b)
  56 {
  57   b->i = 0;
  58 }
  59
  60 /* --- @mpmul_add@ --- *
  61  *
  62  * Arguments:   @mpmul *b@ = pointer to multiplier context
  63  *              @mp *x@ = the next factor to multiply in
  64  *
  65  * Returns:     ---
  66  *
  67  * Use:         Contributes another factor to the mix.  It's important that
  68  *              the integer lasts at least as long as the multiplication
  69  *              context; this sort of rules out @mp_build@ integers.
  70  */
  71
  72 void mpmul_add(mpmul *b, mp *x)
  73 {
  74   size_t i = b->i;
  75
  76   /* --- Now do the reduction step --- */
  77
  78   x = MP_COPY(x);
  79
  80   while (i > 0) {
  81     if (MP_LEN(b->v[i - 1]) > MP_LEN(x))
  82       break;
  83     i--;
  84     x = mp_mul(x, x, b->v[i]);
  85     MP_DROP(b->v[i]);
  86   }
  87
  88   if (i > HWM) {
  89     while (i > LWM) {
  90       i--;
  91       x = mp_mul(x, x, b->v[i]);
  92       MP_DROP(b->v[i]);
  93     }
  94   }
  95
  96   b->v[i] = x;
  97   b->i = i;
  98 }
  99
 100 /* --- @mpmul_done@ --- *
 101  *
 102  * Arguments:   @mpmul *b@ = pointer to big multiplication context
 103  *
 104  * Returns:     The product of all the numbers contributed.
 105  *
 106  * Use:         Returns a (large) product of numbers.  The context is
 107  *              deallocated.
 108  */
 109
 110 mp *mpmul_done(mpmul *b)
 111 {
 112   size_t i = b->i;
 113   mp *x;
 114
 115   if (!i)
 116     return (MP_ONE);
 117   i--;
 118   x = b->v[i];
 119   while (i > 0) {
 120     i--;
 121     x = mp_mul(x, x, b->v[i]);
 122     MP_DROP(b->v[i]);
 123   }
 124   return (x);
 125 }
 126
 127 /* --- @mp_factorial@ --- *
 128  *
 129  * Arguments:   @unsigned long i@ = number whose factorial should be
 130  *                      computed.
 131  *
 132  * Returns:     The requested factorial.
 133  */
 134
 135 mp *mp_factorial(unsigned long i)
 136 {
 137   unsigned long j;
 138   mp *x = MP_NEW;
 139   mpmul b = MPMUL_INIT;
 140
 141   for (j = 1; j <= i; j++) {
 142     x = mp_fromulong(x, j);
 143     mpmul_add(&b, x);
 144   }
 145   mp_drop(x);
 146   return (mpmul_done(&b));
 147 }
 148
 149 /*----- That's all, folks -------------------------------------------------*/