Essentia  2.1-beta5-dev
MersenneTwister.h
Go to the documentation of this file.
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 
70 class MTRand {
71 // Data
72 public:
73  typedef unsigned long uint32; // unsigned integer type, at least 32 bits
74 
75  enum { N = 624 }; // length of state vector
76  enum { SAVE = N + 1 }; // length of array for save()
77 
78 protected:
79  enum { M = 397 }; // period parameter
80 
81  uint32 state[N]; // internal state
82  uint32 *pNext; // next value to get from state
83  int left; // number of values left before reload needed
84 
85 
86 //Methods
87 public:
88  MTRand( const uint32& oneSeed ); // initialize with a simple uint32
89  MTRand( uint32 *const bigSeed, uint32 const seedLength = N ); // or an array
90  MTRand(); // auto-initialize with /dev/urandom or time() and clock()
91 
92  // Do NOT use for CRYPTOGRAPHY without securely hashing several returned
93  // values together, otherwise the generator state can be learned after
94  // reading 624 consecutive values.
95 
96  // Access to 32-bit random numbers
97  double rand(); // real number in [0,1]
98  double rand( const double& n ); // real number in [0,n]
99  double randExc(); // real number in [0,1)
100  double randExc( const double& n ); // real number in [0,n)
101  double randDblExc(); // real number in (0,1)
102  double randDblExc( const double& n ); // real number in (0,n)
103  uint32 randInt(); // integer in [0,2^32-1]
104  uint32 randInt( const uint32& n ); // integer in [0,n] for n < 2^32
105  double operator()() { return rand(); } // same as rand()
106 
107  // Access to 53-bit random numbers (capacity of IEEE double precision)
108  double rand53(); // real number in [0,1)
109 
110  // Access to nonuniform random number distributions
111  double randNorm( const double& mean = 0.0, const double& variance = 0.0 );
112 
113  // Re-seeding functions with same behavior as initializers
114  void seed( const uint32 oneSeed );
115  void seed( uint32 *const bigSeed, const uint32 seedLength = N );
116  void seed();
117 
118  // Saving and loading generator state
119  void save( uint32* saveArray ) const; // to array of size SAVE
120  void load( uint32 *const loadArray ); // from such array
121  friend std::ostream& operator<<( std::ostream& os, const MTRand& mtrand );
122  friend std::istream& operator>>( std::istream& is, MTRand& mtrand );
123 
124 protected:
125  void initialize( const uint32 oneSeed );
126  void reload();
127  uint32 hiBit( const uint32& u ) const { return u & 0x80000000UL; }
128  uint32 loBit( const uint32& u ) const { return u & 0x00000001UL; }
129  uint32 loBits( const uint32& u ) const { return u & 0x7fffffffUL; }
130  uint32 mixBits( const uint32& u, const uint32& v ) const
131  { return hiBit(u) | loBits(v); }
132  uint32 twist( const uint32& m, const uint32& s0, const uint32& s1 ) const
133  { return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL); }
134  static uint32 hash( time_t t, clock_t c );
135 };
136 
137 
138 inline MTRand::MTRand( const uint32& oneSeed )
139  { seed(oneSeed); }
140 
141 inline MTRand::MTRand( uint32 *const bigSeed, const uint32 seedLength )
142  { seed(bigSeed,seedLength); }
143 
145  { seed(); }
146 
147 inline double MTRand::rand()
148  { return double(randInt()) * (1.0/4294967295.0); }
149 
150 inline double MTRand::rand( const double& n )
151  { return rand() * n; }
152 
153 inline double MTRand::randExc()
154  { return double(randInt()) * (1.0/4294967296.0); }
155 
156 inline double MTRand::randExc( const double& n )
157  { return randExc() * n; }
158 
159 inline double MTRand::randDblExc()
160  { return ( double(randInt()) + 0.5 ) * (1.0/4294967296.0); }
161 
162 inline double MTRand::randDblExc( const double& n )
163  { return randDblExc() * n; }
164 
165 inline double MTRand::rand53()
166 {
167  uint32 a = randInt() >> 5, b = randInt() >> 6;
168  return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0); // by Isaku Wada
169 }
170 
171 inline double MTRand::randNorm( const double& mean, const double& variance )
172 {
173  // Return a real number from a normal (Gaussian) distribution with given
174  // mean and variance by Box-Muller method
175  double r = sqrt( -2.0 * log( 1.0-randDblExc()) ) * variance;
176  double phi = 2.0 * 3.14159265358979323846264338328 * randExc();
177  return mean + r * cos(phi);
178 }
179 
181 {
182  // Pull a 32-bit integer from the generator state
183  // Every other access function simply transforms the numbers extracted here
184 
185  if( left == 0 ) reload();
186  --left;
187 
188  register uint32 s1;
189  s1 = *pNext++;
190  s1 ^= (s1 >> 11);
191  s1 ^= (s1 << 7) & 0x9d2c5680UL;
192  s1 ^= (s1 << 15) & 0xefc60000UL;
193  return ( s1 ^ (s1 >> 18) );
194 }
195 
197 {
198  // Find which bits are used in n
199  // Optimized by Magnus Jonsson (magnus@smartelectronix.com)
200  uint32 used = n;
201  used |= used >> 1;
202  used |= used >> 2;
203  used |= used >> 4;
204  used |= used >> 8;
205  used |= used >> 16;
206 
207  // Draw numbers until one is found in [0,n]
208  uint32 i;
209  do
210  i = randInt() & used; // toss unused bits to shorten search
211  while( i > n );
212  return i;
213 }
214 
215 
216 inline void MTRand::seed( const uint32 oneSeed )
217 {
218  // Seed the generator with a simple uint32
219  initialize(oneSeed);
220  reload();
221 }
222 
223 
224 inline void MTRand::seed( uint32 *const bigSeed, const uint32 seedLength )
225 {
226  // Seed the generator with an array of uint32's
227  // There are 2^19937-1 possible initial states. This function allows
228  // all of those to be accessed by providing at least 19937 bits (with a
229  // default seed length of N = 624 uint32's). Any bits above the lower 32
230  // in each element are discarded.
231  // Just call seed() if you want to get array from /dev/urandom
232  initialize(19650218UL);
233  register int i = 1;
234  register uint32 j = 0;
235  register int k = ( N > seedLength ? N : seedLength );
236  for( ; k; --k )
237  {
238  state[i] =
239  state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL );
240  state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
241  state[i] &= 0xffffffffUL;
242  ++i; ++j;
243  if( i >= N ) { state[0] = state[N-1]; i = 1; }
244  if( j >= seedLength ) j = 0;
245  }
246  for( k = N - 1; k; --k )
247  {
248  state[i] =
249  state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL );
250  state[i] -= i;
251  state[i] &= 0xffffffffUL;
252  ++i;
253  if( i >= N ) { state[0] = state[N-1]; i = 1; }
254  }
255  state[0] = 0x80000000UL; // MSB is 1, assuring non-zero initial array
256  reload();
257 }
258 
259 
260 inline void MTRand::seed()
261 {
262  // Seed the generator with an array from /dev/urandom if available
263  // Otherwise use a hash of time() and clock() values
264 
265 #ifdef OS_LINUX
266  // First try getting an array from /dev/urandom
267  FILE* urandom = fopen( "/dev/urandom", "rb" );
268  if( urandom )
269  {
270  uint32 bigSeed[N];
271  register uint32 *s = bigSeed;
272  register int i = N;
273  register bool success = true;
274  while( success && i-- )
275  success = (fread( s++, sizeof(uint32), 1, urandom ) > 0);
276  fclose(urandom);
277  if( success ) { seed( bigSeed, N ); return; }
278  }
279 #endif // OS_LINUX
280 
281  // Was not successful, so use time() and clock() instead
282  seed( hash( time(NULL), clock() ) );
283 }
284 
285 
286 inline void MTRand::initialize( const uint32 seed )
287 {
288  // Initialize generator state with seed
289  // See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
290  // In previous versions, most significant bits (MSBs) of the seed affect
291  // only MSBs of the state array. Modified 9 Jan 2002 by Makoto Matsumoto.
292  register uint32 *s = state;
293  register uint32 *r = state;
294  register int i = 1;
295  *s++ = seed & 0xffffffffUL;
296  for( ; i < N; ++i )
297  {
298  *s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL;
299  r++;
300  }
301 }
302 
303 
304 inline void MTRand::reload()
305 {
306  // Generate N new values in state
307  // Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
308  register uint32 *p = state;
309  register int i;
310  for( i = N - M; i--; ++p )
311  *p = twist( p[M], p[0], p[1] );
312  for( i = M; --i; ++p )
313  *p = twist( p[M-N], p[0], p[1] );
314  *p = twist( p[M-N], p[0], state[0] );
315 
316  left = N, pNext = state;
317 }
318 
319 
320 inline MTRand::uint32 MTRand::hash( time_t t, clock_t c )
321 {
322  // Get a uint32 from t and c
323  // Better than uint32(x) in case x is floating point in [0,1]
324  // Based on code by Lawrence Kirby (fred@genesis.demon.co.uk)
325 
326  static uint32 differ = 0; // guarantee time-based seeds will change
327 
328  uint32 h1 = 0;
329  unsigned char *p = (unsigned char *) &t;
330  for( size_t i = 0; i < sizeof(t); ++i )
331  {
332  h1 *= UCHAR_MAX + 2U;
333  h1 += p[i];
334  }
335  uint32 h2 = 0;
336  p = (unsigned char *) &c;
337  for( size_t j = 0; j < sizeof(c); ++j )
338  {
339  h2 *= UCHAR_MAX + 2U;
340  h2 += p[j];
341  }
342  return ( h1 + differ++ ) ^ h2;
343 }
344 
345 
346 inline void MTRand::save( uint32* saveArray ) const
347 {
348  register uint32 *sa = saveArray;
349  register const uint32 *s = state;
350  register int i = N;
351  for( ; i--; *sa++ = *s++ ) {}
352  *sa = left;
353 }
354 
355 
356 inline void MTRand::load( uint32 *const loadArray )
357 {
358  register uint32 *s = state;
359  register uint32 *la = loadArray;
360  register int i = N;
361  for( ; i--; *s++ = *la++ ) {}
362  left = *la;
363  pNext = &state[N-left];
364 }
365 
366 
367 inline std::ostream& operator<<( std::ostream& os, const MTRand& mtrand )
368 {
369  register const MTRand::uint32 *s = mtrand.state;
370  register int i = mtrand.N;
371  for( ; i--; os << *s++ << "\t" ) {}
372  return os << mtrand.left;
373 }
374 
375 
376 inline std::istream& operator>>( std::istream& is, MTRand& mtrand )
377 {
378  register MTRand::uint32 *s = mtrand.state;
379  register int i = mtrand.N;
380  for( ; i--; is >> *s++ ) {}
381  is >> mtrand.left;
382  mtrand.pNext = &mtrand.state[mtrand.N-mtrand.left];
383  return is;
384 }
385 
386 #endif // MERSENNETWISTER_H
387 
388 // Change log:
389 //
390 // v0.1 - First release on 15 May 2000
391 // - Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
392 // - Translated from C to C++
393 // - Made completely ANSI compliant
394 // - Designed convenient interface for initialization, seeding, and
395 // obtaining numbers in default or user-defined ranges
396 // - Added automatic seeding from /dev/urandom or time() and clock()
397 // - Provided functions for saving and loading generator state
398 //
399 // v0.2 - Fixed bug which reloaded generator one step too late
400 //
401 // v0.3 - Switched to clearer, faster reload() code from Matthew Bellew
402 //
403 // v0.4 - Removed trailing newline in saved generator format to be consistent
404 // with output format of built-in types
405 //
406 // v0.5 - Improved portability by replacing static const int's with enum's and
407 // clarifying return values in seed(); suggested by Eric Heimburg
408 // - Removed MAXINT constant; use 0xffffffffUL instead
409 //
410 // v0.6 - Eliminated seed overflow when uint32 is larger than 32 bits
411 // - Changed integer [0,n] generator to give better uniformity
412 //
413 // v0.7 - Fixed operator precedence ambiguity in reload()
414 // - Added access for real numbers in (0,1) and (0,n)
415 //
416 // v0.8 - Included time.h header to properly support time_t and clock_t
417 //
418 // v1.0 - Revised seeding to match 26 Jan 2002 update of Nishimura and Matsumoto
419 // - Allowed for seeding with arrays of any length
420 // - Added access for real numbers in [0,1) with 53-bit resolution
421 // - Added access for real numbers from normal (Gaussian) distributions
422 // - Increased overall speed by optimizing twist()
423 // - Doubled speed of integer [0,n] generation
424 // - Fixed out-of-range number generation on 64-bit machines
425 // - Improved portability by substituting literal constants for long enum's
426 // - Changed license from GNU LGPL to BSD
void initialize(const uint32 oneSeed)
Definition: MersenneTwister.h:286
T variance(const std::vector< T > &array, const T mean)
Definition: essentiamath.h:447
friend std::istream & operator>>(std::istream &is, MTRand &mtrand)
Definition: MersenneTwister.h:376
double rand()
Definition: MersenneTwister.h:147
static uint32 hash(time_t t, clock_t c)
Definition: MersenneTwister.h:320
uint32 loBits(const uint32 &u) const
Definition: MersenneTwister.h:129
unsigned long uint32
Definition: MersenneTwister.h:73
uint32 mixBits(const uint32 &u, const uint32 &v) const
Definition: MersenneTwister.h:130
#define NULL
Definition: tnt_i_refvec.h:33
int left
Definition: MersenneTwister.h:83
void save(uint32 *saveArray) const
Definition: MersenneTwister.h:346
double randDblExc()
Definition: MersenneTwister.h:159
void load(uint32 *const loadArray)
Definition: MersenneTwister.h:356
friend std::ostream & operator<<(std::ostream &os, const MTRand &mtrand)
Definition: MersenneTwister.h:367
uint32 * pNext
Definition: MersenneTwister.h:82
uint32 randInt()
Definition: MersenneTwister.h:180
MTRand()
Definition: MersenneTwister.h:144
Definition: MersenneTwister.h:75
double rand53()
Definition: MersenneTwister.h:165
double randExc()
Definition: MersenneTwister.h:153
Definition: MersenneTwister.h:79
Definition: MersenneTwister.h:70
void reload()
Definition: MersenneTwister.h:304
double operator()()
Definition: MersenneTwister.h:105
uint32 loBit(const uint32 &u) const
Definition: MersenneTwister.h:128
uint32 state[N]
Definition: MersenneTwister.h:81
T mean(const std::vector< T > &array, int start, int end)
Definition: essentiamath.h:142
void seed()
Definition: MersenneTwister.h:260
uint32 twist(const uint32 &m, const uint32 &s0, const uint32 &s1) const
Definition: MersenneTwister.h:132
uint32 hiBit(const uint32 &u) const
Definition: MersenneTwister.h:127
double randNorm(const double &mean=0.0, const double &variance=0.0)
Definition: MersenneTwister.h:171
Definition: MersenneTwister.h:76