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