#include "k.h"
/* ******************************CONTENTS******************************* */
/*      random number generators */
/*      -   rgamma       - Gamma deviate of parameters a,b */
/*                         (Bratley et al. 1987) */
/*           RGS         - a < 1 */
/*           RGKM3       - a >=1 */
/* ********************************************************************* */
/*     Generate X from gamma(na,nb),na=shape,nb=scale */
/*     E(X)=na*nb and Var(X)= na*nb^2 */
/*     pp.324,325 in Bratley, Fox and Schrage */
/*     (1987, Springer-Verlag) */
/*     na > 1 RGKM3 */
/*     0 < na <= 1 RGS */
/* ********************************************************************* */
void rgamma(int *m,double *na,double *nb,double *x,long *idum)
{
    /* System generated locals */
    int na_dim1, na_offset, nb_dim1, nb_offset, x_dim1, x_offset, i_1;

    /* Subroutine */ int s_stop();

    /* Local variables */
    double work[6];
    extern double rgkm3_();
    int i, j, k;
    extern double rgs_();

    /* Parameter adjustments */
    na_dim1 = *m;
    na_offset = na_dim1 + 1;
    na -= na_offset;
    nb_dim1 = *m;
    nb_offset = nb_dim1 + 1;
    nb -= nb_offset;
    x_dim1 = *m;
    x_offset = x_dim1 + 1;
    x -= x_offset;

    work[0] = -1.;
    i_1 = *m;
/*    printf("in rgamma2 %d %f %f %f %f %ld\n",
    *m,na[1],na[2],nb[1],nb[2],*idum); */
    for (i = 1; i <= i_1; ++i) {
	for (j = 1; j <= 2; ++j) {
	    if (na[i + j * na_dim1] > 0. && nb[i + j * nb_dim1] > 0.) {
L10:
		if (na[i + j * na_dim1] <= 1.) {
		    x[i + j * x_dim1] = rgs_(&na[i + j * na_dim1], idum);
		} else {
		    x[i + j * x_dim1] = rgkm3_(&na[i + j * na_dim1], work, &k,
			     idum);
		}
		if (x[i + j * x_dim1] <= DBL_EPSILON) {
		    goto L10;
		}
	    } else {
		printf(" Error in rgamma routine\n");
		exit(1);
	    }
	    x[i + j * x_dim1] *= nb[i + j * nb_dim1];
	}
    }
    return;
} 
/* ******************************************************************** */
/*     Gamma variate with parameter  0 < alpha <= 1 */
/*      alp  = Gamma distribution parameter */
/*      idum = random number seed */
/*      runif = uniform random number generator */
/*     References: Ahrens,J.H. and U. Dieter (1972), Computer */
/*      methods for sampling from Gamma, Beta, Poisson and */
/*      Binomial distributions, computing, vol. 12, pp. 223-246. */
/*      Tadikamalla, P.R. and M.E. Johnson (1981), A complete */
/*      guide to Gamma variate generation, amer. j. of math. and */
/*      man. sci., vol. 1, pp. 213-236. */
/* ******************************************************************** */
double rgs_(double *alp,long *idum)
{
    /* System generated locals */
    double ret_val;

    /* Builtin functions */
    double log(), exp();

    /* Local variables */
    static double b, p, x, u1, u2;
    extern double runif();
L100:
    u1 = runif(idum);
    b = (*alp + 2.718281828) / 2.718281828;
    p = b * u1;
    u2 = runif(idum);
    if (p > 1.) {
	goto L300;
    }
    x = exp(log(p) / *alp);
    if (u2 > exp(-x)) {
	goto L100;
    }
    ret_val = x;
    return ret_val;
L300:
    x = -log((b - p) / *alp);
    if (log(u2) > (*alp - 1.) * log(x)) {
	goto L100;
    }
    ret_val = x;
    return ret_val;
} 
/* ******************************************************************** */
/*     Gamma variate with parameter alpha > 1 */
/*     inputs: */
/*      alp    = distribution parameter */
/*      work() = vector of work cells; on first call work(1)=-1. */
/*       work() must be preserved by caller between calls. */
/*      k = work cell. K must be preserved by caller between calls. */
/*      idum = random number seed */
/*      runif = uniform random number generator */
/*     References: Cheng,R.C. and G.M. Feast (1979), Some */
/*      simple Gamma variate generators, applied stat. vol 28, */
/*      pp. 290-295. Tadikamalla, P.R. and M.E. Johnson (1981), */
/*      A complete guide to Gamma variate generation, amer. j. */
/*      of math. and man. sci., vol. 1, pp. 213-236. */
/* ******************************************************************* */
double rgkm3_(double *alp,double *work,int *k,long *idum)
{
    /* System generated locals */
    double ret_val;

    /* Builtin functions */
    double sqrt(), log();

    /* Local variables */
    static double u, w, u1, u2;
    extern double runif();

    /* Parameter adjustments */
    --work;

    if (work[1] == *alp) {
	goto L100;
    }
    work[1] = *alp;
    *k = 1;
    if (*alp > 2.5) {
	*k = 2;
    }
    work[2] = *alp - 1.;
    work[3] = (*alp - 1. / (*alp * 6.)) / work[2];
    work[4] = 2. / work[2];
    work[5] = work[4] + 2.;
    switch (*k) {
	case 1:  goto L1;
	case 2:  goto L11;
    }
L1:
    u1 = runif(idum);
    u2 = runif(idum);
    goto L20;
L11:
    work[6] = sqrt(*alp);
L15:
    u1 = runif(idum);
    u = runif(idum);
    u2 = u1 + (1. - u * 1.86) / work[6];
    if (u2 <= 0. || u2 >= 1.) {
	goto L15;
    }
L20:
    w = work[3] * u1 / u2;
    if (work[4] * u2 - work[5] + w + 1. / w <= 0.) {
	goto L200;
    }
    if (work[4] * log(u2) - log(w) + w - 1. < 0.) {
	goto L200;
    }
L100:
    switch (*k) {
	case 1:  goto L1;
	case 2:  goto L15;
    }
L200:
    ret_val = work[2] * w;
    return ret_val;
}