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
|
/* Mersenne Twister pseudo-random number generator functions.
Copyright 2002, 2003 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 the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
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 Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
#include "gmp.h"
#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_srcptr b_orig)
{
mpz_t t, b;
unsigned long e = 0x40118124;
unsigned long bit = 0x20000000;
mpz_init (t);
mpz_init_set (b, b_orig); /* in case r==b_orig */
mpz_set (r, b);
do
{
mpz_mul (r, r, r);
reduce:
for (;;)
{
mpz_tdiv_q_2exp (t, r, 19937L);
if (mpz_sgn (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.
*/
static void
randseed_mt (gmp_randstate_t 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_init (mod);
mpz_init (seed1);
mpz_set_ui (mod, 0L);
mpz_setbit (mod, 19937L);
mpz_sub_ui (mod, mod, 20027L);
mpz_mod (seed1, seed, mod); /* Reduce `seed' modulo `mod'. */
mpz_add_ui (seed1, seed1, 2L); /* seed1 is now ready. */
mangle_seed (seed1, 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);
cnt++;
ASSERT (cnt <= N);
while (cnt < N)
p->mt[cnt++] = 0;
mpz_clear (mod);
mpz_clear (seed1);
/* 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_t rstate)
{
__gmp_randinit_mt_noseed (rstate);
RNG_FNPTR (rstate) = (void *) &Mersenne_Twister_Generator;
}
|
