#include <math.h>
#include <stdlib.h>
#include <float.h>

/*
*  Generates a uniformly distributed r.v. between 0 and 1.
*  Kris Popat  6/85
*  Ref: Applied Statistics, 1982 vol 31 no 2 pp 188-190
*  Based on FORTRAN routine by H. Malvar.
*/

static long ix = 101;
static long iy = 1001;
static long iz = 10001;

double Rand(void)
{
  static float u;
  
  ix = 171*(ix % 177)-2*(ix/177);
  iy = 172*(iy % 176)-2*(iy/176);
  iz = 170*(iz % 178)-2*(iz/178);
  
  if (ix<0) ix = ix + 30269;
  if (iy<0) iy = iy + 30307;
  if (iz<0) iz = iz + 30323;
  
  u = ((float) ix)/30269.0 +
                ((float) iy)/30307.0 + ((float) iz)/30323.0;
  u -= (float)(int)u;
  return(u);
}

/* Returns a sample from Gamma(a, 1).
 * For Gamma(a,b), scale the result by b.
 * This is from BUGS.
 */
double GammaRand(double a)
{
  double c1, c2, c3, c4, c5, c6;
  if(isnan(a)) return a;
  if(a < DBL_EPSILON) return 0;
  /* If a == 1, then gamma is exponential. */
  /* It is important that Rand() is never zero. */
  if(fabs(a - 1) < DBL_EPSILON) return -log(Rand());
  if(a < 1) {
    c6 = exp(log(Rand())/a);
    a = a + 1;
  }
  else c6 = 1;
  c1 = a - 1;
  c2 = (a - 1/(6*a))/c1;
  c3 = 2/c1;
  c4 = c3 + 2;
  c5 = 1/sqrt(a);
  c6 *= c1;
  while(1) {
    double u1, u2, w;
    /* loop until u1 is valid */
    do {
      u1 = Rand();
      u2 = Rand();
      if(a > 2.5) u1 = u2 + c5*(1 - 1.86*u1);
    } while((u1 <= 0) || (u1 >= 1));
    w = c2*u2/u1;
    if((c3*u1 + w + 1/w <= c4) ||
       (c3*log(u1) - log(w) + w < 1))
      return c6*w;
  }
}