86 |
|
|
87 |
|
//Methods |
88 |
|
public: |
89 |
< |
MTRand( const uint32& oneSeed, nstrides = 1, stride = 0); // initialize with a simple uint32 |
90 |
< |
MTRand( uint32 *const bigSeed, uint32 const seedLength = N, nstrides = 1, stride = 0); // or an array |
91 |
< |
MTRand(nstrides = 1, stride = 0); // auto-initialize with /dev/urandom or time() and clock() |
89 |
> |
MTRand( const uint32& oneSeed, int nstrides = 1, int stride = 0); // initialize with a simple uint32 |
90 |
> |
MTRand( uint32 *const bigSeed, uint32 const seedLength = N, int nstrides = 1, int stride = 0); // or an array |
91 |
> |
MTRand(int nstrides = 1, int stride = 0); // auto-initialize with /dev/urandom or time() and clock() |
92 |
|
|
93 |
|
// Do NOT use for CRYPTOGRAPHY without securely hashing several returned |
94 |
|
// values together, otherwise the generator state can be learned after |
101 |
|
double randExc( const double& n ); // real number in [0,n) |
102 |
|
double randDblExc(); // real number in (0,1) |
103 |
|
double randDblExc( const double& n ); // real number in (0,n) |
104 |
< |
uint32 randInt(); // integer in [0,2^32-1] |
104 |
> |
uint32 randInt(); // integer in [0,2^32-1] (modified for striding) |
105 |
> |
uint32 rawRandInt(); // original randInt |
106 |
|
uint32 randInt( const uint32& n ); // integer in [0,n] for n < 2^32 |
107 |
|
double operator()() { return rand(); } // same as rand() |
108 |
|
|
185 |
|
return mean + r * cos(phi); |
186 |
|
} |
187 |
|
|
188 |
< |
inline MTRand::uint32 MTRand::randInt() |
188 |
> |
/** |
189 |
> |
* This function is modified from the original to allow for random |
190 |
> |
* streams on parallel jobs. It now takes numbers from by striding |
191 |
> |
* through the random stream and picking up only one of the random |
192 |
> |
* numbers per nstrides_. The number it picks is the stride_'th |
193 |
> |
* number in the stride sequence. |
194 |
> |
*/ |
195 |
> |
inline MTRand::uint32 MTRand::randInt() { |
196 |
> |
|
197 |
> |
uint32 ranNums[nstrides_]; |
198 |
> |
|
199 |
> |
for (int i = 0; i < nstrides_; ++i) { |
200 |
> |
ranNums[i] = rawRandInt(); |
201 |
> |
} |
202 |
> |
|
203 |
> |
return ranNums[stride_]; |
204 |
> |
} |
205 |
> |
|
206 |
> |
/** |
207 |
> |
* This is the original randInt function which implements the mersenne |
208 |
> |
* twister. |
209 |
> |
*/ |
210 |
> |
inline MTRand::uint32 MTRand::rawRandInt() |
211 |
|
{ |
212 |
|
// Pull a 32-bit integer from the generator state |
213 |
|
// Every other access function simply transforms the numbers extracted here |
214 |
< |
|
215 |
< |
uint32 ranNums[nstrides]; |
216 |
< |
|
217 |
< |
for (int i = 0; i < nstrides; ++i) { |
218 |
< |
if( left == 0 ) { |
219 |
< |
reload(); |
220 |
< |
} |
221 |
< |
|
222 |
< |
--left; |
223 |
< |
|
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]; |
214 |
> |
|
215 |
> |
if( left == 0 ) reload(); |
216 |
> |
--left; |
217 |
> |
|
218 |
> |
register uint32 s1; |
219 |
> |
s1 = *pNext++; |
220 |
> |
s1 ^= (s1 >> 11); |
221 |
> |
s1 ^= (s1 << 7) & 0x9d2c5680UL; |
222 |
> |
s1 ^= (s1 << 15) & 0xefc60000UL; |
223 |
> |
return ( s1 ^ (s1 >> 18) ); |
224 |
|
} |
225 |
|
|
226 |
|
inline MTRand::uint32 MTRand::randInt( const uint32& n ) |
289 |
|
|
290 |
|
inline void MTRand::seed() |
291 |
|
{ |
292 |
< |
// Seed the generator with an array from /dev/urandom if available |
293 |
< |
// Otherwise use a hash of time() and clock() values |
294 |
< |
|
295 |
< |
// First try getting an array from /dev/urandom |
296 |
< |
FILE* urandom = fopen( "/dev/urandom", "rb" ); |
297 |
< |
if( urandom ) |
298 |
< |
{ |
299 |
< |
uint32 bigSeed[N]; |
300 |
< |
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() ) ); |
292 |
> |
vector<uint32> seeds; |
293 |
> |
|
294 |
> |
seeds = generateSeeds(); |
295 |
> |
|
296 |
> |
if (seeds.size() == 1) { |
297 |
> |
seed( seeds[0] ); |
298 |
> |
} else { |
299 |
> |
seed( &seeds[0], seeds.size() ); |
300 |
> |
} |
301 |
|
} |
302 |
|
|
303 |
|
|
304 |
+ |
inline vector<uint32> MTRand::generateSeeds() { |
305 |
+ |
// Seed the generator with an array from /dev/urandom if available |
306 |
+ |
// Otherwise use a hash of time() and clock() values |
307 |
+ |
|
308 |
+ |
vector<uint32> bigSeed; |
309 |
+ |
|
310 |
+ |
// First try getting an array from /dev/urandom |
311 |
+ |
FILE* urandom = fopen( "/dev/urandom", "rb" ); |
312 |
+ |
if( urandom ) |
313 |
+ |
{ |
314 |
+ |
bigSeed.resize(N); |
315 |
+ |
register uint32 *s = &bigSeed[0]; |
316 |
+ |
register int i = N; |
317 |
+ |
register bool success = true; |
318 |
+ |
while( success && i-- ) |
319 |
+ |
success = fread( s++, sizeof(uint32), 1, urandom ); |
320 |
+ |
fclose(urandom); |
321 |
+ |
if( success ) { return bigSeed; } |
322 |
+ |
} |
323 |
+ |
|
324 |
+ |
// Was not successful, so use time() and clock() instead |
325 |
+ |
|
326 |
+ |
bigSeed.push_back(hash( time(NULL), clock())); |
327 |
+ |
return bigSeed; |
328 |
+ |
} |
329 |
+ |
|
330 |
+ |
|
331 |
|
inline void MTRand::initialize( const uint32 seed ) |
332 |
|
{ |
333 |
|
// Initialize generator state with seed |