Basic formula for generating simple random distributions

rand.cpp

double gaussian(double mu, double sigma)
{
  //TODO: This is not thread safe!
  static double lastRandValue = RAND_MAX;
  if (lastRandValue != RAND_MAX) {
    double rt = mu + sigma*lastRandValue;
    lastRandValue = RAND_MAX;
    return rt;
  }
  double s = 0;
  double u = 0;
  double v = 0;
  while ((s == 0) || (s >= 1)) {
    u = uniform(-1, 1);
    v = uniform(-1, 1);
    s = u*u + v*v;
  }
  double factor = sqrt(-2*log(s)/s);
  lastRandValue = u*factor;
  return mu + sigma*v*factor;
}

double uniform(double min, double max)
{
  if(max < min) {
    return RAND_MAX;
  }
  return min + (max-min)*(double)rand() / ((double)RAND_MAX + 1);
}

// Returns randomly sampled differential times for a poissionian process.
// Parameter rate is the mean event frequency. The dimensions of the returned
// time are the inverse of the dimensions of rate
double next_poisson_time(double meanPeriod)
{
  // Get a uniformly distributed psuedo-random number between 0 and RAND_MAX.
  // Be careful to ensure that 0 is excluded because we are passing this to log
  double uniform_rand =
    static_cast<double>(rand() + 1.0)/(static_cast<double>(RAND_MAX) + 1);

  // return the inverted poission CDF
  return -log(uniform_rand)*meanPeriod;
}