math/gfreduce.[ch]: Fix out-of-bounds memory access.

[u/mdw/catacomb] / math / mpmul.c

/* -*-c-*-
 *
 * 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.
 */

/*----- 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.
 */

#define HWM (MPMUL_DEPTH - 20)
#define LWM (MPMUL_DEPTH / 2)

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 > 0 && MP_LEN(b->v[i - 1]) <= MP_LEN(x))) {
      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));
}

/*----- Test rig ----------------------------------------------------------*/

#ifdef TEST_RIG

#include <mLib/testrig.h>

static int vfact(dstr *v)
{
  unsigned long x = *(unsigned long *)v[0].buf;
  mp *fx = *(mp **)v[1].buf;
  mp *y = mp_factorial(x);
  int ok = 1;
  if (!MP_EQ(fx, y)) {
    fprintf(stderr, "factorial failed\n");
    MP_FPRINTF(stderr, (stderr, "%lu! = ", x), fx);
    MP_EPRINT("result", y);
    ok = 0;
  }
  mp_drop(fx);
  mp_drop(y);
  assert(mparena_count(MPARENA_GLOBAL) == 0);
  return (ok);
}

static test_chunk tests[] = {
  { "factorial", vfact, { &type_ulong, &type_mp, 0 } },
  { 0, 0, { 0 } }
};

int main(int argc, char *argv[])
{
  test_run(argc, argv, tests, SRCDIR "/t/mp");
  return (0);
}

#endif

/*----- That's all, folks -------------------------------------------------*/