QGIS API Documentation  2.2.0-Valmiera
mersenne-twister.cpp
Go to the documentation of this file.
1 /*
2  * The Mersenne Twister pseudo-random number generator (PRNG)
3  *
4  * This is an implementation of fast PRNG called MT19937,
5  * meaning it has a period of 2^19937-1, which is a Mersenne
6  * prime.
7  *
8  * This PRNG is fast and suitable for non-cryptographic code.
9  * For instance, it would be perfect for Monte Carlo simulations,
10  * etc.
11  *
12  * For all the details on this algorithm, see the original paper:
13  * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf
14  *
15  * Written by Christian Stigen Larsen
16  * http://csl.sublevel3.org
17  *
19  *
20  * 2012-01-11
21  */
22
23 #include <stdio.h>
24 #include "mersenne-twister.h"
25
26 /*
27  * We have an array of 624 32-bit values, and there are
28  * 31 unused bits, so we have a seed value of
29  * 624*32-31 = 19937 bits.
30  */
31 static const unsigned SIZE = 624;
32 static const unsigned PERIOD = 397;
33 static const unsigned DIFF = SIZE - PERIOD;
34
35 static uint32_t MT[SIZE];
36 static unsigned index = 0;
37
38 #define M32(x) (0x80000000 & x) // 32nd Most Significant Bit
39 #define L31(x) (0x7FFFFFFF & x) // 31 Least Significant Bits
40 #define ODD(x) (x & 1) // Check if number is odd
41
42 #define MD_UINT32_MAX std::numeric_limits<uint32_t>::max()
43
44 static inline void generate_numbers()
45 {
46  /*
47  * Originally, we had one loop with i going from [0, SIZE) and
48  * two modulus operations:
49  *
50  * for ( register unsigned i=0; i<SIZE; ++i ) {
51  * register uint32_t y = M32(MT[i]) | L31(MT[(i+1) % SIZE]);
52  * MT[i] = MT[(i + PERIOD) % SIZE] ^ (y>>1);
53  * if ( ODD(y) ) MT[i] ^= 0x9908b0df;
54  * }
55  *
56  * For performance reasons, we've unrolled the loop three times, thus
57  * mitigating the need for any modulus operations.
58  *
59  * Anyway, it seems this trick is old hat:
61  *
62  */
63
64  static const uint32_t MATRIX[2] = {0, 0x9908b0df};
65  register uint32_t y, i;
66
67  // i = [0 ... 226]
68  for ( i = 0; i < DIFF; ++i )
69  {
70  /*
71  * We're doing 226 = 113*2, an even number of steps, so we can
72  * safely unroll one more step here for speed:
73  */
74  y = M32( MT[i] ) | L31( MT[i+1] );
75  MT[i] = MT[i+PERIOD] ^( y >> 1 ) ^ MATRIX[ODD( y )];
76
77  ++i;
78  y = M32( MT[i] ) | L31( MT[i+1] );
79  MT[i] = MT[i+PERIOD] ^( y >> 1 ) ^ MATRIX[ODD( y )];
80  }
81
82 #define UNROLL \
83  y = M32(MT[i]) | L31(MT[i+1]); \
84  MT[i] = MT[i-DIFF] ^ (y>>1) ^ MATRIX[ODD(y)]; \
85  ++i;
86
87  // i = [227 ... 622]
88  for ( i = DIFF; i < ( SIZE - 1 ); )
89  {
90  /*
91  * 623-227 = 396 = 2*2*3*3*11, so we can unroll this loop in any
92  * number that evenly divides 396 (2, 4, 6, etc).
93  */
94
95  // 11 times
98
100  UNROLL; UNROLL;
101  }
102
103  // i = [623]
104  y = M32( MT[SIZE-1] ) | L31( MT[SIZE-1] );
105  MT[SIZE-1] = MT[PERIOD-1] ^( y >> 1 ) ^ MATRIX[ODD( y )];
106 }
107
108 extern "C" void seed( uint32_t value )
109 {
110  /*
111  * The equation below is a linear congruential generator (LCG),
112  * one of the oldest known pseudo-random number generator
113  * algorithms, in the form X_(n+1) = = (a*X_n + c) (mod m).
114  *
115  * We've implicitly got m=32 (mask + word size of 32 bits), so
116  * there is no need to explicitly use modulus.
117  *
118  * What is interesting is the multiplier a. The one we have
119  * below is 0x6c07865 --- 1812433253 in decimal, and is called
120  * the Borosh-Niederreiter multiplier for modulus 2^32.
121  *
122  * It is mentioned in passing in Knuth's THE ART OF COMPUTER
123  * PROGRAMMING, Volume 2, page 106, Table 1, line 13. LCGs are
124  * treated in the same book, pp. 10-26
125  *
126  * You can read the original paper by Borosh and Niederreiter
127  * as well. It's called OPTIMAL MULTIPLIERS FOR PSEUDO-RANDOM
128  * NUMBER GENERATION BY THE LINEAR CONGRUENTIAL METHOD (1983) at
130  *
132  * http://en.wikipedia.org/wiki/Linear_congruential_generator
133  *
134  * From that page, it says:
135  * "A common Mersenne twister implementation, interestingly
136  * enough, uses an LCG to generate seed data.",
137  *
138  * Since our we're using 32-bits data types for our MT array,
139  * we can skip the masking with 0xFFFFFFFF below.
140  */
141
142  MT[0] = value;
143  index = 0;
144
145  for ( register unsigned i = 1; i < SIZE; ++i )
146  MT[i] = 0x6c078965 * ( MT[i-1] ^ MT[i-1] >> 30 ) + i;
147 }
148
149 extern "C" uint32_t rand_u32()
150 {
151  if ( !index )
153
154  register uint32_t y = MT[index];
155
156  // Tempering
157  y ^= y >> 11;
158  y ^= y << 7 & 0x9d2c5680;
159  y ^= y << 15 & 0xefc60000;
160  y ^= y >> 18;
161
162  if ( ++index == SIZE )
163  index = 0;
164
165  return y;
166 }
167
168 extern "C" int mt_rand()
169 {
170  /*
171  * PORTABILITY WARNING:
172  *
173  * rand_u32() uses all 32-bits for the pseudo-random number,
174  * but rand() must return a number from 0 ... MD_RAND_MAX.
175  *
176  * We'll just assume that rand() only uses 31 bits worth of
177  * data, and that we're on a two's complement system.
178  *
179  * So, to output an integer compatible with rand(), we have
180  * two options: Either mask off the highest (32nd) bit, or
181  * shift right by one bit. Masking with 0x7FFFFFFF will be
182  * compatible with 64-bit MT[], so we'll just use that here.
183  *
184  */
185  return static_cast<int>( 0x7FFFFFFF & rand_u32() );
186 }
187
188 extern "C" void mt_srand( unsigned value )
189 {
190  seed( static_cast<uint32_t>( value ) );
191 }
192
193 extern "C" float randf_cc()
194 {
195  return static_cast<float>( rand_u32() ) / MD_UINT32_MAX;
196 }
197
198 extern "C" float randf_co()
199 {
200  return static_cast<float>( rand_u32() ) / ( MD_UINT32_MAX + 1.0f );
201 }
202
203 extern "C" float randf_oo()
204 {
205  return ( static_cast<float>( rand_u32() ) + 0.5f ) / ( MD_UINT32_MAX + 1.0f );
206 }
207
208 extern "C" double randd_cc()
209 {
210  return static_cast<double>( rand_u32() ) / MD_UINT32_MAX;
211 }
212
213 extern "C" double randd_co()
214 {
215  return static_cast<double>( rand_u32() ) / ( MD_UINT32_MAX + 1.0 );
216 }
217
218 extern "C" double randd_oo()
219 {
220  return ( static_cast<double>( rand_u32() ) + 0.5 ) / ( MD_UINT32_MAX + 1.0 );
221 }
222
223 extern "C" uint64_t rand_u64()
224 {
225  return static_cast<uint64_t>( rand_u32() ) << 32 | rand_u32();
226 }