-/************************** MERSENNE.CPP ******************** AgF 2007-08-01 *\r
-* Random Number generator 'Mersenne Twister' *\r
-* *\r
-* This random number generator is described in the article by *\r
-* M. Matsumoto & T. Nishimura, in: *\r
-* ACM Transactions on Modeling and Computer Simulation, *\r
-* vol. 8, no. 1, 1998, pp. 3-30. *\r
-* Details on the initialization scheme can be found at *\r
-* http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html *\r
-* *\r
-* Experts consider this an excellent random number generator. *\r
-* *\r
-* © 2001 - 2007 A. Fog. Published under the GNU General Public License *\r
-* (www.gnu.org/copyleft/gpl.html) with the further restriction that it *\r
-* cannot be used for gambling applications. *\r
-*****************************************************************************/\r
-\r
-#include "randomc.h"\r
-\r
-void CRandomMersenne::Init0(uint32 seed) {\r
- // Detect computer architecture\r
- union {double f; uint32 i[2];} convert;\r
- convert.f = 1.0;\r
- if (convert.i[1] == 0x3FF00000) Architecture = LITTLE_ENDIAN1;\r
- else if (convert.i[0] == 0x3FF00000) Architecture = BIG_ENDIAN1;\r
- else Architecture = NONIEEE;\r
-\r
- // Seed generator\r
- mt[0]= seed;\r
- for (mti=1; mti < MERS_N; mti++) {\r
- mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);\r
- }\r
-}\r
-\r
-void CRandomMersenne::RandomInit(uint32 seed) {\r
- // Initialize and seed\r
- Init0(seed);\r
-\r
- // Randomize some more\r
- for (int i = 0; i < 37; i++) BRandom();\r
-}\r
-\r
-\r
-void CRandomMersenne::RandomInitByArray(uint32 seeds[], int length) {\r
- // Seed by more than 32 bits\r
- int i, j, k;\r
-\r
- // Initialize\r
- Init0(19650218);\r
-\r
- if (length <= 0) return;\r
-\r
- // Randomize mt[] using whole seeds[] array\r
- i = 1; j = 0;\r
- k = (MERS_N > length ? MERS_N : length);\r
- for (; k; k--) {\r
- mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + seeds[j] + j;\r
- i++; j++;\r
- if (i >= MERS_N) {mt[0] = mt[MERS_N-1]; i=1;}\r
- if (j >= length) j=0;}\r
- for (k = MERS_N-1; k; k--) {\r
- mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) - i;\r
- if (++i >= MERS_N) {mt[0] = mt[MERS_N-1]; i=1;}}\r
- mt[0] = 0x80000000UL; // MSB is 1; assuring non-zero initial array\r
-\r
- // Randomize some more\r
- mti = 0;\r
- for (int i = 0; i <= MERS_N; i++) BRandom();\r
-}\r
-\r
-\r
-uint32 CRandomMersenne::BRandom() {\r
- // Generate 32 random bits\r
- uint32 y;\r
-\r
- if (mti >= MERS_N) {\r
- // Generate MERS_N words at one time\r
- const uint32 LOWER_MASK = (1LU << MERS_R) - 1; // Lower MERS_R bits\r
- const uint32 UPPER_MASK = 0xFFFFFFFF << MERS_R; // Upper (32 - MERS_R) bits\r
- static const uint32 mag01[2] = {0, MERS_A};\r
-\r
- int kk;\r
- for (kk=0; kk < MERS_N-MERS_M; kk++) { \r
- y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);\r
- mt[kk] = mt[kk+MERS_M] ^ (y >> 1) ^ mag01[y & 1];}\r
-\r
- for (; kk < MERS_N-1; kk++) { \r
- y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);\r
- mt[kk] = mt[kk+(MERS_M-MERS_N)] ^ (y >> 1) ^ mag01[y & 1];} \r
-\r
- y = (mt[MERS_N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);\r
- mt[MERS_N-1] = mt[MERS_M-1] ^ (y >> 1) ^ mag01[y & 1];\r
- mti = 0;\r
- }\r
-\r
- y = mt[mti++];\r
-\r
-#if 1\r
- // Tempering (May be omitted):\r
- y ^= y >> MERS_U;\r
- y ^= (y << MERS_S) & MERS_B;\r
- y ^= (y << MERS_T) & MERS_C;\r
- y ^= y >> MERS_L;\r
-#endif\r
-\r
- return y;\r
-}\r
-\r
-\r
-double CRandomMersenne::Random() {\r
- // Output random float number in the interval 0 <= x < 1\r
- union {double f; uint32 i[2];} convert;\r
- uint32 r = BRandom(); // Get 32 random bits\r
- // The fastest way to convert random bits to floating point is as follows:\r
- // Set the binary exponent of a floating point number to 1+bias and set\r
- // the mantissa to random bits. This will give a random number in the \r
- // interval [1,2). Then subtract 1.0 to get a random number in the interval\r
- // [0,1). This procedure requires that we know how floating point numbers\r
- // are stored. The storing method is tested in function RandomInit and saved \r
- // in the variable Architecture.\r
-\r
- // This shortcut allows the compiler to optimize away the following switch\r
- // statement for the most common architectures:\r
-#if defined(_M_IX86) || defined(_M_X64) || defined(__LITTLE_ENDIAN__)\r
- Architecture = LITTLE_ENDIAN1;\r
-#elif defined(__BIG_ENDIAN__)\r
- Architecture = BIG_ENDIAN1;\r
-#endif\r
-\r
- switch (Architecture) {\r
- case LITTLE_ENDIAN1:\r
- convert.i[0] = r << 20;\r
- convert.i[1] = (r >> 12) | 0x3FF00000;\r
- return convert.f - 1.0;\r
- case BIG_ENDIAN1:\r
- convert.i[1] = r << 20;\r
- convert.i[0] = (r >> 12) | 0x3FF00000;\r
- return convert.f - 1.0;\r
- case NONIEEE: default: ;\r
- } \r
- // This somewhat slower method works for all architectures, including \r
- // non-IEEE floating point representation:\r
- return (double)r * (1./((double)(uint32)(-1L)+1.));\r
-}\r
-\r
-\r
-int CRandomMersenne::IRandom(int min, int max) {\r
- // Output random integer in the interval min <= x <= max\r
- // Relative error on frequencies < 2^-32\r
- if (max <= min) {\r
- if (max == min) return min; else return 0x80000000;\r
- }\r
- // Multiply interval with random and truncate\r
- int r = int((max - min + 1) * Random()) + min; \r
- if (r > max) r = max;\r
- return r;\r
-}\r
-\r
-\r
-int CRandomMersenne::IRandomX(int min, int max) {\r
- // Output random integer in the interval min <= x <= max\r
- // Each output value has exactly the same probability.\r
- // This is obtained by rejecting certain bit values so that the number\r
- // of possible bit values is divisible by the interval length\r
- if (max <= min) {\r
- if (max == min) return min; else return 0x80000000;\r
- }\r
-#ifdef INT64_DEFINED\r
- // 64 bit integers available. Use multiply and shift method\r
- uint32 interval; // Length of interval\r
- uint64 longran; // Random bits * interval\r
- uint32 iran; // Longran / 2^32\r
- uint32 remainder; // Longran % 2^32\r
-\r
- interval = uint32(max - min + 1);\r
- if (interval != LastInterval) {\r
- // Interval length has changed. Must calculate rejection limit\r
- // Reject when remainder = 2^32 / interval * interval\r
- // RLimit will be 0 if interval is a power of 2. No rejection then\r
- RLimit = uint32(((uint64)1 << 32) / interval) * interval - 1;\r
- LastInterval = interval;\r
- }\r
- do { // Rejection loop\r
- longran = (uint64)BRandom() * interval;\r
- iran = (uint32)(longran >> 32);\r
- remainder = (uint32)longran;\r
- } while (remainder > RLimit);\r
- // Convert back to signed and return result\r
- return (int32)iran + min;\r
-\r
-#else\r
- // 64 bit integers not available. Use modulo method\r
- uint32 interval; // Length of interval\r
- uint32 bran; // Random bits\r
- uint32 iran; // bran / interval\r
- uint32 remainder; // bran % interval\r
-\r
- interval = uint32(max - min + 1);\r
- if (interval != LastInterval) {\r
- // Interval length has changed. Must calculate rejection limit\r
- // Reject when iran = 2^32 / interval\r
- // We can't make 2^32 so we use 2^32-1 and correct afterwards\r
- RLimit = (uint32)0xFFFFFFFF / interval;\r
- if ((uint32)0xFFFFFFFF % interval == interval - 1) RLimit++;\r
- }\r
- do { // Rejection loop\r
- bran = BRandom();\r
- iran = bran / interval;\r
- remainder = bran % interval;\r
- } while (iran >= RLimit);\r
- // Convert back to signed and return result\r
- return (int32)remainder + min;\r
-\r
-#endif\r
-}\r