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 2.1 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; see the file COPYING.LIB.  If not, write to
  19 the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
  20 MA 02110-1301, USA.  */
  21
  22 #include "gmp.h"
  23 #include "gmp-impl.h"
  24 #include "randmt.h"
  25
  26
  27 /* Calculate (b^e) mod (2^n-k) for e=1074888996, n=19937 and k=20023,
  28    needed by the seeding function below.  */
  29 static void
  30 mangle_seed (mpz_ptr r, mpz_srcptr b_orig)
  31 {
  32   mpz_t          t, b;
  33   unsigned long  e = 0x40118124;
  34   unsigned long  bit = 0x20000000;
  35
  36   mpz_init (t);
  37   mpz_init_set (b, b_orig);  /* in case r==b_orig */
  38
  39   mpz_set (r, b);
  40   do
  41     {
  42       mpz_mul (r, r, r);
  43
  44     reduce:
  45       for (;;)
  46         {
  47           mpz_tdiv_q_2exp (t, r, 19937L);
  48           if (mpz_sgn (t) == 0)
  49             break;
  50           mpz_tdiv_r_2exp (r, r, 19937L);
  51           mpz_addmul_ui (r, t, 20023L);
  52         }
  53
  54       if ((e & bit) != 0)
  55         {
  56           e &= ~bit;
  57           mpz_mul (r, r, b);
  58           goto reduce;
  59         }
  60
  61       bit >>= 1;
  62     }
  63   while (bit != 0);
  64
  65   mpz_clear (t);
  66   mpz_clear (b);
  67 }
  68
  69
  70 /* Seeding function.  Uses powering modulo a non-Mersenne prime to obtain
  71    a permutation of the input seed space.  The modulus is 2^19937-20023,
  72    which is probably prime.  The power is 1074888996.  In order to avoid
  73    seeds 0 and 1 generating invalid or strange output, the input seed is
  74    first manipulated as follows:
  75
  76      seed1 = seed mod (2^19937-20027) + 2
  77
  78    so that seed1 lies between 2 and 2^19937-20026 inclusive. Then the
  79    powering is performed as follows:
  80
  81      seed2 = (seed1^1074888996) mod (2^19937-20023)
  82
  83    and then seed2 is used to bootstrap the buffer.
  84
  85    This method aims to give guarantees that:
  86      a) seed2 will never be zero,
  87      b) seed2 will very seldom have a very low population of ones in its
  88         binary representation, and
  89      c) every seed between 0 and 2^19937-20028 (inclusive) will yield a
  90         different sequence.
  91
  92    CAVEATS:
  93
  94    The period of the seeding function is 2^19937-20027.  This means that
  95    with seeds 2^19937-20027, 2^19937-20026, ... the exact same sequences
  96    are obtained as with seeds 0, 1, etc.; it also means that seed -1
  97    produces the same sequence as seed 2^19937-20028, etc.
  98  */
  99
 100 static void
 101 randseed_mt (gmp_randstate_t rstate, mpz_srcptr seed)
 102 {
 103   int i;
 104   size_t cnt;
 105
 106   gmp_rand_mt_struct *p;
 107   mpz_t mod;    /* Modulus.  */
 108   mpz_t seed1;  /* Intermediate result.  */
 109
 110   p = (gmp_rand_mt_struct *) RNG_STATE (rstate);
 111
 112   mpz_init (mod);
 113   mpz_init (seed1);
 114
 115   mpz_set_ui (mod, 0L);
 116   mpz_setbit (mod, 19937L);
 117   mpz_sub_ui (mod, mod, 20027L);
 118   mpz_mod (seed1, seed, mod);   /* Reduce `seed' modulo `mod'.  */
 119   mpz_add_ui (seed1, seed1, 2L);        /* seed1 is now ready.  */
 120   mangle_seed (seed1, seed1);   /* Perform the mangling by powering.  */
 121
 122   /* Copy the last bit into bit 31 of mt[0] and clear it.  */
 123   p->mt[0] = (mpz_tstbit (seed1, 19936L) != 0) ? 0x80000000 : 0;
 124   mpz_clrbit (seed1, 19936L);
 125
 126   /* Split seed1 into N-1 32-bit chunks.  */
 127   mpz_export (&p->mt[1], &cnt, -1, sizeof (p->mt[1]), 0,
 128               8 * sizeof (p->mt[1]) - 32, seed1);
 129   cnt++;
 130   ASSERT (cnt <= N);
 131   while (cnt < N)
 132     p->mt[cnt++] = 0;
 133
 134   mpz_clear (mod);
 135   mpz_clear (seed1);
 136
 137   /* Warm the generator up if necessary.  */
 138   if (WARM_UP != 0)
 139     for (i = 0; i < WARM_UP / N; i++)
 140       __gmp_mt_recalc_buffer (p->mt);
 141
 142   p->mti = WARM_UP % N;
 143 }
 144
 145
 146 static const gmp_randfnptr_t Mersenne_Twister_Generator = {
 147   randseed_mt,
 148   __gmp_randget_mt,
 149   __gmp_randclear_mt,
 150   __gmp_randiset_mt
 151 };
 152
 153 /* Initialize MT-specific data.  */
 154 void
 155 gmp_randinit_mt (gmp_randstate_t rstate)
 156 {
 157   __gmp_randinit_mt_noseed (rstate);
 158   RNG_FNPTR (rstate) = (void *) &Mersenne_Twister_Generator;
 159 }