66 |
|
#include <stdio.h> |
67 |
|
#include <time.h> |
68 |
|
#include <math.h> |
69 |
+ |
#include <vector> |
70 |
+ |
namespace oopse { |
71 |
|
|
72 |
|
class MTRand { |
73 |
|
// Data |
88 |
|
|
89 |
|
//Methods |
90 |
|
public: |
91 |
< |
MTRand( const uint32& oneSeed, int nstrides = 1, int stride = 0); // initialize with a simple uint32 |
92 |
< |
MTRand( uint32 *const bigSeed, uint32 const seedLength = N, int nstrides = 1, int stride = 0); // or an array |
93 |
< |
MTRand(int nstrides = 1, int stride = 0); // auto-initialize with /dev/urandom or time() and clock() |
91 |
> |
MTRand( const uint32& oneSeed, int nstrides, int stride); // initialize with a simple uint32 |
92 |
> |
MTRand( uint32 *const bigSeed, uint32 const seedLength, int nstrides, int stride); // or an array |
93 |
> |
MTRand(int nstrides, int stride); // auto-initialize with /dev/urandom or time() and clock() |
94 |
|
|
95 |
|
// Do NOT use for CRYPTOGRAPHY without securely hashing several returned |
96 |
|
// values together, otherwise the generator state can be learned after |
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 |
> |
std::vector<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: |