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 #ifndef MERSENNETWISTER_H
00044 #define MERSENNETWISTER_H
00045
00046
00047
00048 #include <stdio.h>
00049 #include <time.h>
00050 #include <limits.h>
00051 #include <iostream.h>
00052
00053 class MTRand {
00054
00055 public:
00056 typedef unsigned long uint32;
00057
00058 enum { N = 624 };
00059 enum { SAVE = N + 1 };
00060
00061 protected:
00062 enum { M = 397 };
00063 enum { MAGIC = 0x9908b0dfU };
00064
00065 uint32 state[N];
00066 uint32 *pNext;
00067 int left;
00068
00069
00070
00071 public:
00072 MTRand( const uint32& oneSeed );
00073 MTRand( uint32 *const bigSeed );
00074 MTRand();
00075
00076
00077
00078
00079
00080 double rand();
00081 double rand( const double& n );
00082 double randExc();
00083 double randExc( const double& n );
00084 uint32 randInt();
00085 uint32 randInt( const uint32& n );
00086 double operator()() { return rand(); }
00087
00088
00089 void seed( uint32 oneSeed );
00090 void seed( uint32 *const bigSeed );
00091 void seed();
00092
00093
00094 void save( uint32* saveArray ) const;
00095 void load( uint32 *const loadArray );
00096 friend ostream& operator<<( ostream& os, const MTRand& mtrand );
00097 friend istream& operator>>( istream& is, MTRand& mtrand );
00098
00099 protected:
00100 void reload();
00101 uint32 hiBit( const uint32& u ) const { return u & 0x80000000U; }
00102 uint32 loBit( const uint32& u ) const { return u & 0x00000001U; }
00103 uint32 loBits( const uint32& u ) const { return u & 0x7fffffffU; }
00104 uint32 mixBits( const uint32& u, const uint32& v ) const
00105 { return hiBit(u) | loBits(v); }
00106 uint32 twist( const uint32& m, const uint32& s0, const uint32& s1 ) const
00107 { return m ^ (mixBits(s0,s1)>>1) ^ (loBit(s1) ? MAGIC : 0U); }
00108 static uint32 hash( time_t t, clock_t c );
00109 };
00110
00111
00112 inline MTRand::MTRand( const uint32& oneSeed )
00113 { seed(oneSeed); }
00114
00115 inline MTRand::MTRand( uint32 *const bigSeed )
00116 { seed(bigSeed); }
00117
00118 inline MTRand::MTRand()
00119 { seed(); }
00120
00121 inline double MTRand::rand()
00122 { return double(randInt()) * 2.3283064370807974e-10; }
00123
00124 inline double MTRand::rand( const double& n )
00125 { return rand() * n; }
00126
00127 inline double MTRand::randExc()
00128 { return double(randInt()) * 2.3283064365386963e-10; }
00129
00130 inline double MTRand::randExc( const double& n )
00131 { return randExc() * n; }
00132
00133 inline MTRand::uint32 MTRand::randInt()
00134 {
00135 if( left == 0 ) reload();
00136 --left;
00137
00138 register uint32 s1;
00139 s1 = *pNext++;
00140 s1 ^= (s1 >> 11);
00141 s1 ^= (s1 << 7) & 0x9d2c5680U;
00142 s1 ^= (s1 << 15) & 0xefc60000U;
00143 return ( s1 ^ (s1 >> 18) );
00144 }
00145
00146
00147 inline MTRand::uint32 MTRand::randInt( const uint32& n )
00148 { return int( randExc() * (n+1) ); }
00149
00150
00151
00152 inline void MTRand::seed( uint32 oneSeed )
00153 {
00154
00155 register uint32 *s;
00156 register int i;
00157 for( i = N, s = state;
00158 i--;
00159 *s = oneSeed & 0xffff0000,
00160 *s++ |= ( (oneSeed *= 69069U)++ & 0xffff0000 ) >> 16,
00161 (oneSeed *= 69069U)++ ) {}
00162 reload();
00163 }
00164
00165
00166 inline void MTRand::seed( uint32 *const bigSeed )
00167 {
00168
00169
00170
00171
00172
00173 register uint32 *s = state, *b = bigSeed;
00174 register int i = N;
00175 for( ; i--; *s++ = *b++ ) {}
00176 reload();
00177 }
00178
00179
00180 inline void MTRand::seed()
00181 {
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204 seed( hash( time((long*)NULL), clock() ) );
00205 }
00206
00207
00208 inline void MTRand::reload()
00209 {
00210
00211
00212 register uint32 *p = state;
00213 register int i;
00214 for( i = N - M; i--; )
00215 *p++ = twist( p[M], p[0], p[1] );
00216 for( i = M; --i; )
00217 *p++ = twist( p[M-N], p[0], p[1] );
00218 *p = twist( p[M-N], p[0], state[0] );
00219
00220 left = N, pNext = state;
00221 }
00222
00223
00224 inline MTRand::uint32 MTRand::hash( time_t t, clock_t c )
00225 {
00226
00227
00228
00229
00230 static uint32 differ = 0;
00231
00232 uint32 h1 = 0;
00233 unsigned char *p = (unsigned char *) &t;
00234 for( size_t i = 0; i < sizeof(t); ++i )
00235 {
00236 h1 *= UCHAR_MAX + 2U;
00237 h1 += p[i];
00238 }
00239 uint32 h2 = 0;
00240 p = (unsigned char *) &c;
00241 for( size_t j = 0; j < sizeof(c); ++j )
00242 {
00243 h2 *= UCHAR_MAX + 2U;
00244 h2 += p[j];
00245 }
00246 return ( h1 + differ++ ) ^ h2;
00247 }
00248
00249
00250 inline void MTRand::save( uint32* saveArray ) const
00251 {
00252 register uint32 *sa = saveArray;
00253 register const uint32 *s = state;
00254 register int i = N;
00255 for( ; i--; *sa++ = *s++ ) {}
00256 *sa = left;
00257 }
00258
00259
00260 inline void MTRand::load( uint32 *const loadArray )
00261 {
00262 register uint32 *s = state;
00263 register uint32 *la = loadArray;
00264 register int i = N;
00265 for( ; i--; *s++ = *la++ ) {}
00266 left = *la;
00267 pNext = &state[N-left];
00268 }
00269
00270
00271 inline ostream& operator<<( ostream& os, const MTRand& mtrand )
00272 {
00273 register const MTRand::uint32 *s = mtrand.state;
00274 register int i = mtrand.N;
00275 for( ; i--; os << *s++ << "\t" ) {}
00276 return os << mtrand.left;
00277 }
00278
00279
00280 inline istream& operator>>( istream& is, MTRand& mtrand )
00281 {
00282 register MTRand::uint32 *s = mtrand.state;
00283 register int i = mtrand.N;
00284 for( ; i--; is >> *s++ ) {}
00285 is >> mtrand.left;
00286 mtrand.pNext = &mtrand.state[mtrand.N-mtrand.left];
00287 return is;
00288 }
00289
00290 #endif //MERSENNETWISTER
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312