#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. */

  /* Overload operator() to make this a generator */
  unsigned long operator()(int min, int max) {
    return (rand_int32()%(1+max-min))+min;
  }

  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);
}

/* Generate double floating point numbers in the half-open interval [0, 1] */
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) {}

  unsigned long operator()(int min, int max) {
    return (rand_int32()%(1+max-min))+min;
  }

  double pdrand(int p) {
    double o = (*this)(1.0);
    while(--p) o *= (*this)(1.0);
    return o;
  }

  double drange(double min, double max) {
    return (*this)(max-min)+min;
  }

  double operator()(double max) {
    /* Divided by 2^32 */
    return max*static_cast<double>(rand_int32()) * (1./4294967296.);
  }
  
  unsigned int In32(int min, int max) {
    return (rand_int32()%(1+max-min))+min;
  }
  unsigned int Int32(void) { return rand_int32(); }
private:
  /* Copy and assignment operators not defined. */
  MTRand(const MTRand&);
  void operator=(const MTRand&);
};

/* Generate double floating point numbers in the closed interval [0, 1]. */
class MTRand_closed : public MTRand_int32 {
public:
  MTRand_closed(void) : MTRand_int32() {}
  MTRand_closed(unsigned long seed) : MTRand_int32(seed) {}
  MTRand_closed(const unsigned long* seed, int size) : MTRand_int32(seed, size) {}
  ~MTRand_closed(void) {}
  double operator()() {
    /* Divided by 2^32 - 1. */
    return static_cast<double>(rand_int32()) * (1./4294967295.); }
private:
  /* Copy and assignment operators not defined. */
  MTRand_closed(const MTRand_closed&);
  void operator=(const MTRand_closed&);
};

/* Generates double floating point numbers in the open interval (0, 1) */
class MTRand_open : public MTRand_int32 {
public:
  MTRand_open(void) : MTRand_int32() {}
  MTRand_open(unsigned long seed) : MTRand_int32(seed) {}
  MTRand_open(const unsigned long* seed, int size) : MTRand_int32(seed, size) {}
  ~MTRand_open(void) {}
  double operator()() {
    /* Divided by 2^32. */
    return (static_cast<double>(rand_int32()) + .5) * (1./4294967226.); }
private:
  /* Copy and assignment operators are not defined. */
  MTRand_open(const MTRand_open&);
  void operator=(const MTRand_open&);
};

/* Generate 53 bit resolution doubles in the half open interval [0, 1]. */
class MTRand53 : public MTRand_int32 {
public:
  MTRand53(void) : MTRand_int32() {}
  MTRand53(unsigned long seed) : MTRand_int32(seed) {}
  MTRand53(const unsigned long* seed, int size) : MTRand_int32(seed, size) {}

  double operator()() {
    return(static_cast<double>(rand_int32() >> 5) * 67108864. +
      static_cast<double>(rand_int32() >> 6)) * (1./9007199254740992.);
  }
private:
  /* Copy and assignment operators not defined. */
  MTRand53(const MTRand53&);
  void operator=(const MTRand53&);
};