randmts.c

/* Mersenne Twister pseudo-random number generator functions.

Copyright 2002, 2003 Free Software Foundation, Inc.

This file is part of the GNU MP Library.

The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

The GNU MP Library 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 Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */

#include "gmp.h"
#include "gmp-impl.h"
#include "randmt.h"


/* Calculate (b^e) mod (2^n-k) for e=1074888996, n=19937 and k=20023,
   needed by the seeding function below.  */
static void
mangle_seed (mpz_ptr r, mpz_srcptr b_orig)
{
  mpz_t          t, b;
  unsigned long  e = 0x40118124;
  unsigned long  bit = 0x20000000;

  mpz_init (t);
  mpz_init_set (b, b_orig);  /* in case r==b_orig */

  mpz_set (r, b);
  do
    {
      mpz_mul (r, r, r);

    reduce:
      for (;;)
        {
          mpz_tdiv_q_2exp (t, r, 19937L);
          if (mpz_sgn (t) == 0)
            break;
          mpz_tdiv_r_2exp (r, r, 19937L);
          mpz_addmul_ui (r, t, 20023L);
        }

      if ((e & bit) != 0)
        {
          e &= ~bit;
          mpz_mul (r, r, b);
          goto reduce;
        }

      bit >>= 1;
    }
  while (bit != 0);

  mpz_clear (t);
  mpz_clear (b);
}


/* Seeding function.  Uses powering modulo a non-Mersenne prime to obtain
   a permutation of the input seed space.  The modulus is 2^19937-20023,
   which is probably prime.  The power is 1074888996.  In order to avoid
   seeds 0 and 1 generating invalid or strange output, the input seed is
   first manipulated as follows:

     seed1 = seed mod (2^19937-20027) + 2

   so that seed1 lies between 2 and 2^19937-20026 inclusive. Then the
   powering is performed as follows:

     seed2 = (seed1^1074888996) mod (2^19937-20023)

   and then seed2 is used to bootstrap the buffer.

   This method aims to give guarantees that:
     a) seed2 will never be zero,
     b) seed2 will very seldom have a very low population of ones in its
	binary representation, and
     c) every seed between 0 and 2^19937-20028 (inclusive) will yield a
	different sequence.

   CAVEATS:

   The period of the seeding function is 2^19937-20027.  This means that
   with seeds 2^19937-20027, 2^19937-20026, ... the exact same sequences
   are obtained as with seeds 0, 1, etc.; it also means that seed -1
   produces the same sequence as seed 2^19937-20028, etc.
 */

static void
randseed_mt (gmp_randstate_t rstate, mpz_srcptr seed)
{
  int i;
  size_t cnt;

  gmp_rand_mt_struct *p;
  mpz_t mod;    /* Modulus.  */
  mpz_t seed1;  /* Intermediate result.  */

  p = (gmp_rand_mt_struct *) RNG_STATE (rstate);

  mpz_init (mod);
  mpz_init (seed1);

  mpz_set_ui (mod, 0L);
  mpz_setbit (mod, 19937L);
  mpz_sub_ui (mod, mod, 20027L);
  mpz_mod (seed1, seed, mod);	/* Reduce `seed' modulo `mod'.  */
  mpz_add_ui (seed1, seed1, 2L);	/* seed1 is now ready.  */
  mangle_seed (seed1, seed1);	/* Perform the mangling by powering.  */

  /* Copy the last bit into bit 31 of mt[0] and clear it.  */
  p->mt[0] = (mpz_tstbit (seed1, 19936L) != 0) ? 0x80000000 : 0;
  mpz_clrbit (seed1, 19936L);

  /* Split seed1 into N-1 32-bit chunks.  */
  mpz_export (&p->mt[1], &cnt, -1, sizeof (p->mt[1]), 0,
              8 * sizeof (p->mt[1]) - 32, seed1);
  cnt++;
  ASSERT (cnt <= N);
  while (cnt < N)
    p->mt[cnt++] = 0;

  mpz_clear (mod);
  mpz_clear (seed1);

  /* Warm the generator up if necessary.  */
  if (WARM_UP != 0)
    for (i = 0; i < WARM_UP / N; i++)
      __gmp_mt_recalc_buffer (p->mt);

  p->mti = WARM_UP % N;
}


static const gmp_randfnptr_t Mersenne_Twister_Generator = {
  randseed_mt,
  __gmp_randget_mt,
  __gmp_randclear_mt,
  __gmp_randiset_mt
};

/* Initialize MT-specific data.  */
void
gmp_randinit_mt (gmp_randstate_t rstate)
{
  __gmp_randinit_mt_noseed (rstate);
  RNG_FNPTR (rstate) = (void *) &Mersenne_Twister_Generator;
}