/*  RandomGenerator.cc    Implementation of the RandomGenerator class    Copyright (c) 2005, 2006              by Hyuk Cho    Copyright (c) 2003, 2004    	      by Hyuk Cho, Yuqiang Guan, and Suvrit Sra                {hyukcho, yguan, suvrit}@cs.utexas.edu*/#include <cmath>#include "RandomGenerator.h"#define PI 3.14159265358979323846264338328// Construction/DestructionRandomGenerator::RandomGenerator(){}RandomGenerator::~RandomGenerator(){}// Implementation of base classdouble RandomGenerator::GetUniformPos(){  double x ;  do {    x = GetUniform();  } while (x == 0) ;  return x ;}unsigned long int RandomGenerator::GetUniformInt(unsigned long int n){  unsigned long int offset = min;  unsigned long int range = max - offset;  unsigned long int scale = range / n;  unsigned long int k;  if (n > range)     n = range;  do {    k = (Get() - offset) / scale;  } while (k >= n);  return k;}double RandomGenerator::GetGaussian(const double sigma){#if 1 /* Polar (Box-Mueller) method; See Knuth v2, 3rd ed, p122 */  double x, y, r2;  do {      /* choose x,y in uniform square (-1,-1) to (+1,+1) */    x = -1 + 2 * GetUniform();    y = -1 + 2 * GetUniform();      /* see if it is in the unit circle */    r2 = x * x + y * y;  } while (r2 > 1.0 || r2 == 0);  /* Box-Muller transform */  return sigma * y * sqrt (-2.0 * log (r2) / r2);   /* only one random deviate is produced, theother is discarded. Because saving      the other in a static variable would screw up re-entrant or theaded code. */#endif#if 0 /* Ratio method (Kinderman-Monahan); see Knuth v2, 3rd ed, p130 */      /* K+M, ACM Trans Math Software 3 (1977) 257-260. */  double u, v, x, xx;  do {    v = GetUniform(r);    do {      u = GetUniform(r);    } while (u == 0);      /* Const 1.715... = sqrt(8/e) */    x = 1.71552776992141359295*(v-0.5)/u;  } while (x * x > -4.0 * log(u));  return sigma * x;#endif}double RandomGenerator::GetGaussianPDF(const double x, const double sigma){  double u = x / fabs(sigma) ;  double p = (1 / (sqrt (2 * PI) * fabs(sigma)) ) * exp (-u * u / 2);  return p;}// Implementation of MT_19937 generatorRandomGeneratorMT19937::RandomGeneratorMT19937() : UPPER_MASK(0x80000000UL), LOWER_MASK(0x7fffffffUL){  name = "mt19937";			//??  max = 0xffffffffUL;			/* RAND_MAX  */  min = 0;				/* RAND_MIN  */  size = sizeof(mti) + sizeof(mt);	// ??}RandomGeneratorMT19937::~RandomGeneratorMT19937(){}unsigned long RandomGeneratorMT19937::Get(){  unsigned long k ;#define MAGIC(y) (((y)&0x1) ? 0x9908b0dfUL : 0)  if (mti >= MT_N){	/* generate N words at one time */    int kk;    for (kk = 0; kk < MT_N - MT_M; kk++){      unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);      mt[kk] = mt[kk + MT_M] ^ (y >> 1) ^ MAGIC(y);    }    for (; kk < MT_N - 1; kk++){      unsigned long y = (mt[kk] & UPPER_MASK) | (mt[kk + 1] & LOWER_MASK);      mt[kk] = mt[kk + (MT_M - MT_N)] ^ (y >> 1) ^ MAGIC(y);    }    {      unsigned long y = (mt[MT_N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);      mt[MT_N - 1] = mt[MT_M - 1] ^ (y >> 1) ^ MAGIC(y);    }    mti = 0;  }  /* Tempering */  k = mt[mti];  k ^= (k >> 11);  k ^= (k << 7) & 0x9d2c5680UL;  k ^= (k << 15) & 0xefc60000UL;  k ^= (k >> 18);  mti++;  return k;}double RandomGeneratorMT19937::GetUniform(){  return Get() / 4294967296.0 ;}void RandomGeneratorMT19937::Set(unsigned long int seed){  int i;  if (seed == 0)    seed = 4357;	/* the default seed is 4357 */  mt[0] = seed & 0xffffffffUL;  /* We use the congruence s_{n+1} = (69069*s_n) mod 2^32 to     initialize the state. This works because ANSI-C unsigned long     integer arithmetic is automatically modulo 2^32 (or a higher     power of two), so we can safely ignore overflow. */#define LCG(n) ((69069 * n) & 0xffffffffUL)  for (i = 1; i < MT_N; i++)    mt[i] = LCG (mt[i - 1]);  mti = i;}// Implementation of Taus generatorRandomGeneratorTaus::RandomGeneratorTaus(){  name = "taus";			/* name */  max = 0xffffffffUL;			/* RAND_MAX */  min = 0;			        /* RAND_MIN */  size = sizeof(unsigned long int) * 3;}RandomGeneratorTaus::~RandomGeneratorTaus(){}unsigned long int RandomGeneratorTaus::Get(){#define MASK 0xffffffffUL#define TAUSWORTHE(s,a,b,c,d) (((s &c) <<d) &MASK) ^ ((((s <<a) &MASK)^s) >>b)  s1 = TAUSWORTHE (s1, 13, 19, 4294967294UL, 12);  s2 = TAUSWORTHE (s2, 2, 25, 4294967288UL, 4);  s3 = TAUSWORTHE (s3, 3, 11, 4294967280UL, 17);  return (s1 ^ s2 ^ s3);}double RandomGeneratorTaus::GetUniform(){  return Get() / 4294967296.0 ;}void RandomGeneratorTaus::Set(unsigned long seed){  if (seed == 0)    seed = 1;	/* default seed is 1 */#define LCG(n) ((69069 * n) & 0xffffffffUL)  s1 = LCG (seed);  s2 = LCG (s1);  s3 = LCG (s2);  /* "warm it up" */  Get();  Get();  Get();  Get();  Get();  Get();  return;}// Implementation of TT800 generatorRandomGeneratorTT800::RandomGeneratorTT800(){  name = "tt800";			/* name */  max = 0xffffffffUL;			/* RAND_MAX */  min = 0;			        /* RAND_MIN */  size = sizeof(n) + sizeof(x);}RandomGeneratorTT800::~RandomGeneratorTT800(){}unsigned long int RandomGeneratorTT800::Get(){  /* this is the magic vector, a */  const unsigned long mag01[2] = {0x00000000, 0x8ebfd028UL};  unsigned long int y;  if (n >= TT_N){    int i;    for (i = 0; i < TT_N - TT_M; i++)      x[i] = x[i + TT_M] ^ (x[i] >> 1) ^ mag01[x[i] % 2];    for (; i < TT_N; i++)      x[i] = x[i + (TT_M - TT_N)] ^ (x[i] >> 1) ^ mag01[x[i] % 2];    n = 0;  }  y = x[n];  y ^= (y << 7) & 0x2b5b2500UL;		/* s and b, magic vectors */  y ^= (y << 15) & 0xdb8b0000UL;	/* t and c, magic vectors */  y &= 0xffffffffUL;			/* you may delete this line if word size = 32 */  /* The following line was added by Makoto Matsumoto in the 1996     version to improve lower bit's correlation.  Delete this line     to use the code published in 1994.  */  y ^= (y >> 16);	/* added to the 1994 version */  n = n + 1;  return y;}double RandomGeneratorTT800::GetUniform(){  return Get() / 4294967296.0 ;}void RandomGeneratorTT800::Set(unsigned long int seed){  const int init_n = 0;  const unsigned long int init_x[TT_N] =		{0x95f24dabUL, 0x0b685215UL, 0xe76ccae7UL,		 0xaf3ec239UL, 0x715fad23UL, 0x24a590adUL,		 0x69e4b5efUL, 0xbf456141UL, 0x96bc1b7bUL,		 0xa7bdf825UL, 0xc1de75b7UL, 0x8858a9c9UL,		 0x2da87693UL, 0xb657f9ddUL, 0xffdc8a9fUL,		 0x8121da71UL, 0x8b823ecbUL, 0x885d05f5UL,		 0x4e20cd47UL, 0x5a9ad5d9UL, 0x512c0c03UL,		 0xea857ccdUL, 0x4cc1d30fUL, 0x8891a8a1UL,		 0xa6b7aadbUL};  if (seed == 0){	/* default seed is given explicitly in the original code */    n = init_n;    for (int i=0; i<TT_N; i++)       x[i] = init_x[i];  } else {    n = 0;    x[0] = seed & 0xffffffffUL;    for (int i = 1; i < TT_N; i++)      x[i] = (69069 * x[i - 1]) & 0xffffffffUL;  }  return;}// Implementation of R250 generatorRandomGeneratorR250::RandomGeneratorR250(){  name = "r250";			/* name */  max = 0xffffffffUL;			/* RAND_MAX */  min = 0;			        /* RAND_MIN */  size = sizeof(i) + sizeof(x);}RandomGeneratorR250::~RandomGeneratorR250(){}unsigned long int RandomGeneratorR250::Get(){  unsigned long int k;  int j;  if (i >= 147)    j = i - 147;  else    j = i + 103;  k = x[i] ^ x[j];  x[i] = k;  if (i >= 249)    i = 0;  else    i = i + 1;  return k;}double RandomGeneratorR250::GetUniform(){  return Get() /  4294967296.0 ;}void RandomGeneratorR250::Set(unsigned long int seed){  int j;  if (seed == 0)    seed = 1;	/* default seed is 1 */  i = 0;#define LCG(n) ((69069 * n) & 0xffffffffUL)  for (j = 0; j < 250; j++){	/* Fill the buffer  */    seed = LCG (seed);    x[j] = seed;  }  {    /* Masks for turning on the diagonal bit and turning off the leftmost bits */    unsigned long int msb = 0x80000000UL;    unsigned long int mask = 0xffffffffUL;    for (j = 0; j < 32; j++){      int k = 7 * j + 3;	/* Select a word to operate on        */      x[k] &= mask;	/* Turn off bits left of the diagonal */      x[k] |= msb;	/* Turn on the diagonal bit           */      mask >>= 1;      msb >>= 1;    }  }  return;}