#include "randlib.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)***********************************************************************/{/* JJV changed expmax (log(1.0E38)==87.49823), and added minlog */#define expmax 87.49823#define infnty 1.0E38#define minlog 1.0E-37static float olda = -1.0E37;static float oldb = -1.0E37;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 < minlog || bb < minlog)) goto S10;    fputs(" AA or BB < 1.0E-37 in GENBET - Abort!\n",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));/* JJV altered this */    if(v > expmax) goto S55;/* * JJV added checker to see if a*exp(v) will overflow * JJV S50 _was_ w = a*exp(v); also note here a > 1.0 */    w = exp(v);    if(w > infnty/a) goto S55;    w *= a;    goto S60;S55:    w = infnty;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 * *    JJV added checker to see if log(alpha/(b+w)) will  *    JJV overflow.  If so, we count the log as -INF, and *    JJV consequently evaluate conditional as true, i.e. *    JJV the algorithm rejects the trial and starts over *    JJV May not need this here since alpha > 2.0 */    if(alpha/(b+w) < minlog) goto S40;    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));/* *    JJV instead of checking v > expmax at top, I will check *    JJV if a < 1, then check the appropriate values */    if(a > 1.0) goto S135;/*   JJV a < 1 so it can help out if exp(v) would overflow */    if(v > expmax) goto S132;    w = a*exp(v);    goto S200;S132:    w = v + log(a);    if(w > expmax) goto S140;    w = exp(w);    goto S200;S135:/*   JJV in this case a > 1 */    if(v > expmax) goto S140;    w = exp(v);    if(w > infnty/a) goto S140;    w *= a;    goto S200;S140:    w = infnty;    goto S200;/* * JJV old code *    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));/*   JJV same kind of checking as above */    if(a > 1.0) goto S175;/* JJV a < 1 so it can help out if exp(v) would overflow */    if(v > expmax) goto S172;    w = a*exp(v);    goto S190;S172:    w = v + log(a);    if(w > expmax) goto S180;    w = exp(w);    goto S190;S175:/* JJV in this case a > 1.0 */    if(v > expmax) goto S180;    w = exp(v);    if(w > infnty/a) goto S180;    w *= a;    goto S190;S180:    w = infnty;/* *   JJV old code *    if(!(v > expmax)) goto S180; *    w = infnty; *    goto S190; *S180: *    w = a*exp(v); */S190:/* * JJV here we also check to see if log overlows; if so, we treat it * JJV as -INF, which means condition is true, i.e. restart */    if(alpha/(b+w) < minlog) goto S120;    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#undef minlog}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\n",stderr);    fprintf(stderr," Value of DF: %16.6E\n",df);    exit(1);S10:/* * JJV changed the code to call SGAMMA directly *    genchi = 2.0*gengam(1.0,df/2.0); <- OLD */    genchi = 2.0*sgamma(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.        JJV (av >= 0)                              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;/* JJV added check that av >= 0 */    if(av >= 0.0) goto S10;    fputs(" AV < 0 in GENEXP - ABORT\n",stderr);    fprintf(stderr," Value of AV: %16.6E\n",av);    exit(1);S10:    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!\n",stderr);    fprintf(stderr," DFN value: %16.6E DFD value: %16.6E\n",dfn,dfd);    exit(1);S10:/* * JJV changed this to call SGAMMA directly * *     GENF = ( GENCHI( DFN ) / DFN ) / ( GENCHI( DFD ) / DFD ) *   xnum = genchi(dfn)/dfn; <- OLD *   xden = genchi(dfd)/dfd; <- OLD */    xnum = 2.0*sgamma(dfn/2.0)/dfn;    xden = 2.0*sgamma(dfd/2.0)/dfd;/* * JJV changed constant to prevent underflow at compile time. *   if(!(xden <= 9.999999999998E-39*xnum)) goto S20; */    if(!(xden <= 1.0E-37*xnum)) goto S20;    fputs(" GENF - generated numbers would cause overflow\n",stderr);    fprintf(stderr," Numerator %16.6E Denominator %16.6E\n",xnum,xden);/* * JJV changed next 2 lines to reflect constant change above in the * JJV truncated value returned. *   fputs(" GENF returning 1.0E38\n",stderr); *   genf = 1.0E38; */    fputs(" GENF returning 1.0E37\n",stderr);    genf = 1.0E37;    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     JJV   (a > 0)     r --> Shape parameter of Gamma distribution     JJV   (r > 0)                              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     JJV altered following to reflect argument ranges               (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;/* JJV added argument checker */    if(a > 0.0 && r > 0.0) goto S10;    fputs(" A or R nonpositive in GENGAM - abort!\n",stderr);    fprintf(stderr," A value: %16.6E R value: %16.6E\n",a,r);    exit(1);S10:    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