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