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
172
173
174
175
176
177
178
179
180
181
182
183
|
/************************** mersenne.cpp **********************************
* Author: Agner Fog
* Date created: 2001
* Last modified: 2008-11-16
* Project: randomc.h
* Platform: Any C++
* Description:
* Random Number generator of type 'Mersenne Twister'
*
* This random number generator is described in the article by
* M. Matsumoto & T. Nishimura, in:
* ACM Transactions on Modeling and Computer Simulation,
* vol. 8, no. 1, 1998, pp. 3-30.
* Details on the initialization scheme can be found at
* http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
*
* Further documentation:
* The file ran-instructions.pdf contains further documentation and
* instructions.
*
* Copyright 2001-2008 by Agner Fog.
* GNU General Public License http://www.gnu.org/licenses/gpl.html
*******************************************************************************/
#include "randomc.h"
void CRandomMersenne::Init0(int seed) {
// Seed generator
const uint32_t factor = 1812433253UL;
mt[0]= seed;
for (mti=1; mti < MERS_N; mti++) {
mt[mti] = (factor * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
}
}
void CRandomMersenne::RandomInit(int seed) {
// Initialize and seed
Init0(seed);
// Randomize some more
for (int i = 0; i < 37; i++) BRandom();
}
void CRandomMersenne::RandomInitByArray(int const seeds[], int NumSeeds) {
// Seed by more than 32 bits
int i, j, k;
// Initialize
Init0(19650218);
if (NumSeeds <= 0) return;
// Randomize mt[] using whole seeds[] array
i = 1; j = 0;
k = (MERS_N > NumSeeds ? MERS_N : NumSeeds);
for (; k; k--) {
mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + (uint32_t)seeds[j] + j;
i++; j++;
if (i >= MERS_N) {mt[0] = mt[MERS_N-1]; i=1;}
if (j >= NumSeeds) j=0;}
for (k = MERS_N-1; k; k--) {
mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) - i;
if (++i >= MERS_N) {mt[0] = mt[MERS_N-1]; i=1;}}
mt[0] = 0x80000000UL; // MSB is 1; assuring non-zero initial array
// Randomize some more
mti = 0;
for (int i = 0; i <= MERS_N; i++) BRandom();
}
uint32_t CRandomMersenne::BRandom() {
// Generate 32 random bits
uint32_t y;
if (mti >= MERS_N) {
// Generate MERS_N words at one time
const uint32_t LOWER_MASK = (1LU << MERS_R) - 1; // Lower MERS_R bits
const uint32_t UPPER_MASK = 0xFFFFFFFF << MERS_R; // Upper (32 - MERS_R) bits
static const uint32_t mag01[2] = {0, MERS_A};
int kk;
for (kk=0; kk < MERS_N-MERS_M; kk++) {
y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
mt[kk] = mt[kk+MERS_M] ^ (y >> 1) ^ mag01[y & 1];}
for (; kk < MERS_N-1; kk++) {
y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
mt[kk] = mt[kk+(MERS_M-MERS_N)] ^ (y >> 1) ^ mag01[y & 1];}
y = (mt[MERS_N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
mt[MERS_N-1] = mt[MERS_M-1] ^ (y >> 1) ^ mag01[y & 1];
mti = 0;
}
y = mt[mti++];
// Tempering (May be omitted):
y ^= y >> MERS_U;
y ^= (y << MERS_S) & MERS_B;
y ^= (y << MERS_T) & MERS_C;
y ^= y >> MERS_L;
return y;
}
double CRandomMersenne::Random() {
// Output random float number in the interval 0 <= x < 1
// Multiply by 2^(-32)
return (double)BRandom() * (1./(65536.*65536.));
}
int CRandomMersenne::IRandom(int min, int max) {
// Output random integer in the interval min <= x <= max
// Relative error on frequencies < 2^-32
if (max <= min) {
if (max == min) return min; else return 0x80000000;
}
// Multiply interval with random and truncate
int r = int((double)(uint32_t)(max - min + 1) * Random() + min);
if (r > max) r = max;
return r;
}
int CRandomMersenne::IRandomX(int min, int max) {
// Output random integer in the interval min <= x <= max
// Each output value has exactly the same probability.
// This is obtained by rejecting certain bit values so that the number
// of possible bit values is divisible by the interval length
if (max <= min) {
if (max == min) return min; else return 0x80000000;
}
#ifdef INT64_SUPPORTED
// 64 bit integers available. Use multiply and shift method
uint32_t interval; // Length of interval
uint64_t longran; // Random bits * interval
uint32_t iran; // Longran / 2^32
uint32_t remainder; // Longran % 2^32
interval = uint32_t(max - min + 1);
if (interval != LastInterval) {
// Interval length has changed. Must calculate rejection limit
// Reject when remainder >= 2^32 / interval * interval
// RLimit will be 0 if interval is a power of 2. No rejection then
RLimit = uint32_t(((uint64_t)1 << 32) / interval) * interval - 1;
LastInterval = interval;
}
do { // Rejection loop
longran = (uint64_t)BRandom() * interval;
iran = (uint32_t)(longran >> 32);
remainder = (uint32_t)longran;
} while (remainder > RLimit);
// Convert back to signed and return result
return (int32_t)iran + min;
#else
// 64 bit integers not available. Use modulo method
uint32_t interval; // Length of interval
uint32_t bran; // Random bits
uint32_t iran; // bran / interval
uint32_t remainder; // bran % interval
interval = uint32_t(max - min + 1);
if (interval != LastInterval) {
// Interval length has changed. Must calculate rejection limit
// Reject when iran = 2^32 / interval
// We can't make 2^32 so we use 2^32-1 and correct afterwards
RLimit = (uint32_t)0xFFFFFFFF / interval;
if ((uint32_t)0xFFFFFFFF % interval == interval - 1) RLimit++;
}
do { // Rejection loop
bran = BRandom();
iran = bran / interval;
remainder = bran % interval;
} while (iran >= RLimit);
// Convert back to signed and return result
return (int32_t)remainder + min;
#endif
}
|
