Statistics
| Branch: | Revision:

root / include / mersennetwister.h @ fbe00e73

History | View | Annotate | Download (14.6 KB)

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
NAMESPACE_BEGIN
70

    
71
class MTRand {
72
// Data
73
public:
74
    typedef unsigned long uint32;  // unsigned integer type, at least 32 bits
75

    
76
    enum Dummy1 { N = 624 };       // length of state vector
77
    enum Dummy2 { SAVE = N + 1 };  // length of array for save()
78
    // Note: DummyX names needed by buggy gcc 4.0.1 on OS/X (Andras)
79

    
80
protected:
81
    enum Dummy3 { M = 397 };  // period parameter
82

    
83
    uint32 state[N];   // internal state
84
    uint32 *pNext;     // next value to get from state
85
    int left;          // number of values left before reload needed
86

    
87

    
88
//Methods
89
public:
90
    MTRand( const uint32& oneSeed );  // initialize with a simple uint32
91
    MTRand( uint32 *const bigSeed, uint32 const seedLength = N );  // or an array
92
    MTRand();  // auto-initialize with /dev/urandom or time() and clock()
93

    
94
    // Do NOT use for CRYPTOGRAPHY without securely hashing several returned
95
    // values together, otherwise the generator state can be learned after
96
    // reading 624 consecutive values.
97

    
98
    // Access to 32-bit random numbers
99
    double rand();                          // real number in [0,1]
100
    double rand( const double& n );         // real number in [0,n]
101
    double randExc();                       // real number in [0,1)
102
    double randExc( const double& n );      // real number in [0,n)
103
    double randDblExc();                    // real number in (0,1)
104
    double randDblExc( const double& n );   // real number in (0,n)
105
    uint32 randInt();                       // integer in [0,2^32-1]
106
    uint32 randInt( const uint32& n );      // integer in [0,n] for n < 2^32
107
    double operator()() { return rand(); }  // same as rand()
108

    
109
    // Access to 53-bit random numbers (capacity of IEEE double precision)
110
    double rand53();  // real number in [0,1)
111

    
112
    // Access to nonuniform random number distributions
113
    double randNorm( const double& mean = 0.0, const double& variance = 0.0 );
114

    
115
    // Re-seeding functions with same behavior as initializers
116
    void seed( const uint32 oneSeed );
117
    void seed( uint32 *const bigSeed, const uint32 seedLength = N );
118
    void seed();
119

    
120
    // Saving and loading generator state
121
    void save( uint32* saveArray ) const;  // to array of size SAVE
122
    void load( uint32 *const loadArray );  // from such array
123
    friend std::ostream& operator<<( std::ostream& os, const MTRand& mtrand );
124
    friend std::istream& operator>>( std::istream& is, MTRand& mtrand );
125

    
126
protected:
127
    void initialize( const uint32 oneSeed );
128
    void reload();
129
    uint32 hiBit( const uint32& u ) const { return u & 0x80000000UL; }
130
    uint32 loBit( const uint32& u ) const { return u & 0x00000001UL; }
131
    uint32 loBits( const uint32& u ) const { return u & 0x7fffffffUL; }
132
    uint32 mixBits( const uint32& u, const uint32& v ) const
133
        { return hiBit(u) | loBits(v); }
134
    uint32 twist( const uint32& m, const uint32& s0, const uint32& s1 ) const
135
        { return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL); }
136
    static uint32 hash( time_t t, clock_t c );
137
};
138

    
139

    
140
inline MTRand::MTRand( const uint32& oneSeed )
141
    { seed(oneSeed); }
142

    
143
inline MTRand::MTRand( uint32 *const bigSeed, const uint32 seedLength )
144
    { seed(bigSeed,seedLength); }
145

    
146
inline MTRand::MTRand()
147
    { seed(); }
148

    
149
inline double MTRand::rand()
150
    { return double(randInt()) * (1.0/4294967295.0); }
151

    
152
inline double MTRand::rand( const double& n )
153
    { return rand() * n; }
154

    
155
inline double MTRand::randExc()
156
    { return double(randInt()) * (1.0/4294967296.0); }
157

    
158
inline double MTRand::randExc( const double& n )
159
    { return randExc() * n; }
160

    
161
inline double MTRand::randDblExc()
162
    { return ( double(randInt()) + 0.5 ) * (1.0/4294967296.0); }
163

    
164
inline double MTRand::randDblExc( const double& n )
165
    { return randDblExc() * n; }
166

    
167
inline double MTRand::rand53()
168
{
169
    uint32 a = randInt() >> 5, b = randInt() >> 6;
170
    return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0);  // by Isaku Wada
