1 /* Mersenne Twister pseudo-random number generator functions.
2
3 Copyright 2002, 2003, 2013, 2014 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 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
14 or
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
20 or both in parallel, as here.
21
22 The GNU MP Library is distributed in the hope that it will be useful, but
23 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
24 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
25 for more details.
26
27 You should have received copies of the GNU General Public License and the
28 GNU Lesser General Public License along with the GNU MP Library. If not,
29 see 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. */
37 static void
38 mangle_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 Moreover, c) is not guaranted, there are many seeds yielding to the
107 same sequence, because gcd (1074888996, 2^19937 - 20023 - 1) = 12.
108 E.g. x and x'=x*19^((2^19937-20023-1) / 12) mod (2^19937-20023), if
109 chosen as seed1, generate the same seed2, for every x.
110 Similarly x" can be obtained from x', obtaining 12 different
111 values.
112 */
113
114 static void
115 randseed_mt (gmp_randstate_ptr rstate, mpz_srcptr seed)
116 {
117 int i;
118 size_t cnt;
119
120 gmp_rand_mt_struct *p;
121 mpz_t mod; /* Modulus. */
122 mpz_t seed1; /* Intermediate result. */
123
124 p = (gmp_rand_mt_struct *) RNG_STATE (rstate);
125
126 mpz_init2 (mod, 19938L);
127 mpz_init2 (seed1, 19937L);
128
129 mpz_setbit (mod, 19937L);
130 mpz_sub_ui (mod, mod, 20027L);
131 mpz_mod (seed1, seed, mod); /* Reduce `seed' modulo `mod'. */
132 mpz_clear (mod);
133 mpz_add_ui (seed1, seed1, 2L); /* seed1 is now ready. */
134 mangle_seed (seed1); /* Perform the mangling by powering. */
135
136 /* Copy the last bit into bit 31 of mt[0] and clear it. */
137 p->mt[0] = (mpz_tstbit (seed1, 19936L) != 0) ? 0x80000000 : 0;
138 mpz_clrbit (seed1, 19936L);
139
140 /* Split seed1 into N-1 32-bit chunks. */
141 mpz_export (&p->mt[1], &cnt, -1, sizeof (p->mt[1]), 0,
142 8 * sizeof (p->mt[1]) - 32, seed1);
143 mpz_clear (seed1);
144 cnt++;
145 ASSERT (cnt <= N);
146 while (cnt < N)
147 p->mt[cnt++] = 0;
148
149 /* Warm the generator up if necessary. */
150 if (WARM_UP != 0)
151 for (i = 0; i < WARM_UP / N; i++)
152 __gmp_mt_recalc_buffer (p->mt);
153
154 p->mti = WARM_UP % N;
155 }
156
157
158 static const gmp_randfnptr_t Mersenne_Twister_Generator = {
159 randseed_mt,
160 __gmp_randget_mt,
161 __gmp_randclear_mt,
162 __gmp_randiset_mt
163 };
164
165 /* Initialize MT-specific data. */
166 void
167 gmp_randinit_mt (gmp_randstate_ptr rstate)
168 {
169 __gmp_randinit_mt_noseed (rstate);
170 RNG_FNPTR (rstate) = (void *) &Mersenne_Twister_Generator;
171 }