contrib/gmp/randmts.c

   1 /* Mersenne Twister pseudo-random number generator functions.
   2
   3 Copyright 2002, 2003 Free Software Foundation, Inc.
   4
   5 This file is part of the GNU MP Library.
   6
   7 The GNU MP Library is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU Lesser General Public License as published by
   9 the Free Software Foundation; either version 3 of the License, or (at your
  10 option) any later version.
  11
  12 The GNU MP Library is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  19
  20 #include "gmp.h"
  21 #include "gmp-impl.h"
  22 #include "randmt.h"
  23
  24
  25 /* Calculate (b^e) mod (2^n-k) for e=1074888996, n=19937 and k=20023,
  26    needed by the seeding function below.  */
  27 static void
  28 mangle_seed (mpz_ptr r, mpz_srcptr b_orig)
  29 {
  30   mpz_t          t, b;
  31   unsigned long  e = 0x40118124;
  32   unsigned long  bit = 0x20000000;
  33
  34   mpz_init (t);
  35   mpz_init_set (b, b_orig);  /* in case r==b_orig */
  36
  37   mpz_set (r, b);
  38   do
  39     {
  40       mpz_mul (r, r, r);
  41
  42     reduce:
  43       for (;;)
  44         {
  45           mpz_tdiv_q_2exp (t, r, 19937L);
  46           if (mpz_sgn (t) == 0)
  47             break;
  48           mpz_tdiv_r_2exp (r, r, 19937L);
  49           mpz_addmul_ui (r, t, 20023L);
  50         }
  51
  52       if ((e & bit) != 0)
  53         {
  54           e &= ~bit;
  55           mpz_mul (r, r, b);
  56           goto reduce;
  57         }
  58
  59       bit >>= 1;
  60     }
  61   while (bit != 0);
  62
  63   mpz_clear (t);
  64   mpz_clear (b);
  65 }
  66
  67
  68 /* Seeding function.  Uses powering modulo a non-Mersenne prime to obtain
  69    a permutation of the input seed space.  The modulus is 2^19937-20023,
  70    which is probably prime.  The power is 1074888996.  In order to avoid
  71    seeds 0 and 1 generating invalid or strange output, the input seed is
  72    first manipulated as follows:
  73
  74      seed1 = seed mod (2^19937-20027) + 2
  75
  76    so that seed1 lies between 2 and 2^19937-20026 inclusive. Then the
  77    powering is performed as follows:
  78
  79      seed2 = (seed1^1074888996) mod (2^19937-20023)
  80
  81    and then seed2 is used to bootstrap the buffer.
  82
  83    This method aims to give guarantees that:
  84      a) seed2 will never be zero,
  85      b) seed2 will very seldom have a very low population of ones in its
  86         binary representation, and
  87      c) every seed between 0 and 2^19937-20028 (inclusive) will yield a
  88         different sequence.
  89
  90    CAVEATS:
  91
  92    The period of the seeding function is 2^19937-20027.  This means that
  93    with seeds 2^19937-20027, 2^19937-20026, ... the exact same sequences
  94    are obtained as with seeds 0, 1, etc.; it also means that seed -1
  95    produces the same sequence as seed 2^19937-20028, etc.
  96  */
  97
  98 static void
  99 randseed_mt (gmp_randstate_t rstate, mpz_srcptr seed)
 100 {
 101   int i;
 102   size_t cnt;
 103
 104   gmp_rand_mt_struct *p;
 105   mpz_t mod;    /* Modulus.  */
 106   mpz_t seed1;  /* Intermediate result.  */
 107
 108   p = (gmp_rand_mt_struct *) RNG_STATE (rstate);
 109
 110   mpz_init (mod);
 111   mpz_init (seed1);
 112
 113   mpz_set_ui (mod, 0L);
 114   mpz_setbit (mod, 19937L);
 115   mpz_sub_ui (mod, mod, 20027L);
 116   mpz_mod (seed1, seed, mod);   /* Reduce `seed' modulo `mod'.  */
 117   mpz_add_ui (seed1, seed1, 2L);        /* seed1 is now ready.  */
 118   mangle_seed (seed1, seed1);   /* Perform the mangling by powering.  */
 119
 120   /* Copy the last bit into bit 31 of mt[0] and clear it.  */
 121   p->mt[0] = (mpz_tstbit (seed1, 19936L) != 0) ? 0x80000000 : 0;
 122   mpz_clrbit (seed1, 19936L);
 123
 124   /* Split seed1 into N-1 32-bit chunks.  */
 125   mpz_export (&p->mt[1], &cnt, -1, sizeof (p->mt[1]), 0,
 126               8 * sizeof (p->mt[1]) - 32, seed1);
 127   cnt++;
 128   ASSERT (cnt <= N);
 129   while (cnt < N)
 130     p->mt[cnt++] = 0;
 131
 132   mpz_clear (mod);
 133   mpz_clear (seed1);
 134
 135   /* Warm the generator up if necessary.  */
 136   if (WARM_UP != 0)
 137     for (i = 0; i < WARM_UP / N; i++)
 138       __gmp_mt_recalc_buffer (p->mt);
 139
 140   p->mti = WARM_UP % N;
 141 }
 142
 143
 144 static const gmp_randfnptr_t Mersenne_Twister_Generator = {
 145   randseed_mt,
 146   __gmp_randget_mt,
 147   __gmp_randclear_mt,
 148   __gmp_randiset_mt
 149 };
 150
 151 /* Initialize MT-specific data.  */
 152 void
 153 gmp_randinit_mt (gmp_randstate_t rstate)
 154 {
 155   __gmp_randinit_mt_noseed (rstate);
 156   RNG_FNPTR (rstate) = (void *) &Mersenne_Twister_Generator;
 157 }