1 /************************** MERSENNE.CPP ******************** AgF 2007-08-01 *
\r
2 * Random Number generator 'Mersenne Twister' *
\r
4 * This random number generator is described in the article by *
\r
5 * M. Matsumoto & T. Nishimura, in: *
\r
6 * ACM Transactions on Modeling and Computer Simulation, *
\r
7 * vol. 8, no. 1, 1998, pp. 3-30. *
\r
8 * Details on the initialization scheme can be found at *
\r
9 * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html *
\r
11 * Experts consider this an excellent random number generator. *
\r
13 * © 2001 - 2007 A. Fog. Published under the GNU General Public License *
\r
14 * (www.gnu.org/copyleft/gpl.html) with the further restriction that it *
\r
15 * cannot be used for gambling applications. *
\r
16 *****************************************************************************/
\r
18 #include "randomc.h"
\r
20 void CRandomMersenne::Init0(uint32 seed) {
\r
21 // Detect computer architecture
\r
22 union {double f; uint32 i[2];} convert;
\r
24 if (convert.i[1] == 0x3FF00000) Architecture = LITTLE_ENDIAN1;
\r
25 else if (convert.i[0] == 0x3FF00000) Architecture = BIG_ENDIAN1;
\r
26 else Architecture = NONIEEE;
\r
30 for (mti=1; mti < MERS_N; mti++) {
\r
31 mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
\r
35 void CRandomMersenne::RandomInit(uint32 seed) {
\r
36 // Initialize and seed
\r
39 // Randomize some more
\r
40 for (int i = 0; i < 37; i++) BRandom();
\r
44 void CRandomMersenne::RandomInitByArray(uint32 seeds[], int length) {
\r
45 // Seed by more than 32 bits
\r
51 if (length <= 0) return;
\r
53 // Randomize mt[] using whole seeds[] array
\r
55 k = (MERS_N > length ? MERS_N : length);
\r
57 mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + seeds[j] + j;
\r
59 if (i >= MERS_N) {mt[0] = mt[MERS_N-1]; i=1;}
\r
60 if (j >= length) j=0;}
\r
61 for (k = MERS_N-1; k; k--) {
\r
62 mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) - i;
\r
63 if (++i >= MERS_N) {mt[0] = mt[MERS_N-1]; i=1;}}
\r
64 mt[0] = 0x80000000UL; // MSB is 1; assuring non-zero initial array
\r
66 // Randomize some more
\r
68 for (int i = 0; i <= MERS_N; i++) BRandom();
\r
72 uint32 CRandomMersenne::BRandom() {
\r
73 // Generate 32 random bits
\r
76 if (mti >= MERS_N) {
\r
77 // Generate MERS_N words at one time
\r
78 const uint32 LOWER_MASK = (1LU << MERS_R) - 1; // Lower MERS_R bits
\r
79 const uint32 UPPER_MASK = 0xFFFFFFFF << MERS_R; // Upper (32 - MERS_R) bits
\r
80 static const uint32 mag01[2] = {0, MERS_A};
\r
83 for (kk=0; kk < MERS_N-MERS_M; kk++) {
\r
84 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
\r
85 mt[kk] = mt[kk+MERS_M] ^ (y >> 1) ^ mag01[y & 1];}
\r
87 for (; kk < MERS_N-1; kk++) {
\r
88 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
\r
89 mt[kk] = mt[kk+(MERS_M-MERS_N)] ^ (y >> 1) ^ mag01[y & 1];}
\r
91 y = (mt[MERS_N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
\r
92 mt[MERS_N-1] = mt[MERS_M-1] ^ (y >> 1) ^ mag01[y & 1];
\r
99 // Tempering (May be omitted):
\r
101 y ^= (y << MERS_S) & MERS_B;
\r
102 y ^= (y << MERS_T) & MERS_C;
\r
110 double CRandomMersenne::Random() {
\r
111 // Output random float number in the interval 0 <= x < 1
\r
112 union {double f; uint32 i[2];} convert;
\r
113 uint32 r = BRandom(); // Get 32 random bits
\r
114 // The fastest way to convert random bits to floating point is as follows:
\r
115 // Set the binary exponent of a floating point number to 1+bias and set
\r
116 // the mantissa to random bits. This will give a random number in the
\r
117 // interval [1,2). Then subtract 1.0 to get a random number in the interval
\r
118 // [0,1). This procedure requires that we know how floating point numbers
\r
119 // are stored. The storing method is tested in function RandomInit and saved
\r
120 // in the variable Architecture.
\r
122 // This shortcut allows the compiler to optimize away the following switch
\r
123 // statement for the most common architectures:
\r
124 #if defined(_M_IX86) || defined(_M_X64) || defined(__LITTLE_ENDIAN__)
\r
125 Architecture = LITTLE_ENDIAN1;
\r
126 #elif defined(__BIG_ENDIAN__)
\r
127 Architecture = BIG_ENDIAN1;
\r
130 switch (Architecture) {
\r
131 case LITTLE_ENDIAN1:
\r
132 convert.i[0] = r << 20;
\r
133 convert.i[1] = (r >> 12) | 0x3FF00000;
\r
134 return convert.f - 1.0;
\r
136 convert.i[1] = r << 20;
\r
137 convert.i[0] = (r >> 12) | 0x3FF00000;
\r
138 return convert.f - 1.0;
\r
139 case NONIEEE: default: ;
\r
141 // This somewhat slower method works for all architectures, including
\r
142 // non-IEEE floating point representation:
\r
143 return (double)r * (1./((double)(uint32)(-1L)+1.));
\r
147 int CRandomMersenne::IRandom(int min, int max) {
\r
148 // Output random integer in the interval min <= x <= max
\r
149 // Relative error on frequencies < 2^-32
\r
151 if (max == min) return min; else return 0x80000000;
\r
153 // Multiply interval with random and truncate
\r
154 int r = int((max - min + 1) * Random()) + min;
\r
155 if (r > max) r = max;
\r
160 int CRandomMersenne::IRandomX(int min, int max) {
\r
161 // Output random integer in the interval min <= x <= max
\r
162 // Each output value has exactly the same probability.
\r
163 // This is obtained by rejecting certain bit values so that the number
\r
164 // of possible bit values is divisible by the interval length
\r
166 if (max == min) return min; else return 0x80000000;
\r
168 #ifdef INT64_DEFINED
\r
169 // 64 bit integers available. Use multiply and shift method
\r
170 uint32 interval; // Length of interval
\r
171 uint64 longran; // Random bits * interval
\r
172 uint32 iran; // Longran / 2^32
\r
173 uint32 remainder; // Longran % 2^32
\r
175 interval = uint32(max - min + 1);
\r
176 if (interval != LastInterval) {
\r
177 // Interval length has changed. Must calculate rejection limit
\r
178 // Reject when remainder = 2^32 / interval * interval
\r
179 // RLimit will be 0 if interval is a power of 2. No rejection then
\r
180 RLimit = uint32(((uint64)1 << 32) / interval) * interval - 1;
\r
181 LastInterval = interval;
\r
183 do { // Rejection loop
\r
184 longran = (uint64)BRandom() * interval;
\r
185 iran = (uint32)(longran >> 32);
\r
186 remainder = (uint32)longran;
\r
187 } while (remainder > RLimit);
\r
188 // Convert back to signed and return result
\r
189 return (int32)iran + min;
\r
192 // 64 bit integers not available. Use modulo method
\r
193 uint32 interval; // Length of interval
\r
194 uint32 bran; // Random bits
\r
195 uint32 iran; // bran / interval
\r
196 uint32 remainder; // bran % interval
\r
198 interval = uint32(max - min + 1);
\r
199 if (interval != LastInterval) {
\r
200 // Interval length has changed. Must calculate rejection limit
\r
201 // Reject when iran = 2^32 / interval
\r
202 // We can't make 2^32 so we use 2^32-1 and correct afterwards
\r
203 RLimit = (uint32)0xFFFFFFFF / interval;
\r
204 if ((uint32)0xFFFFFFFF % interval == interval - 1) RLimit++;
\r
206 do { // Rejection loop
\r
208 iran = bran / interval;
\r
209 remainder = bran % interval;
\r
210 } while (iran >= RLimit);
\r
211 // Convert back to signed and return result
\r
212 return (int32)remainder + min;
\r