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)17 void 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()25 unsigned 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()54 double 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)67 long 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