#include "ranlib.h"#include <stdio.h>#include <math.h>#include <stdlib.h>#define ABS(x) ((x) >= 0 ? (x) : -(x))#define min(a,b) ((a) <= (b) ? (a) : (b))#define max(a,b) ((a) >= (b) ? (a) : (b))void ftnstop(char*);float genbet(float aa,float bb)/***********************************************************************     float genbet(float aa,float bb)               GeNerate BETa random deviate                              Function     Returns a single random deviate from the beta distribution with     parameters A and B.  The density of the beta is               x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1                              Arguments     aa --> First parameter of the beta distribution            bb --> Second parameter of the beta distribution                                     Method     R. C. H. Cheng     Generating Beta Variatew with Nonintegral Shape Parameters     Communications of the ACM, 21:317-322  (1978)     (Algorithms BB and BC)***********************************************************************/{#define expmax 89.0#define infnty 1.0E38static float olda = -1.0;static float oldb = -1.0;static float genbet,a,alpha,b,beta,delta,gamma,k1,k2,r,s,t,u1,u2,v,w,y,z;static long qsame;    qsame = olda == aa && oldb == bb;    if(qsame) goto S20;    if(!(aa <= 0.0 || bb <= 0.0)) goto S10;    fputs(" AA or BB <= 0 in GENBET - Abort!",stderr);    fprintf(stderr," AA: %16.6E BB %16.6E\n",aa,bb);    exit(1);S10:    olda = aa;    oldb = bb;S20:    if(!(min(aa,bb) > 1.0)) goto S100;/*     Alborithm BB     Initialize*/    if(qsame) goto S30;    a = min(aa,bb);    b = max(aa,bb);    alpha = a+b;    beta = sqrt((alpha-2.0)/(2.0*a*b-alpha));    gamma = a+1.0/beta;S30:S40:    u1 = ranf();/*     Step 1*/    u2 = ranf();    v = beta*log(u1/(1.0-u1));    if(!(v > expmax)) goto S50;    w = infnty;    goto S60;S50:    w = a*exp(v);S60:    z = pow(u1,2.0)*u2;    r = gamma*v-1.3862944;    s = a+r-w;/*     Step 2*/    if(s+2.609438 >= 5.0*z) goto S70;/*     Step 3*/    t = log(z);    if(s > t) goto S70;/*     Step 4*/    if(r+alpha*log(alpha/(b+w)) < t) goto S40;S70:/*     Step 5*/    if(!(aa == a)) goto S80;    genbet = w/(b+w);    goto S90;S80:    genbet = b/(b+w);S90:    goto S230;S100:/*     Algorithm BC     Initialize*/    if(qsame) goto S110;    a = max(aa,bb);    b = min(aa,bb);    alpha = a+b;    beta = 1.0/b;    delta = 1.0+a-b;    k1 = delta*(1.38889E-2+4.16667E-2*b)/(a*beta-0.777778);    k2 = 0.25+(0.5+0.25/delta)*b;S110:S120:    u1 = ranf();/*     Step 1*/    u2 = ranf();    if(u1 >= 0.5) goto S130;/*     Step 2*/    y = u1*u2;    z = u1*y;    if(0.25*u2+z-y >= k1) goto S120;    goto S170;S130:/*     Step 3*/    z = pow(u1,2.0)*u2;    if(!(z <= 0.25)) goto S160;    v = beta*log(u1/(1.0-u1));    if(!(v > expmax)) goto S140;    w = infnty;    goto S150;S140:    w = a*exp(v);S150:    goto S200;S160:    if(z >= k2) goto S120;S170:/*     Step 4     Step 5*/    v = beta*log(u1/(1.0-u1));    if(!(v > expmax)) goto S180;    w = infnty;    goto S190;S180:    w = a*exp(v);S190:    if(alpha*(log(alpha/(b+w))+v)-1.3862944 < log(z)) goto S120;S200:/*     Step 6*/    if(!(a == aa)) goto S210;    genbet = w/(b+w);    goto S220;S210:    genbet = b/(b+w);S230:S220:    return genbet;#undef expmax#undef infnty}float genchi(float df)/***********************************************************************     float genchi(float df)                Generate random value of CHIsquare variable                              Function     Generates random deviate from the distribution of a chisquare     with DF degrees of freedom random variable.                              Arguments     df --> Degrees of freedom of the chisquare            (Must be positive)                                     Method     Uses relation between chisquare and gamma.***********************************************************************/{static float genchi;    if(!(df <= 0.0)) goto S10;    fputs("DF <= 0 in GENCHI - ABORT",stderr);    fprintf(stderr,"Value of DF: %16.6E\n",df);    exit(1);S10:    genchi = 2.0*gengam(1.0,df/2.0);    return genchi;}float genexp(float av)/***********************************************************************     float genexp(float av)                    GENerate EXPonential random deviate                              Function     Generates a single random deviate from an exponential     distribution with mean AV.                              Arguments     av --> The mean of the exponential distribution from which            a random deviate is to be generated.                              Method     Renames SEXPO from TOMS as slightly modified by BWB to use RANF     instead of SUNIF.     For details see:               Ahrens, J.H. and Dieter, U.               Computer Methods for Sampling From the               Exponential and Normal Distributions.               Comm. ACM, 15,10 (Oct. 1972), 873 - 882.***********************************************************************/{static float genexp;    genexp = sexpo()*av;    return genexp;}float genf(float dfn,float dfd)/***********************************************************************     float genf(float dfn,float dfd)                GENerate random deviate from the F distribution                              Function     Generates a random deviate from the F (variance ratio)     distribution with DFN degrees of freedom in the numerator     and DFD degrees of freedom in the denominator.                              Arguments     dfn --> Numerator degrees of freedom             (Must be positive)     dfd --> Denominator degrees of freedom             (Must be positive)                              Method     Directly generates ratio of chisquare variates***********************************************************************/{static float genf,xden,xnum;    if(!(dfn <= 0.0 || dfd <= 0.0)) goto S10;    fputs("Degrees of freedom nonpositive in GENF - abort!",stderr);    fprintf(stderr,"DFN value: %16.6EDFD value: %16.6E\n",dfn,dfd);    exit(1);S10:    xnum = genchi(dfn)/dfn;/*      GENF = ( GENCHI( DFN ) / DFN ) / ( GENCHI( DFD ) / DFD )*/    xden = genchi(dfd)/dfd;    if(!(xden <= 9.999999999998E-39*xnum)) goto S20;    fputs(" GENF - generated numbers would cause overflow",stderr);    fprintf(stderr," Numerator %16.6E Denominator %16.6E\n",xnum,xden);    fputs(" GENF returning 1.0E38",stderr);    genf = 1.0E38;    goto S30;S20:    genf = xnum/xden;S30:    return genf;}float gengam(float a,float r)/***********************************************************************     float gengam(float a,float r)           GENerates random deviates from GAMma distribution                              Function     Generates random deviates from the gamma distribution whose     density is          (A**R)/Gamma(R) * X**(R-1) * Exp(-A*X)                              Arguments     a --> Location parameter of Gamma distribution     r --> Shape parameter of Gamma distribution                              Method     Renames SGAMMA from TOMS as slightly modified by BWB to use RANF     instead of SUNIF.     For details see:               (Case R >= 1.0)               Ahrens, J.H. and Dieter, U.               Generating Gamma Variates by a               Modified Rejection Technique.               Comm. ACM, 25,1 (Jan. 1982), 47 - 54.     Algorithm GD               (Case 0.0 <= R <= 1.0)               Ahrens, J.H. and Dieter, U.               Computer Methods for Sampling from Gamma,               Beta, Poisson and Binomial Distributions.               Computing, 12 (1974), 223-246/     Adapted algorithm GS.***********************************************************************/{static float gengam;    gengam = sgamma(r);    gengam /= a;    return gengam;}void genmn(float *parm,float *x,float *work)/***********************************************************************     void genmn(float *parm,float *x,float *work)              GENerate Multivariate Normal random deviate                              Arguments     parm --> Parameters needed to generate multivariate normal               deviates (MEANV and Cholesky decomposition of               COVM). Set by a previous call to SETGMN.               1 : 1                - size of deviate, P               2 : P + 1            - mean vector               P+2 : P*(P+3)/2 + 1  - upper half of cholesky                                       decomposition of cov matrix     x    <-- Vector deviate generated.     work <--> Scratch array                              Method     1) Generate P independent standard normal deviates - Ei ~ N(0,1)     2) Using Cholesky decomposition find A s.t. trans(A)*A = COVM     3) trans(A)E + MEANV ~ N(MEANV,COVM)***********************************************************************/{static long i,icount,j,p,D1,D2,D3,D4;static float ae;    p = (long) (*parm);/*     Generate P independent normal deviates - WORK ~ N(0,1)*/    for(i=1; i<=p; i++) *(work+i-1) = snorm();    for(i=1,D3=1,D4=(p-i+D3)/D3; D4>0; D4--,i+=D3) {/*     PARM (P+2 : P*(P+3)/2 + 1) contains A, the Cholesky      decomposition of the desired covariance matrix.          trans(A)(1,1) = PARM(P+2)          trans(A)(2,1) = PARM(P+3)          trans(A)(2,2) = PARM(P+2+P)          trans(A)(3,1) = PARM(P+4)          trans(A)(3,2) = PARM(P+3+P)          trans(A)(3,3) = PARM(P+2-1+2P)  ...     trans(A)*WORK + MEANV ~ N(MEANV,COVM)*/        icount = 0;        ae = 0.0;        for(j=1,D1=1,D2=(i-j+D1)/D1; D2>0; D2--,j+=D1) {            icount += (j-1);            ae += (*(parm+i+(j-1)*p-icount+p)**(work+j-1));        }        *(x+i-1) = ae+*(parm+i);    }}void genmul(long n,float *p,long ncat,long *ix)/***********************************************************************      void genmul(int n,float *p,int ncat,int *ix)     GENerate an observation from the MULtinomial distribution                              Arguments     N --> Number of events that will be classified into one of           the categories 1..NCAT     P --> Vector of probabilities.  P(i) is the probability that           an event will be classified into category i.  Thus, P(i)           must be [0,1]. Only the first NCAT-1 P(i) must be defined           since P(NCAT) is 1.0 minus the sum of the first           NCAT-1 P(i).     NCAT --> Number of categories.  Length of P and IX.     IX <-- Observation from multinomial distribution.  All IX(i)            will be nonnegative and their sum will be N.                              Method     Algorithm from page 559 of      Devroye, Luc      Non-Uniform Random Variate Generation.  Springer-Verlag,     New York, 1986. ***********************************************************************/{static float prob,ptot,sum;static long i,icat,ntot;    if(n < 0) ftnstop("N < 0 in GENMUL");    if(ncat <= 1) ftnstop("NCAT <= 1 in GENMUL");    ptot = 0.0F;    for(i=0; i<ncat-1; i++) {        if(*(p+i) < 0.0F) ftnstop("Some P(i) < 0 in GENMUL");        if(*(p+i) > 1.0F) ftnstop("Some P(i) > 1 in GENMUL");        ptot += *(p+i);    }    if(ptot > 0.99999F) ftnstop("Sum of P(i) > 1 in GENMUL");/*     Initialize variables*/    ntot = n;    sum = 1.0F;    for(i=0; i<ncat; i++) ix[i] = 0;/*     Generate the observation*/    for(icat=0; icat<ncat-1; icat++) {        prob = *(p+icat)/sum;        *(ix+icat) = ignbin(ntot,prob);        ntot -= *(ix+icat);	if(ntot <= 0) return;        sum -= *(p+icat);    }    *(ix+ncat-1) = ntot;/*     Finished*/    return;}float gennch(float df,float xnonc)/***********************************************************************     float gennch(float df,float xnonc)           Generate random value of Noncentral CHIsquare variable                              Function     Generates random deviate  from the  distribution  of a  noncentral     chisquare with DF degrees  of freedom and noncentrality  parameter     xnonc.                              Arguments     df --> Degrees of freedom of the chisquare            (Must be > 1.0)     xnonc --> Noncentrality parameter of the chisquare               (Must be >= 0.0)                              Method     Uses fact that  noncentral chisquare  is  the  sum of a  chisquare     deviate with DF-1  degrees of freedom plus the  square of a normal     deviate with mean XNONC and standard deviation 1.***********************************************************************/{static float gennch;    if(!(df <= 1.0 || xnonc < 0.0)) goto S10;    fputs("DF <= 1 or XNONC < 0 in GENNCH - ABORT",stderr);    fprintf(stderr,"Value of DF: %16.6E Value of XNONC%16.6E\n",df,xnonc);    exit(1);S10:    gennch = genchi(df-1.0)+pow(gennor(sqrt(xnonc),1.0),2.0);    return gennch;}float gennf(float dfn,float dfd,float xnonc)/***********************************************************************     float gennf(float dfn,float dfd,float xnonc)           GENerate random deviate from the Noncentral F distribution                              Function     Generates a random deviate from the  noncentral F (variance ratio)     distribution with DFN degrees of freedom in the numerator, and DFD     degrees of freedom in the denominator, and noncentrality parameter     XNONC.                              Arguments     dfn --> Numerator degrees of freedom             (Must be >= 1.0)     dfd --> Denominator degrees of freedom             (Must be positive)     xnonc --> Noncentrality parameter               (Must be nonnegative)                              Method     Directly generates ratio of noncentral numerator chisquare variate     to central denominator chisquare variate.***********************************************************************/{static float gennf,xden,xnum;static long qcond;    qcond = dfn <= 1.0 || dfd <= 0.0 || xnonc < 0.0;    if(!qcond) goto S10;