/* udrnrt.c uniformly disributed random numbers in range 0.0 to 1.0 */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

double udrnrt(void)
{
  static double modulus = 2147483647.0; /* 2**31-1 */
  static double multiplier = 16807.0; /* prime to prime power */
  static double state = 30000000.0; /* arbitrary initial value */
  double t;
  int i;

  t = state * multiplier;
  i = t / modulus;
  state = t - (i * modulus);
  if(state < 0.0) state = state + modulus;
  return (double)state / (double)modulus;
}

double gauss(double mean, double sigma) /* Gaussian distribution, about 
                                           mean, with sigma = std deviation */
{
  int i;
  double a = 0.0;
  for(i=0; i<12; i++) a += udrnrt();
  return (a-6.0)*sigma+mean;  /* max 3 sigma about mean */
}

double gausse(double mean, double sigma) /* using logs */
{
  static int which = 0;
  static double y1, y2;
  double w, x1, x2;
  if(which==0)
  {
    do
    {
      x1 = 2.0*udrnrt()-1.0;
      x2 = 2.0*udrnrt()-1.0;
      w = x1*x1 + x2*x2;
    } while (w>=1.0);
    w = sqrt((-2.0*log(w))/w);
    y1 = x1*w*sigma+mean;
    y2 = x2*w*sigma+mean;
    which = 1;
    return y1;
  }
  else
  {
    which = 0;
    return y2;
  }
}