00001
00025 #ifndef EO_RANDOM_NUMBER_GENERATOR
00026 #define EO_RANDOM_NUMBER_GENERATOR
00027
00028
00029
00030
00031
00032
00033
00034
00035 # if (defined _MSC_VER)
00036 typedef unsigned long uint32_t;
00037 #else
00038 #if (! defined __sun)
00039
00040
00041 #include <stdint.h>
00042 #else
00043 #include <inttypes.h>
00044 #endif
00045 #endif
00046
00047 #include <cmath>
00048 #include <vector>
00049 #include "eoPersistent.h"
00050 #include "eoObject.h"
00051
00052
00111 class eoRng : public eoObject, public eoPersistent
00112 {
00113 public :
00114
00121 eoRng(uint32_t s)
00122 : state(0), next(0), left(-1), cached(false)
00123 {
00124 state = new uint32_t[N+1];
00125 initialize(2*s);
00126 }
00127
00129 ~eoRng()
00130 {
00131 delete [] state;
00132 }
00133
00144 void reseed(uint32_t s)
00145 {
00146 initialize(2*s);
00147 }
00148
00159 void oldReseed(uint32_t s)
00160 {
00161 initialize(s);
00162 }
00163
00169 double uniform(double m = 1.0)
00170 {
00171 return m * double(rand()) / double(1.0 + rand_max());
00172 }
00173
00179 uint32_t random(uint32_t m)
00180 {
00181
00182
00183
00184
00185
00186 return uint32_t(uniform() * double(m));
00187 }
00188
00196 bool flip(double bias=0.5)
00197 {
00198 return uniform() < bias;
00199 }
00200
00207 double normal();
00208
00216 double normal(double stdev)
00217 { return stdev * normal(); }
00218
00227 double normal(double mean, double stdev)
00228 { return mean + normal(stdev); }
00229
00235 double negexp(double mean)
00236 {
00237 return -mean*log(double(rand()) / rand_max());
00238 }
00239
00243 uint32_t rand();
00244
00248 uint32_t rand_max() const { return uint32_t(0xffffffff); }
00249
00256 template <typename TYPE>
00257 int roulette_wheel(const std::vector<TYPE>& vec, TYPE total = 0)
00258 {
00259 if (total == 0)
00260 {
00261 for (unsigned i = 0; i < vec.size(); ++i)
00262 total += vec[i];
00263 }
00264 double fortune = uniform() * total;
00265 int i = 0;
00266 while (fortune >= 0)
00267 {
00268 fortune -= vec[i++];
00269 }
00270 return --i;
00271 }
00272
00273
00278 template <typename TYPE>
00279 const TYPE& choice(const std::vector<TYPE>& vec)
00280 { return vec[random(vec.size())]; }
00281
00282
00293 template <typename TYPE>
00294 TYPE& choice(std::vector<TYPE>& vec)
00295 { return vec[random(vec.size())]; }
00296
00301 void printOn(std::ostream& _os) const
00302 {
00303 for (int i = 0; i < N; ++i)
00304 {
00305 _os << state[i] << ' ';
00306 }
00307 _os << int(next - state) << ' ';
00308 _os << left << ' ' << cached << ' ' << cacheValue;
00309 }
00310
00315 void readFrom(std::istream& _is)
00316 {
00317 for (int i = 0; i < N; ++i)
00318 {
00319 _is >> state[i];
00320 }
00321
00322 int n;
00323 _is >> n;
00324 next = state + n;
00325
00326 _is >> left;
00327 _is >> cached;
00328 _is >> cacheValue;
00329 }
00330
00331 std::string className() const { return "Mersenne-Twister"; }
00332
00333 private:
00334
00335 uint32_t restart();
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383 void initialize(uint32_t seed);
00384
00386 uint32_t *state;
00387
00389 uint32_t *next;
00390
00392 int left;
00393
00395 bool cached;
00396
00398 double cacheValue;
00399
00401 static const int N;
00402
00404 static const int M;
00405
00407 static const uint32_t K;
00408
00409
00418 eoRng(const eoRng&);
00419
00424 eoRng& operator=(const eoRng&);
00425 };
00426
00427
00428
00429 namespace eo
00430 {
00432 extern eoRng rng;
00433 }
00434 using eo::rng;
00435
00436
00437
00438
00439
00440
00441
00442 #define hiBit(u) ((u) & 0x80000000U) // mask all but highest bit of u
00443 #define loBit(u) ((u) & 0x00000001U) // mask all but lowest bit of u
00444 #define loBits(u) ((u) & 0x7FFFFFFFU) // mask the highest bit of u
00445 #define mixBits(u, v) (hiBit(u)|loBits(v)) // move hi bit of u to hi bit of v
00446
00447
00448 inline void eoRng::initialize(uint32_t seed)
00449 {
00450 left = -1;
00451
00452 register uint32_t x = (seed | 1U) & 0xFFFFFFFFU, *s = state;
00453 register int j;
00454
00455 for(left=0, *s++=x, j=N; --j;
00456 *s++ = (x*=69069U) & 0xFFFFFFFFU);
00457 }
00458
00459
00460
00461 inline uint32_t eoRng::restart()
00462 {
00463 register uint32_t *p0=state, *p2=state+2, *pM=state+M, s0, s1;
00464 register int j;
00465
00466 left=N-1, next=state+1;
00467
00468 for(s0=state[0], s1=state[1], j=N-M+1; --j; s0=s1, s1=*p2++)
00469 *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);
00470
00471 for(pM=state, j=M; --j; s0=s1, s1=*p2++)
00472 *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);
00473
00474 s1=state[0], *p0 = *pM ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);
00475 s1 ^= (s1 >> 11);
00476 s1 ^= (s1 << 7) & 0x9D2C5680U;
00477 s1 ^= (s1 << 15) & 0xEFC60000U;
00478 return(s1 ^ (s1 >> 18));
00479 }
00480
00481
00482
00483 inline uint32_t eoRng::rand()
00484 {
00485 if(--left < 0)
00486 return(restart());
00487 uint32_t y = *next++;
00488 y ^= (y >> 11);
00489 y ^= (y << 7) & 0x9D2C5680U;
00490 y ^= (y << 15) & 0xEFC60000U;
00491 return(y ^ (y >> 18));
00492 }
00493
00494
00495
00496 inline double eoRng::normal()
00497 {
00498 if (cached) {
00499 cached = false;
00500 return cacheValue;
00501 }
00502 double rSquare, var1, var2;
00503 do {
00504 var1 = 2.0 * uniform() - 1.0;
00505 var2 = 2.0 * uniform() - 1.0;
00506 rSquare = var1 * var1 + var2 * var2;
00507 } while (rSquare >= 1.0 || rSquare == 0.0);
00508 double factor = sqrt(-2.0 * log(rSquare) / rSquare);
00509 cacheValue = var1 * factor;
00510 cached = true;
00511 return (var2 * factor);
00512 }
00513
00514
00515
00516 namespace eo
00517 {
00534 template <typename T>
00535 inline T random(const T& min, const T& max) {
00536 return static_cast<T>(rng.uniform() * (max-min)) + min; }
00537
00555 template <typename T>
00556 inline T random(const T& max) {
00557 return static_cast<T>(rng.uniform() * max); }
00558
00566 inline double normal() { return rng.normal(); }
00567 }
00568
00569
00570 #endif
00571
00572
00573
00574
00575
00576
00577
00578
00579