Mercurial > hg > pymt
comparison MersenneTwister.h @ 0:0fa810263337
import
author | Peter Meerwald <pmeerw@cosy.sbg.ac.at> |
---|---|
date | Thu, 06 Sep 2007 13:58:46 +0200 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:0fa810263337 |
---|---|
1 // MersenneTwister.h | |
2 // Mersenne Twister random number generator -- a C++ class MTRand | |
3 // Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus | |
4 // Richard J. Wagner v1.0 15 May 2003 rjwagner@writeme.com | |
5 | |
6 // The Mersenne Twister is an algorithm for generating random numbers. It | |
7 // was designed with consideration of the flaws in various other generators. | |
8 // The period, 2^19937-1, and the order of equidistribution, 623 dimensions, | |
9 // are far greater. The generator is also fast; it avoids multiplication and | |
10 // division, and it benefits from caches and pipelines. For more information | |
11 // see the inventors' web page at http://www.math.keio.ac.jp/~matumoto/emt.html | |
12 | |
13 // Reference | |
14 // M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-Dimensionally | |
15 // Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions on | |
16 // Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3-30. | |
17 | |
18 // Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, | |
19 // Copyright (C) 2000 - 2003, Richard J. Wagner | |
20 // All rights reserved. | |
21 // | |
22 // Redistribution and use in source and binary forms, with or without | |
23 // modification, are permitted provided that the following conditions | |
24 // are met: | |
25 // | |
26 // 1. Redistributions of source code must retain the above copyright | |
27 // notice, this list of conditions and the following disclaimer. | |
28 // | |
29 // 2. Redistributions in binary form must reproduce the above copyright | |
30 // notice, this list of conditions and the following disclaimer in the | |
31 // documentation and/or other materials provided with the distribution. | |
32 // | |
33 // 3. The names of its contributors may not be used to endorse or promote | |
34 // products derived from this software without specific prior written | |
35 // permission. | |
36 // | |
37 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS | |
38 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT | |
39 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR | |
40 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR | |
41 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, | |
42 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, | |
43 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR | |
44 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF | |
45 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING | |
46 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS | |
47 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. | |
48 | |
49 // The original code included the following notice: | |
50 // | |
51 // When you use this, send an email to: matumoto@math.keio.ac.jp | |
52 // with an appropriate reference to your work. | |
53 // | |
54 // It would be nice to CC: rjwagner@writeme.com and Cokus@math.washington.edu | |
55 // when you write. | |
56 | |
57 #ifndef MERSENNETWISTER_H | |
58 #define MERSENNETWISTER_H | |
59 | |
60 // Not thread safe (unless auto-initialization is avoided and each thread has | |
61 // its own MTRand object) | |
62 | |
63 #include <iostream> | |
64 #include <limits.h> | |
65 #include <stdio.h> | |
66 #include <time.h> | |
67 #include <math.h> | |
68 | |
69 class MTRand { | |
70 // Data | |
71 public: | |
72 typedef unsigned long uint32; // unsigned integer type, at least 32 bits | |
73 | |
74 enum { N = 624 }; // length of state vector | |
75 enum { SAVE = N + 1 }; // length of array for save() | |
76 | |
77 protected: | |
78 enum { M = 397 }; // period parameter | |
79 | |
80 uint32 state[N]; // internal state | |
81 uint32 *pNext; // next value to get from state | |
82 int left; // number of values left before reload needed | |
83 | |
84 | |
85 //Methods | |
86 public: | |
87 MTRand( const uint32& oneSeed ); // initialize with a simple uint32 | |
88 MTRand( uint32 *const bigSeed, uint32 const seedLength = N ); // or an array | |
89 MTRand(); // auto-initialize with /dev/urandom or time() and clock() | |
90 | |
91 // Do NOT use for CRYPTOGRAPHY without securely hashing several returned | |
92 // values together, otherwise the generator state can be learned after | |
93 // reading 624 consecutive values. | |
94 | |
95 // Access to 32-bit random numbers | |
96 double rand(); // real number in [0,1] | |
97 double rand( const double& n ); // real number in [0,n] | |
98 double randExc(); // real number in [0,1) | |
99 double randExc( const double& n ); // real number in [0,n) | |
100 double randDblExc(); // real number in (0,1) | |
101 double randDblExc( const double& n ); // real number in (0,n) | |
102 uint32 randInt(); // integer in [0,2^32-1] | |
103 uint32 randInt( const uint32& n ); // integer in [0,n] for n < 2^32 | |
104 double operator()() { return rand(); } // same as rand() | |
105 | |
106 // Access to 53-bit random numbers (capacity of IEEE double precision) | |
107 double rand53(); // real number in [0,1) | |
108 | |
109 // Access to nonuniform random number distributions | |
110 double randNorm( const double& mean = 0.0, const double& variance = 0.0 ); | |
111 | |
112 // Re-seeding functions with same behavior as initializers | |
113 void seed( const uint32 oneSeed ); | |
114 void seed( uint32 *const bigSeed, const uint32 seedLength = N ); | |
115 void seed(); | |
116 | |
117 // Saving and loading generator state | |
118 void save( uint32* saveArray ) const; // to array of size SAVE | |
119 void load( uint32 *const loadArray ); // from such array | |
120 friend std::ostream& operator<<( std::ostream& os, const MTRand& mtrand ); | |
121 friend std::istream& operator>>( std::istream& is, MTRand& mtrand ); | |
122 | |
123 protected: | |
124 void initialize( const uint32 oneSeed ); | |
125 void reload(); | |
126 uint32 hiBit( const uint32& u ) const { return u & 0x80000000UL; } | |
127 uint32 loBit( const uint32& u ) const { return u & 0x00000001UL; } | |
128 uint32 loBits( const uint32& u ) const { return u & 0x7fffffffUL; } | |
129 uint32 mixBits( const uint32& u, const uint32& v ) const | |
130 { return hiBit(u) | loBits(v); } | |
131 uint32 twist( const uint32& m, const uint32& s0, const uint32& s1 ) const | |
132 { return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL); } | |
133 static uint32 hash( time_t t, clock_t c ); | |
134 }; | |
135 | |
136 | |
137 inline MTRand::MTRand( const uint32& oneSeed ) | |
138 { seed(oneSeed); } | |
139 | |
140 inline MTRand::MTRand( uint32 *const bigSeed, const uint32 seedLength ) | |
141 { seed(bigSeed,seedLength); } | |
142 | |
143 inline MTRand::MTRand() | |
144 { seed(); } | |
145 | |
146 inline double MTRand::rand() | |
147 { return double(randInt()) * (1.0/4294967295.0); } | |
148 | |
149 inline double MTRand::rand( const double& n ) | |
150 { return rand() * n; } | |
151 | |
152 inline double MTRand::randExc() | |
153 { return double(randInt()) * (1.0/4294967296.0); } | |
154 | |
155 inline double MTRand::randExc( const double& n ) | |
156 { return randExc() * n; } | |
157 | |
158 inline double MTRand::randDblExc() | |
159 { return ( double(randInt()) + 0.5 ) * (1.0/4294967296.0); } | |
160 | |
161 inline double MTRand::randDblExc( const double& n ) | |
162 { return randDblExc() * n; } | |
163 | |
164 inline double MTRand::rand53() | |
165 { | |
166 uint32 a = randInt() >> 5, b = randInt() >> 6; | |
167 return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0); // by Isaku Wada | |
168 } | |
169 | |
170 inline double MTRand::randNorm( const double& mean, const double& variance ) | |
171 { | |
172 // Return a real number from a normal (Gaussian) distribution with given | |
173 // mean and variance by Box-Muller method | |
174 double r = sqrt( -2.0 * log( 1.0-randDblExc()) ) * variance; | |
175 double phi = 2.0 * 3.14159265358979323846264338328 * randExc(); | |
176 return mean + r * cos(phi); | |
177 } | |
178 | |
179 inline MTRand::uint32 MTRand::randInt() | |
180 { | |
181 // Pull a 32-bit integer from the generator state | |
182 // Every other access function simply transforms the numbers extracted here | |
183 | |
184 if( left == 0 ) reload(); | |
185 --left; | |
186 | |
187 register uint32 s1; | |
188 s1 = *pNext++; | |
189 s1 ^= (s1 >> 11); | |
190 s1 ^= (s1 << 7) & 0x9d2c5680UL; | |
191 s1 ^= (s1 << 15) & 0xefc60000UL; | |
192 return ( s1 ^ (s1 >> 18) ); | |
193 } | |
194 | |
195 inline MTRand::uint32 MTRand::randInt( const uint32& n ) | |
196 { | |
197 // Find which bits are used in n | |
198 // Optimized by Magnus Jonsson (magnus@smartelectronix.com) | |
199 uint32 used = n; | |
200 used |= used >> 1; | |
201 used |= used >> 2; | |
202 used |= used >> 4; | |
203 used |= used >> 8; | |
204 used |= used >> 16; | |
205 | |
206 // Draw numbers until one is found in [0,n] | |
207 uint32 i; | |
208 do | |
209 i = randInt() & used; // toss unused bits to shorten search | |
210 while( i > n ); | |
211 return i; | |
212 } | |
213 | |
214 | |
215 inline void MTRand::seed( const uint32 oneSeed ) | |
216 { | |
217 // Seed the generator with a simple uint32 | |
218 initialize(oneSeed); | |
219 reload(); | |
220 } | |
221 | |
222 | |
223 inline void MTRand::seed( uint32 *const bigSeed, const uint32 seedLength ) | |
224 { | |
225 // Seed the generator with an array of uint32's | |
226 // There are 2^19937-1 possible initial states. This function allows | |
227 // all of those to be accessed by providing at least 19937 bits (with a | |
228 // default seed length of N = 624 uint32's). Any bits above the lower 32 | |
229 // in each element are discarded. | |
230 // Just call seed() if you want to get array from /dev/urandom | |
231 initialize(19650218UL); | |
232 register int i = 1; | |
233 register uint32 j = 0; | |
234 register int k = ( N > seedLength ? N : seedLength ); | |
235 for( ; k; --k ) | |
236 { | |
237 state[i] = | |
238 state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL ); | |
239 state[i] += ( bigSeed[j] & 0xffffffffUL ) + j; | |
240 state[i] &= 0xffffffffUL; | |
241 ++i; ++j; | |
242 if( i >= N ) { state[0] = state[N-1]; i = 1; } | |
243 if( j >= seedLength ) j = 0; | |
244 } | |
245 for( k = N - 1; k; --k ) | |
246 { | |
247 state[i] = | |
248 state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL ); | |
249 state[i] -= i; | |
250 state[i] &= 0xffffffffUL; | |
251 ++i; | |
252 if( i >= N ) { state[0] = state[N-1]; i = 1; } | |
253 } | |
254 state[0] = 0x80000000UL; // MSB is 1, assuring non-zero initial array | |
255 reload(); | |
256 } | |
257 | |
258 | |
259 inline void MTRand::seed() | |
260 { | |
261 // Seed the generator with an array from /dev/urandom if available | |
262 // Otherwise use a hash of time() and clock() values | |
263 | |
264 // First try getting an array from /dev/urandom | |
265 FILE* urandom = fopen( "/dev/urandom", "rb" ); | |
266 if( urandom ) | |
267 { | |
268 uint32 bigSeed[N]; | |
269 register uint32 *s = bigSeed; | |
270 register int i = N; | |
271 register bool success = true; | |
272 while( success && i-- ) | |
273 success = fread( s++, sizeof(uint32), 1, urandom ); | |
274 fclose(urandom); | |
275 if( success ) { seed( bigSeed, N ); return; } | |
276 } | |
277 | |
278 // Was not successful, so use time() and clock() instead | |
279 seed( hash( time(NULL), clock() ) ); | |
280 } | |
281 | |
282 | |
283 inline void MTRand::initialize( const uint32 seed ) | |
284 { | |
285 // Initialize generator state with seed | |
286 // See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier. | |
287 // In previous versions, most significant bits (MSBs) of the seed affect | |
288 // only MSBs of the state array. Modified 9 Jan 2002 by Makoto Matsumoto. | |
289 register uint32 *s = state; | |
290 register uint32 *r = state; | |
291 register int i = 1; | |
292 *s++ = seed & 0xffffffffUL; | |
293 for( ; i < N; ++i ) | |
294 { | |
295 *s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL; | |
296 r++; | |
297 } | |
298 } | |
299 | |
300 | |
301 inline void MTRand::reload() | |
302 { | |
303 // Generate N new values in state | |
304 // Made clearer and faster by Matthew Bellew (matthew.bellew@home.com) | |
305 register uint32 *p = state; | |
306 register int i; | |
307 for( i = N - M; i--; ++p ) | |
308 *p = twist( p[M], p[0], p[1] ); | |
309 for( i = M; --i; ++p ) | |
310 *p = twist( p[M-N], p[0], p[1] ); | |
311 *p = twist( p[M-N], p[0], state[0] ); | |
312 | |
313 left = N, pNext = state; | |
314 } | |
315 | |
316 | |
317 inline MTRand::uint32 MTRand::hash( time_t t, clock_t c ) | |
318 { | |
319 // Get a uint32 from t and c | |
320 // Better than uint32(x) in case x is floating point in [0,1] | |
321 // Based on code by Lawrence Kirby (fred@genesis.demon.co.uk) | |
322 | |
323 static uint32 differ = 0; // guarantee time-based seeds will change | |
324 | |
325 uint32 h1 = 0; | |
326 unsigned char *p = (unsigned char *) &t; | |
327 for( size_t i = 0; i < sizeof(t); ++i ) | |
328 { | |
329 h1 *= UCHAR_MAX + 2U; | |
330 h1 += p[i]; | |
331 } | |
332 uint32 h2 = 0; | |
333 p = (unsigned char *) &c; | |
334 for( size_t j = 0; j < sizeof(c); ++j ) | |
335 { | |
336 h2 *= UCHAR_MAX + 2U; | |
337 h2 += p[j]; | |
338 } | |
339 return ( h1 + differ++ ) ^ h2; | |
340 } | |
341 | |
342 | |
343 inline void MTRand::save( uint32* saveArray ) const | |
344 { | |
345 register uint32 *sa = saveArray; | |
346 register const uint32 *s = state; | |
347 register int i = N; | |
348 for( ; i--; *sa++ = *s++ ) {} | |
349 *sa = left; | |
350 } | |
351 | |
352 | |
353 inline void MTRand::load( uint32 *const loadArray ) | |
354 { | |
355 register uint32 *s = state; | |
356 register uint32 *la = loadArray; | |
357 register int i = N; | |
358 for( ; i--; *s++ = *la++ ) {} | |
359 left = *la; | |
360 pNext = &state[N-left]; | |
361 } | |
362 | |
363 | |
364 inline std::ostream& operator<<( std::ostream& os, const MTRand& mtrand ) | |
365 { | |
366 register const MTRand::uint32 *s = mtrand.state; | |
367 register int i = mtrand.N; | |
368 for( ; i--; os << *s++ << "\t" ) {} | |
369 return os << mtrand.left; | |
370 } | |
371 | |
372 | |
373 inline std::istream& operator>>( std::istream& is, MTRand& mtrand ) | |
374 { | |
375 register MTRand::uint32 *s = mtrand.state; | |
376 register int i = mtrand.N; | |
377 for( ; i--; is >> *s++ ) {} | |
378 is >> mtrand.left; | |
379 mtrand.pNext = &mtrand.state[mtrand.N-mtrand.left]; | |
380 return is; | |
381 } | |
382 | |
383 #endif // MERSENNETWISTER_H | |
384 | |
385 // Change log: | |
386 // | |
387 // v0.1 - First release on 15 May 2000 | |
388 // - Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus | |
389 // - Translated from C to C++ | |
390 // - Made completely ANSI compliant | |
391 // - Designed convenient interface for initialization, seeding, and | |
392 // obtaining numbers in default or user-defined ranges | |
393 // - Added automatic seeding from /dev/urandom or time() and clock() | |
394 // - Provided functions for saving and loading generator state | |
395 // | |
396 // v0.2 - Fixed bug which reloaded generator one step too late | |
397 // | |
398 // v0.3 - Switched to clearer, faster reload() code from Matthew Bellew | |
399 // | |
400 // v0.4 - Removed trailing newline in saved generator format to be consistent | |
401 // with output format of built-in types | |
402 // | |
403 // v0.5 - Improved portability by replacing static const int's with enum's and | |
404 // clarifying return values in seed(); suggested by Eric Heimburg | |
405 // - Removed MAXINT constant; use 0xffffffffUL instead | |
406 // | |
407 // v0.6 - Eliminated seed overflow when uint32 is larger than 32 bits | |
408 // - Changed integer [0,n] generator to give better uniformity | |
409 // | |
410 // v0.7 - Fixed operator precedence ambiguity in reload() | |
411 // - Added access for real numbers in (0,1) and (0,n) | |
412 // | |
413 // v0.8 - Included time.h header to properly support time_t and clock_t | |
414 // | |
415 // v1.0 - Revised seeding to match 26 Jan 2002 update of Nishimura and Matsumoto | |
416 // - Allowed for seeding with arrays of any length | |
417 // - Added access for real numbers in [0,1) with 53-bit resolution | |
418 // - Added access for real numbers from normal (Gaussian) distributions | |
419 // - Increased overall speed by optimizing twist() | |
420 // - Doubled speed of integer [0,n] generation | |
421 // - Fixed out-of-range number generation on 64-bit machines | |
422 // - Improved portability by substituting literal constants for long enum's | |
423 // - Changed license from GNU LGPL to BSD |