mersenne.cpp

   1 /************************** MERSENNE.CPP ******************** AgF 2007-08-01 *\r
   2 *  Random Number generator 'Mersenne Twister'                                *\r
   3 *                                                                            *\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
  10 *                                                                            *\r
  11 *  Experts consider this an excellent random number generator.               *\r
  12 *                                                                            *\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
  17 \r
  18 #include "randomc.h"\r
  19 \r
  20 void CRandomMersenne::Init0(uint32 seed) {\r
  21    // Detect computer architecture\r
  22    union {double f; uint32 i[2];} convert;\r
  23    convert.f = 1.0;\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
  27 \r
  28    // Seed generator\r
  29    mt[0]= seed;\r
  30    for (mti=1; mti < MERS_N; mti++) {\r
  31       mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);\r
  32    }\r
  33 }\r
  34 \r
  35 void CRandomMersenne::RandomInit(uint32 seed) {\r
  36    // Initialize and seed\r
  37    Init0(seed);\r
  38 \r
  39    // Randomize some more\r
  40    for (int i = 0; i < 37; i++) BRandom();\r
  41 }\r
  42 \r
  43 \r
  44 void CRandomMersenne::RandomInitByArray(uint32 seeds[], int length) {\r
  45    // Seed by more than 32 bits\r
  46    int i, j, k;\r
  47 \r
  48    // Initialize\r
  49    Init0(19650218);\r
  50 \r
  51    if (length <= 0) return;\r
  52 \r
  53    // Randomize mt[] using whole seeds[] array\r
  54    i = 1;  j = 0;\r
  55    k = (MERS_N > length ? MERS_N : length);\r
  56    for (; k; k--) {\r
  57       mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + seeds[j] + j;\r
  58       i++; 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
  65 \r
  66    // Randomize some more\r
  67    mti = 0;\r
  68    for (int i = 0; i <= MERS_N; i++) BRandom();\r
  69 }\r
  70 \r
  71 \r
  72 uint32 CRandomMersenne::BRandom() {\r
  73    // Generate 32 random bits\r
  74    uint32 y;\r
  75 \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
  81 \r
  82       int kk;\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
  86 \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
  90 \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
  93       mti = 0;\r
  94    }\r
  95 \r
  96    y = mt[mti++];\r
  97 \r
  98 #if 1\r
  99    // Tempering (May be omitted):\r
 100    y ^=  y >> MERS_U;\r
 101    y ^= (y << MERS_S) & MERS_B;\r
 102    y ^= (y << MERS_T) & MERS_C;\r
 103    y ^=  y >> MERS_L;\r
 104 #endif\r
 105 \r
 106    return y;\r
 107 }\r
 108 \r
 109 \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
 121 \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
 128 #endif\r
 129 \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
 135    case BIG_ENDIAN1:\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
 140    } \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
 144 }\r
 145 \r
 146 \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
 150    if (max <= min) {\r
 151       if (max == min) return min; else return 0x80000000;\r
 152    }\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
 156    return r;\r
 157 }\r
 158 \r
 159 \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
 165    if (max <= min) {\r
 166       if (max == min) return min; else return 0x80000000;\r
 167    }\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
 174 \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
 182    }\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
 190 \r
 191 #else\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
 197 \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
 205    }\r
 206    do { // Rejection loop\r
 207       bran = BRandom();\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
 213 \r
 214 #endif\r
 215 }\r