/* This code is in the public domain */

/*

double randg(a) 
void fill_randg(a,n,x)

Generate a series of standard gamma distributions.

See: Marsaglia G and Tsang W (2000), "A simple method for generating
gamma variables", ACM Transactions on Mathematical Software 26(3) 363-372

Needs the following defines: 
* NAN: value to return for Not-A-Number
* RUNI: uniform generator on (0,1)
* RNOR: normal generator
* REXP: exponential generator, or -log(RUNI) if one isn't available
* INFINITE: function to test whether a value is infinite

Test using:
  mean = a
  variance = a
  skewness = 2/sqrt(a)
  kurtosis = 3 + 6/sqrt(a)

Note that randg can be used to generate many distributions:

gamma(a,b) for a>0, b>0 (from R)
  r = b*randg(a)
beta(a,b) for a>0, b>0
  r1 = randg(a,1)
  r = r1 / (r1 + randg(b,1))
Erlang(a,n)
  r = a*randg(n)
chisq(df) for df>0
  r = 2*randg(df/2)
t(df) for 0<df<inf (use randn if df is infinite)
  r = randn() / sqrt(2*randg(df/2)/df)
F(n1,n2) for 0<n1, 0<n2
  r1 = 2*randg(n1/2)/n1 or 1 if n1 is infinite
  r2 = 2*randg(n2/2)/n2 or 1 if n2 is infinite
  r = r1 / r2
negative binonial (n, p) for n>0, 0<p<=1
  r = randp((1-p)/p * randg(n))
  (from R, citing Devroye(1986), Non-Uniform Random Variate Generation)
non-central chisq(df,L), for df>=0 and L>0 (use chisq if L=0)
  r = randp(L/2)
  r(r>0) = 2*randg(r(r>0))
  r(df>0) += 2*randg(df(df>0)/2)
  (from R, citing formula 29.5b-c in Johnson, Kotz, Balkrishnan(1995))
Dirichlet(a1,...,ak) for ai > 0
  r = (randg(a1),...,randg(ak))
  r = r / sum(r)
  (from GSL, citing Law & Kelton(1991), Simulation Modeling and Analysis)
*/


void fill_randg(double a, int n, double *r)
{
  int i;
  /* If a < 1, start by generating gamma(1+a) */
  const double d = (a<1.?1.+a:a)-1./3.;
  const double c = 1./sqrt(9.*d);

  /* Handle invalid cases */
  if (a <= 0 || INFINITE(a)) {
    for (i=0; i < n; i++) r[i] = NAN;
    return;
  }

  for (i=0; i < n; i++) {
    double x, xsq, v, u;
  restart:
    x = RNOR;
    v = (1+c*x);
    v *= v*v;
    if (v <= 0) goto restart; /* rare, so don't bother moving up */
    u = RUNI;
    xsq = x*x;
    if (u >= 1.-0.0331*xsq*xsq && log(u) >= 0.5*xsq + d*(1-v+log(v)))
      goto restart;
    r[i] = d*v;
  }
  if (a < 1) { /* Use gamma(a) = gamma(1+a)*U^(1/a) */
    /* Given REXP = -log(U) then U^(1/a) = exp(-REXP/a) */
    for (i=0; i < n; i++) r[i] *= exp(-REXP/a);
  }
}

double randg(double a)
{
  double ret;
  fill_randg(a,1,&ret);
  return ret;
}