From: mdw Date: Sat, 1 Jul 2000 11:21:39 +0000 (+0000) Subject: New interface for computing products of many (small) integers. X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/commitdiff_plain/da38281c333dcb1fa319ce300fe459dae193b99f New interface for computing products of many (small) integers. --- diff --git a/mpmul.c b/mpmul.c new file mode 100644 index 0000000..46f7b72 --- /dev/null +++ b/mpmul.c @@ -0,0 +1,149 @@ +/* -*-c-*- + * + * $Id: mpmul.c,v 1.1 2000/07/01 11:21:39 mdw Exp $ + * + * Multiply many small numbers together + * + * (c) 2000 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mpmul.c,v $ + * Revision 1.1 2000/07/01 11:21:39 mdw + * New interface for computing products of many (small) integers. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include "mp.h" +#include "mpint.h" +#include "mpmul.h" + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @mpmul_init@ --- * + * + * Arguments: @mpmul *b@ = pointer to multiplier context to initialize + * + * Returns: --- + * + * Use: Initializes a big multiplier context for use. + */ + +void mpmul_init(mpmul *b) +{ + b->i = 0; +} + +/* --- @mpmul_add@ --- * + * + * Arguments: @mpmul *b@ = pointer to multiplier context + * @mp *x@ = the next factor to multiply in + * + * Returns: --- + * + * Use: Contributes another factor to the mix. It's important that + * the integer lasts at least as long as the multiplication + * context; this sort of rules out @mp_build@ integers. + */ + +void mpmul_add(mpmul *b, mp *x) +{ + size_t i = b->i; + + /* --- Now do the reduction step --- */ + + x = MP_COPY(x); + + while (i > 0) { + if (MP_LEN(b->v[i - 1]) > MP_LEN(x)) + break; + i--; + x = mp_mul(x, x, b->v[i]); + MP_DROP(b->v[i]); + } + + if (i > HWM) { + while (i > LWM) { + i--; + x = mp_mul(x, x, b->v[i]); + MP_DROP(b->v[i]); + } + } + + b->v[i] = x; + b->i = i; +} + +/* --- @mpmul_done@ --- * + * + * Arguments: @mpmul *b@ = pointer to big multiplication context + * + * Returns: The product of all the numbers contributed. + * + * Use: Returns a (large) product of numbers. The context is + * deallocated. + */ + +mp *mpmul_done(mpmul *b) +{ + size_t i = b->i; + mp *x; + + if (!i) + return (MP_ONE); + i--; + x = b->v[i]; + while (i > 0) { + i--; + x = mp_mul(x, x, b->v[i]); + MP_DROP(b->v[i]); + } + return (x); +} + +/* --- @mp_factorial@ --- * + * + * Arguments: @unsigned long i@ = number whose factorial should be + * computed. + * + * Returns: The requested factorial. + */ + +mp *mp_factorial(unsigned long i) +{ + unsigned long j; + mp *x = MP_NEW; + mpmul b = MPMUL_INIT; + + for (j = 1; j <= i; j++) { + x = mp_fromulong(x, j); + mpmul_add(&b, x); + } + mp_drop(x); + return (mpmul_done(&b)); +} + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/mpmul.h b/mpmul.h new file mode 100644 index 0000000..76c9559 --- /dev/null +++ b/mpmul.h @@ -0,0 +1,152 @@ +/* -*-c-*- + * + * $Id: mpmul.h,v 1.1 2000/07/01 11:21:39 mdw Exp $ + * + * Multiply many small numbers together + * + * (c) 2000 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mpmul.h,v $ + * Revision 1.1 2000/07/01 11:21:39 mdw + * New interface for computing products of many (small) integers. + * + */ + +#ifndef CATACOMB_MPMUL_H +#define CATACOMB_MPMUL_H + +#ifdef __cplusplus + extern "C" { +#endif + +/*----- Header files ------------------------------------------------------*/ + +#ifndef CATACOMB_MP_H +# include "mp.h" +#endif + +/*----- Magic numbers -----------------------------------------------------*/ + +/* --- How the algorithm works --- * + * + * Multiplication on large integers is least wasteful when the numbers + * multiplied are approximately the same size. When a new multiplier is + * added to the system, we push it onto a stack. Then we `reduce' the stack: + * while the value on the top of the stack is not shorter than the value + * below it, replace the top two elements by their product. + * + * Let %$b$% be the radix of our multiprecision integers, and let %$Z$% be + * the maximum number of digits. Then the largest integer we can represent + * is %$M - 1 = b^Z - 1$%. We could assume that all of the integers we're + * given are about the same size. This would give us the same upper bound as + * that derived in `mptext.c'. + * + * However, we're in less control over our inputs. In particular, if a + * sequence of integers with strictly decreasing lengths is input then we're + * sunk. Suppose that the stack contains, from top to bottom, %$b^i$%, + * %$b^{i+1}$%, ..., %$b^n$%. The final product will therefore be + * %$p = b^{(n+i)(n-i+1)/2}$%. We must now find the maximum stack depth + * %$d = n - i$% such that %$p > M$%. + * + * Taking logs of both sides gives that %$(d + 2 i)(d + 1) > 2 Z$%. We can + * maximize %$d$% by taking %$i = 0$%, which gives that %$d^2 + d > 2 Z$%, so + * %$d$% must be approximately %$(\sqrt{8 Z + 1} - 1)/2$%, which is + * uncomfortably large. + * + * We compromise by choosing double the `mptext' bound and imposing high- and + * low-water marks for forced reduction. + */ + +#define MPMUL_DEPTH (2 * (CHAR_BIT * sizeof(size_t) + 10)) + +#define HWM (MPMUL_DEPTH - 20) +#define LWM (MPMUL_DEPTH / 2) + +/*----- Data structures ---------------------------------------------------*/ + +typedef struct mpmul { + size_t i; + mp *v[MPMUL_DEPTH]; +} mpmul; + +#define MPMUL_INIT { 0 } + +/*----- Functions provided ------------------------------------------------*/ + +/* --- @mpmul_init@ --- * + * + * Arguments: @mpmul *b@ = pointer to multiplier context to initialize + * + * Returns: --- + * + * Use: Initializes a big multiplier context for use. + */ + +extern void mpmul_init(mpmul */*b*/); + +/* --- @mpmul_add@ --- * + * + * Arguments: @mpmul *b@ = pointer to multiplier context + * @mp *x@ = the next factor to multiply in + * + * Returns: --- + * + * Use: Contributes another factor to the mix. It's important that + * the integer lasts at least as long as the multiplication + * context; this sort of rules out @mp_build@ integers. + */ + +extern void mpmul_add(mpmul */*b*/, mp */*x*/); + +/* --- @mpmul_done@ --- * + * + * Arguments: @mpmul *b@ = pointer to big multiplication context + * + * Returns: The product of all the numbers contributed. + * + * Use: Returns a (large) product of numbers. The context is + * deallocated. + */ + +extern mp *mpmul_done(mpmul */*b*/); + +/* --- @mp_factorial@ --- * + * + * Arguments: @unsigned long i@ = number whose factorial should be + * computed. + * + * Returns: The requested factorial. + */ + +extern mp *mp_factorial(unsigned long /*i*/); + +/*----- That's all, folks -------------------------------------------------*/ + +#ifdef __cplusplus + } +#endif + +#endif