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
mangle_seed(mpz_ptr r)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
107 static void
randseed_mt(gmp_randstate_t rstate,mpz_srcptr seed)108 randseed_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
151 static 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. */
159 void
gmp_randinit_mt(gmp_randstate_t rstate)160 gmp_randinit_mt (gmp_randstate_t rstate)
161 {
162 __gmp_randinit_mt_noseed (rstate);
163 RNG_FNPTR (rstate) = (void *) &Mersenne_Twister_Generator;
164 }
165