60 #ifndef MERSENNETWISTER_WFMATH_H 61 #define MERSENNETWISTER_WFMATH_H 77 typedef unsigned long uint32;
79 static const uint32 N = 624;
80 static const uint32 SAVE = N + 1;
83 static const uint32 M = 397;
92 MTRand(
const uint32& oneSeed );
93 MTRand( uint32 *
const bigSeed, uint32
const seedLength = N );
101 template<
typename FloatT>
104 float randf(
const float& n );
106 double rand(
const double& n );
108 double randExc(
const double& n );
110 double randDblExc(
const double& n );
112 uint32 randInt(
const uint32& n );
113 double operator()() {
return rand(); }
119 double randNorm(
const double& mean = 0.0,
const double& variance = 0.0 );
122 void seed(
const uint32 oneSeed );
123 void seed( uint32 *
const bigSeed,
const uint32 seedLength = N );
127 void save( uint32* saveArray )
const;
128 void load( uint32 *
const loadArray );
129 friend std::ostream& operator<<( std::ostream& os,
const MTRand& mtrand );
130 friend std::istream& operator>>( std::istream& is, MTRand& mtrand );
132 static MTRand instance;
135 void initialize(
const uint32 oneSeed );
137 uint32 hiBit(
const uint32& u )
const {
return u & 0x80000000UL; }
138 uint32 loBit(
const uint32& u )
const {
return u & 0x00000001UL; }
139 uint32 loBits(
const uint32& u )
const {
return u & 0x7fffffffUL; }
140 uint32 mixBits(
const uint32& u,
const uint32& v )
const 141 {
return hiBit(u) | loBits(v); }
142 uint32 twist(
const uint32& m,
const uint32& s0,
const uint32& s1 )
const 143 {
return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL); }
147 inline MTRand::MTRand(
const uint32& oneSeed ) : pNext(0), left(0)
150 inline MTRand::MTRand( uint32 *
const bigSeed,
const uint32 seedLength ) : pNext(0), left(0)
151 { seed(bigSeed,seedLength); }
153 inline MTRand::MTRand() : pNext(0), left(0)
157 inline float MTRand::rand<float>()
158 {
return float(randInt()) * (1.0f/4294967295.0f); }
161 inline double MTRand::rand<double>()
162 {
return double(randInt()) * (1.0/4294967295.0); }
164 inline float MTRand::randf()
165 {
return float(randInt()) * (1.0f/4294967295.0f); }
167 inline float MTRand::randf(
const float& n )
168 {
return randf() * n; }
170 inline double MTRand::rand()
171 {
return double(randInt()) * (1.0/4294967295.0); }
173 inline double MTRand::rand(
const double& n )
174 {
return rand() * n; }
176 inline double MTRand::randExc()
177 {
return double(randInt()) * (1.0/4294967296.0); }
179 inline double MTRand::randExc(
const double& n )
180 {
return randExc() * n; }
182 inline double MTRand::randDblExc()
183 {
return (
double(randInt()) + 0.5 ) * (1.0/4294967296.0); }
185 inline double MTRand::randDblExc(
const double& n )
186 {
return randDblExc() * n; }
188 inline double MTRand::rand53()
190 uint32 a = randInt() >> 5, b = randInt() >> 6;
191 return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0);
194 inline double MTRand::randNorm(
const double& mean,
const double& variance )
198 double r = std::sqrt( -2.0 * std::log( 1.0-randDblExc()) ) * variance;
199 double phi = 2.0 * 3.14159265358979323846264338328 * randExc();
200 return mean + r * std::cos(phi);
203 inline MTRand::uint32 MTRand::randInt()
208 if( left == 0 ) reload();
214 s1 ^= (s1 << 7) & 0x9d2c5680UL;
215 s1 ^= (s1 << 15) & 0xefc60000UL;
216 return ( s1 ^ (s1 >> 18) );
219 inline MTRand::uint32 MTRand::randInt(
const uint32& n )
233 i = randInt() & used;
239 inline void MTRand::seed(
const uint32 oneSeed )
247 inline void MTRand::seed( uint32 *
const bigSeed,
const uint32 seedLength )
255 initialize(19650218UL);
256 register unsigned i = 1;
257 register uint32 j = 0;
258 register unsigned k = ( N > seedLength ? N : seedLength );
262 state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL );
263 state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
264 state[i] &= 0xffffffffUL;
266 if( i >= N ) { state[0] = state[N-1]; i = 1; }
267 if( j >= seedLength ) j = 0;
269 for( k = N - 1; k; --k )
272 state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL );
274 state[i] &= 0xffffffffUL;
276 if( i >= N ) { state[0] = state[N-1]; i = 1; }
278 state[0] = 0x80000000UL;
283 inline void MTRand::initialize(
const uint32 seed )
289 register uint32 *s = state;
290 register uint32 *r = state;
291 register unsigned i = 1;
292 *s++ = seed & 0xffffffffUL;
295 *s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL;
301 inline void MTRand::reload()
305 register uint32 *p = state;
307 for( i = N - M; i--; ++p )
308 *p = twist( p[M], p[0], p[1] );
309 for( i = M; --i; ++p )
310 *p = twist( p[M-N], p[0], p[1] );
311 *p = twist( p[M-N], p[0], state[0] );
313 left = N, pNext = state;
318 inline void MTRand::save( uint32* saveArray )
const 320 register uint32 *sa = saveArray;
321 register const uint32 *s = state;
323 for( ; i--; *sa++ = *s++ ) {}
328 inline void MTRand::load( uint32 *
const loadArray )
330 register uint32 *s = state;
331 register uint32 *la = loadArray;
333 for( ; i--; *s++ = *la++ ) {}
335 pNext = &state[N-left];
341 #endif // MERSENNETWISTER_H Generic library namespace.
Definition: atlasconv.h:45