testing/random.h

/* -*- Mode: C++ -*-  */
/* This is public-domain Mersenne Twister code,
 * attributed to Michael Brundage.  Thanks!
 * http://www.qbrundage.com/michaelb/pubs/essays/random_number_generation.html
 */
#undef MT_LEN
#undef MT_IA
class MTRandom {
 public:
  enum Constants { 
    MT_LEN = 624,
    MT_IA = 397
  };

  static const uint32_t TEST_SEED1;
  static const uint32_t UPPER_MASK;
  static const uint32_t LOWER_MASK;
  static const uint32_t MATRIX_A;

  MTRandom() {
    Init(TEST_SEED1);
  }

  MTRandom(uint32_t seed) {
    Init(seed);
  }

  uint32_t Rand32 () {
    uint32_t y;
    static unsigned long mag01[2] = { 
      0 , MATRIX_A
    };

    if (mt_index_ >= MT_LEN) {
      int kk;

      for (kk = 0; kk < MT_LEN - MT_IA; kk++) {
	y = (mt_buffer_[kk] & UPPER_MASK) | (mt_buffer_[kk + 1] & LOWER_MASK);
	mt_buffer_[kk] = mt_buffer_[kk + MT_IA] ^ (y >> 1) ^ mag01[y & 0x1UL];
      }
      for (;kk < MT_LEN - 1; kk++) {
	y = (mt_buffer_[kk] & UPPER_MASK) | (mt_buffer_[kk + 1] & LOWER_MASK);
	mt_buffer_[kk] = mt_buffer_[kk + (MT_IA - MT_LEN)] ^ (y >> 1) ^ mag01[y & 0x1UL];
      }
      y = (mt_buffer_[MT_LEN - 1] & UPPER_MASK) | (mt_buffer_[0] & LOWER_MASK);
      mt_buffer_[MT_LEN - 1] = mt_buffer_[MT_IA - 1] ^ (y >> 1) ^ mag01[y & 0x1UL];

      mt_index_ = 0;
    }
  
    y = mt_buffer_[mt_index_++];

    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680UL;
    y ^= (y << 15) & 0xefc60000UL;
    y ^= (y >> 18);

    return y;
  }

  uint32_t ExpRand32(uint32_t mean) {
    double mean_d = mean;
    double erand  = log (1.0 / (Rand32() / (double)UINT32_MAX));
    uint32_t x = (uint32_t) (mean_d * erand + 0.5);
    return x;
  }

  uint64_t Rand64() {
    return ((uint64_t)Rand32() << 32) | Rand32();
  }

  uint64_t ExpRand64(uint64_t mean) {
    double mean_d = mean;
    double erand  = log (1.0 / (Rand64() / (double)UINT32_MAX));
    uint64_t x = (uint64_t) (mean_d * erand + 0.5);
    return x;
  }

  template <typename T>
  T Rand() {
    switch (sizeof(T)) {
    case sizeof(uint32_t):
      return Rand32();
    case sizeof(uint64_t):
      return Rand64();
    default:
      cerr << "Invalid sizeof T" << endl;
      abort();
    }
  }

  template <typename T>
  T ExpRand(T mean) {
    switch (sizeof(T)) {
    case sizeof(uint32_t):
      return ExpRand32(mean);
    case sizeof(uint64_t):
      return ExpRand64(mean);
    default:
      cerr << "Invalid sizeof T" << endl;
      abort();
    }
  }

 private:
  void Init(uint32_t seed) {
    mt_buffer_[0] = seed;
    mt_index_ = MT_LEN;
    for (int i = 1; i < MT_LEN; i++) {
      /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
      /* In the previous versions, MSBs of the seed affect   */
      /* only MSBs of the array mt[].                        */
      /* 2002/01/09 modified by Makoto Matsumoto             */
      mt_buffer_[i] = 
	(1812433253UL * (mt_buffer_[i-1] ^ (mt_buffer_[i-1] >> 30)) + i);
    }
  }

  int mt_index_;
  uint32_t mt_buffer_[MT_LEN];
};

const uint32_t MTRandom::TEST_SEED1 = 5489UL;
const uint32_t MTRandom::UPPER_MASK = 0x80000000;
const uint32_t MTRandom::LOWER_MASK = 0x7FFFFFFF;
const uint32_t MTRandom::MATRIX_A = 0x9908B0DF;

class MTRandom8 {
public:
  MTRandom8(MTRandom *rand)
    : rand_(rand) {
  }

  uint8_t Rand8() {
    uint32_t r = rand_->Rand32();

    // TODO: make this use a single byte at a time?
    return (r & 0xff) ^ (r >> 7) ^ (r >> 15) ^ (r >> 21);
  }

private:
  MTRandom *rand_;
};