| 66 |
|
#include <stdio.h> |
| 67 |
|
#include <time.h> |
| 68 |
|
#include <math.h> |
| 69 |
+ |
#include <vector> |
| 70 |
+ |
namespace oopse { |
| 71 |
|
|
| 72 |
|
class MTRand { |
| 73 |
|
// Data |
| 103 |
|
double randExc( const double& n ); // real number in [0,n) |
| 104 |
|
double randDblExc(); // real number in (0,1) |
| 105 |
|
double randDblExc( const double& n ); // real number in (0,n) |
| 106 |
< |
uint32 randInt(); // integer in [0,2^32-1] |
| 106 |
> |
uint32 randInt(); // integer in [0,2^32-1] (modified for striding) |
| 107 |
> |
uint32 rawRandInt(); // original randInt |
| 108 |
|
uint32 randInt( const uint32& n ); // integer in [0,n] for n < 2^32 |
| 109 |
|
double operator()() { return rand(); } // same as rand() |
| 110 |
|
|
| 118 |
|
void seed( const uint32 oneSeed ); |
| 119 |
|
void seed( uint32 *const bigSeed, const uint32 seedLength = N ); |
| 120 |
|
void seed(); |
| 121 |
< |
|
| 121 |
> |
|
| 122 |
> |
std::vector<uint32>generateSeeds(); |
| 123 |
> |
|
| 124 |
|
// Saving and loading generator state |
| 125 |
|
void save( uint32* saveArray ) const; // to array of size SAVE |
| 126 |
|
void load( uint32 *const loadArray ); // from such array |
| 189 |
|
return mean + r * cos(phi); |
| 190 |
|
} |
| 191 |
|
|
| 192 |
< |
inline MTRand::uint32 MTRand::randInt() |
| 192 |
> |
/** |
| 193 |
> |
* This function is modified from the original to allow for random |
| 194 |
> |
* streams on parallel jobs. It now takes numbers from by striding |
| 195 |
> |
* through the random stream and picking up only one of the random |
| 196 |
> |
* numbers per nstrides_. The number it picks is the stride_'th |
| 197 |
> |
* number in the stride sequence. |
| 198 |
> |
*/ |
| 199 |
> |
inline MTRand::uint32 MTRand::randInt() { |
| 200 |
> |
|
| 201 |
> |
uint32 ranNums[nstrides_]; |
| 202 |
> |
|
| 203 |
> |
for (int i = 0; i < nstrides_; ++i) { |
| 204 |
> |
ranNums[i] = rawRandInt(); |
| 205 |
> |
} |
| 206 |
> |
|
| 207 |
> |
return ranNums[stride_]; |
| 208 |
> |
} |
| 209 |
> |
|
| 210 |
> |
/** |
| 211 |
> |
* This is the original randInt function which implements the mersenne |
| 212 |
> |
* twister. |
| 213 |
> |
*/ |
| 214 |
> |
inline MTRand::uint32 MTRand::rawRandInt() |
| 215 |
|
{ |
| 216 |
|
// Pull a 32-bit integer from the generator state |
| 217 |
|
// Every other access function simply transforms the numbers extracted here |
| 218 |
< |
|
| 219 |
< |
uint32 ranNums[nstrides_]; |
| 220 |
< |
|
| 221 |
< |
for (int i = 0; i < nstrides_; ++i) { |
| 222 |
< |
if( left == 0 ) { |
| 223 |
< |
reload(); |
| 224 |
< |
} |
| 225 |
< |
|
| 226 |
< |
--left; |
| 227 |
< |
|
| 201 |
< |
register uint32 s1; |
| 202 |
< |
s1 = *pNext++; |
| 203 |
< |
s1 ^= (s1 >> 11); |
| 204 |
< |
s1 ^= (s1 << 7) & 0x9d2c5680UL; |
| 205 |
< |
s1 ^= (s1 << 15) & 0xefc60000UL; |
| 206 |
< |
ranNums[i] = s1 ^ (s1 >> 18); |
| 207 |
< |
} |
| 208 |
< |
|
| 209 |
< |
return ranNums[stride_]; |
| 218 |
> |
|
| 219 |
> |
if( left == 0 ) reload(); |
| 220 |
> |
--left; |
| 221 |
> |
|
| 222 |
> |
register uint32 s1; |
| 223 |
> |
s1 = *pNext++; |
| 224 |
> |
s1 ^= (s1 >> 11); |
| 225 |
> |
s1 ^= (s1 << 7) & 0x9d2c5680UL; |
| 226 |
> |
s1 ^= (s1 << 15) & 0xefc60000UL; |
| 227 |
> |
return ( s1 ^ (s1 >> 18) ); |
| 228 |
|
} |
| 229 |
|
|
| 230 |
|
inline MTRand::uint32 MTRand::randInt( const uint32& n ) |
| 293 |
|
|
| 294 |
|
inline void MTRand::seed() |
| 295 |
|
{ |
| 296 |
< |
// Seed the generator with an array from /dev/urandom if available |
| 297 |
< |
// Otherwise use a hash of time() and clock() values |
| 298 |
< |
|
| 299 |
< |
// First try getting an array from /dev/urandom |
| 300 |
< |
FILE* urandom = fopen( "/dev/urandom", "rb" ); |
| 301 |
< |
if( urandom ) |
| 302 |
< |
{ |
| 303 |
< |
uint32 bigSeed[N]; |
| 304 |
< |
register uint32 *s = bigSeed; |
| 287 |
< |
register int i = N; |
| 288 |
< |
register bool success = true; |
| 289 |
< |
while( success && i-- ) |
| 290 |
< |
success = fread( s++, sizeof(uint32), 1, urandom ); |
| 291 |
< |
fclose(urandom); |
| 292 |
< |
if( success ) { seed( bigSeed, N ); return; } |
| 293 |
< |
} |
| 294 |
< |
|
| 295 |
< |
// Was not successful, so use time() and clock() instead |
| 296 |
< |
seed( hash( time(NULL), clock() ) ); |
| 296 |
> |
std::vector<uint32> seeds; |
| 297 |
> |
|
| 298 |
> |
seeds = generateSeeds(); |
| 299 |
> |
|
| 300 |
> |
if (seeds.size() == 1) { |
| 301 |
> |
seed( seeds[0] ); |
| 302 |
> |
} else { |
| 303 |
> |
seed( &seeds[0], seeds.size() ); |
| 304 |
> |
} |
| 305 |
|
} |
| 306 |
|
|
| 307 |
|
|
| 308 |
+ |
inline std::vector<MTRand::uint32> MTRand::generateSeeds() { |
| 309 |
+ |
// Seed the generator with an array from /dev/urandom if available |
| 310 |
+ |
// Otherwise use a hash of time() and clock() values |
| 311 |
+ |
|
| 312 |
+ |
std::vector<uint32> bigSeed; |
| 313 |
+ |
|
| 314 |
+ |
// First try getting an array from /dev/urandom |
| 315 |
+ |
FILE* urandom = fopen( "/dev/urandom", "rb" ); |
| 316 |
+ |
if( urandom ) |
| 317 |
+ |
{ |
| 318 |
+ |
bigSeed.resize(N); |
| 319 |
+ |
register uint32 *s = &bigSeed[0]; |
| 320 |
+ |
register int i = N; |
| 321 |
+ |
register bool success = true; |
| 322 |
+ |
while( success && i-- ) |
| 323 |
+ |
success = fread( s++, sizeof(uint32), 1, urandom ); |
| 324 |
+ |
fclose(urandom); |
| 325 |
+ |
if( success ) { return bigSeed; } |
| 326 |
+ |
} |
| 327 |
+ |
|
| 328 |
+ |
// Was not successful, so use time() and clock() instead |
| 329 |
+ |
|
| 330 |
+ |
bigSeed.push_back(hash( time(NULL), clock())); |
| 331 |
+ |
return bigSeed; |
| 332 |
+ |
} |
| 333 |
+ |
|
| 334 |
+ |
|
| 335 |
|
inline void MTRand::initialize( const uint32 seed ) |
| 336 |
|
{ |
| 337 |
|
// Initialize generator state with seed |
| 432 |
|
return is; |
| 433 |
|
} |
| 434 |
|
|
| 435 |
+ |
} |
| 436 |
|
#endif // MERSENNETWISTER_H |
| 437 |
|
|
| 438 |
|
// Change log: |
