testing/random.h

   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
  20   MTRandom() {
  21     Init(TEST_SEED1);
  22   }
  23
  24   explicit MTRandom(uint32_t seed) {
  25     Init(seed);
  26   }
  27
  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
  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
  68   uint64_t Rand64() {
  69     return ((uint64_t)Rand32() << 32) | Rand32();
  70   }
  71
  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>
  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>
  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:
 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:
 130   MTRandom8(MTRandom *rand)
 131     : rand_(rand) {
 132   }
 133
 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 };