MersenneTwister.h
00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060 #ifndef MERSENNETWISTER_WFMATH_H
00061 #define MERSENNETWISTER_WFMATH_H
00062
00063
00064
00065
00066 #include <iosfwd>
00067 #include <climits>
00068 #include <cmath>
00069
00070
00071
00072 namespace WFMath {
00073
00074 class MTRand {
00075
00076 public:
00077 typedef unsigned long uint32;
00078
00079 enum { N = 624 };
00080 enum { SAVE = N + 1 };
00081
00082 protected:
00083 enum { M = 397 };
00084
00085 uint32 state[N];
00086 uint32 *pNext;
00087 int left;
00088
00089
00090
00091 public:
00092 MTRand( const uint32& oneSeed );
00093 MTRand( uint32 *const bigSeed, uint32 const seedLength = N );
00094 MTRand();
00095
00096
00097
00098
00099
00100
00101 double rand();
00102 double rand( const double& n );
00103 double randExc();
00104 double randExc( const double& n );
00105 double randDblExc();
00106 double randDblExc( const double& n );
00107 uint32 randInt();
00108 uint32 randInt( const uint32& n );
00109 double operator()() { return rand(); }
00110
00111
00112 double rand53();
00113
00114
00115 double randNorm( const double& mean = 0.0, const double& variance = 0.0 );
00116
00117
00118 void seed( const uint32 oneSeed );
00119 void seed( uint32 *const bigSeed, const uint32 seedLength = N );
00120 void seed();
00121
00122
00123 void save( uint32* saveArray ) const;
00124 void load( uint32 *const loadArray );
00125 friend std::ostream& operator<<( std::ostream& os, const MTRand& mtrand );
00126 friend std::istream& operator>>( std::istream& is, MTRand& mtrand );
00127
00128 static MTRand instance;
00129
00130 protected:
00131 void initialize( const uint32 oneSeed );
00132 void reload();
00133 uint32 hiBit( const uint32& u ) const { return u & 0x80000000UL; }
00134 uint32 loBit( const uint32& u ) const { return u & 0x00000001UL; }
00135 uint32 loBits( const uint32& u ) const { return u & 0x7fffffffUL; }
00136 uint32 mixBits( const uint32& u, const uint32& v ) const
00137 { return hiBit(u) | loBits(v); }
00138 uint32 twist( const uint32& m, const uint32& s0, const uint32& s1 ) const
00139 { return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL); }
00140 };
00141
00142
00143 inline MTRand::MTRand( const uint32& oneSeed )
00144 { seed(oneSeed); }
00145
00146 inline MTRand::MTRand( uint32 *const bigSeed, const uint32 seedLength )
00147 { seed(bigSeed,seedLength); }
00148
00149 inline MTRand::MTRand()
00150 { seed(); }
00151
00152 inline double MTRand::rand()
00153 { return double(randInt()) * (1.0/4294967295.0); }
00154
00155 inline double MTRand::rand( const double& n )
00156 { return rand() * n; }
00157
00158 inline double MTRand::randExc()
00159 { return double(randInt()) * (1.0/4294967296.0); }
00160
00161 inline double MTRand::randExc( const double& n )
00162 { return randExc() * n; }
00163
00164 inline double MTRand::randDblExc()
00165 { return ( double(randInt()) + 0.5 ) * (1.0/4294967296.0); }
00166
00167 inline double MTRand::randDblExc( const double& n )
00168 { return randDblExc() * n; }
00169
00170 inline double MTRand::rand53()
00171 {
00172 uint32 a = randInt() >> 5, b = randInt() >> 6;
00173 return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0);
00174 }
00175
00176 inline double MTRand::randNorm( const double& mean, const double& variance )
00177 {
00178
00179
00180 double r = sqrt( -2.0 * log( 1.0-randDblExc()) ) * variance;
00181 double phi = 2.0 * 3.14159265358979323846264338328 * randExc();
00182 return mean + r * cos(phi);
00183 }
00184
00185 inline MTRand::uint32 MTRand::randInt()
00186 {
00187
00188
00189
00190 if( left == 0 ) reload();
00191 --left;
00192
00193 register uint32 s1;
00194 s1 = *pNext++;
00195 s1 ^= (s1 >> 11);
00196 s1 ^= (s1 << 7) & 0x9d2c5680UL;
00197 s1 ^= (s1 << 15) & 0xefc60000UL;
00198 return ( s1 ^ (s1 >> 18) );
00199 }
00200
00201 inline MTRand::uint32 MTRand::randInt( const uint32& n )
00202 {
00203
00204
00205 uint32 used = n;
00206 used |= used >> 1;
00207 used |= used >> 2;
00208 used |= used >> 4;
00209 used |= used >> 8;
00210 used |= used >> 16;
00211
00212
00213 uint32 i;
00214 do
00215 i = randInt() & used;
00216 while( i > n );
00217 return i;
00218 }
00219
00220
00221 inline void MTRand::seed( const uint32 oneSeed )
00222 {
00223
00224 initialize(oneSeed);
00225 reload();
00226 }
00227
00228
00229 inline void MTRand::seed( uint32 *const bigSeed, const uint32 seedLength )
00230 {
00231
00232
00233
00234
00235
00236
00237 initialize(19650218UL);
00238 register int i = 1;
00239 register uint32 j = 0;
00240 register int k = ( N > seedLength ? N : seedLength );
00241 for( ; k; --k )
00242 {
00243 state[i] =
00244 state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL );
00245 state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
00246 state[i] &= 0xffffffffUL;
00247 ++i; ++j;
00248 if( i >= N ) { state[0] = state[N-1]; i = 1; }
00249 if( j >= seedLength ) j = 0;
00250 }
00251 for( k = N - 1; k; --k )
00252 {
00253 state[i] =
00254 state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL );
00255 state[i] -= i;
00256 state[i] &= 0xffffffffUL;
00257 ++i;
00258 if( i >= N ) { state[0] = state[N-1]; i = 1; }
00259 }
00260 state[0] = 0x80000000UL;
00261 reload();
00262 }
00263
00264
00265 inline void MTRand::initialize( const uint32 seed )
00266 {
00267
00268
00269
00270
00271 register uint32 *s = state;
00272 register uint32 *r = state;
00273 register int i = 1;
00274 *s++ = seed & 0xffffffffUL;
00275 for( ; i < N; ++i )
00276 {
00277 *s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL;
00278 r++;
00279 }
00280 }
00281
00282
00283 inline void MTRand::reload()
00284 {
00285
00286
00287 register uint32 *p = state;
00288 register int i;
00289 for( i = N - M; i--; ++p )
00290 *p = twist( p[M], p[0], p[1] );
00291 for( i = M; --i; ++p )
00292 *p = twist( p[M-N], p[0], p[1] );
00293 *p = twist( p[M-N], p[0], state[0] );
00294
00295 left = N, pNext = state;
00296 }
00297
00298
00299
00300 inline void MTRand::save( uint32* saveArray ) const
00301 {
00302 register uint32 *sa = saveArray;
00303 register const uint32 *s = state;
00304 register int i = N;
00305 for( ; i--; *sa++ = *s++ ) {}
00306 *sa = left;
00307 }
00308
00309
00310 inline void MTRand::load( uint32 *const loadArray )
00311 {
00312 register uint32 *s = state;
00313 register uint32 *la = loadArray;
00314 register int i = N;
00315 for( ; i--; *s++ = *la++ ) {}
00316 left = *la;
00317 pNext = &state[N-left];
00318 }
00319
00320
00321 }
00322
00323 #endif // MERSENNETWISTER_H
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
Generated on Tue Jul 27 21:41:56 2004 for WFMath by
1.3.7