#include <stdio.h>#include <stdlib.h>#include <math.h>#include "cdflib.h"/* * A comment about ints and longs - whether ints or longs are used should * make no difference, but where double r-values are assigned to ints the * r-value is cast converted to a long, which is then assigned to the int * to be compatible with the operation of fifidint. *//*-----------------------------------------------------------------------      COMPUTATION OF LN(GAMMA(B)/GAMMA(A+B)) WHEN B .GE. 8                          --------      IN THIS ALGORITHM, DEL(X) IS THE FUNCTION DEFINED BY     LN(GAMMA(X)) = (X - 0.5)*LN(X) - X + 0.5*LN(2*PI) + DEL(X). -----------------------------------------------------------------------*/double algdiv(double *a,double *b){static double c0 = .833333333333333e-01;static double c1 = -.277777777760991e-02;static double c2 = .793650666825390e-03;static double c3 = -.595202931351870e-03;static double c4 = .837308034031215e-03;static double c5 = -.165322962780713e-02;static double algdiv,c,d,h,s11,s3,s5,s7,s9,t,u,v,w,x,x2,T1;/*     ..     .. Executable Statements ..*/    if(*a <= *b) goto S10;    h = *b/ *a;    c = 1.0e0/(1.0e0+h);    x = h/(1.0e0+h);    d = *a+(*b-0.5e0);    goto S20;S10:    h = *a/ *b;    c = h/(1.0e0+h);    x = 1.0e0/(1.0e0+h);    d = *b+(*a-0.5e0);S20:/*                SET SN = (1 - X**N)/(1 - X)*/    x2 = x*x;    s3 = 1.0e0+(x+x2);    s5 = 1.0e0+(x+x2*s3);    s7 = 1.0e0+(x+x2*s5);    s9 = 1.0e0+(x+x2*s7);    s11 = 1.0e0+(x+x2*s9);/*                SET W = DEL(B) - DEL(A + B)*/    t = pow(1.0e0/ *b,2.0);    w = ((((c5*s11*t+c4*s9)*t+c3*s7)*t+c2*s5)*t+c1*s3)*t+c0;    w *= (c/ *b);/*                    COMBINE THE RESULTS*/    T1 = *a/ *b;    u = d*alnrel(&T1);    v = *a*(log(*b)-1.0e0);    if(u <= v) goto S30;    algdiv = w-v-u;    return algdiv;S30:    algdiv = w-u-v;    return algdiv;}double alngam(double *x)/***********************************************************************      double alngam(double *x)                 double precision LN of the GAMma function                                Function       Returns the natural logarithm of GAMMA(X).                                Arguments       X --> value at which scaled log gamma is to be returned                    X is DOUBLE PRECISION                                Method       If X .le. 6.0, then use recursion to get X below 3     then apply rational approximation number 5236 of     Hart et al, Computer Approximations, John Wiley and     Sons, NY, 1968.      If X .gt. 6.0, then use recursion to get X to at least 12 and     then use formula 5423 of the same source. ***********************************************************************/{#define hln2pi 0.91893853320467274178e0static double coef[5] = {    0.83333333333333023564e-1,-0.27777777768818808e-2,0.79365006754279e-3,    -0.594997310889e-3,0.8065880899e-3};static double scoefd[4] = {    0.62003838007126989331e2,0.9822521104713994894e1,-0.8906016659497461257e1,    0.1000000000000000000e1};static double scoefn[9] = {    0.62003838007127258804e2,0.36036772530024836321e2,0.20782472531792126786e2,    0.6338067999387272343e1,0.215994312846059073e1,0.3980671310203570498e0,    0.1093115956710439502e0,0.92381945590275995e-2,0.29737866448101651e-2};static int K1 = 9;static int K3 = 4;static int K5 = 5;static double alngam,offset,prod,xx;static int i,n;static double T2,T4,T6;/*     ..     .. Executable Statements ..*/    if(!(*x <= 6.0e0)) goto S70;    prod = 1.0e0;    xx = *x;    if(!(*x > 3.0e0)) goto S30;S10:    if(!(xx > 3.0e0)) goto S20;    xx -= 1.0e0;    prod *= xx;    goto S10;S30:S20:    if(!(*x < 2.0e0)) goto S60;S40:    if(!(xx < 2.0e0)) goto S50;    prod /= xx;    xx += 1.0e0;    goto S40;S60:S50:    T2 = xx-2.0e0;    T4 = xx-2.0e0;    alngam = devlpl(scoefn,&K1,&T2)/devlpl(scoefd,&K3,&T4);/*     COMPUTE RATIONAL APPROXIMATION TO GAMMA(X)*/    alngam *= prod;    alngam = log(alngam);    goto S110;S70:    offset = hln2pi;/*     IF NECESSARY MAKE X AT LEAST 12 AND CARRY CORRECTION IN OFFSET*/    n = fifidint(12.0e0-*x);    if(!(n > 0)) goto S90;    prod = 1.0e0;    for(i=1; i<=n; i++) prod *= (*x+(double)(i-1));    offset -= log(prod);    xx = *x+(double)n;    goto S100;S90:    xx = *x;S100:/*     COMPUTE POWER SERIES*/    T6 = 1.0e0/pow(xx,2.0);    alngam = devlpl(coef,&K5,&T6)/xx;    alngam += (offset+(xx-0.5e0)*log(xx)-xx);S110:    return alngam;#undef hln2pi}double alnrel(double *a)/*-----------------------------------------------------------------------            EVALUATION OF THE FUNCTION LN(1 + A)-----------------------------------------------------------------------*/{static double p1 = -.129418923021993e+01;static double p2 = .405303492862024e+00;static double p3 = -.178874546012214e-01;static double q1 = -.162752256355323e+01;static double q2 = .747811014037616e+00;static double q3 = -.845104217945565e-01;static double alnrel,t,t2,w,x;/*     ..     .. Executable Statements ..*/    if(fabs(*a) > 0.375e0) goto S10;    t = *a/(*a+2.0e0);    t2 = t*t;    w = (((p3*t2+p2)*t2+p1)*t2+1.0e0)/(((q3*t2+q2)*t2+q1)*t2+1.0e0);    alnrel = 2.0e0*t*w;    return alnrel;S10:    x = 1.e0+*a;    alnrel = log(x);    return alnrel;}double apser(double *a,double *b,double *x,double *eps)/*-----------------------------------------------------------------------     APSER YIELDS THE INCOMPLETE BETA RATIO I(SUB(1-X))(B,A) FOR     A .LE. MIN(EPS,EPS*B), B*X .LE. 1, AND X .LE. 0.5. USED WHEN     A IS VERY SMALL. USE ONLY IF ABOVE INEQUALITIES ARE SATISFIED.-----------------------------------------------------------------------*/{static double g = .577215664901533e0;static double apser,aj,bx,c,j,s,t,tol;/*     ..     .. Executable Statements ..*/    bx = *b**x;    t = *x-bx;    if(*b**eps > 2.e-2) goto S10;    c = log(*x)+psi(b)+g+t;    goto S20;S10:    c = log(bx)+g+t;S20:    tol = 5.0e0**eps*fabs(c);    j = 1.0e0;    s = 0.0e0;S30:    j += 1.0e0;    t *= (*x-bx/j);    aj = t/j;    s += aj;    if(fabs(aj) > tol) goto S30;    apser = -(*a*(c+s));    return apser;}double basym(double *a,double *b,double *lambda,double *eps)/*-----------------------------------------------------------------------     ASYMPTOTIC EXPANSION FOR IX(A,B) FOR LARGE A AND B.     LAMBDA = (A + B)*Y - B  AND EPS IS THE TOLERANCE USED.     IT IS ASSUMED THAT LAMBDA IS NONNEGATIVE AND THAT     A AND B ARE GREATER THAN OR EQUAL TO 15.-----------------------------------------------------------------------*/{static double e0 = 1.12837916709551e0;static double e1 = .353553390593274e0;static int num = 20;/*------------------------     ****** NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP            ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN.            THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1.------------------------     E0 = 2/SQRT(PI)     E1 = 2**(-3/2)------------------------*/static int K3 = 1;static double basym,bsum,dsum,f,h,h2,hn,j0,j1,r,r0,r1,s,sum,t,t0,t1,u,w,w0,z,z0,    z2,zn,znm1;static int i,im1,imj,j,m,mm1,mmj,n,np1;static double a0[21],b0[21],c[21],d[21],T1,T2;/*     ..     .. Executable Statements ..*/    basym = 0.0e0;    if(*a >= *b) goto S10;    h = *a/ *b;    r0 = 1.0e0/(1.0e0+h);    r1 = (*b-*a)/ *b;    w0 = 1.0e0/sqrt(*a*(1.0e0+h));    goto S20;S10:    h = *b/ *a;    r0 = 1.0e0/(1.0e0+h);    r1 = (*b-*a)/ *a;    w0 = 1.0e0/sqrt(*b*(1.0e0+h));S20:    T1 = -(*lambda/ *a);    T2 = *lambda/ *b;    f = *a*rlog1(&T1)+*b*rlog1(&T2);    t = exp(-f);    if(t == 0.0e0) return basym;    z0 = sqrt(f);    z = 0.5e0*(z0/e1);    z2 = f+f;    a0[0] = 2.0e0/3.0e0*r1;    c[0] = -(0.5e0*a0[0]);    d[0] = -c[0];    j0 = 0.5e0/e0*erfc1(&K3,&z0);    j1 = e1;    sum = j0+d[0]*w0*j1;    s = 1.0e0;    h2 = h*h;    hn = 1.0e0;    w = w0;    znm1 = z;    zn = z2;    for(n=2; n<=num; n+=2) {        hn = h2*hn;        a0[n-1] = 2.0e0*r0*(1.0e0+h*hn)/((double)n+2.0e0);        np1 = n+1;        s += hn;        a0[np1-1] = 2.0e0*r1*s/((double)n+3.0e0);        for(i=n; i<=np1; i++) {            r = -(0.5e0*((double)i+1.0e0));            b0[0] = r*a0[0];            for(m=2; m<=i; m++) {                bsum = 0.0e0;                mm1 = m-1;                for(j=1; j<=mm1; j++) {                    mmj = m-j;                    bsum += (((double)j*r-(double)mmj)*a0[j-1]*b0[mmj-1]);                }                b0[m-1] = r*a0[m-1]+bsum/(double)m;            }            c[i-1] = b0[i-1]/((double)i+1.0e0);            dsum = 0.0e0;            im1 = i-1;            for(j=1; j<=im1; j++) {                imj = i-j;                dsum += (d[imj-1]*c[j-1]);            }            d[i-1] = -(dsum+c[i-1]);        }        j0 = e1*znm1+((double)n-1.0e0)*j0;        j1 = e1*zn+(double)n*j1;        znm1 = z2*znm1;        zn = z2*zn;        w = w0*w;        t0 = d[n-1]*w*j0;        w = w0*w;        t1 = d[np1-1]*w*j1;        sum += (t0+t1);        if(fabs(t0)+fabs(t1) <= *eps*sum) goto S80;    }S80:    u = exp(-bcorr(a,b));    basym = e0*t*u*sum;    return basym;}double bcorr(double *a0,double *b0)/*-----------------------------------------------------------------------      EVALUATION OF  DEL(A0) + DEL(B0) - DEL(A0 + B0)  WHERE     LN(GAMMA(A)) = (A - 0.5)*LN(A) - A + 0.5*LN(2*PI) + DEL(A).     IT IS ASSUMED THAT A0 .GE. 8 AND B0 .GE. 8. -----------------------------------------------------------------------*/{static double c0 = .833333333333333e-01;static double c1 = -.277777777760991e-02;static double c2 = .793650666825390e-03;static double c3 = -.595202931351870e-03;static double c4 = .837308034031215e-03;static double c5 = -.165322962780713e-02;static double bcorr,a,b,c,h,s11,s3,s5,s7,s9,t,w,x,x2;/*     ..     .. Executable Statements ..*/    a = fifdmin1(*a0,*b0);    b = fifdmax1(*a0,*b0);    h = a/b;    c = h/(1.0e0+h);    x = 1.0e0/(1.0e0+h);    x2 = x*x;/*                SET SN = (1 - X**N)/(1 - X)*/    s3 = 1.0e0+(x+x2);    s5 = 1.0e0+(x+x2*s3);    s7 = 1.0e0+(x+x2*s5);    s9 = 1.0e0+(x+x2*s7);    s11 = 1.0e0+(x+x2*s9);/*                SET W = DEL(B) - DEL(A + B)*/    t = pow(1.0e0/b,2.0);    w = ((((c5*s11*t+c4*s9)*t+c3*s7)*t+c2*s5)*t+c1*s3)*t+c0;    w *= (c/b);/*                   COMPUTE  DEL(A) + W*/    t = pow(1.0e0/a,2.0);    bcorr = (((((c5*t+c4)*t+c3)*t+c2)*t+c1)*t+c0)/a+w;    return bcorr;}double betaln(double *a0,double *b0)/*-----------------------------------------------------------------------     EVALUATION OF THE LOGARITHM OF THE BETA FUNCTION-----------------------------------------------------------------------     E = 0.5*LN(2*PI)--------------------------*/{static double e = .918938533204673e0;static double betaln,a,b,c,h,u,v,w,z;static int i,n;static double T1;/*     ..     .. Executable Statements ..*/    a = fifdmin1(*a0,*b0);    b = fifdmax1(*a0,*b0);    if(a >= 8.0e0) goto S100;    if(a >= 1.0e0) goto S20;/*-----------------------------------------------------------------------                   PROCEDURE WHEN A .LT. 1-----------------------------------------------------------------------*/    if(b >= 8.0e0) goto S10;    T1 = a+b;    betaln = gamln(&a)+(gamln(&b)-gamln(&T1));    return betaln;S10:    betaln = gamln(&a)+algdiv(&a,&b);    return betaln;S20:/*-----------------------------------------------------------------------                PROCEDURE WHEN 1 .LE. A .LT. 8-----------------------------------------------------------------------*/    if(a > 2.0e0) goto S40;    if(b > 2.0e0) goto S30;    betaln = gamln(&a)+gamln(&b)-gsumln(&a,&b);    return betaln;S30:    w = 0.0e0;    if(b < 8.0e0) goto S60;    betaln = gamln(&a)+algdiv(&a,&b);    return betaln;S40:/*                REDUCTION OF A WHEN B .LE. 1000*/    if(b > 1000.0e0) goto S80;    n = (long)(a - 1.0e0);    w = 1.0e0;    for(i=1; i<=n; i++) {        a -= 1.0e0;        h = a/b;        w *= (h/(1.0e0+h));    }    w = log(w);    if(b < 8.0e0) goto S60;    betaln = w+gamln(&a)+algdiv(&a,&b);    return betaln;S60:/*                 REDUCTION OF B WHEN B .LT. 8*/    n = (long)(b - 1.0e0);    z = 1.0e0;    for(i=1; i<=n; i++) {        b -= 1.0e0;        z *= (b/(a+b));    }    betaln = w+log(z)+(gamln(&a)+(gamln(&b)-gsumln(&a,&b)));    return betaln;S80:/*                REDUCTION OF A WHEN B .GT. 1000*/    n = (long)(a - 1.0e0);    w = 1.0e0;    for(i=1; i<=n; i++) {        a -= 1.0e0;        w *= (a/(1.0e0+a/b));    }    betaln = log(w)-(double)n*log(b)+(gamln(&a)+algdiv(&a,&b));    return betaln;S100:/*-----------------------------------------------------------------------                   PROCEDURE WHEN A .GE. 8-----------------------------------------------------------------------*/    w = bcorr(&a,&b);    h = a/b;    c = h/(1.0e0+h);    u = -((a-0.5e0)*log(c));    v = b*alnrel(&h);    if(u <= v) goto S110;    betaln = -(0.5e0*log(b))+e+w-v-u;    return betaln;S110:    betaln = -(0.5e0*log(b))+e+w-u-v;    return betaln;}double bfrac(double *a,double *b,double *x,double *y,double *lambda,	     double *eps)/*-----------------------------------------------------------------------     CONTINUED FRACTION EXPANSION FOR IX(A,B) WHEN A,B .GT. 1.     IT IS ASSUMED THAT  LAMBDA = (A + B)*Y - B.-----------------------------------------------------------------------*/{static double bfrac,alpha,an,anp1,beta,bn,bnp1,c,c0,c1,e,n,p,r,r0,s,t,w,yp1;/*     ..     .. Executable Statements ..*/    bfrac = brcomp(a,b,x,y);    if(bfrac == 0.0e0) return bfrac;    c = 1.0e0+*lambda;    c0 = *b/ *a;    c1 = 1.0e0+1.0e0/ *a;    yp1 = *y+1.0e0;    n = 0.0e0;    p = 1.0e0;    s = *a+1.0e0;    an = 0.0e0;    bn = anp1 = 1.0e0;    bnp1 = c/c1;    r = c1/c;S10:/*        CONTINUED FRACTION CALCULATION*/    n += 1.0e0;    t = n/ *a;    w = n*(*b-n)**x;    e = *a/s;    alpha = p*(p+c0)*e*e*(w**x);    e = (1.0e0+t)/(c1+t+t);    beta = n+w/s+e*(c+n*yp1);    p = 1.0e0+t;    s += 2.0e0;/*        UPDATE AN, BN, ANP1, AND BNP1*/    t = alpha*an+beta*anp1;    an = anp1;    anp1 = t;    t = alpha*bn+beta*bnp1;    bn = bnp1;    bnp1 = t;    r0 = r;