1 /* -*- Mode: C++ -*- */ 2 /* This is public-domain Mersenne Twister code, 3 * attributed to Michael Brundage. Thanks! 4 * http://www.qbrundage.com/michaelb/pubs/essays/random_number_generation.html 5 */ 6 #undef MT_LEN 7 #undef MT_IA 8 class MTRandom { 9 public: 10 enum Constants { 11 MT_LEN = 624, 12 MT_IA = 397 13 }; 14 15 static const uint32_t TEST_SEED1; 16 static const uint32_t UPPER_MASK; 17 static const uint32_t LOWER_MASK; 18 static const uint32_t MATRIX_A; 19 MTRandom()20 MTRandom() { 21 Init(TEST_SEED1); 22 } 23 MTRandom(uint32_t seed)24 explicit MTRandom(uint32_t seed) { 25 Init(seed); 26 } 27 Rand32()28 uint32_t Rand32 () { 29 uint32_t y; 30 static unsigned long mag01[2] = { 31 0 , MATRIX_A 32 }; 33 34 if (mt_index_ >= MT_LEN) { 35 int kk; 36 37 for (kk = 0; kk < MT_LEN - MT_IA; kk++) { 38 y = (mt_buffer_[kk] & UPPER_MASK) | (mt_buffer_[kk + 1] & LOWER_MASK); 39 mt_buffer_[kk] = mt_buffer_[kk + MT_IA] ^ (y >> 1) ^ mag01[y & 0x1UL]; 40 } 41 for (;kk < MT_LEN - 1; kk++) { 42 y = (mt_buffer_[kk] & UPPER_MASK) | (mt_buffer_[kk + 1] & LOWER_MASK); 43 mt_buffer_[kk] = mt_buffer_[kk + (MT_IA - MT_LEN)] ^ (y >> 1) ^ mag01[y & 0x1UL]; 44 } 45 y = (mt_buffer_[MT_LEN - 1] & UPPER_MASK) | (mt_buffer_[0] & LOWER_MASK); 46 mt_buffer_[MT_LEN - 1] = mt_buffer_[MT_IA - 1] ^ (y >> 1) ^ mag01[y & 0x1UL]; 47 48 mt_index_ = 0; 49 } 50 51 y = mt_buffer_[mt_index_++]; 52 53 y ^= (y >> 11); 54 y ^= (y << 7) & 0x9d2c5680UL; 55 y ^= (y << 15) & 0xefc60000UL; 56 y ^= (y >> 18); 57 58 return y; 59 } 60 ExpRand32(uint32_t mean)61 uint32_t ExpRand32(uint32_t mean) { 62 double mean_d = mean; 63 double erand = log (1.0 / (Rand32() / (double)UINT32_MAX)); 64 uint32_t x = (uint32_t) (mean_d * erand + 0.5); 65 return x; 66 } 67 Rand64()68 uint64_t Rand64() { 69 return ((uint64_t)Rand32() << 32) | Rand32(); 70 } 71 ExpRand64(uint64_t mean)72 uint64_t ExpRand64(uint64_t mean) { 73 double mean_d = mean; 74 double erand = log (1.0 / (Rand64() / (double)UINT32_MAX)); 75 uint64_t x = (uint64_t) (mean_d * erand + 0.5); 76 return x; 77 } 78 79 template <typename T> Rand()80 T Rand() { 81 switch (sizeof(T)) { 82 case sizeof(uint32_t): 83 return Rand32(); 84 case sizeof(uint64_t): 85 return Rand64(); 86 default: 87 cerr << "Invalid sizeof T" << endl; 88 abort(); 89 } 90 } 91 92 template <typename T> ExpRand(T mean)93 T ExpRand(T mean) { 94 switch (sizeof(T)) { 95 case sizeof(uint32_t): 96 return ExpRand32(mean); 97 case sizeof(uint64_t): 98 return ExpRand64(mean); 99 default: 100 cerr << "Invalid sizeof T" << endl; 101 abort(); 102 } 103 } 104 105 private: Init(uint32_t seed)106 void Init(uint32_t seed) { 107 mt_buffer_[0] = seed; 108 mt_index_ = MT_LEN; 109 for (int i = 1; i < MT_LEN; i++) { 110 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ 111 /* In the previous versions, MSBs of the seed affect */ 112 /* only MSBs of the array mt[]. */ 113 /* 2002/01/09 modified by Makoto Matsumoto */ 114 mt_buffer_[i] = 115 (1812433253UL * (mt_buffer_[i-1] ^ (mt_buffer_[i-1] >> 30)) + i); 116 } 117 } 118 119 int mt_index_; 120 uint32_t mt_buffer_[MT_LEN]; 121 }; 122 123 const uint32_t MTRandom::TEST_SEED1 = 5489UL; 124 const uint32_t MTRandom::UPPER_MASK = 0x80000000; 125 const uint32_t MTRandom::LOWER_MASK = 0x7FFFFFFF; 126 const uint32_t MTRandom::MATRIX_A = 0x9908B0DF; 127 128 class MTRandom8 { 129 public: MTRandom8(MTRandom * rand)130 MTRandom8(MTRandom *rand) 131 : rand_(rand) { 132 } 133 Rand8()134 uint8_t Rand8() { 135 uint32_t r = rand_->Rand32(); 136 137 // TODO: make this use a single byte at a time? 138 return (r & 0xff) ^ (r >> 7) ^ (r >> 15) ^ (r >> 21); 139 } 140 141 private: 142 MTRandom *rand_; 143 }; 144