//
// Implementation of a random number generator that draws from a
// gaussian distribution.  The underlying (uniform deviate)
// pseudorandom number generator is the Park and Miller "Minimum
// standard" generator with Bayes-Durham shuffle
//
// Copyright (c) 2001 Wesley H. Huang.  All rights reserved.
//
#include <math.h>
#include "rand.h"

long pmRand();
void pmRandSeed(long s);
long pmRand();
long pmRandShuffle();
float myRand();

//
// returns a number drawn from a gaussian distribution with the
// specified mean and standard deviation
// 
float gaussian(float mean, float stdev)
{
  static bool cached = false;
  static float extra;  	// the extra number from a 0 mean unit stdev gaussian

  if (cached) {
    cached = false;
    return extra*stdev + mean;
  }
  else {
    // pick a random point in the unit circle (excluding the origin)
    float a,b,c;
    do {
      a = 2*myRand()-1.0;
      b = 2*myRand()-1.0;
      c = a*a + b*b;
    }
    while (c >= 1.0 || c == 0.0);
    // transform it into two values drawn from a gaussian
    float t = sqrt(-2*log(c)/c);
    extra = t*a;
    cached = true;
    return t*b*stdev + mean;
  }
}

// default random generator seed (nothing special about this number
// except that it's not 0)
const long defaultSeed = 18846224; 
long pmRandGenSeed = defaultSeed;

// set the seed (optional)
void pmRandSeed(long s)
{
  if (s == 0) // can't use 0 as a seed
    pmRandGenSeed = defaultSeed;	// use the default
  else
    pmRandGenSeed = s;
}
 
// parameters for the random number generator
const long a = 16807;
const long m = 2147483647;
const long q = m/a;		// 127773
const long r = m % a;		// 2836
const long pmRandMax = m;

// Park and Miller's "Minimal Standard" generator using Schrage's
// method to do the computation in 32 bit arithmetic without overflow.
//
// This is a multiplicative congruential generator: I_{j+1} = a I_j mod m
// using the above constants (which are very carefully picked).
//
// The result is an integer between 1 and m-1 inclusive.
//
// Reference: Numerical Recipes in C, 2nd ed, pg 278.
long pmRand()
{

  long k = pmRandGenSeed/q;
  pmRandGenSeed = a*(pmRandGenSeed - k*q) - r*k;
  if (pmRandGenSeed < 0)
    pmRandGenSeed += m;
  return pmRandGenSeed;
}

// the Bayes-Durham shuffle
long pmRandShuffle()
{
  const int size = 32;
  const long div = 1 + (pmRandMax-1)/size;
  static bool initialized = false;
  static long rn[size];
  static long next;
  int index;

  if (!initialized) {
    for (int k= 0; k < size; k++) // warm up
      pmRand();
    for (int k= 0; k < size; k++) // fill table
      rn[k] =pmRand();
    next = pmRand();
    initialized = true;
  }
  index = next/div;		// find out which table entry to use next
  next = rn[index];		// that random nember determines the next index
  rn[index] = pmRand();		// refill the table slot
  return next;			// return the random number
}
  
// 
// This is a simple but pretty good random number generator.  I am
// providing it so that everyone will have the same sequence of
// pseudorandom numbers which will probably be useful for me for
// testing.
//
// Returns a floating point number in (0.0, 1.0)
float myRand()
{
  return pmRandShuffle()/(pmRandMax + 1.0);
}

void myRandSeed(long s)
{
  pmRandSeed(s);
}