58 #ifndef NPSTAT_MERSENNETWISTERIMPL_HH_
59 #define NPSTAT_MERSENNETWISTERIMPL_HH_
76 typedef unsigned long uint32;
79 enum { SAVE = N + 1 };
90 MTRand(
const uint32 oneSeed );
91 MTRand( uint32 *
const bigSeed, uint32
const seedLength = N );
101 uint32 randInt(
const uint32 n );
103 double rand(
const double n );
105 double randExc(
const double n );
107 double randDblExc(
const double n );
114 double randNorm(
const double mean = 0.0,
const double stddev = 1.0 );
117 void seed(
const uint32 oneSeed );
118 void seed( uint32 *
const bigSeed,
const uint32 seedLength = N );
122 void save( uint32* saveArray )
const;
123 void load( uint32 *
const loadArray );
124 friend std::ostream& operator<<( std::ostream& os,
const MTRand& mtrand );
125 friend std::istream& operator>>( std::istream& is,
MTRand& mtrand );
129 void initialize(
const uint32 oneSeed );
131 uint32 hiBit(
const uint32 u )
const {
return u & 0x80000000UL; }
132 uint32 loBit(
const uint32 u )
const {
return u & 0x00000001UL; }
133 uint32 loBits(
const uint32 u )
const {
return u & 0x7fffffffUL; }
134 uint32 mixBits(
const uint32 u,
const uint32 v )
const
135 {
return hiBit(u) | loBits(v); }
136 uint32 magic(
const uint32 u )
const
137 {
return loBit(u) ? 0x9908b0dfUL : 0x0UL; }
138 uint32 twist(
const uint32 m,
const uint32 s0,
const uint32 s1 )
const
139 {
return m ^ (mixBits(s0,s1)>>1) ^ magic(s1); }
140 static uint32 hash( time_t t, clock_t c );
145 inline MTRand::uint32 MTRand::hash( time_t t, clock_t c )
151 static uint32 differ = 0;
154 unsigned char *p = (
unsigned char *) &t;
155 for(
size_t i = 0; i <
sizeof(t); ++i )
157 h1 *= UCHAR_MAX + 2U;
161 p = (
unsigned char *) &c;
162 for(
size_t j = 0; j <
sizeof(c); ++j )
164 h2 *= UCHAR_MAX + 2U;
167 return ( h1 + differ++ ) ^ h2;
170 inline void MTRand::initialize(
const uint32 seed )
176 register uint32 *s = state;
177 register uint32 *r = state;
179 *s++ = seed & 0xffffffffUL;
182 *s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL;
187 inline void MTRand::reload()
191 static const int MmN = int(M) - int(N);
192 register uint32 *p = state;
194 for( i = N - M; i--; ++p )
195 *p = twist( p[M], p[0], p[1] );
196 for( i = M; --i; ++p )
197 *p = twist( p[MmN], p[0], p[1] );
198 *p = twist( p[MmN], p[0], state[0] );
200 left = N, pNext = state;
203 inline void MTRand::seed(
const uint32 oneSeed )
210 inline void MTRand::seed( uint32 *
const bigSeed,
const uint32 seedLength )
218 initialize(19650218UL);
220 register uint32 j = 0;
221 register int k = ( N > seedLength ?
static_cast<uint32
>(N) : seedLength );
225 state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL );
226 state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
227 state[i] &= 0xffffffffUL;
229 if( i >= N ) { state[0] = state[N-1]; i = 1; }
230 if( j >= seedLength ) j = 0;
232 for( k = N - 1; k; --k )
235 state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL );
237 state[i] &= 0xffffffffUL;
239 if( i >= N ) { state[0] = state[N-1]; i = 1; }
241 state[0] = 0x80000000UL;
245 inline void MTRand::seed()
251 FILE* urandom = fopen(
"/dev/urandom",
"rb" );
255 register uint32 *s = bigSeed;
257 register bool success =
true;
258 while( success && i-- )
259 success = fread( s++,
sizeof(uint32), 1, urandom );
261 if( success ) { seed( bigSeed, N );
return; }
265 seed( hash( time(NULL), clock() ) );
268 inline MTRand::MTRand(
const uint32 oneSeed )
271 inline MTRand::MTRand( uint32 *
const bigSeed,
const uint32 seedLength )
272 { seed(bigSeed,seedLength); }
274 inline MTRand::MTRand()
277 inline MTRand::MTRand(
const MTRand& o )
279 register const uint32 *t = o.state;
280 register uint32 *s = state;
282 for( ; i--; *s++ = *t++ ) {}
284 pNext = &state[N-left];
287 inline MTRand::uint32 MTRand::randInt()
292 if( left == 0 ) reload();
298 s1 ^= (s1 << 7) & 0x9d2c5680UL;
299 s1 ^= (s1 << 15) & 0xefc60000UL;
300 return ( s1 ^ (s1 >> 18) );
303 inline MTRand::uint32 MTRand::randInt(
const uint32 n )
317 i = randInt() & used;
322 inline double MTRand::rand()
323 {
return double(randInt()) * (1.0/4294967295.0); }
325 inline double MTRand::rand(
const double n )
326 {
return rand() * n; }
328 inline double MTRand::randExc()
329 {
return double(randInt()) * (1.0/4294967296.0); }
331 inline double MTRand::randExc(
const double n )
332 {
return randExc() * n; }
334 inline double MTRand::randDblExc()
335 {
return (
double(randInt()) + 0.5 ) * (1.0/4294967296.0); }
337 inline double MTRand::randDblExc(
const double n )
338 {
return randDblExc() * n; }
340 inline double MTRand::rand53()
342 uint32 a = randInt() >> 5, b = randInt() >> 6;
343 return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0);
346 inline double MTRand::randNorm(
const double mean,
const double stddev )
353 x = 2.0 * rand() - 1.0;
354 y = 2.0 * rand() - 1.0;
357 while ( r >= 1.0 || r == 0.0 );
358 double s = sqrt( -2.0 * log(r) / r );
359 return mean + x * s * stddev;
362 inline double MTRand::operator()()
367 inline void MTRand::save( uint32* saveArray )
const
369 register const uint32 *s = state;
370 register uint32 *sa = saveArray;
372 for( ; i--; *sa++ = *s++ ) {}
376 inline void MTRand::load( uint32 *
const loadArray )
378 register uint32 *s = state;
379 register uint32 *la = loadArray;
381 for( ; i--; *s++ = *la++ ) {}
383 pNext = &state[N-left];
386 inline std::ostream& operator<<( std::ostream& os,
const MTRand& mtrand )
388 register const MTRand::uint32 *s = mtrand.state;
389 register int i = mtrand.N;
390 for( ; i--; os << *s++ <<
"\t" ) {}
391 return os << mtrand.left;
394 inline std::istream& operator>>( std::istream& is, MTRand& mtrand )
396 register MTRand::uint32 *s = mtrand.state;
397 register int i = mtrand.N;
398 for( ; i--; is >> *s++ ) {}
400 mtrand.pNext = &mtrand.state[mtrand.N-mtrand.left];
404 inline MTRand& MTRand::operator=(
const MTRand& o )
406 if(
this == &o )
return (*
this);
407 register const uint32 *t = o.state;
408 register uint32 *s = state;
410 for( ; i--; *s++ = *t++ ) {}
412 pNext = &state[N-left];
Definition: MersenneTwisterImpl.hh:73
Definition: AbsArrayProjector.hh:14