1/* Mersenne Twister pseudo-random number generator functions.
2
3Copyright 2002, 2003 Free Software Foundation, Inc.
4
5This file is part of the GNU MP Library.
6
7The GNU MP Library is free software; you can redistribute it and/or modify
8it under the terms of the GNU Lesser General Public License as published by
9the Free Software Foundation; either version 3 of the License, or (at your
10option) any later version.
11
12The GNU MP Library is distributed in the hope that it will be useful, but
13WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
15License for more details.
16
17You should have received a copy of the GNU Lesser General Public License
18along 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.  */
27static void
28mangle_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
98static void
99randseed_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
144static 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.  */
152void
153gmp_randinit_mt (gmp_randstate_t rstate)
154{
155  __gmp_randinit_mt_noseed (rstate);
156  RNG_FNPTR (rstate) = (void *) &Mersenne_Twister_Generator;
157}
158