1 /************************** MERSENNE.CPP ******************** AgF 2001-10-18 * 2 * Random Number generator 'Mersenne Twister' * 3 * * 4 * This random number generator is described in the article by * 5 * M. Matsumoto & T. Nishimura, in: * 6 * ACM Transactions on Modeling and Computer Simulation, * 7 * vol. 8, no. 1, 1998, pp. 3-30. * 8 * * 9 * Experts consider this an excellent random number generator. * 10 * * 11 * � 2002 A. Fog. GNU General Public License www.gnu.org/copyleft/gpl.html * 12 *****************************************************************************/ 13 14 #include "randomc.h" 15 16 RandomInit(long int seed)17void TRandomMersenne::RandomInit(long int seed) { 18 // re-seed generator 19 unsigned long s = (unsigned long)seed; 20 for (mti=0; mti<N; mti++) { 21 s = s * 29943829 - 1; 22 mt[mti] = s;}} 23 24 BRandom()25unsigned long TRandomMersenne::BRandom() { 26 // generate 32 random bits 27 unsigned long y; 28 29 if (mti >= N) { 30 // generate N words at one time 31 const unsigned long LOWER_MASK = (1u << R) - 1; // lower R bits 32 const unsigned long UPPER_MASK = -1 << R; // upper 32-R bits 33 int kk, km; 34 for (kk=0, km=M; kk < N-1; kk++) { 35 y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK); 36 mt[kk] = mt[km] ^ (y >> 1) ^ (-(signed long)(y & 1) & MATRIX_A); 37 if (++km >= N) km = 0;} 38 39 y = (mt[N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK); 40 mt[N-1] = mt[M-1] ^ (y >> 1) ^ (-(signed long)(y & 1) & MATRIX_A); 41 mti = 0;} 42 43 y = mt[mti++]; 44 45 // Tempering (May be omitted): 46 y ^= y >> TEMU; 47 y ^= (y << TEMS) & TEMB; 48 y ^= (y << TEMT) & TEMC; 49 y ^= y >> TEML; 50 51 return y;} 52 53 Random()54double TRandomMersenne::Random() { 55 // output random float number in the interval 0 <= x < 1 56 union { 57 double f; 58 unsigned long i[2];} 59 convert; 60 // get 32 random bits and convert to float 61 unsigned long r = BRandom(); 62 convert.i[0] = r << 20; 63 convert.i[1] = (r >> 12) | 0x3FF00000; 64 return convert.f - 1.0;} 65 66 IRandom(long min,long max)67long TRandomMersenne::IRandom(long min, long max) { 68 // output random integer in the interval min <= x <= max 69 long r; 70 r = long((max - min + 1) * Random()) + min; // multiply interval with random and truncate 71 if (r > max) r = max; 72 if (max < min) return 0x80000000; 73 return r;} 74 75