r = anp1/bnp1;    if(fabs(r-r0) <= *eps*r) goto S20;/*        RESCALE AN, BN, ANP1, AND BNP1*/    an /= bnp1;    bn /= bnp1;    anp1 = r;    bnp1 = 1.0e0;    goto S10;S20:/*                 TERMINATION*/    bfrac *= r;    return bfrac;}void bgrat(double *a,double *b,double *x,double *y,double *w,	   double *eps,int *ierr)/*-----------------------------------------------------------------------     ASYMPTOTIC EXPANSION FOR IX(A,B) WHEN A IS LARGER THAN B.     THE RESULT OF THE EXPANSION IS ADDED TO W. IT IS ASSUMED     THAT A .GE. 15 AND B .LE. 1.  EPS IS THE TOLERANCE USED.     IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.-----------------------------------------------------------------------*/{static double bm1,bp2n,cn,coef,dj,j,l,lnx,n2,nu,p,q,r,s,sum,t,t2,u,v,z;static int i,n,nm1;static double c[30],d[30],T1;/*     ..     .. Executable Statements ..*/    bm1 = *b-0.5e0-0.5e0;    nu = *a+0.5e0*bm1;    if(*y > 0.375e0) goto S10;    T1 = -*y;    lnx = alnrel(&T1);    goto S20;S10:    lnx = log(*x);S20:    z = -(nu*lnx);    if(*b*z == 0.0e0) goto S70;/*                 COMPUTATION OF THE EXPANSION                 SET R = EXP(-Z)*Z**B/GAMMA(B)*/    r = *b*(1.0e0+gam1(b))*exp(*b*log(z));    r *= (exp(*a*lnx)*exp(0.5e0*bm1*lnx));    u = algdiv(b,a)+*b*log(nu);    u = r*exp(-u);    if(u == 0.0e0) goto S70;    grat1(b,&z,&r,&p,&q,eps);    v = 0.25e0*pow(1.0e0/nu,2.0);    t2 = 0.25e0*lnx*lnx;    l = *w/u;    j = q/r;    sum = j;    t = cn = 1.0e0;    n2 = 0.0e0;    for(n=1; n<=30; n++) {        bp2n = *b+n2;        j = (bp2n*(bp2n+1.0e0)*j+(z+bp2n+1.0e0)*t)*v;        n2 += 2.0e0;        t *= t2;        cn /= (n2*(n2+1.0e0));        c[n-1] = cn;        s = 0.0e0;        if(n == 1) goto S40;        nm1 = n-1;        coef = *b-(double)n;        for(i=1; i<=nm1; i++) {            s += (coef*c[i-1]*d[n-i-1]);            coef += *b;        }S40:        d[n-1] = bm1*cn+s/(double)n;        dj = d[n-1]*j;        sum += dj;        if(sum <= 0.0e0) goto S70;        if(fabs(dj) <= *eps*(sum+l)) goto S60;    }S60:/*                    ADD THE RESULTS TO W*/    *ierr = 0;    *w += (u*sum);    return;S70:/*               THE EXPANSION CANNOT BE COMPUTED*/    *ierr = 1;    return;}double bpser(double *a,double *b,double *x,double *eps)/*-----------------------------------------------------------------------     POWER SERIES EXPANSION FOR EVALUATING IX(A,B) WHEN B .LE. 1     OR B*X .LE. 0.7.  EPS IS THE TOLERANCE USED.-----------------------------------------------------------------------*/{static double bpser,a0,apb,b0,c,n,sum,t,tol,u,w,z;static int i,m;/*     ..     .. Executable Statements ..*/    bpser = 0.0e0;    if(*x == 0.0e0) return bpser;/*-----------------------------------------------------------------------            COMPUTE THE FACTOR X**A/(A*BETA(A,B))-----------------------------------------------------------------------*/    a0 = fifdmin1(*a,*b);    if(a0 < 1.0e0) goto S10;    z = *a*log(*x)-betaln(a,b);    bpser = exp(z)/ *a;    goto S100;S10:    b0 = fifdmax1(*a,*b);    if(b0 >= 8.0e0) goto S90;    if(b0 > 1.0e0) goto S40;/*            PROCEDURE FOR A0 .LT. 1 AND B0 .LE. 1*/    bpser = pow(*x,*a);    if(bpser == 0.0e0) return bpser;    apb = *a+*b;    if(apb > 1.0e0) goto S20;    z = 1.0e0+gam1(&apb);    goto S30;S20:    u = *a+*b-1.e0;    z = (1.0e0+gam1(&u))/apb;S30:    c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;    bpser *= (c*(*b/apb));    goto S100;S40:/*         PROCEDURE FOR A0 .LT. 1 AND 1 .LT. B0 .LT. 8*/    u = gamln1(&a0);    m = (long)(b0 - 1.0e0);    if(m < 1) goto S60;    c = 1.0e0;    for(i=1; i<=m; i++) {        b0 -= 1.0e0;        c *= (b0/(a0+b0));    }    u = log(c)+u;S60:    z = *a*log(*x)-u;    b0 -= 1.0e0;    apb = a0+b0;    if(apb > 1.0e0) goto S70;    t = 1.0e0+gam1(&apb);    goto S80;S70:    u = a0+b0-1.e0;    t = (1.0e0+gam1(&u))/apb;S80:    bpser = exp(z)*(a0/ *a)*(1.0e0+gam1(&b0))/t;    goto S100;S90:/*            PROCEDURE FOR A0 .LT. 1 AND B0 .GE. 8*/    u = gamln1(&a0)+algdiv(&a0,&b0);    z = *a*log(*x)-u;    bpser = a0/ *a*exp(z);S100:    if(bpser == 0.0e0 || *a <= 0.1e0**eps) return bpser;/*-----------------------------------------------------------------------                     COMPUTE THE SERIES-----------------------------------------------------------------------*/    sum = n = 0.0e0;    c = 1.0e0;    tol = *eps/ *a;S110:    n += 1.0e0;    c *= ((0.5e0+(0.5e0-*b/n))**x);    w = c/(*a+n);    sum += w;    if(fabs(w) > tol) goto S110;    bpser *= (1.0e0+*a*sum);    return bpser;}void bratio(double *a,double *b,double *x,double *y,double *w,	    double *w1,int *ierr)/*-----------------------------------------------------------------------             EVALUATION OF THE INCOMPLETE BETA FUNCTION IX(A,B)                      --------------------      IT IS ASSUMED THAT A AND B ARE NONNEGATIVE, AND THAT X .LE. 1     AND Y = 1 - X.  BRATIO ASSIGNS W AND W1 THE VALUES                       W  = IX(A,B)                      W1 = 1 - IX(A,B)      IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.     IF NO INPUT ERRORS ARE DETECTED THEN IERR IS SET TO 0 AND     W AND W1 ARE COMPUTED. OTHERWISE, IF AN ERROR IS DETECTED,     THEN W AND W1 ARE ASSIGNED THE VALUE 0 AND IERR IS SET TO     ONE OF THE FOLLOWING VALUES ...         IERR = 1  IF A OR B IS NEGATIVE        IERR = 2  IF A = B = 0        IERR = 3  IF X .LT. 0 OR X .GT. 1        IERR = 4  IF Y .LT. 0 OR Y .GT. 1        IERR = 5  IF X + Y .NE. 1        IERR = 6  IF X = A = 0        IERR = 7  IF Y = B = 0 --------------------     WRITTEN BY ALFRED H. MORRIS, JR.        NAVAL SURFACE WARFARE CENTER        DAHLGREN, VIRGINIA     REVISED ... NOV 1991-----------------------------------------------------------------------*/{static int K1 = 1;static double a0,b0,eps,lambda,t,x0,y0,z;static int ierr1,ind,n;static double T2,T3,T4,T5;/*     ..     .. Executable Statements ..*//*     ****** EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE SMALLEST            FLOATING POINT NUMBER FOR WHICH 1.0 + EPS .GT. 1.0*/    eps = spmpar(&K1);    *w = *w1 = 0.0e0;    if(*a < 0.0e0 || *b < 0.0e0) goto S270;    if(*a == 0.0e0 && *b == 0.0e0) goto S280;    if(*x < 0.0e0 || *x > 1.0e0) goto S290;    if(*y < 0.0e0 || *y > 1.0e0) goto S300;    z = *x+*y-0.5e0-0.5e0;    if(fabs(z) > 3.0e0*eps) goto S310;    *ierr = 0;    if(*x == 0.0e0) goto S210;    if(*y == 0.0e0) goto S230;    if(*a == 0.0e0) goto S240;    if(*b == 0.0e0) goto S220;    eps = fifdmax1(eps,1.e-15);    if(fifdmax1(*a,*b) < 1.e-3*eps) goto S260;    ind = 0;    a0 = *a;    b0 = *b;    x0 = *x;    y0 = *y;    if(fifdmin1(a0,b0) > 1.0e0) goto S40;/*             PROCEDURE FOR A0 .LE. 1 OR B0 .LE. 1*/    if(*x <= 0.5e0) goto S10;    ind = 1;    a0 = *b;    b0 = *a;    x0 = *y;    y0 = *x;S10:    if(b0 < fifdmin1(eps,eps*a0)) goto S90;    if(a0 < fifdmin1(eps,eps*b0) && b0*x0 <= 1.0e0) goto S100;    if(fifdmax1(a0,b0) > 1.0e0) goto S20;    if(a0 >= fifdmin1(0.2e0,b0)) goto S110;    if(pow(x0,a0) <= 0.9e0) goto S110;    if(x0 >= 0.3e0) goto S120;    n = 20;    goto S140;S20:    if(b0 <= 1.0e0) goto S110;    if(x0 >= 0.3e0) goto S120;    if(x0 >= 0.1e0) goto S30;    if(pow(x0*b0,a0) <= 0.7e0) goto S110;S30:    if(b0 > 15.0e0) goto S150;    n = 20;    goto S140;S40:/*             PROCEDURE FOR A0 .GT. 1 AND B0 .GT. 1*/    if(*a > *b) goto S50;    lambda = *a-(*a+*b)**x;    goto S60;S50:    lambda = (*a+*b)**y-*b;S60:    if(lambda >= 0.0e0) goto S70;    ind = 1;    a0 = *b;    b0 = *a;    x0 = *y;    y0 = *x;    lambda = fabs(lambda);S70:    if(b0 < 40.0e0 && b0*x0 <= 0.7e0) goto S110;    if(b0 < 40.0e0) goto S160;    if(a0 > b0) goto S80;    if(a0 <= 100.0e0) goto S130;    if(lambda > 0.03e0*a0) goto S130;    goto S200;S80:    if(b0 <= 100.0e0) goto S130;    if(lambda > 0.03e0*b0) goto S130;    goto S200;S90:/*            EVALUATION OF THE APPROPRIATE ALGORITHM*/    *w = fpser(&a0,&b0,&x0,&eps);    *w1 = 0.5e0+(0.5e0-*w);    goto S250;S100:    *w1 = apser(&a0,&b0,&x0,&eps);    *w = 0.5e0+(0.5e0-*w1);    goto S250;S110:    *w = bpser(&a0,&b0,&x0,&eps);    *w1 = 0.5e0+(0.5e0-*w);    goto S250;S120:    *w1 = bpser(&b0,&a0,&y0,&eps);    *w = 0.5e0+(0.5e0-*w1);    goto S250;S130:    T2 = 15.0e0*eps;    *w = bfrac(&a0,&b0,&x0,&y0,&lambda,&T2);    *w1 = 0.5e0+(0.5e0-*w);    goto S250;S140:    *w1 = bup(&b0,&a0,&y0,&x0,&n,&eps);    b0 += (double)n;S150:    T3 = 15.0e0*eps;    bgrat(&b0,&a0,&y0,&x0,w1,&T3,&ierr1);    *w = 0.5e0+(0.5e0-*w1);    goto S250;S160:    n = (long)(b0);    b0 -= (double)n;    if(b0 != 0.0e0) goto S170;    n -= 1;    b0 = 1.0e0;S170:    *w = bup(&b0,&a0,&y0,&x0,&n,&eps);    if(x0 > 0.7e0) goto S180;    *w += bpser(&a0,&b0,&x0,&eps);    *w1 = 0.5e0+(0.5e0-*w);    goto S250;S180:    if(a0 > 15.0e0) goto S190;    n = 20;    *w += bup(&a0,&b0,&x0,&y0,&n,&eps);    a0 += (double)n;S190:    T4 = 15.0e0*eps;    bgrat(&a0,&b0,&x0,&y0,w,&T4,&ierr1);    *w1 = 0.5e0+(0.5e0-*w);    goto S250;S200:    T5 = 100.0e0*eps;    *w = basym(&a0,&b0,&lambda,&T5);    *w1 = 0.5e0+(0.5e0-*w);    goto S250;S210:/*               TERMINATION OF THE PROCEDURE*/    if(*a == 0.0e0) goto S320;S220:    *w = 0.0e0;    *w1 = 1.0e0;    return;S230:    if(*b == 0.0e0) goto S330;S240:    *w = 1.0e0;    *w1 = 0.0e0;    return;S250:    if(ind == 0) return;    t = *w;    *w = *w1;    *w1 = t;    return;S260:/*           PROCEDURE FOR A AND B .LT. 1.E-3*EPS*/    *w = *b/(*a+*b);    *w1 = *a/(*a+*b);    return;S270:/*                       ERROR RETURN*/    *ierr = 1;    return;S280:    *ierr = 2;    return;S290:    *ierr = 3;    return;S300:    *ierr = 4;    return;S310:    *ierr = 5;    return;S320:    *ierr = 6;    return;S330:    *ierr = 7;    return;}double brcmp1(int *mu,double *a,double *b,double *x,double *y)/*-----------------------------------------------------------------------          EVALUATION OF  EXP(MU) * (X**A*Y**B/BETA(A,B))-----------------------------------------------------------------------*/{static double Const = .398942280401433e0;static double brcmp1,a0,apb,b0,c,e,h,lambda,lnx,lny,t,u,v,x0,y0,z;static int i,n;/*-----------------     CONST = 1/SQRT(2*PI)-----------------*/static double T1,T2,T3,T4;/*     ..     .. Executable Statements ..*/    a0 = fifdmin1(*a,*b);    if(a0 >= 8.0e0) goto S130;    if(*x > 0.375e0) goto S10;    lnx = log(*x);    T1 = -*x;    lny = alnrel(&T1);    goto S30;S10:    if(*y > 0.375e0) goto S20;    T2 = -*y;    lnx = alnrel(&T2);    lny = log(*y);    goto S30;S20:    lnx = log(*x);    lny = log(*y);S30:    z = *a*lnx+*b*lny;    if(a0 < 1.0e0) goto S40;    z -= betaln(a,b);    brcmp1 = esum(mu,&z);    return brcmp1;S40:/*-----------------------------------------------------------------------              PROCEDURE FOR A .LT. 1 OR B .LT. 1-----------------------------------------------------------------------*/    b0 = fifdmax1(*a,*b);    if(b0 >= 8.0e0) goto S120;    if(b0 > 1.0e0) goto S70;/*                   ALGORITHM FOR B0 .LE. 1*/    brcmp1 = esum(mu,&z);    if(brcmp1 == 0.0e0) return brcmp1;    apb = *a+*b;    if(apb > 1.0e0) goto S50;    z = 1.0e0+gam1(&apb);    goto S60;S50:    u = *a+*b-1.e0;    z = (1.0e0+gam1(&u))/apb;S60:    c = (1.0e0+gam1(a))*(1.0e0+gam1(b))/z;    brcmp1 = brcmp1*(a0*c)/(1.0e0+a0/b0);    return brcmp1;S70:/*                ALGORITHM FOR 1 .LT. B0 .LT. 8*/    u = gamln1(&a0);    n = (long)(b0 - 1.0e0);    if(n < 1) goto S90;    c = 1.0e0;    for(i=1; i<=n; i++) {        b0 -= 1.0e0;        c *= (b0/(a0+b0));    }    u = log(c)+u;S90:    z -= u;    b0 -= 1.0e0;    apb = a0+b0;    if(apb > 1.0e0) goto S100;    t = 1.0e0+gam1(&apb);    goto S110;S100:    u = a0+b0-1.e0;    t = (1.0e0+gam1(&u))/apb;S110:    brcmp1 = a0*esum(mu,&z)*(1.0e0+gam1(&b0))/t;    return brcmp1;S120:/*                   ALGORITHM FOR B0 .GE. 8*/    u = gamln1(&a0)+algdiv(&a0,&b0);    T3 = z-u;    brcmp1 = a0*esum(mu,&T3);    return brcmp1;S130:/*-----------------------------------------------------------------------              PROCEDURE FOR A .GE. 8 AND B .GE. 8-----------------------------------------------------------------------*/    if(*a > *b) goto S140;    h = *a/ *b;    x0 = h/(1.0e0+h);    y0 = 1.0e0/(1.0e0+h);    lambda = *a-(*a+*b)**x;    goto S150;S140:    h = *b/ *a;    x0 = 1.0e0/(1.0e0+h);    y0 = h/(1.0e0+h);    lambda = (*a+*b)**y-*b;S150:    e = -(lambda/ *a);    if(fabs(e) > 0.6e0) goto S160;    u = rlog1(&e);    goto S170;S160:    u = e-log(*x/x0);S170:    e = lambda/ *b;    if(fabs(e) > 0.6e0) goto S180;    v = rlog1(&e);    goto S190;