QGIS API Documentation  3.8.0-Zanzibar (11aff65)
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.name
17  *
18  * Distributed under the modified BSD license.
19  *
20  * 2015-02-17
21  */
22 
23 #include <cstdio>
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 UNROLL(expr) \
43  y = M32(MT[i]) | L31(MT[i+1]); \
44  MT[i] = MT[expr] ^ (y>>1) ^ MATRIX[ODD(y)]; \
45  ++i;
46 
47 #define MD_UINT32_MAX std::numeric_limits<uint32_t>::max()
48 
49 static inline void generate_numbers()
50 {
51  /*
52  * Originally, we had one loop with i going from [0, SIZE) and
53  * two modulus operations:
54  *
55  * for ( register unsigned i=0; i<SIZE; ++i ) {
56  * register uint32_t y = M32(MT[i]) | L31(MT[(i+1) % SIZE]);
57  * MT[i] = MT[(i + PERIOD) % SIZE] ^ (y>>1);
58  * if ( ODD(y) ) MT[i] ^= 0x9908b0df;
59  * }
60  *
61  * For performance reasons, we've unrolled the loop three times, thus
62  * mitigating the need for any modulus operations.
63  *
64  * Anyway, it seems this trick is old hat:
65  * http://www.quadibloc.com/crypto/co4814.htm
66  *
67  */
68 
69  static const uint32_t MATRIX[2] = {0, 0x9908b0df};
70  uint32_t y, i = 0;
71 
72  // i = [0 ... 226]
73  while ( i < ( DIFF - 1 ) )
74  {
75  /*
76  * We're doing 226 = 113*2, an even number of steps, so we can
77  * safely unroll one more step here for speed:
78  */
79  UNROLL( i + PERIOD );
80  UNROLL( i + PERIOD );
81  }
82 
83  // i = 226
84  UNROLL( ( i + PERIOD ) % SIZE );
85 
86  // i = [227 ... 622]
87  while ( i < ( SIZE - 1 ) )
88  {
89  /*
90  * 623-227 = 396 = 2*2*3*3*11, so we can unroll this loop in any number
91  * that evenly divides 396 (2, 4, 6, etc). Here we'll unroll 11 times.
92  */
93  UNROLL( i - DIFF );
94  UNROLL( i - DIFF );
95  UNROLL( i - DIFF );
96  UNROLL( i - DIFF );
97  UNROLL( i - DIFF );
98  UNROLL( i - DIFF );
99  UNROLL( i - DIFF );
100  UNROLL( i - DIFF );
101  UNROLL( i - DIFF );
102  UNROLL( i - DIFF );
103  UNROLL( i - DIFF );
104  }
105 
106  // i = 623
107  y = M32( MT[SIZE - 1] ) | L31( MT[0] );
108  MT[SIZE - 1] = MT[PERIOD - 1] ^ ( y >> 1 ) ^ MATRIX[ODD( y )];
109 }
110 
111 extern "C" void seed( uint32_t value )
112 {
113  /*
114  * The equation below is a linear congruential generator (LCG),
115  * one of the oldest known pseudo-random number generator
116  * algorithms, in the form X_(n+1) = = (a*X_n + c) (mod m).
117  *
118  * We've implicitly got m=32 (mask + word size of 32 bits), so
119  * there is no need to explicitly use modulus.
120  *
121  * What is interesting is the multiplier a. The one we have
122  * below is 0x6c07865 --- 1812433253 in decimal, and is called
123  * the Borosh-Niederreiter multiplier for modulus 2^32.
124  *
125  * It is mentioned in passing in Knuth's THE ART OF COMPUTER
126  * PROGRAMMING, Volume 2, page 106, Table 1, line 13. LCGs are
127  * treated in the same book, pp. 10-26
128  *
129  * You can read the original paper by Borosh and Niederreiter
130  * as well. It's called OPTIMAL MULTIPLIERS FOR PSEUDO-RANDOM
131  * NUMBER GENERATION BY THE LINEAR CONGRUENTIAL METHOD (1983) at
132  * http://www.springerlink.com/content/n7765ku70w8857l7/
133  *
134  * You can read about LCGs at:
135  * http://en.wikipedia.org/wiki/Linear_congruential_generator
136  *
137  * From that page, it says:
138  * "A common Mersenne twister implementation, interestingly
139  * enough, uses an LCG to generate seed data.",
140  *
141  * Since we're using 32-bits data types for our MT array, we can skip the
142  * masking with 0xFFFFFFFF below.
143  */
144 
145  MT[0] = value;
146  index = 0;
147 
148  for ( unsigned i = 1; i < SIZE; ++i )
149  MT[i] = 0x6c078965 * ( MT[i - 1] ^ MT[i - 1] >> 30 ) + i;
150 }
151 
152 extern "C" uint32_t rand_u32()
153 {
154  if ( !index )
155  generate_numbers();
156 
157  uint32_t y = MT[index];
158 
159  // Tempering
160  y ^= y >> 11;
161  y ^= y << 7 & 0x9d2c5680;
162  y ^= y << 15 & 0xefc60000;
163  y ^= y >> 18;
164 
165  if ( ++index == SIZE )
166  index = 0;
167 
168  return y;
169 }
170 
171 extern "C" int mt_rand()
172 {
173  /*
174  * PORTABILITY WARNING:
175  *
176  * rand_u32() uses all 32-bits for the pseudo-random number,
177  * but rand() must return a number from 0 ... MD_RAND_MAX.
178  *
179  * We'll just assume that rand() only uses 31 bits worth of
180  * data, and that we're on a two's complement system.
181  *
182  * So, to output an integer compatible with rand(), we have
183  * two options: Either mask off the highest (32nd) bit, or
184  * shift right by one bit. Masking with 0x7FFFFFFF will be
185  * compatible with 64-bit MT[], so we'll just use that here.
186  *
187  */
188  return static_cast<int>( 0x7FFFFFFF & rand_u32() );
189 }
190 
191 extern "C" void mt_srand( unsigned value )
192 {
193  seed( static_cast<uint32_t>( value ) );
194 }
195 
196 extern "C" float randf_cc()
197 {
198  return static_cast<float>( rand_u32() ) / MD_UINT32_MAX;
199 }
200 
201 extern "C" float randf_co()
202 {
203  return static_cast<float>( rand_u32() ) / ( MD_UINT32_MAX + 1.0f );
204 }
205 
206 extern "C" float randf_oo()
207 {
208  return ( static_cast<float>( rand_u32() ) + 0.5f ) / ( MD_UINT32_MAX + 1.0f );
209 }
210 
211 extern "C" double randd_cc()
212 {
213  return static_cast<double>( rand_u32() ) / MD_UINT32_MAX;
214 }
215 
216 extern "C" double randd_co()
217 {
218  return static_cast<double>( rand_u32() ) / ( MD_UINT32_MAX + 1.0 );
219 }
220 
221 extern "C" double randd_oo()
222 {
223  return ( static_cast<double>( rand_u32() ) + 0.5 ) / ( MD_UINT32_MAX + 1.0 );
224 }
225 
226 extern "C" uint64_t rand_u64()
227 {
228  return static_cast<uint64_t>( rand_u32() ) << 32 | rand_u32();
229 }
double randd_cc()
#define MD_UINT32_MAX
uint64_t rand_u64()
#define UNROLL(expr)
float randf_cc()
int mt_rand()
float randf_co()
double randd_oo()
void seed(uint32_t value)
void mt_srand(unsigned value)
#define M32(x)
double randd_co()
#define L31(x)
#define ODD(x)
float randf_oo()
uint32_t rand_u32()