#pragma once

/* Mersenne Twister rng. */
class MTRand_int32 {
public:
  /* The default constructor uses default seed only if this is the first instance. */
  MTRand_int32(void) { seed(5489UL); }
  /* Constructor with 32 bit int as seed. */
  MTRand_int32(unsigned long s) { seed(s); }
  /* Consstructor with array of size 32 bit ints as seed. */
  MTRand_int32(const unsigned long* array, int size) { seed(array, size); }

  void seed(unsigned long); /* Seed with 32 bit integer. */
  void seed(const unsigned long*, int size); /* Seed with array. */
  
  virtual ~MTRand_int32(void) { }
protected: /* Used by derived classes, otherwise not accessible; use the ()-operator. */
  unsigned long rand_int32(void);
private:
  unsigned long operator()() { return rand_int32(); }
  static const int n = 624, m = 397; /* Compile time constants. */

  /* Static: */
  unsigned long state[n]; /* State vector array. */
  int p; /* Position in array state. */

  /* Generate the pseudo random numbers. */
  unsigned long twiddle(unsigned long, unsigned long);
  void gen_state(void);

  /* Doesn't make much sense having a copy or assignment operator laying around. */
  MTRand_int32(const MTRand_int32&);
  void operator=(const MTRand_int32&);
};

/* Inline for speed, must therefore reside in header. */
inline unsigned long MTRand_int32::twiddle(unsigned long u, unsigned long v) {
  return (((u & 0x80000000UL) | (v & 0x7FFFFFFFUL)) >> 1)
    ^ ((v & 1UL) ? 0x9908B0DFUL : 0x0UL);
}

/* Generate 32 bit random int. */
inline unsigned long MTRand_int32::rand_int32(void) {
  if(p == n) gen_state(); /* New state vector needed. */
  /* gen_state() is split off to be non-inline, because it is only called once
   * in every 624 calls, otherwise irand() would become too big to get inlined.
   */
  unsigned long x = state[p++];
  x ^= (x >> 11);
  x ^= (x << 7)  & 0x9D2C5680UL;
  x ^= (x << 15) & 0xEFC60000UL;
  return x ^ (x>>18);
}

class MTRand : public MTRand_int32 {
public:
  MTRand(void) : MTRand_int32() {}
  MTRand(unsigned int long seed) : MTRand_int32(seed) {}
  MTRand(const unsigned long* seed, int size) : MTRand_int32(seed, size) {}
  ~MTRand(void) {}

  double NDouble(int p) {
    double o = Double(1.0);
    while(--p) o *= Double(1.0);
    return o;
  }

  double Double(double min, double max) {
    return Double(max-min)+min;
  }

  /* Interval [0, 1]. */
  double Double(double max) {
    /* Divided by 2^32 */
    return max*static_cast<double>(rand_int32()) * (1./4294967296.);
  }

  /* Interval [0, 1]. */
  double Double(void) {
    /* Divided by 2^32. */
    return static_cast<double>(rand_int32()) * (1./4294967296.);
  }
  
  /* [0, 1] */
  double Double_closed(void) {
    /* Divided by 2^32 - 1. */
    return static_cast<double>(rand_int32()) * (1./4294967295.);
  }

  /* [0, 1]. */
  double Double_open(void) {
    /* Divided by 2^32. */
    return (static_cast<double>(rand_int32()) + .5) * (1./4294967296.);
  }

  /* Generates 53 bit resolution doubles in the half-open interval [0, 1]. */
  double Double53(void) {
    return (static_cast<double>(rand_int32() >> 5) * 67108864. +
        static_cast<double>(rand_int32() >> 6)) * (1./9007199254740992.);
  }
  /* [min,max] */
  unsigned int Int32(int min, int max) {
    return (rand_int32()%(1+max-min))+min;
  }

  /* [0,max] */
  unsigned int Int32(int max) { return rand_int32()%max; }
  unsigned int Int32(void) { return rand_int32(); }
private:
  /* Copy and assignment operators not defined. */
  MTRand(const MTRand&);
  void operator=(const MTRand&);
};