171
}
172

    
173
inline double MTRand::randNorm( const double& mean, const double& variance )
174
{
175
    // Return a real number from a normal (Gaussian) distribution with given
176
    // mean and variance by Box-Muller method
177
    double r = sqrt( -2.0 * log( 1.0-randDblExc()) ) * variance;
178
    double phi = 2.0 * 3.14159265358979323846264338328 * randExc();
179
    return mean + r * cos(phi);
180
}
181

    
182
inline MTRand::uint32 MTRand::randInt()
183
{
184
    // Pull a 32-bit integer from the generator state
185
    // Every other access function simply transforms the numbers extracted here
186

    
187
    if( left == 0 ) reload();
188
    --left;
189

    
190
    register uint32 s1;
191
    s1 = *pNext++;
192
    s1 ^= (s1 >> 11);
193
    s1 ^= (s1 <<  7) & 0x9d2c5680UL;
194
    s1 ^= (s1 << 15) & 0xefc60000UL;
195
    return ( s1 ^ (s1 >> 18) );
196
}
197

    
198
inline MTRand::uint32 MTRand::randInt( const uint32& n )
199
{
200
    // Find which bits are used in n
201
    // Optimized by Magnus Jonsson (magnus@smartelectronix.com)
202
    uint32 used = n;
203
    used |= used >> 1;
204
    used |= used >> 2;
205
    used |= used >> 4;
206
    used |= used >> 8;
207
    used |= used >> 16;
208

    
209
    // Draw numbers until one is found in [0,n]
210
    uint32 i;
211
    do
212
        i = randInt() & used;  // toss unused bits to shorten search
213
    while( i > n );
214
    return i;
215
}
216

    
217

    
218
inline void MTRand::seed( const uint32 oneSeed )
219
{
220
    // Seed the generator with a simple uint32
221
    initialize(oneSeed);
222
    reload();
223
}
224

    
225

    
226
inline void MTRand::seed( uint32 *const bigSeed, const uint32 seedLength )
227
{
228
    // Seed the generator with an array of uint32's
229
    // There are 2^19937-1 possible initial states.  This function allows
230
    // all of those to be accessed by providing at least 19937 bits (with a
231
    // default seed length of N = 624 uint32's).  Any bits above the lower 32
232
    // in each element are discarded.
233
    // Just call seed() if you want to get array from /dev/urandom
234
    initialize(19650218UL);
235
    register int i = 1;
236
    register uint32 j = 0;
237
    register int k = ( N > seedLength ? N : seedLength );
238
    for( ; k; --k )
239
    {
240
        state[i] =
241
            state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL );
242
        state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
243
        state[i] &= 0xffffffffUL;
244
        ++i;  ++j;
245
        if( i >= N ) { state[0] = state[N-1];  i = 1; }
246
        if( j >= seedLength ) j = 0;
247
    }
248
    for( k = N - 1; k; --k )
249
    {
250
        state[i] =
251
            state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL );
252
        state[i] -= i;
253
        state[i] &= 0xffffffffUL;
254
        ++i;
255
        if( i >= N ) { state[0] = state[N-1];  i = 1; }
256
    }
257
    state[0] = 0x80000000UL;  // MSB is 1, assuring non-zero initial array
258
    reload();
259
}
260

    
261

    
262
inline void MTRand::seed()
263
{
264
    // Seed the generator with an array from /dev/urandom if available
265
    // Otherwise use a hash of time() and clock() values
266

    
267
    // First try getting an array from /dev/urandom
268
    FILE* urandom = fopen( "/dev/urandom", "rb" );
269
    if( urandom )
270
    {
271
        uint32 bigSeed[N];
272
        register uint32 *s = bigSeed;
273
        register int i = N;
274
        register bool success = true;
275
        while( success && i-- )
276
            success = fread( s++, sizeof(uint32), 1, urandom );
277
        fclose(urandom);
278
        if( success ) { seed( bigSeed, N );  return; }
279
    }
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
NAMESPACE_END
387

    
388

    
389
#endif  // MERSENNETWISTER_H
390

    
391
// Change log:
392
//
393
// v0.1 - First release on 15 May 2000
394
//      - Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
395
//      - Translated from C to C++
396
//      - Made completely ANSI compliant
397
//      - Designed convenient interface for initialization, seeding, and
398
//        obtaining numbers in default or user-defined ranges
399
//      - Added automatic seeding from /dev/urandom or time() and clock()
400
//      - Provided functions for saving and loading generator state
401
//
402
// v0.2 - Fixed bug which reloaded generator one step too late
403
//
404
// v0.3 - Switched to clearer, faster reload() code from Matthew Bellew
405
//
406
// v0.4 - Removed trailing newline in saved generator format to be consistent
407
//        with output format of built-in types
408
//
409
// v0.5 - Improved portability by replacing static const int's with enum's and
410
//        clarifying return values in seed(); suggested by Eric Heimburg
411
//      - Removed MAXINT constant; use 0xffffffffUL instead
412
//
413
// v0.6 - Eliminated seed overflow when uint32 is larger than 32 bits
414
//      - Changed integer [0,n] generator to give better uniformity
415
//
416
// v0.7 - Fixed operator precedence ambiguity in reload()
417
//      - Added access for real numbers in (0,1) and (0,n)
418
//
419
// v0.8 - Included time.h header to properly support time_t and clock_t
420
//
421
// v1.0 - Revised seeding to match 26 Jan 2002 update of Nishimura and Matsumoto
422
//      - Allowed for seeding with arrays of any length
423
//      - Added access for real numbers in [0,1) with 53-bit resolution
424
//      - Added access for real numbers from normal (Gaussian) distributions
425
//      - Increased overall speed by optimizing twist()
426
//      - Doubled speed of integer [0,n] generation
427
//      - Fixed out-of-range number generation on 64-bit machines
428
//      - Improved portability by substituting literal constants for long enum's
429
//      - Changed license from GNU LGPL to BSD