/** ------------------------------------------------------------------------ 
* These are C lang routines for generating random variables from eight     
* continuous distributions.                                               
*                                                                        
* There are generators for 8 continuous distributions:                     
*                                                                          
*      Generator        Range (x)    Mean           Variance               
*                                                                          
*      Uniform(a, b)    a < x < b    (a+b)/2        (b-a)*(b-a)/12          
*      Exponential(m)   x > 0        m              m*m                     
*      Erlang(n, b)     x > 0        n*b            n*b*b                   
*      Normal           all x        0              1                       
*      Gauss(m, s)      all x        m              s*s                    
*      Lognormal(a, b)  x > 0           see below                           
*      Chisquare(n)     x > 0        n              2*n                     
*      Student(n)       all x        0  (n > 1)     n/(n - 2)   (n > 2)     
*                                                                           
* For the Lognormal(a, b), the mean and variance are:                       
*                                                                           
*                        mean = Exp(a + 0.5*b*b)                            
*                    variance = (Exp(b*b) - 1)*Exp(2*a + b*b)               
*                                                                           
*                                                                           
* Purpose         : Continuous Generator Routines                           
* Author          : Steve Park                                              
* Latest Revision : 10-24-90                                                
* Converted to C  : David W. Geyer  09-03-91                                
* ------------------------------------------------------------------------- 
*/

/* include files */
#include <math.h>                          /* for log(), exp() */

/** ========================================================================= 
*                             double Uniform()                              
* ========================================================================= 
* Returns a uniformly distributed real between a and b.                     
* NOTE: use a < b                                                           
* ========================================================================= 
*/
double Uniform(double a, double b)
{ double Random();

  return(a + (b - a) * Random());
}  /* Uniform */

/** ========================================================================= 
*                             double Exponential()                          
* ========================================================================= 
* Returns an exponentially distributed positive real.                       
* NOTE: use m > 0                                                           
* ========================================================================= 
*/
double Exponential(double m)
{ double Random();

  return(- m * log(1.0 - Random()));
} /* Exponential */

/** ========================================================================= 
*                             double Erlang()                               
* ========================================================================= 
* Returns an Erlang distributed positive real.                              
* NOTE: use n > 0 and b > 0                                                 
* ========================================================================= 
*/
double Erlang(long n, double b)
{ long   i;
  double x, Exponential(double);

  x = 0.0;
  for (i=0;i < n;i++)
    x += Exponential(b);
  return(x);
} /* Erlang */

/** ========================================================================= 
*                             double Normal()                               
* ========================================================================= 
* Returns a standard normal distributed real.                               
* Uses a very accurate approximation of the Normal idf due to Odeh & Evans, 
* J. Applied Statistics, 1974, vol 23, pp 96-97.                            
* ========================================================================= 
*/
double Normal()
{ double p0 = 0.322232431088;        double q0 = 0.099348462606;
  double p1 = 1.0;                   double q1 = 0.588581570495;
  double p2 = 0.342242088547;        double q2 = 0.531103462366;
  double p3 = 0.204231210245e-1;     double q3 = 0.103537752850;
  double p4 = 0.453642210148e-4;     double q4 = 0.385607006340e-2;
  double u, t, p, q, Random();

  u = Random();
  if (u < 0.5)
    t = sqrt(- 2 * log(u));
  else
    t = sqrt(- 2 * log(1.0 - u));
  p = p0 + t * (p1 + t * (p2 + t * (p3 + t * p4)));
  q = q0 + t * (q1 + t * (q2 + t * (q3 + t * q4)));
  if (u < 0.5)
    return( (p / q) - t );
  else
    return( t - (p / q) );
}  /* Normal */

/** ========================================================================= 
*                             double Gauss()                                
* ========================================================================= 
* Returns a Gaussian distributed real.                                      
* NOTE: use s > 0                                                           
* ========================================================================= 
*/
double Gauss(double m, double s)
{ double Normal();

  return(m + s * Normal());
}  /* Gauss */

/** ======================================================================= 
*                             double Lognormal()                            
* ========================================================================= 
* Returns a lognormal distributed positive real.                            
* NOTE: use b > 0                                                           
* ========================================================================= 
*/
double Lognormal(double a, double b)
{ double Normal();

  return( exp(a + b * Normal()) );
}  /* Lognormal */

/** ======================================================================= 
*                             double Chisquare()                            
* ========================================================================= 
* Returns a chi-square distributed positive real.                           
* NOTE: use n > 0                                                           
* ========================================================================= 
*/
double Chisquare(long n)
{ long    i;
  double  z, x, Normal();

  x = 0.0;
  for(i=0;i < n;i++)  {
    z = Normal();
    x += z * z;
  }
  return(x);
} /* Chisquare */

/** ======================================================================= 
*                             double Student()                              
* ========================================================================= 
* Returns a Student-t distributed real.                                     
* NOTE: use n > 0                                                           
* ========================================================================= 
*/
double Student(long n)
{ double Chisquare(long), Normal();

  return( Normal() / sqrt(Chisquare(n) / (double)n) );
} /* Student */