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 #ifndef MERSENNETWISTER_H
00058 #define MERSENNETWISTER_H
00059
00060
00061
00062
00063 #include <iostream>
00064 #include <limits.h>
00065 #include <stdio.h>
00066 #include <time.h>
00067 #include <math.h>
00068
00069 NAMESPACE_BEGIN
00070
00071 class MTRand {
00072
00073 public:
00074 typedef unsigned long uint32;
00075
00076 enum Dummy1 { N = 624 };
00077 enum Dummy2 { SAVE = N + 1 };
00078
00079
00080 protected:
00081 enum Dummy3 { M = 397 };
00082
00083 uint32 state[N];
00084 uint32 *pNext;
00085 int left;
00086
00087
00088
00089 public:
00090 MTRand( const uint32& oneSeed );
00091 MTRand( uint32 *const bigSeed, uint32 const seedLength = N );
00092 MTRand();
00093
00094
00095
00096
00097
00098
00099 double rand();
00100 double rand( const double& n );
00101 double randExc();
00102 double randExc( const double& n );
00103 double randDblExc();
00104 double randDblExc( const double& n );
00105 uint32 randInt();
00106 uint32 randInt( const uint32& n );
00107 double operator()() { return rand(); }
00108
00109
00110 double rand53();
00111
00112
00113 double randNorm( const double& mean = 0.0, const double& variance = 0.0 );
00114
00115
00116 void seed( const uint32 oneSeed );
00117 void seed( uint32 *const bigSeed, const uint32 seedLength = N );
00118 void seed();
00119
00120
00121 void save( uint32* saveArray ) const;
00122 void load( uint32 *const loadArray );
00123 friend std::ostream& operator<<( std::ostream& os, const MTRand& mtrand );
00124 friend std::istream& operator>>( std::istream& is, MTRand& mtrand );
00125
00126 protected:
00127 void initialize( const uint32 oneSeed );
00128 void reload();
00129 uint32 hiBit( const uint32& u ) const { return u & 0x80000000UL; }
00130 uint32 loBit( const uint32& u ) const { return u & 0x00000001UL; }
00131 uint32 loBits( const uint32& u ) const { return u & 0x7fffffffUL; }
00132 uint32 mixBits( const uint32& u, const uint32& v ) const
00133 { return hiBit(u) | loBits(v); }
00134 uint32 twist( const uint32& m, const uint32& s0, const uint32& s1 ) const
00135 { return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL); }
00136 static uint32 hash( time_t t, clock_t c );
00137 };
00138
00139
00140 inline MTRand::MTRand( const uint32& oneSeed )
00141 { seed(oneSeed); }
00142
00143 inline MTRand::MTRand( uint32 *const bigSeed, const uint32 seedLength )
00144 { seed(bigSeed,seedLength); }
00145
00146 inline MTRand::MTRand()
00147 { seed(); }
00148
00149 inline double MTRand::rand()
00150 { return double(randInt()) * (1.0/4294967295.0); }
00151
00152 inline double MTRand::rand( const double& n )
00153 { return rand() * n; }
00154
00155 inline double MTRand::randExc()
00156 { return double(randInt()) * (1.0/4294967296.0); }
00157
00158 inline double MTRand::randExc( const double& n )
00159 { return randExc() * n; }
00160
00161 inline double MTRand::randDblExc()
00162 { return ( double(randInt()) + 0.5 ) * (1.0/4294967296.0); }
00163
00164 inline double MTRand::randDblExc( const double& n )
00165 { return randDblExc() * n; }
00166
00167 inline double MTRand::rand53()
00168 {
00169 uint32 a = randInt() >> 5, b = randInt() >> 6;
00170 return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0);
00171 }
00172
00173 inline double MTRand::randNorm( const double& mean, const double& variance )
00174 {
00175
00176
00177 double r = sqrt( -2.0 * log( 1.0-randDblExc()) ) * variance;
00178 double phi = 2.0 * 3.14159265358979323846264338328 * randExc();
00179 return mean + r * cos(phi);
00180 }
00181
00182 inline MTRand::uint32 MTRand::randInt()
00183 {
00184
00185
00186
00187 if( left == 0 ) reload();
00188 --left;
00189
00190 uint32 s1;
00191 s1 = *pNext++;
00192 s1 ^= (s1 >> 11);
00193 s1 ^= (s1 << 7) & 0x9d2c5680UL;
00194 s1 ^= (s1 << 15) & 0xefc60000UL;
00195 return ( s1 ^ (s1 >> 18) );
00196 }
00197
00198 inline MTRand::uint32 MTRand::randInt( const uint32& n )
00199 {
00200
00201
00202 uint32 used = n;
00203 used |= used >> 1;
00204 used |= used >> 2;
00205 used |= used >> 4;
00206 used |= used >> 8;
00207 used |= used >> 16;
00208
00209
00210 uint32 i;
00211 do
00212 i = randInt() & used;
00213 while( i > n );
00214 return i;
00215 }
00216
00217
00218 inline void MTRand::seed( const uint32 oneSeed )
00219 {
00220
00221 initialize(oneSeed);
00222 reload();
00223 }
00224
00225
00226 inline void MTRand::seed( uint32 *const bigSeed, const uint32 seedLength )
00227 {
00228
00229
00230
00231
00232
00233
00234 initialize(19650218UL);
00235 int i = 1;
00236 uint32 j = 0;
00237 int k = ( N > seedLength ? N : seedLength );
00238 for( ; k; --k )
00239 {
00240 state[i] =
00241 state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL );
00242 state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
00243 state[i] &= 0xffffffffUL;
00244 ++i; ++j;
00245 if( i >= N ) { state[0] = state[N-1]; i = 1; }
00246 if( j >= seedLength ) j = 0;
00247 }
00248 for( k = N - 1; k; --k )
00249 {
00250 state[i] =
00251 state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL );
00252 state[i] -= i;
00253 state[i] &= 0xffffffffUL;
00254 ++i;
00255 if( i >= N ) { state[0] = state[N-1]; i = 1; }
00256 }
00257 state[0] = 0x80000000UL;
00258 reload();
00259 }
00260
00261
00262 inline void MTRand::seed()
00263 {
00264
00265
00266
00267
00268 FILE* urandom = fopen( "/dev/urandom", "rb" );
00269 if( urandom )
00270 {
00271 uint32 bigSeed[N];
00272 uint32 *s = bigSeed;
00273 int i = N;
00274 bool success = true;
00275 while( success && i-- )
00276 success = fread( s++, sizeof(uint32), 1, urandom );
00277 fclose(urandom);
00278 if( success ) { seed( bigSeed, N ); return; }
00279 }
00280
00281
00282 seed( hash( time(NULL), clock() ) );
00283 }
00284
00285
00286 inline void MTRand::initialize( const uint32 seed )
00287 {
00288
00289
00290
00291
00292 uint32 *s = state;
00293 uint32 *r = state;
00294 int i = 1;
00295 *s++ = seed & 0xffffffffUL;
00296 for( ; i < N; ++i )
00297 {
00298 *s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL;
00299 r++;
00300 }
00301 }
00302
00303
00304 inline void MTRand::reload()
00305 {
00306
00307
00308 uint32 *p = state;
00309 int i;
00310 for( i = N - M; i--; ++p )
00311 *p = twist( p[M], p[0], p[1] );
00312 for( i = M; --i; ++p )
00313 *p = twist( p[M-N], p[0], p[1] );
00314 *p = twist( p[M-N], p[0], state[0] );
00315
00316 left = N, pNext = state;
00317 }
00318
00319
00320 inline MTRand::uint32 MTRand::hash( time_t t, clock_t c )
00321 {
00322
00323
00324
00325
00326 static uint32 differ = 0;
00327
00328 uint32 h1 = 0;
00329 unsigned char *p = (unsigned char *) &t;
00330 for( size_t i = 0; i < sizeof(t); ++i )
00331 {
00332 h1 *= UCHAR_MAX + 2U;
00333 h1 += p[i];
00334 }
00335 uint32 h2 = 0;
00336 p = (unsigned char *) &c;
00337 for( size_t j = 0; j < sizeof(c); ++j )
00338 {
00339 h2 *= UCHAR_MAX + 2U;
00340 h2 += p[j];
00341 }
00342 return ( h1 + differ++ ) ^ h2;
00343 }
00344
00345
00346 inline void MTRand::save( uint32* saveArray ) const
00347 {
00348 uint32 *sa = saveArray;
00349 const uint32 *s = state;
00350 int i = N;
00351 for( ; i--; *sa++ = *s++ ) {}
00352 *sa = left;
00353 }
00354
00355
00356 inline void MTRand::load( uint32 *const loadArray )
00357 {
00358 uint32 *s = state;
00359 uint32 *la = loadArray;
00360 int i = N;
00361 for( ; i--; *s++ = *la++ ) {}
00362 left = *la;
00363 pNext = &state[N-left];
00364 }
00365
00366
00367 inline std::ostream& operator<<( std::ostream& os, const MTRand& mtrand )
00368 {
00369 const MTRand::uint32 *s = mtrand.state;
00370 int i = mtrand.N;
00371 for( ; i--; os << *s++ << "\t" ) {}
00372 return os << mtrand.left;
00373 }
00374
00375
00376 inline std::istream& operator>>( std::istream& is, MTRand& mtrand )
00377 {
00378 MTRand::uint32 *s = mtrand.state;
00379 int i = mtrand.N;
00380 for( ; i--; is >> *s++ ) {}
00381 is >> mtrand.left;
00382 mtrand.pNext = &mtrand.state[mtrand.N-mtrand.left];
00383 return is;
00384 }
00385
00386 NAMESPACE_END
00387
00388
00389 #endif // MERSENNETWISTER_H
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429