60 #ifndef MERSENNETWISTER_WFMATH_H
61 #define MERSENNETWISTER_WFMATH_H
77 typedef unsigned long uint32;
80 enum { SAVE = N + 1 };
92 MTRand(
const uint32& oneSeed );
93 MTRand( uint32 *
const bigSeed, uint32
const seedLength = N );
102 double rand(
const double& n );
104 double randExc(
const double& n );
106 double randDblExc(
const double& n );
108 uint32 randInt(
const uint32& n );
109 double operator()() {
return rand(); }
115 double randNorm(
const double& mean = 0.0,
const double& variance = 0.0 );
118 void seed(
const uint32 oneSeed );
119 void seed( uint32 *
const bigSeed,
const uint32 seedLength = N );
123 void save( uint32* saveArray )
const;
124 void load( uint32 *
const loadArray );
125 friend std::ostream& operator<<( std::ostream& os,
const MTRand& mtrand );
126 friend std::istream& operator>>( std::istream& is, MTRand& mtrand );
128 static MTRand instance;
131 void initialize(
const uint32 oneSeed );
133 uint32 hiBit(
const uint32& u )
const {
return u & 0x80000000UL; }
134 uint32 loBit(
const uint32& u )
const {
return u & 0x00000001UL; }
135 uint32 loBits(
const uint32& u )
const {
return u & 0x7fffffffUL; }
136 uint32 mixBits(
const uint32& u,
const uint32& v )
const
137 {
return hiBit(u) | loBits(v); }
138 uint32 twist(
const uint32& m,
const uint32& s0,
const uint32& s1 )
const
139 {
return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL); }
143 inline MTRand::MTRand(
const uint32& oneSeed ) : pNext(0), left(0)
146 inline MTRand::MTRand( uint32 *
const bigSeed,
const uint32 seedLength ) : pNext(0), left(0)
147 { seed(bigSeed,seedLength); }
149 inline MTRand::MTRand() : pNext(0), left(0)
152 inline double MTRand::rand()
153 {
return double(randInt()) * (1.0/4294967295.0); }
155 inline double MTRand::rand(
const double& n )
156 {
return rand() * n; }
158 inline double MTRand::randExc()
159 {
return double(randInt()) * (1.0/4294967296.0); }
161 inline double MTRand::randExc(
const double& n )
162 {
return randExc() * n; }
164 inline double MTRand::randDblExc()
165 {
return (
double(randInt()) + 0.5 ) * (1.0/4294967296.0); }
167 inline double MTRand::randDblExc(
const double& n )
168 {
return randDblExc() * n; }
170 inline double MTRand::rand53()
172 uint32 a = randInt() >> 5, b = randInt() >> 6;
173 return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0);
176 inline double MTRand::randNorm(
const double& mean,
const double& variance )
180 double r = sqrt( -2.0 * log( 1.0-randDblExc()) ) * variance;
181 double phi = 2.0 * 3.14159265358979323846264338328 * randExc();
182 return mean + r * cos(phi);
185 inline MTRand::uint32 MTRand::randInt()
190 if( left == 0 ) reload();
196 s1 ^= (s1 << 7) & 0x9d2c5680UL;
197 s1 ^= (s1 << 15) & 0xefc60000UL;
198 return ( s1 ^ (s1 >> 18) );
201 inline MTRand::uint32 MTRand::randInt(
const uint32& n )
215 i = randInt() & used;
221 inline void MTRand::seed(
const uint32 oneSeed )
229 inline void MTRand::seed( uint32 *
const bigSeed,
const uint32 seedLength )
237 initialize(19650218UL);
239 register uint32 j = 0;
240 register int k = ( N > seedLength ? N : seedLength );
244 state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL );
245 state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
246 state[i] &= 0xffffffffUL;
248 if( i >= N ) { state[0] = state[N-1]; i = 1; }
249 if( j >= seedLength ) j = 0;
251 for( k = N - 1; k; --k )
254 state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL );
256 state[i] &= 0xffffffffUL;
258 if( i >= N ) { state[0] = state[N-1]; i = 1; }
260 state[0] = 0x80000000UL;
265 inline void MTRand::initialize(
const uint32 seed )
271 register uint32 *s = state;
272 register uint32 *r = state;
274 *s++ = seed & 0xffffffffUL;
277 *s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL;
283 inline void MTRand::reload()
287 register uint32 *p = state;
289 for( i = N - M; i--; ++p )
290 *p = twist( p[M], p[0], p[1] );
291 for( i = M; --i; ++p )
292 *p = twist( p[M-N], p[0], p[1] );
293 *p = twist( p[M-N], p[0], state[0] );
295 left = N, pNext = state;
300 inline void MTRand::save( uint32* saveArray )
const
302 register uint32 *sa = saveArray;
303 register const uint32 *s = state;
305 for( ; i--; *sa++ = *s++ ) {}
310 inline void MTRand::load( uint32 *
const loadArray )
312 register uint32 *s = state;
313 register uint32 *la = loadArray;
315 for( ; i--; *s++ = *la++ ) {}
317 pNext = &state[N-left];
323 #endif // MERSENNETWISTER_H