# include <cstdlib>
# include <iostream>
# include <iomanip>
# include <cmath>
# include <ctime>

using namespace std;

# include "dcdflib.H"

//****************************************************************************80

double algdiv ( double *a, double *b )

//****************************************************************************80
//
//  Purpose:
// 
//    ALGDIV computes ln ( Gamma ( B ) / Gamma ( A + B ) ) when 8 <= B.
//
//  Discussion:
//
//    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).
//
//  Parameters:
//
//    Input, double *A, *B, define the arguments.
//
//    Output, double ALGDIV, the value of ln(Gamma(B)/Gamma(A+B)).
//
{
  static double algdiv;
  static double c;
  static double c0 =  0.833333333333333e-01;
  static double c1 = -0.277777777760991e-02;
  static double c2 =  0.793650666825390e-03;
  static double c3 = -0.595202931351870e-03;
  static double c4 =  0.837308034031215e-03;
  static double c5 = -0.165322962780713e-02;
  static double d;
  static double h;
  static double s11;
  static double s3;
  static double s5;
  static double s7;
  static double s9;
  static double t;
  static double T1;
  static double u;
  static double v;
  static double w;
  static double x;
  static double x2;

  if ( *b <= *a )
  {
    h = *b / *a;
    c = 1.0e0 / ( 1.0e0 + h );
    x = h / ( 1.0e0 + h );
    d = *a + ( *b - 0.5e0 );
  }
  else
  {
    h = *a / *b;
    c = h / ( 1.0e0 + h );
    x = 1.0e0 / ( 1.0e0 + h );
    d = *b + ( *a - 0.5e0 );
  }
//
//  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 ( v < u )
  {
    algdiv = w - v - u;
  }
  else
  {
    algdiv = w - u - v;
  }
  return algdiv;
}
//****************************************************************************80

double alnrel ( double *a )

//****************************************************************************80
//
//  Purpose:
// 
//    ALNREL evaluates the function ln ( 1 + A ).
//
//  Modified:
//
//    17 November 2006
//
//  Reference:
//
//    Armido DiDinato, Alfred Morris,
//    Algorithm 708: 
//    Significant Digit Computation of the Incomplete Beta Function Ratios,
//    ACM Transactions on Mathematical Software,
//    Volume 18, 1993, pages 360-373.
//
//  Parameters:
//
//    Input, double *A, the argument.
//
//    Output, double ALNREL, the value of ln ( 1 + A ).
//
{
  double alnrel;
  static double p1 = -0.129418923021993e+01;
  static double p2 =  0.405303492862024e+00;
  static double p3 = -0.178874546012214e-01;
  static double q1 = -0.162752256355323e+01;
  static double q2 =  0.747811014037616e+00;
  static double q3 = -0.845104217945565e-01;
  double t;
  double t2;
  double w;
  double x;

  if ( fabs ( *a ) <= 0.375e0 )
  {
    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;
  }
  else
  {
    x = 1.0e0 + *a;
    alnrel = log ( x );
  }
  return alnrel;
}
//****************************************************************************80

double apser ( double *a, double *b, double *x, double *eps )

//****************************************************************************80
//
//  Purpose:
// 
//    APSER computes the incomplete beta ratio I(SUB(1-X))(B,A).
//
//  Discussion:
//
//    APSER is used only for cases where
//
//      A <= min ( EPS, EPS * B ), 
//      B * X <= 1, and 
//      X <= 0.5.
//
//  Parameters:
//
//    Input, double *A, *B, *X, the parameters of teh
//    incomplete beta ratio.
//
//    Input, double *EPS, a tolerance.
//
//    Output, double APSER, the computed value of the
//    incomplete beta ratio.
//
{
  static double g = 0.577215664901533e0;
  static double apser,aj,bx,c,j,s,t,tol;

    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 = j + 1.0e0;
    t = t * (*x-bx/j);
    aj = t/j;
    s = s + aj;
    if(fabs(aj) > tol) goto S30;
    apser = -(*a*(c+s));
    return apser;
}
//****************************************************************************80

double bcorr ( double *a0, double *b0 )

//****************************************************************************80
//
//  Purpose:
// 
//    BCORR evaluates DEL(A0) + DEL(B0) - DEL(A0 + B0).
//
//  Discussion:
//
//    The function DEL(A) is a remainder term that is used in the expression:
//
//      ln ( Gamma ( A ) ) = ( A - 0.5 ) * ln ( A ) 
//        - A + 0.5 * ln ( 2 * PI ) + DEL ( A ),
//
//    or, in other words, DEL ( A ) is defined as:
//
//      DEL ( A ) = ln ( Gamma ( A ) ) - ( A - 0.5 ) * ln ( A ) 
//        + A + 0.5 * ln ( 2 * PI ).
//
//  Parameters:
//
//    Input, double *A0, *B0, the arguments.
//    It is assumed that 8 <= A0 and 8 <= B0.
//
//    Output, double *BCORR, the value of the function.
//
{
  static double c0 =  0.833333333333333e-01;
  static double c1 = -0.277777777760991e-02;
  static double c2 =  0.793650666825390e-03;
  static double c3 = -0.595202931351870e-03;
  static double c4 =  0.837308034031215e-03;
  static double c5 = -0.165322962780713e-02;
  static double bcorr,a,b,c,h,s11,s3,s5,s7,s9,t,w,x,x2;

  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;
}
//****************************************************************************80**

double beta ( double a, double b )

//****************************************************************************80**
//
//  Purpose:
//
//    BETA evaluates the beta function.
//
//  Modified:
//
//    03 December 1999
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double A, B, the arguments of the beta function.
//
//    Output, double BETA, the value of the beta function.
//
{
  return ( exp ( beta_log ( &a, &b ) ) );
}
//****************************************************************************80

double beta_asym ( double *a, double *b, double *lambda, double *eps )

//****************************************************************************80
//
//  Purpose:
//
//    BETA_ASYM computes an asymptotic expansion for IX(A,B), for large A and B.
//
//  Parameters:
//
//    Input, double *A, *B, the parameters of the function.
//    A and B should be nonnegative.  It is assumed that both A and B
//    are greater than or equal to 15.
//
//    Input, double *LAMBDA, the value of ( A + B ) * Y - B.
//    It is assumed that 0 <= LAMBDA.
//
//    Input, double *EPS, the tolerance.
//
{
  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 value;
  static double 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;

    value = 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 value;
    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 * error_fc ( &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));
    value = e0*t*u*sum;
    return value;
}
//****************************************************************************80

double beta_frac ( double *a, double *b, double *x, double *y, double *lambda,
  double *eps )

//****************************************************************************80
//
//  Purpose:
// 
//    BETA_FRAC evaluates a continued fraction expansion for IX(A,B).
//
//  Parameters:
//
//    Input, double *A, *B, the parameters of the function.
//    A and B should be nonnegative.  It is assumed that both A and
//    B are greater than 1.
//
//    Input, double *X, *Y.  X is the argument of the
//    function, and should satisy 0 <= X <= 1.  Y should equal 1 - X.
//
//    Input, double *LAMBDA, the value of ( A + B ) * Y - B.
//
//    Input, double *EPS, a tolerance.
//
//    Output, double BETA_FRAC, the value of the continued
//    fraction approximation for IX(A,B).
//
{
  static double bfrac,alpha,an,anp1,beta,bn,bnp1,c,c0,c1,e,n,p,r,r0,s,t,w,yp1;

  bfrac = beta_rcomp ( 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;
//
//  CONTINUED FRACTION CALCULATION
//
S10:
  n = 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;
  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;
//
//  TERMINATION
//
S20:
  bfrac = bfrac * r;
  return bfrac;
}
//****************************************************************************80

void beta_grat ( double *a, double *b, double *x, double *y, double *w,
  double *eps,int *ierr )

//****************************************************************************80
//
//  Purpose:
// 
//    BETA_GRAT evaluates an asymptotic expansion for IX(A,B).
//
//  Parameters:
//
//    Input, double *A, *B, the parameters of the function.
//    A and B should be nonnegative.  It is assumed that 15 <= A 
//    and B <= 1, and that B is less than A.
//
//    Input, double *X, *Y.  X is the argument of the
//    function, and should satisy 0 <= X <= 1.  Y should equal 1 - X.
//
//    Input/output, double *W, a quantity to which the
//    result of the computation is to be added on output.
//
//    Input, double *EPS, a tolerance.
//
//    Output, int *IERR, an error flag, which is 0 if no error
//    was detected.
//
{
  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;

    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;
    gamma_rat1 ( 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 = 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 = s + (coef*c[i-1]*d[n-i-1]);
            coef = coef + *b;
        }
S40:
        d[n-1] = bm1*cn+s/(double)n;
        dj = d[n-1]*j;
        sum = 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 = *w + (u*sum);
    return;
S70:
//
//  THE EXPANSION CANNOT BE COMPUTED
//
    *ierr = 1;
    return;
}
//****************************************************************************80

void beta_inc ( double *a, double *b, double *x, double *y, double *w,
  double *w1, int *ierr )

//****************************************************************************80
//
//  Purpose:
// 
//    BETA_INC evaluates the incomplete beta function IX(A,B).
//
//  Author:
//
//    Alfred H Morris, Jr,
//    Naval Surface Weapons Center,
//    Dahlgren, Virginia.
//
//  Parameters:
//
//    Input, double *A, *B, the parameters of the function.
//    A and B should be nonnegative.
//
//    Input, double *X, *Y.  X is the argument of the
//    function, and should satisy 0 <= X <= 1.  Y should equal 1 - X.
//
//    Output, double *W, *W1, the values of IX(A,B) and
//    1-IX(A,B).
//
//    Output, int *IERR, the error flag.
//    0, no error was detected.
//    1, A or B is negative;
//    2, A = B = 0;
//    3, X < 0 or 1 < X;
//    4, Y < 0 or 1 < Y;
//    5, X + Y /= 1;
//    6, X = A = 0;
//    7, Y = B = 0.
//
{
  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;
//
//  EPS IS A MACHINE DEPENDENT CONSTANT. EPS IS THE SMALLEST
//  NUMBER FOR WHICH 1.0 + EPS .GT. 1.0
//
    eps = dpmpar ( &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 = beta_pser(&a0,&b0,&x0,&eps);
    *w1 = 0.5e0+(0.5e0-*w);
    goto S250;
S120:
    *w1 = beta_pser(&b0,&a0,&y0,&eps);
    *w = 0.5e0+(0.5e0-*w1);
    goto S250;
S130:
    T2 = 15.0e0*eps;
    *w = beta_frac ( &a0,&b0,&x0,&y0,&lambda,&T2 );
    *w1 = 0.5e0+(0.5e0-*w);
    goto S250;
S140:
    *w1 = beta_up ( &b0, &a0, &y0, &x0, &n, &eps );
    b0 = b0 + (double)n;
S150:
    T3 = 15.0e0*eps;
    beta_grat (&b0,&a0,&y0,&x0,w1,&T3,&ierr1);
    *w = 0.5e0+(0.5e0-*w1);
    goto S250;
S160:
    n = ( int ) b0;
    b0 -= (double)n;
    if(b0 != 0.0e0) goto S170;
    n -= 1;
    b0 = 1.0e0;
S170:
    *w = beta_up ( &b0, &a0, &y0, &x0, &n, &eps );
    if(x0 > 0.7e0) goto S180;
    *w = *w + beta_pser(&a0,&b0,&x0,&eps);
    *w1 = 0.5e0+(0.5e0-*w);
    goto S250;
S180:
    if(a0 > 15.0e0) goto S190;
    n = 20;
    *w = *w + beta_up ( &a0, &b0, &x0, &y0, &n, &eps );
    a0 = a0 + (double)n;
S190:
    T4 = 15.0e0*eps;
    beta_grat ( &a0, &b0, &x0, &y0, w, &T4, &ierr1 );
    *w1 = 0.5e0+(0.5e0-*w);
    goto S250;
S200:
    T5 = 100.0e0*eps;
    *w = beta_asym ( &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;
}
//****************************************************************************80***

void beta_inc_values ( int *n_data, double *a, double *b, double *x, 
  double *fx )

//****************************************************************************80***
//
//  Purpose: 
//
//    BETA_INC_VALUES returns some values of the incomplete Beta function.
//
//  Discussion:
//
//    The incomplete Beta function may be written
//
//      BETA_INC(A,B,X) = Integral (0 to X) T**(A-1) * (1-T)**(B-1) dT
//                      / Integral (0 to 1) T**(A-1) * (1-T)**(B-1) dT
//
//    Thus,
//
//      BETA_INC(A,B,0.0) = 0.0
//      BETA_INC(A,B,1.0) = 1.0
//
//    Note that in Mathematica, the expressions:
//
//      BETA[A,B]   = Integral (0 to 1) T**(A-1) * (1-T)**(B-1) dT
//      BETA[X,A,B] = Integral (0 to X) T**(A-1) * (1-T)**(B-1) dT
//
//    and thus, to evaluate the incomplete Beta function requires:
//
//      BETA_INC(A,B,X) = BETA[X,A,B] / BETA[A,B]
//
//  Modified:
//
//    09 June 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz and Irene Stegun,
//    Handbook of Mathematical Functions,
//    US Department of Commerce, 1964.
//
//    Karl Pearson,
//    Tables of the Incomplete Beta Function,
//    Cambridge University Press, 1968.
//
//  Parameters:
//
//    Input/output, int *N_DATA.  The user sets N_DATA to 0 before the
//    first call.  On each call, the routine increments N_DATA by 1, and
//    returns the corresponding data; when there is no more data, the
//    output value of N_DATA will be 0 again.
//
//    Output, double *A, *B, the parameters of the function.
//
//    Output, double *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 30

  double a_vec[N_MAX] = {
     0.5E+00,  0.5E+00,  0.5E+00,  1.0E+00, 
     1.0E+00,  1.0E+00,  1.0E+00,  1.0E+00, 
     2.0E+00,  2.0E+00,  2.0E+00,  2.0E+00, 
     2.0E+00,  2.0E+00,  2.0E+00,  2.0E+00, 
     2.0E+00,  5.5E+00, 10.0E+00, 10.0E+00, 
    10.0E+00, 10.0E+00, 20.0E+00, 20.0E+00, 
    20.0E+00, 20.0E+00, 20.0E+00, 30.0E+00, 
    30.0E+00, 40.0E+00 };
  double b_vec[N_MAX] = { 
     0.5E+00,  0.5E+00,  0.5E+00,  0.5E+00, 
     0.5E+00,  0.5E+00,  0.5E+00,  1.0E+00, 
     2.0E+00,  2.0E+00,  2.0E+00,  2.0E+00, 
     2.0E+00,  2.0E+00,  2.0E+00,  2.0E+00, 
     2.0E+00,  5.0E+00,  0.5E+00,  5.0E+00, 
     5.0E+00, 10.0E+00,  5.0E+00, 10.0E+00, 
    10.0E+00, 20.0E+00, 20.0E+00, 10.0E+00, 
    10.0E+00, 20.0E+00 };
  double fx_vec[N_MAX] = { 
    0.0637686E+00, 0.2048328E+00, 1.0000000E+00, 0.0E+00,       
    0.0050126E+00, 0.0513167E+00, 0.2928932E+00, 0.5000000E+00, 
    0.028E+00,     0.104E+00,     0.216E+00,     0.352E+00,     
    0.500E+00,     0.648E+00,     0.784E+00,     0.896E+00,     
    0.972E+00,     0.4361909E+00, 0.1516409E+00, 0.0897827E+00, 
    1.0000000E+00, 0.5000000E+00, 0.4598773E+00, 0.2146816E+00, 
    0.9507365E+00, 0.5000000E+00, 0.8979414E+00, 0.2241297E+00, 
    0.7586405E+00, 0.7001783E+00 };
  double x_vec[N_MAX] = { 
    0.01E+00, 0.10E+00, 1.00E+00, 0.0E+00,  
    0.01E+00, 0.10E+00, 0.50E+00, 0.50E+00, 
    0.1E+00,  0.2E+00,  0.3E+00,  0.4E+00,  
    0.5E+00,  0.6E+00,  0.7E+00,  0.8E+00,  
    0.9E+00,  0.50E+00, 0.90E+00, 0.50E+00, 
    1.00E+00, 0.50E+00, 0.80E+00, 0.60E+00, 
    0.80E+00, 0.50E+00, 0.60E+00, 0.70E+00, 
    0.80E+00, 0.70E+00 };

  if ( *n_data < 0 )
  {
    *n_data = 0;
  }

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *a = 0.0E+00;
    *b = 0.0E+00;
    *x = 0.0E+00;
    *fx = 0.0E+00;
  }
  else
  {
    *a = a_vec[*n_data-1];
    *b = b_vec[*n_data-1];
    *x = x_vec[*n_data-1];
    *fx = fx_vec[*n_data-1];
  }
  return;
# undef N_MAX
}
//****************************************************************************80

double beta_log ( double *a0, double *b0 )

//****************************************************************************80
//
//  Purpose:
//
//    BETA_LOG evaluates the logarithm of the beta function.
//
//  Reference:
//
//    Armido DiDinato and Alfred Morris,
//    Algorithm 708: 
//    Significant Digit Computation of the Incomplete Beta Function Ratios,
//    ACM Transactions on Mathematical Software, 
//    Volume 18, 1993, pages 360-373.
//
//  Parameters:
//
//    Input, double *A0, *B0, the parameters of the function.
//    A0 and B0 should be nonnegative.
//
//    Output, double *BETA_LOG, the value of the logarithm
//    of the Beta function.
//
{
  static double e = .918938533204673e0;
  static double value,a,b,c,h,u,v,w,z;
  static int i,n;
  static double T1;

    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;
    value = gamma_log ( &a )+( gamma_log ( &b )- gamma_log ( &T1 ));
    return value;
S10:
    value = gamma_log ( &a )+algdiv(&a,&b);
    return value;
S20:
//
//  PROCEDURE WHEN 1 .LE. A .LT. 8
//
    if(a > 2.0e0) goto S40;
    if(b > 2.0e0) goto S30;
    value = gamma_log ( &a )+ gamma_log ( &b )-gsumln(&a,&b);
    return value;
S30:
    w = 0.0e0;
    if(b < 8.0e0) goto S60;
    value = gamma_log ( &a )+algdiv(&a,&b);
    return value;
S40:
//
//  REDUCTION OF A WHEN B .LE. 1000
//
    if(b > 1000.0e0) goto S80;
    n = ( int ) ( 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;
    value = w+ gamma_log ( &a )+algdiv(&a,&b);
    return value;
S60:
//
//  REDUCTION OF B WHEN B .LT. 8
//
    n = ( int ) ( b - 1.0e0 );
    z = 1.0e0;
    for ( i = 1; i <= n; i++ )
    {
        b -= 1.0e0;
        z *= (b/(a+b));
    }
    value = w+log(z)+( gamma_log ( &a )+( gamma_log ( &b )-gsumln(&a,&b)));
    return value;
S80:
//
//  REDUCTION OF A WHEN B .GT. 1000
//
    n = ( int ) ( a - 1.0e0 );
    w = 1.0e0;
    for ( i = 1; i <= n; i++ )
    {
        a -= 1.0e0;
        w *= (a/(1.0e0+a/b));
    }
    value = log(w)-(double)n*log(b)+( gamma_log ( &a )+algdiv(&a,&b));
    return value;
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;
    value = -(0.5e0*log(b))+e+w-v-u;
    return value;
S110:
    value = -(0.5e0*log(b))+e+w-u-v;
    return value;
}
//****************************************************************************80

double beta_pser ( double *a, double *b, double *x, double *eps )

//****************************************************************************80
//
//  Purpose:
// 
//    BETA_PSER uses a power series expansion to evaluate IX(A,B)(X).
//
//  Discussion:
//
//    BETA_PSER is used when B <= 1 or B*X <= 0.7.
//
//  Parameters:
//
//    Input, double *A, *B, the parameters.
//
//    Input, double *X, the point where the function
//    is to be evaluated.
//
//    Input, double *EPS, the tolerance.
//
//    Output, double BETA_PSER, the approximate value of IX(A,B)(X).
//
{
  static double bpser,a0,apb,b0,c,n,sum,t,tol,u,w,z;
  static int i,m;

    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)-beta_log(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 = gamma_ln1 ( &a0 );
    m = ( int ) ( 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 = gamma_ln1 ( &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 = n + 1.0e0;
    c *= ((0.5e0+(0.5e0-*b/n))**x);
    w = c/(*a+n);
    sum = sum + w;
    if(fabs(w) > tol) goto S110;
    bpser *= (1.0e0+*a*sum);
    return bpser;
}
//****************************************************************************80

double beta_rcomp ( double *a, double *b, double *x, double *y )

//****************************************************************************80
//
//  Purpose:
// 
//    BETA_RCOMP evaluates X**A * Y**B / Beta(A,B).
//
//  Parameters:
//
//    Input, double *A, *B, the parameters of the Beta function.
//    A and B should be nonnegative.
//
//    Input, double *X, *Y, define the numerator of the fraction.
//
//    Output, double BETA_RCOMP, the value of X**A * Y**B / Beta(A,B).
//
{
  static double Const = .398942280401433e0;
  static double brcomp,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;

    brcomp = 0.0e0;
    if(*x == 0.0e0 || *y == 0.0e0) return brcomp;
    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 -= beta_log(a,b);
    brcomp = exp(z);
    return brcomp;
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
//
    brcomp = exp(z);
    if(brcomp == 0.0e0) return brcomp;
    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;
    brcomp = brcomp*(a0*c)/(1.0e0+a0/b0);
    return brcomp;
S70:
//
//  ALGORITHM FOR 1 .LT. B0 .LT. 8
//
    u = gamma_ln1 ( &a0 );
    n = ( int ) ( 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:
    brcomp = a0*exp(z)*(1.0e0+gam1(&b0))/t;
    return brcomp;
S120:
//
//  ALGORITHM FOR B0 .GE. 8
//
    u = gamma_ln1 ( &a0 ) + algdiv ( &a0, &b0 );
    brcomp = a0*exp(z-u);
    return brcomp;
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;
S180:
    v = e-log(*y/y0);
S190:
    z = exp(-(*a*u+*b*v));
    brcomp = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b));
    return brcomp;
}
//****************************************************************************80

double beta_rcomp1 ( int *mu, double *a, double *b, double *x, double *y )

//****************************************************************************80
//
//  Purpose:
// 
//    BETA_RCOMP1 evaluates exp(MU) * X**A * Y**B / Beta(A,B).
//
//  Parameters:
//
//    Input, int MU, ?
//
//    Input, double A, B, the parameters of the Beta function.
//    A and B should be nonnegative.
//
//    Input, double X, Y, ?
//
//    Output, double BETA_RCOMP1, the value 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;

    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 -= beta_log(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 = gamma_ln1 ( &a0 );
    n = ( int ) ( 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 = gamma_ln1 ( &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;
S180:
    v = e-log(*y/y0);
S190:
    T4 = -(*a*u+*b*v);
    z = esum(mu,&T4);
    brcmp1 = Const*sqrt(*b*x0)*z*exp(-bcorr(a,b));
    return brcmp1;
}
//****************************************************************************80

double beta_up ( double *a, double *b, double *x, double *y, int *n, 
  double *eps )

//****************************************************************************80
//
//  Purpose:
// 
//    BETA_UP evaluates IX(A,B) - IX(A+N,B) where N is a positive integer.
//
//  Parameters:
//
//    Input, double *A, *B, the parameters of the function.
//    A and B should be nonnegative.
//
//    Input, double *X, *Y, ?
//
//    Input, int *N, ?
//
//    Input, double *EPS, the tolerance.
//
//    Output, double BETA_UP, the value of IX(A,B) - IX(A+N,B).
//
{
  static int K1 = 1;
  static int K2 = 0;
  static double bup,ap1,apb,d,l,r,t,w;
  static int i,k,kp1,mu,nm1;
//
//  OBTAIN THE SCALING FACTOR EXP(-MU) AND
//  EXP(MU)*(X**A*Y**B/BETA(A,B))/A
//
    apb = *a+*b;
    ap1 = *a+1.0e0;
    mu = 0;
    d = 1.0e0;
    if(*n == 1 || *a < 1.0e0) goto S10;
    if(apb < 1.1e0*ap1) goto S10;
    mu = ( int ) fabs ( exparg(&K1) );
    k = ( int ) exparg ( &K2 );
    if(k < mu) mu = k;
    t = mu;
    d = exp(-t);
S10:
    bup = beta_rcomp1 ( &mu, a, b, x, y ) / *a;
    if(*n == 1 || bup == 0.0e0) return bup;
    nm1 = *n-1;
    w = d;
//
//  LET K BE THE INDEX OF THE MAXIMUM TERM
//
    k = 0;
    if(*b <= 1.0e0) goto S50;
    if(*y > 1.e-4) goto S20;
    k = nm1;
    goto S30;
S20:
    r = (*b-1.0e0)**x/ *y-*a;
    if(r < 1.0e0) goto S50;
    t = ( double ) nm1;
    k = nm1;
    if ( r < t ) k = ( int ) r;
S30:
//
//          ADD THE INCREASING TERMS OF THE SERIES
//
    for ( i = 1; i <= k; i++ )
    {
        l = i-1;
        d = (apb+l)/(ap1+l)**x*d;
        w = w + d;
    }
    if(k == nm1) goto S70;
S50:
//
//          ADD THE REMAINING TERMS OF THE SERIES
//
    kp1 = k+1;
    for ( i = kp1; i <= nm1; i++ )
    {
        l = i-1;
        d = (apb+l)/(ap1+l)**x*d;
        w = w + d;
        if(d <= *eps*w) goto S70;
    }
S70:
//
//  TERMINATE THE PROCEDURE
//
    bup *= w;
    return bup;
}
//****************************************************************************80***

void binomial_cdf_values ( int *n_data, int *a, double *b, int *x, double *fx )

//****************************************************************************80***
//
//  Purpose: 
//
//    BINOMIAL_CDF_VALUES returns some values of the binomial CDF.
//
//  Discussion:
//
//    CDF(X)(A,B) is the probability of at most X successes in A trials,
//    given that the probability of success on a single trial is B.
//
//  Modified:
//
//    31 May 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz and Irene Stegun,
//    Handbook of Mathematical Functions,
//    US Department of Commerce, 1964.
//
//    Daniel Zwillinger,
//    CRC Standard Mathematical Tables and Formulae,
//    30th Edition, CRC Press, 1996, pages 651-652.
//
//  Parameters:
//
//    Input/output, int *N_DATA.  The user sets N_DATA to 0 before the
//    first call.  On each call, the routine increments N_DATA by 1, and
//    returns the corresponding data; when there is no more data, the
//    output value of N_DATA will be 0 again.
//
//    Output, int *A, double *B, the parameters of the function.
//
//    Output, int *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 17

  int a_vec[N_MAX] = { 
     2,  2,  2,  2, 
     2,  4,  4,  4, 
     4, 10, 10, 10, 
    10, 10, 10, 10, 
    10 };
  double b_vec[N_MAX] = { 
    0.05E+00, 0.05E+00, 0.05E+00, 0.50E+00, 
    0.50E+00, 0.25E+00, 0.25E+00, 0.25E+00, 
    0.25E+00, 0.05E+00, 0.10E+00, 0.15E+00, 
    0.20E+00, 0.25E+00, 0.30E+00, 0.40E+00, 
    0.50E+00 };
  double fx_vec[N_MAX] = { 
    0.9025E+00, 0.9975E+00, 1.0000E+00, 0.2500E+00, 
    0.7500E+00, 0.3164E+00, 0.7383E+00, 0.9492E+00, 
    0.9961E+00, 0.9999E+00, 0.9984E+00, 0.9901E+00, 
    0.9672E+00, 0.9219E+00, 0.8497E+00, 0.6331E+00, 
    0.3770E+00 };
  int x_vec[N_MAX] = { 
     0, 1, 2, 0, 
     1, 0, 1, 2, 
     3, 4, 4, 4, 
     4, 4, 4, 4, 
     4 };

  if ( *n_data < 0 )
  {
    *n_data = 0;
  }

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *a = 0;
    *b = 0.0E+00;
    *x = 0;
    *fx = 0.0E+00;
  }
  else
  {
    *a = a_vec[*n_data-1];
    *b = b_vec[*n_data-1];
    *x = x_vec[*n_data-1];
    *fx = fx_vec[*n_data-1];
  }
  return;
# undef N_MAX
}
//****************************************************************************80

void cdfbet ( int *which, double *p, double *q, double *x, double *y,
  double *a, double *b, int *status, double *bound )

//****************************************************************************80
//
//  Purpose:
// 
//    CDFBET evaluates the CDF of the Beta Distribution.
//
//  Discussion:
//
//    This routine calculates any one parameter of the beta distribution 
//    given the others.
//
//    The value P of the cumulative distribution function is calculated 
//    directly by code associated with the reference.
//
//    Computation of the other parameters involves a seach for a value that
//    produces the desired value of P.  The search relies on the
//    monotonicity of P with respect to the other parameters.
//
//    The beta density is proportional to t^(A-1) * (1-t)^(B-1).
//
//  Modified:
//
//    09 June 2004
//
//  Reference:
//
//    Armido DiDinato and Alfred Morris,
//    Algorithm 708: 
//    Significant Digit Computation of the Incomplete Beta Function Ratios,
//    ACM Transactions on Mathematical Software, 
//    Volume 18, 1993, pages 360-373.
//
//  Parameters:
//
//    Input, int *WHICH, indicates which of the next four argument
//    values is to be calculated from the others.
//    1: Calculate P and Q from X, Y, A and B;
//    2: Calculate X and Y from P, Q, A and B;
//    3: Calculate A from P, Q, X, Y and B;
//    4: Calculate B from P, Q, X, Y and A.
//
//    Input/output, double *P, the integral from 0 to X of the
//    chi-square distribution.  Input range: [0, 1].
//
//    Input/output, double *Q, equals 1-P.  Input range: [0, 1].
//
//    Input/output, double *X, the upper limit of integration 
//    of the beta density.  If it is an input value, it should lie in
//    the range [0,1].  If it is an output value, it will be searched for
//    in the range [0,1].
//
//    Input/output, double *Y, equal to 1-X.  If it is an input
//    value, it should lie in the range [0,1].  If it is an output value,
//    it will be searched for in the range [0,1].
//
//    Input/output, double *A, the first parameter of the beta
//    density.  If it is an input value, it should lie in the range
//    (0, +infinity).  If it is an output value, it will be searched
//    for in the range [1D-300,1D300].
//
//    Input/output, double *B, the second parameter of the beta
//    density.  If it is an input value, it should lie in the range
//    (0, +infinity).  If it is an output value, it will be searched
//    for in the range [1D-300,1D300].
//
//    Output, int *STATUS, reports the status of the computation.
//     0, if the calculation completed correctly;
//    -I, if the input parameter number I is out of range;
//    +1, if the answer appears to be lower than lowest search bound;
//    +2, if the answer appears to be higher than greatest search bound;
//    +3, if P + Q /= 1;
//    +4, if X + Y /= 1.
//
//    Output, double *BOUND, is only defined if STATUS is nonzero.
//    If STATUS is negative, then this is the value exceeded by parameter I.
//    if STATUS is 1 or 2, this is the search bound that was exceeded.
//
{
# define tol (1.0e-8)
# define atol (1.0e-50)
# define zero (1.0e-300)
# define inf 1.0e300
# define one 1.0e0

  static int K1 = 1;
  static double K2 = 0.0e0;
  static double K3 = 1.0e0;
  static double K8 = 0.5e0;
  static double K9 = 5.0e0;
  static double fx,xhi,xlo,cum,ccum,xy,pq;
  static unsigned long qhi,qleft,qporq;
  static double T4,T5,T6,T7,T10,T11,T12,T13,T14,T15;

  *status = 0;
  *bound = 0.0;
//
//     Check arguments
//
    if(!(*which < 1 || *which > 4)) goto S30;
    if(!(*which < 1)) goto S10;
    *bound = 1.0e0;
    goto S20;
S10:
    *bound = 4.0e0;
S20:
    *status = -1;
    return;
S30:
    if(*which == 1) goto S70;
//
//     P
//
    if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
    if(!(*p < 0.0e0)) goto S40;
    *bound = 0.0e0;
    goto S50;
S40:
    *bound = 1.0e0;
S50:
    *status = -2;
    return;
S70:
S60:
    if(*which == 1) goto S110;
//
//     Q
//
    if(!(*q < 0.0e0 || *q > 1.0e0)) goto S100;
    if(!(*q < 0.0e0)) goto S80;
    *bound = 0.0e0;
    goto S90;
S80:
    *bound = 1.0e0;
S90:
    *status = -3;
    return;
S110:
S100:
    if(*which == 2) goto S150;
//
//     X
//
    if(!(*x < 0.0e0 || *x > 1.0e0)) goto S140;
    if(!(*x < 0.0e0)) goto S120;
    *bound = 0.0e0;
    goto S130;
S120:
    *bound = 1.0e0;
S130:
    *status = -4;
    return;
S150:
S140:
    if(*which == 2) goto S190;
//
//     Y
//
    if(!(*y < 0.0e0 || *y > 1.0e0)) goto S180;
    if(!(*y < 0.0e0)) goto S160;
    *bound = 0.0e0;
    goto S170;
S160:
    *bound = 1.0e0;
S170:
    *status = -5;
    return;
S190:
S180:
    if(*which == 3) goto S210;
//
//     A
//
    if(!(*a <= 0.0e0)) goto S200;
    *bound = 0.0e0;
    *status = -6;
    return;
S210:
S200:
    if(*which == 4) goto S230;
//
//     B
//
    if(!(*b <= 0.0e0)) goto S220;
    *bound = 0.0e0;
    *status = -7;
    return;
S230:
S220:
    if(*which == 1) goto S270;
//
//     P + Q
//
    pq = *p+*q;
    if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0 * dpmpar ( &K1 ) ) ) goto S260;
    if(!(pq < 0.0e0)) goto S240;
    *bound = 0.0e0;
    goto S250;
S240:
    *bound = 1.0e0;
S250:
    *status = 3;
    return;
S270:
S260:
    if(*which == 2) goto S310;
//
//     X + Y
//
    xy = *x+*y;
    if(!(fabs(xy-0.5e0-0.5e0) > 3.0e0 * dpmpar ( &K1 ) ) ) goto S300;
    if(!(xy < 0.0e0)) goto S280;
    *bound = 0.0e0;
    goto S290;
S280:
    *bound = 1.0e0;
S290:
    *status = 4;
    return;
S310:
S300:
    if(!(*which == 1)) qporq = *p <= *q;
//
//     Select the minimum of P or Q
//     Calculate ANSWERS
//
    if(1 == *which) {
//
//     Calculating P and Q
//
        cumbet(x,y,a,b,p,q);
        *status = 0;
    }
    else if(2 == *which) {
//
//     Calculating X and Y
//
        T4 = atol;
        T5 = tol;
        dstzr(&K2,&K3,&T4,&T5);
        if(!qporq) goto S340;
        *status = 0;
        dzror(status,x,&fx,&xlo,&xhi,&qleft,&qhi);
        *y = one-*x;
S320:
        if(!(*status == 1)) goto S330;
        cumbet(x,y,a,b,&cum,&ccum);
        fx = cum-*p;
        dzror(status,x,&fx,&xlo,&xhi,&qleft,&qhi);
        *y = one-*x;
        goto S320;
S330:
        goto S370;
S340:
        *status = 0;
        dzror(status,y,&fx,&xlo,&xhi,&qleft,&qhi);
        *x = one-*y;
S350:
        if(!(*status == 1)) goto S360;
        cumbet(x,y,a,b,&cum,&ccum);
        fx = ccum-*q;
        dzror(status,y,&fx,&xlo,&xhi,&qleft,&qhi);
        *x = one-*y;
        goto S350;
S370:
S360:
        if(!(*status == -1)) goto S400;
        if(!qleft) goto S380;
        *status = 1;
        *bound = 0.0e0;
        goto S390;
S380:
        *status = 2;
        *bound = 1.0e0;
S400:
S390:
        ;
    }
    else if(3 == *which) {
//
//     Computing A
//
        *a = 5.0e0;
        T6 = zero;
        T7 = inf;
        T10 = atol;
        T11 = tol;
        dstinv(&T6,&T7,&K8,&K8,&K9,&T10,&T11);
        *status = 0;
        dinvr(status,a,&fx,&qleft,&qhi);
S410:
        if(!(*status == 1)) goto S440;
        cumbet(x,y,a,b,&cum,&ccum);
        if(!qporq) goto S420;
        fx = cum-*p;
        goto S430;
S420:
        fx = ccum-*q;
S430:
        dinvr(status,a,&fx,&qleft,&qhi);
        goto S410;
S440:
        if(!(*status == -1)) goto S470;
        if(!qleft) goto S450;
        *status = 1;
        *bound = zero;
        goto S460;
S450:
        *status = 2;
        *bound = inf;
S470:
S460:
        ;
    }
    else if(4 == *which) {
//
//     Computing B
//
        *b = 5.0e0;
        T12 = zero;
        T13 = inf;
        T14 = atol;
        T15 = tol;
        dstinv(&T12,&T13,&K8,&K8,&K9,&T14,&T15);
        *status = 0;
        dinvr(status,b,&fx,&qleft,&qhi);
S480:
        if(!(*status == 1)) goto S510;
        cumbet(x,y,a,b,&cum,&ccum);
        if(!qporq) goto S490;
        fx = cum-*p;
        goto S500;
S490:
        fx = ccum-*q;
S500:
        dinvr(status,b,&fx,&qleft,&qhi);
        goto S480;
S510:
        if(!(*status == -1)) goto S540;
        if(!qleft) goto S520;
        *status = 1;
        *bound = zero;
        goto S530;
S520:
        *status = 2;
        *bound = inf;
S530:
        ;
    }
S540:
    return;
# undef tol
# undef atol
# undef zero
# undef inf
# undef one
}
//****************************************************************************80

void cdfbin ( int *which, double *p, double *q, double *s, double *xn,
  double *pr, double *ompr, int *status, double *bound )

//****************************************************************************80
//
//  Purpose:
// 
//    CDFBIN evaluates the CDF of the Binomial distribution.
//
//  Discussion:
//
//    This routine calculates any one parameter of the binomial distribution 
//    given the others.
//
//    The value P of the cumulative distribution function is calculated 
//    directly.
//
//    Computation of the other parameters involves a seach for a value that
//    produces the desired value of P.  The search relies on the
//    monotonicity of P with respect to the other parameters.
//
//    P is the probablility of S or fewer successes in XN binomial trials,
//    each trial having an individual probability of success of PR.  
//
//  Modified:
//
//    09 June 2004
//
//  Reference:
//
//    Milton Abramowitz and Irene Stegun,
//    Handbook of Mathematical Functions 
//    1966, Formula 26.5.24.
//
//  Parameters:
//
//    Input, int *WHICH, indicates which of argument values is to 
//    be calculated from the others.
//    1: Calculate P and Q from S, XN, PR and OMPR;
//    2: Calculate S from P, Q, XN, PR and OMPR;
//    3: Calculate XN from P, Q, S, PR and OMPR;
//    4: Calculate PR and OMPR from P, Q, S and XN.
//
//    Input/output, double *P, the cumulation, from 0 to S,
//    of the binomial distribution.  If P is an input value, it should 
//    lie in the range [0,1].
//
//    Input/output, double *Q, equal to 1-P.  If Q is an input
//    value, it should lie in the range [0,1].  If Q is an output value,
//    it will lie in the range [0,1].
//
//    Input/output, double *S, the number of successes observed.
//    Whether this is an input or output value, it should lie in the
//    range [0,XN].
//
//    Input/output, double *XN, the number of binomial trials.
//    If this is an input value it should lie in the range: (0, +infinity).
//    If it is an output value it will be searched for in the
//    range [1.0D-300, 1.0D+300].
//
//    Input/output, double *PR, the probability of success in each 
//    binomial trial.  Whether this is an input or output value, it should
//    lie in the range: [0,1].
//
//    Input/output, double *OMPR, equal to 1-PR.  Whether this is an
//    input or output value, it should lie in the range [0,1].  Also, it should
//    be the case that PR + OMPR = 1.
//
//    Output, int *STATUS, reports the status of the computation.
//     0, if the calculation completed correctly;
//    -I, if the input parameter number I is out of range;
//    +1, if the answer appears to be lower than lowest search bound;
//    +2, if the answer appears to be higher than greatest search bound;
//    +3, if P + Q /= 1;
//    +4, if PR + OMPR /= 1.
//
//    Output, double *BOUND, is only defined if STATUS is nonzero.
//    If STATUS is negative, then this is the value exceeded by parameter I.
//    if STATUS is 1 or 2, this is the search bound that was exceeded.
//
{
# define atol (1.0e-50)
# define tol (1.0e-8)
# define zero (1.0e-300)
# define inf 1.0e300
# define one 1.0e0

  static int K1 = 1;
  static double K2 = 0.0e0;
  static double K3 = 0.5e0;
  static double K4 = 5.0e0;
  static double K11 = 1.0e0;
  static double fx,xhi,xlo,cum,ccum,pq,prompr;
  static unsigned long qhi,qleft,qporq;
  static double T5,T6,T7,T8,T9,T10,T12,T13;

  *status = 0;
  *bound = 0.0;
//
//     Check arguments
//
    if(!(*which < 1 && *which > 4)) goto S30;
    if(!(*which < 1)) goto S10;
    *bound = 1.0e0;
    goto S20;
S10:
    *bound = 4.0e0;
S20:
    *status = -1;
    return;
S30:
    if(*which == 1) goto S70;
//
//     P
//
    if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
    if(!(*p < 0.0e0)) goto S40;
    *bound = 0.0e0;
    goto S50;
S40:
    *bound = 1.0e0;
S50:
    *status = -2;
    return;
S70:
S60:
    if(*which == 1) goto S110;
//
//     Q
//
    if(!(*q < 0.0e0 || *q > 1.0e0)) goto S100;
    if(!(*q < 0.0e0)) goto S80;
    *bound = 0.0e0;
    goto S90;
S80:
    *bound = 1.0e0;
S90:
    *status = -3;
    return;
S110:
S100:
    if(*which == 3) goto S130;
//
//     XN
//
    if(!(*xn <= 0.0e0)) goto S120;
    *bound = 0.0e0;
    *status = -5;
    return;
S130:
S120:
    if(*which == 2) goto S170;
//
//     S
//
    if(!(*s < 0.0e0 || *which != 3 && *s > *xn)) goto S160;
    if(!(*s < 0.0e0)) goto S140;
    *bound = 0.0e0;
    goto S150;
S140:
    *bound = *xn;
S150:
    *status = -4;
    return;
S170:
S160:
    if(*which == 4) goto S210;
//
//     PR
//
    if(!(*pr < 0.0e0 || *pr > 1.0e0)) goto S200;
    if(!(*pr < 0.0e0)) goto S180;
    *bound = 0.0e0;
    goto S190;
S180:
    *bound = 1.0e0;
S190:
    *status = -6;
    return;
S210:
S200:
    if(*which == 4) goto S250;
//
//     OMPR
//
    if(!(*ompr < 0.0e0 || *ompr > 1.0e0)) goto S240;
    if(!(*ompr < 0.0e0)) goto S220;
    *bound = 0.0e0;
    goto S230;
S220:
    *bound = 1.0e0;
S230:
    *status = -7;
    return;
S250:
S240:
    if(*which == 1) goto S290;
//
//     P + Q
//
    pq = *p+*q;
    if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0 * dpmpar ( &K1 ) ) ) goto S280;
    if(!(pq < 0.0e0)) goto S260;
    *bound = 0.0e0;
    goto S270;
S260:
    *bound = 1.0e0;
S270:
    *status = 3;
    return;
S290:
S280:
    if(*which == 4) goto S330;
//
//     PR + OMPR
//
    prompr = *pr+*ompr;
    if(!(fabs(prompr-0.5e0-0.5e0) > 3.0e0 * dpmpar ( &K1 ) ) ) goto S320;
    if(!(prompr < 0.0e0)) goto S300;
    *bound = 0.0e0;
    goto S310;
S300:
    *bound = 1.0e0;
S310:
    *status = 4;
    return;
S330:
S320:
    if(!(*which == 1)) qporq = *p <= *q;
//
//     Select the minimum of P or Q
//     Calculate ANSWERS
//
    if(1 == *which) {
//
//     Calculating P
//
        cumbin(s,xn,pr,ompr,p,q);
        *status = 0;
    }
    else if(2 == *which) {
//
//     Calculating S
//
        *s = 5.0e0;
        T5 = atol;
        T6 = tol;
        dstinv(&K2,xn,&K3,&K3,&K4,&T5,&T6);
        *status = 0;
        dinvr(status,s,&fx,&qleft,&qhi);
S340:
        if(!(*status == 1)) goto S370;
        cumbin(s,xn,pr,ompr,&cum,&ccum);
        if(!qporq) goto S350;
        fx = cum-*p;
        goto S360;
S350:
        fx = ccum-*q;
S360:
        dinvr(status,s,&fx,&qleft,&qhi);
        goto S340;
S370:
        if(!(*status == -1)) goto S400;
        if(!qleft) goto S380;
        *status = 1;
        *bound = 0.0e0;
        goto S390;
S380:
        *status = 2;
        *bound = *xn;
S400:
S390:
        ;
    }
    else if(3 == *which) {
//
//     Calculating XN
//
        *xn = 5.0e0;
        T7 = zero;
        T8 = inf;
        T9 = atol;
        T10 = tol;
        dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
        *status = 0;
        dinvr(status,xn,&fx,&qleft,&qhi);
S410:
        if(!(*status == 1)) goto S440;
        cumbin(s,xn,pr,ompr,&cum,&ccum);
        if(!qporq) goto S420;
        fx = cum-*p;
        goto S430;
S420:
        fx = ccum-*q;
S430:
        dinvr(status,xn,&fx,&qleft,&qhi);
        goto S410;
S440:
        if(!(*status == -1)) goto S470;
        if(!qleft) goto S450;
        *status = 1;
        *bound = zero;
        goto S460;
S450:
        *status = 2;
        *bound = inf;
S470:
S460:
        ;
    }
    else if(4 == *which) {
//
//     Calculating PR and OMPR
//
        T12 = atol;
        T13 = tol;
        dstzr(&K2,&K11,&T12,&T13);
        if(!qporq) goto S500;
        *status = 0;
        dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
        *ompr = one-*pr;
S480:
        if(!(*status == 1)) goto S490;
        cumbin(s,xn,pr,ompr,&cum,&ccum);
        fx = cum-*p;
        dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
        *ompr = one-*pr;
        goto S480;
S490:
        goto S530;
S500:
        *status = 0;
        dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
        *pr = one-*ompr;
S510:
        if(!(*status == 1)) goto S520;
        cumbin(s,xn,pr,ompr,&cum,&ccum);
        fx = ccum-*q;
        dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
        *pr = one-*ompr;
        goto S510;
S530:
S520:
        if(!(*status == -1)) goto S560;
        if(!qleft) goto S540;
        *status = 1;
        *bound = 0.0e0;
        goto S550;
S540:
        *status = 2;
        *bound = 1.0e0;
S550:
        ;
    }
S560:
    return;
# undef atol
# undef tol
# undef zero
# undef inf
# undef one
}
//****************************************************************************80

void cdfchi ( int *which, double *p, double *q, double *x, double *df,
  int *status, double *bound )

//****************************************************************************80
//
//  Purpose:
// 
//    CDFCHI evaluates the CDF of the chi square distribution.
//
//  Discussion:
//
//    This routine calculates any one parameter of the chi square distribution 
//    given the others.
//
//    The value P of the cumulative distribution function is calculated 
//    directly.
//
//    Computation of the other parameters involves a seach for a value that
//    produces the desired value of P.  The search relies on the
//    monotonicity of P with respect to the other parameters.
//
//    The CDF of the chi square distribution can be evaluated 
//    within Mathematica by commands such as:
//
//      Needs["Statistics`ContinuousDistributions`"]
//      CDF [ ChiSquareDistribution [ DF ], X ]
//
//  Reference:
//
//    Milton Abramowitz and Irene Stegun,
//    Handbook of Mathematical Functions 
//    1966, Formula 26.4.19.
//
//    Stephen Wolfram,
//    The Mathematica Book,
//    Fourth Edition,
//    Wolfram Media / Cambridge University Press, 1999.
//
//  Parameters:
//
//    Input, int *WHICH, indicates which argument is to be calculated
//    from the others.
//    1: Calculate P and Q from X and DF;
//    2: Calculate X from P, Q and DF;
//    3: Calculate DF from P, Q and X.
//
//    Input/output, double *P, the integral from 0 to X of 
//    the chi-square distribution.  If this is an input value, it should
//    lie in the range [0,1].
//
//    Input/output, double *Q, equal to 1-P.  If Q is an input
//    value, it should lie in the range [0,1].  If Q is an output value,
//    it will lie in the range [0,1].
//
//    Input/output, double *X, the upper limit of integration 
//    of the chi-square distribution.  If this is an input 
//    value, it should lie in the range: [0, +infinity).  If it is an output
//    value, it will be searched for in the range: [0,1.0D+300].
//
//    Input/output, double *DF, the degrees of freedom of the
//    chi-square distribution.  If this is an input value, it should lie
//    in the range: (0, +infinity).  If it is an output value, it will be
//    searched for in the range: [ 1.0D-300, 1.0D+300].
//
//    Output, int *STATUS, reports the status of the computation.
//     0, if the calculation completed correctly;
//    -I, if the input parameter number I is out of range;
//    +1, if the answer appears to be lower than lowest search bound;
//    +2, if the answer appears to be higher than greatest search bound;
//    +3, if P + Q /= 1;
//    +10, an error was returned from CUMGAM.
//
//    Output, double *BOUND, is only defined if STATUS is nonzero.
//    If STATUS is negative, then this is the value exceeded by parameter I.
//    if STATUS is 1 or 2, this is the search bound that was exceeded.
//
{
# define tol (1.0e-8)
# define atol (1.0e-50)
# define zero (1.0e-300)
# define inf 1.0e300

  static int K1 = 1;
  static double K2 = 0.0e0;
  static double K4 = 0.5e0;
  static double K5 = 5.0e0;
  static double fx,cum,ccum,pq,porq;
  static unsigned long qhi,qleft,qporq;
  static double T3,T6,T7,T8,T9,T10,T11;

  *status = 0;
  *bound = 0.0;
//
//     Check arguments
//
    if(!(*which < 1 || *which > 3)) goto S30;
    if(!(*which < 1)) goto S10;
    *bound = 1.0e0;
    goto S20;
S10:
    *bound = 3.0e0;
S20:
    *status = -1;
    return;
S30:
    if(*which == 1) goto S70;
//
//     P
//
    if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
    if(!(*p < 0.0e0)) goto S40;
    *bound = 0.0e0;
    goto S50;
S40:
    *bound = 1.0e0;
S50:
    *status = -2;
    return;
S70:
S60:
    if(*which == 1) goto S110;
//
//     Q
//
    if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
    if(!(*q <= 0.0e0)) goto S80;
    *bound = 0.0e0;
    goto S90;
S80:
    *bound = 1.0e0;
S90:
    *status = -3;
    return;
S110:
S100:
    if(*which == 2) goto S130;
//
//     X
//
    if(!(*x < 0.0e0)) goto S120;
    *bound = 0.0e0;
    *status = -4;
    return;
S130:
S120:
    if(*which == 3) goto S150;
//
//     DF
//
    if(!(*df <= 0.0e0)) goto S140;
    *bound = 0.0e0;
    *status = -5;
    return;
S150:
S140:
    if(*which == 1) goto S190;
//
//     P + Q
//
    pq = *p+*q;
    if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0 * dpmpar ( &K1 ) ) ) goto S180;
    if(!(pq < 0.0e0)) goto S160;
    *bound = 0.0e0;
    goto S170;
S160:
    *bound = 1.0e0;
S170:
    *status = 3;
    return;
S190:
S180:
    if(*which == 1) goto S220;
//
//     Select the minimum of P or Q
//
    qporq = *p <= *q;
    if(!qporq) goto S200;
    porq = *p;
    goto S210;
S200:
    porq = *q;
S220:
S210:
//
//     Calculate ANSWERS
//
    if(1 == *which) {
//
//     Calculating P and Q
//
        *status = 0;
        cumchi(x,df,p,q);
        if(porq > 1.5e0) {
            *status = 10;
            return;
        }
    }
    else if(2 == *which) {
//
//     Calculating X
//
        *x = 5.0e0;
        T3 = inf;
        T6 = atol;
        T7 = tol;
        dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
        *status = 0;
        dinvr(status,x,&fx,&qleft,&qhi);
S230:
        if(!(*status == 1)) goto S270;
        cumchi(x,df,&cum,&ccum);
        if(!qporq) goto S240;
        fx = cum-*p;
        goto S250;
S240:
        fx = ccum-*q;
S250:
        if(!(fx+porq > 1.5e0)) goto S260;
        *status = 10;
        return;
S260:
        dinvr(status,x,&fx,&qleft,&qhi);
        goto S230;
S270:
        if(!(*status == -1)) goto S300;
        if(!qleft) goto S280;
        *status = 1;
        *bound = 0.0e0;
        goto S290;
S280:
        *status = 2;
        *bound = inf;
S300:
S290:
        ;
    }
    else if(3 == *which) {
//
//  Calculating DF
//
        *df = 5.0e0;
        T8 = zero;
        T9 = inf;
        T10 = atol;
        T11 = tol;
        dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
        *status = 0;
        dinvr(status,df,&fx,&qleft,&qhi);
S310:
        if(!(*status == 1)) goto S350;
        cumchi(x,df,&cum,&ccum);
        if(!qporq) goto S320;
        fx = cum-*p;
        goto S330;
S320:
        fx = ccum-*q;
S330:
        if(!(fx+porq > 1.5e0)) goto S340;
        *status = 10;
        return;
S340:
        dinvr(status,df,&fx,&qleft,&qhi);
        goto S310;
S350:
        if(!(*status == -1)) goto S380;
        if(!qleft) goto S360;
        *status = 1;
        *bound = zero;
        goto S370;
S360:
        *status = 2;
        *bound = inf;
S370:
        ;
    }
S380:
    return;
# undef tol
# undef atol
# undef zero
# undef inf
}
//****************************************************************************80

void cdfchn ( int *which, double *p, double *q, double *x, double *df,
  double *pnonc, int *status, double *bound )

//****************************************************************************80
//
//  Purpose:
// 
//    CDFCHN evaluates the CDF of the Noncentral Chi-Square.
//
//  Discussion:
//
//    This routine calculates any one parameter of the noncentral chi-square
//    distribution given values for the others.
//
//    The value P of the cumulative distribution function is calculated 
//    directly.
//
//    Computation of the other parameters involves a seach for a value that
//    produces the desired value of P.  The search relies on the
//    monotonicity of P with respect to the other parameters.
//
//    The computation time required for this routine is proportional
//    to the noncentrality parameter (PNONC).  Very large values of
//    this parameter can consume immense computer resources.  This is
//    why the search range is bounded by 10,000.
//
//    The CDF of the noncentral chi square distribution can be evaluated 
//    within Mathematica by commands such as:
//
//      Needs["Statistics`ContinuousDistributions`"]
//      CDF[ NoncentralChiSquareDistribution [ DF, LAMBDA ], X ]
//
//  Reference:
//
//    Milton Abramowitz and Irene Stegun,
//    Handbook of Mathematical Functions 
//    1966, Formula 26.5.25.
//
//    Stephen Wolfram,
//    The Mathematica Book,
//    Fourth Edition,
//    Wolfram Media / Cambridge University Press, 1999.
//
//  Parameters:
//
//    Input, int *WHICH, indicates which argument is to be calculated
//    from the others.
//    1: Calculate P and Q from X, DF and PNONC;
//    2: Calculate X from P, DF and PNONC;
//    3: Calculate DF from P, X and PNONC;
//    4: Calculate PNONC from P, X and DF.
//
//    Input/output, double *P, the integral from 0 to X of 
//    the noncentral chi-square distribution.  If this is an input
//    value, it should lie in the range: [0, 1.0-1.0D-16).
//
//    Input/output, double *Q, is generally not used by this 
//    subroutine and is only included for similarity with other routines.
//    However, if P is to be computed, then a value will also be computed
//    for Q.
//
//    Input, double *X, the upper limit of integration of the 
//    noncentral chi-square distribution.  If this is an input value, it
//    should lie in the range: [0, +infinity).  If it is an output value,
//    it will be sought in the range: [0,1.0D+300].
//
//    Input/output, double *DF, the number of degrees of freedom 
//    of the noncentral chi-square distribution.  If this is an input value,
//    it should lie in the range: (0, +infinity).  If it is an output value,
//    it will be searched for in the range: [ 1.0D-300, 1.0D+300].
//
//    Input/output, double *PNONC, the noncentrality parameter of 
//    the noncentral chi-square distribution.  If this is an input value, it
//    should lie in the range: [0, +infinity).  If it is an output value,
//    it will be searched for in the range: [0,1.0D+4]
//
//    Output, int *STATUS, reports on the calculation.
//    0, if calculation completed correctly;
//    -I, if input parameter number I is out of range;
//    1, if the answer appears to be lower than the lowest search bound;
//    2, if the answer appears to be higher than the greatest search bound.
//
//    Output, double *BOUND, is only defined if STATUS is nonzero.
//    If STATUS is negative, then this is the value exceeded by parameter I.
//    if STATUS is 1 or 2, this is the search bound that was exceeded.
//
{
# define tent4 1.0e4
# define tol (1.0e-8)
# define atol (1.0e-50)
# define zero (1.0e-300)
# define one (1.0e0-1.0e-16)
# define inf 1.0e300

  static double K1 = 0.0e0;
  static double K3 = 0.5e0;
  static double K4 = 5.0e0;
  static double fx,cum,ccum;
  static unsigned long qhi,qleft;
  static double T2,T5,T6,T7,T8,T9,T10,T11,T12,T13;

  *status = 0;
  *bound = 0.0;
//
//     Check arguments
//
    if(!(*which < 1 || *which > 4)) goto S30;
    if(!(*which < 1)) goto S10;
    *bound = 1.0e0;
    goto S20;
S10:
    *bound = 4.0e0;
S20:
    *status = -1;
    return;
S30:
    if(*which == 1) goto S70;
//
//     P
//
    if(!(*p < 0.0e0 || *p > one)) goto S60;
    if(!(*p < 0.0e0)) goto S40;
    *bound = 0.0e0;
    goto S50;
S40:
    *bound = one;
S50:
    *status = -2;
    return;
S70:
S60:
    if(*which == 2) goto S90;
//
//     X
//
    if(!(*x < 0.0e0)) goto S80;
    *bound = 0.0e0;
    *status = -4;
    return;
S90:
S80:
    if(*which == 3) goto S110;
//
//     DF
//
    if(!(*df <= 0.0e0)) goto S100;
    *bound = 0.0e0;
    *status = -5;
    return;
S110:
S100:
    if(*which == 4) goto S130;
//
//     PNONC
//
    if(!(*pnonc < 0.0e0)) goto S120;
    *bound = 0.0e0;
    *status = -6;
    return;
S130:
S120:
//
//     Calculate ANSWERS
//
    if(1 == *which) {
//
//     Calculating P and Q
//
        cumchn(x,df,pnonc,p,q);
        *status = 0;
    }
    else if(2 == *which) {
//
//     Calculating X
//
        *x = 5.0e0;
        T2 = inf;
        T5 = atol;
        T6 = tol;
        dstinv(&K1,&T2,&K3,&K3,&K4,&T5,&T6);
        *status = 0;
        dinvr(status,x,&fx,&qleft,&qhi);
S140:
        if(!(*status == 1)) goto S150;
        cumchn(x,df,pnonc,&cum,&ccum);
        fx = cum-*p;
        dinvr(status,x,&fx,&qleft,&qhi);
        goto S140;
S150:
        if(!(*status == -1)) goto S180;
        if(!qleft) goto S160;
        *status = 1;
        *bound = 0.0e0;
        goto S170;
S160:
        *status = 2;
        *bound = inf;
S180:
S170:
        ;
    }
    else if(3 == *which) {
//
//     Calculating DF
//
        *df = 5.0e0;
        T7 = zero;
        T8 = inf;
        T9 = atol;
        T10 = tol;
        dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
        *status = 0;
        dinvr(status,df,&fx,&qleft,&qhi);
S190:
        if(!(*status == 1)) goto S200;
        cumchn(x,df,pnonc,&cum,&ccum);
        fx = cum-*p;
        dinvr(status,df,&fx,&qleft,&qhi);
        goto S190;
S200:
        if(!(*status == -1)) goto S230;
        if(!qleft) goto S210;
        *status = 1;
        *bound = zero;
        goto S220;
S210:
        *status = 2;
        *bound = inf;
S230:
S220:
        ;
    }
    else if(4 == *which) {
//
//     Calculating PNONC
//
        *pnonc = 5.0e0;
        T11 = tent4;
        T12 = atol;
        T13 = tol;
        dstinv(&K1,&T11,&K3,&K3,&K4,&T12,&T13);
        *status = 0;
        dinvr(status,pnonc,&fx,&qleft,&qhi);
S240:
        if(!(*status == 1)) goto S250;
        cumchn(x,df,pnonc,&cum,&ccum);
        fx = cum-*p;
        dinvr(status,pnonc,&fx,&qleft,&qhi);
        goto S240;
S250:
        if(!(*status == -1)) goto S280;
        if(!qleft) goto S260;
        *status = 1;
        *bound = zero;
        goto S270;
S260:
        *status = 2;
        *bound = tent4;
S270:
        ;
    }
S280:
    return;
# undef tent4
# undef tol
# undef atol
# undef zero
# undef one
# undef inf
}
//****************************************************************************80

void cdff ( int *which, double *p, double *q, double *f, double *dfn,
  double *dfd, int *status, double *bound )

//****************************************************************************80
//
//  Purpose:
// 
//    CDFF evaluates the CDF of the F distribution.
//
//  Discussion:
//
//    This routine calculates any one parameter of the F distribution 
//    given the others.
//
//    The value P of the cumulative distribution function is calculated 
//    directly.
//
//    Computation of the other parameters involves a seach for a value that
//    produces the desired value of P.  The search relies on the
//    monotonicity of P with respect to the other parameters.
//
//    The value of the cumulative F distribution is not necessarily
//    monotone in either degree of freedom.  There thus may be two
//    values that provide a given CDF value.  This routine assumes
//    monotonicity and will find an arbitrary one of the two values.
//
//  Modified:
//
//    14 April 2007
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions 
//    1966, Formula 26.6.2.
//
//  Parameters:
//
//    Input, int *WHICH, indicates which argument is to be calculated
//    from the others.
//    1: Calculate P and Q from F, DFN and DFD;
//    2: Calculate F from P, Q, DFN and DFD;
//    3: Calculate DFN from P, Q, F and DFD;
//    4: Calculate DFD from P, Q, F and DFN.
//
//    Input/output, double *P, the integral from 0 to F of 
//    the F-density.  If it is an input value, it should lie in the
//    range [0,1].
//
//    Input/output, double *Q, equal to 1-P.  If Q is an input
//    value, it should lie in the range [0,1].  If Q is an output value,
//    it will lie in the range [0,1].
//
//    Input/output, double *F, the upper limit of integration 
//    of the F-density.  If this is an input value, it should lie in the
//    range [0, +infinity).  If it is an output value, it will be searched
//    for in the range [0,1.0D+300].
//
//    Input/output, double *DFN, the number of degrees of 
//    freedom of the numerator sum of squares.  If this is an input value,
//    it should lie in the range: (0, +infinity).  If it is an output value,
//    it will be searched for in the range: [ 1.0D-300, 1.0D+300].
//
//    Input/output, double *DFD, the number of degrees of freedom 
//    of the denominator sum of squares.  If this is an input value, it should
//    lie in the range: (0, +infinity).  If it is an output value, it will
//    be searched for in the  range: [ 1.0D-300, 1.0D+300].
//
//    Output, int *STATUS, reports the status of the computation.
//     0, if the calculation completed correctly;
//    -I, if the input parameter number I is out of range;
//    +1, if the answer appears to be lower than lowest search bound;
//    +2, if the answer appears to be higher than greatest search bound;
//    +3, if P + Q /= 1.
//
//    Output, double *BOUND, is only defined if STATUS is nonzero.
//    If STATUS is negative, then this is the value exceeded by parameter I.
//    if STATUS is 1 or 2, this is the search bound that was exceeded.
//
{
# define tol (1.0e-8)
# define atol (1.0e-50)
# define zero (1.0e-300)
# define inf 1.0e300

  static int K1 = 1;
  static double K2 = 0.0e0;
  static double K4 = 0.5e0;
  static double K5 = 5.0e0;
  static double pq,fx,cum,ccum;
  static unsigned long qhi,qleft,qporq;
  static double T3,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15;

  *status = 0;
  *bound = 0.0;
//
//  Check arguments
//
    if(!(*which < 1 || *which > 4)) goto S30;
    if(!(*which < 1)) goto S10;
    *bound = 1.0e0;
    goto S20;
S10:
    *bound = 4.0e0;
S20:
    *status = -1;
    return;
S30:
    if(*which == 1) goto S70;
//
//     P
//
    if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
    if(!(*p < 0.0e0)) goto S40;
    *bound = 0.0e0;
    goto S50;
S40:
    *bound = 1.0e0;
S50:
    *status = -2;
    return;
S70:
S60:
    if(*which == 1) goto S110;
//
//     Q
//
    if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
    if(!(*q <= 0.0e0)) goto S80;
    *bound = 0.0e0;
    goto S90;
S80:
    *bound = 1.0e0;
S90:
    *status = -3;
    return;
S110:
S100:
    if(*which == 2) goto S130;
//
//     F
//
    if(!(*f < 0.0e0)) goto S120;
    *bound = 0.0e0;
    *status = -4;
    return;
S130:
S120:
    if(*which == 3) goto S150;
//
//     DFN
//
    if(!(*dfn <= 0.0e0)) goto S140;
    *bound = 0.0e0;
    *status = -5;
    return;
S150:
S140:
    if(*which == 4) goto S170;
//
//     DFD
//
    if(!(*dfd <= 0.0e0)) goto S160;
    *bound = 0.0e0;
    *status = -6;
    return;
S170:
S160:
    if(*which == 1) goto S210;
//
//     P + Q
//
    pq = *p+*q;
    if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0 * dpmpar ( &K1 ) ) ) goto S200;
    if(!(pq < 0.0e0)) goto S180;
    *bound = 0.0e0;
    goto S190;
S180:
    *bound = 1.0e0;
S190:
    *status = 3;
    return;
S210:
S200:
    if(!(*which == 1)) qporq = *p <= *q;
//
//     Select the minimum of P or Q
//     Calculate ANSWERS
//
    if(1 == *which) {
//
//     Calculating P
//
        cumf(f,dfn,dfd,p,q);
        *status = 0;
    }
    else if(2 == *which) {
//
//     Calculating F
//
        *f = 5.0e0;
        T3 = inf;
        T6 = atol;
        T7 = tol;
        dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
        *status = 0;
        dinvr(status,f,&fx,&qleft,&qhi);
S220:
        if(!(*status == 1)) goto S250;
        cumf(f,dfn,dfd,&cum,&ccum);
        if(!qporq) goto S230;
        fx = cum-*p;
        goto S240;
S230:
        fx = ccum-*q;
S240:
        dinvr(status,f,&fx,&qleft,&qhi);
        goto S220;
S250:
        if(!(*status == -1)) goto S280;
        if(!qleft) goto S260;
        *status = 1;
        *bound = 0.0e0;
        goto S270;
S260:
        *status = 2;
        *bound = inf;
S280:
S270:
        ;
    }
//
//  Calculate DFN.
//
//  Note that, in the original calculation, the lower bound for DFN was 0.
//  Using DFN = 0 causes an error in CUMF when it calls BETA_INC.
//  The lower bound was set to the more reasonable value of 1.
//  JVB, 14 April 2007.
//
  else if ( 3 == *which ) 
  {

    T8 = 1.0;
    T9 = inf;
    T10 = atol;
    T11 = tol;
    dstinv ( &T8, &T9, &K4, &K4, &K5, &T10, &T11 );

    *status = 0;
    *dfn = 5.0;
    fx = 0.0;

    dinvr ( status, dfn, &fx, &qleft, &qhi );

    while ( *status == 1 )
    {
      cumf ( f, dfn, dfd, &cum, &ccum );

      if ( *p <= *q ) 
      {
        fx = cum - *p;
      }
      else
      {
        fx = ccum - *q;
      }
      dinvr ( status, dfn, &fx, &qleft, &qhi );
    }

    if ( *status == -1 )
    {
      if ( qleft )
      {
        *status = 1;
        *bound = 1.0;
      }
      else
      {
        *status = 2;
        *bound = inf;
      }
    }
  }
//
//  Calculate DFD.
//
//  Note that, in the original calculation, the lower bound for DFD was 0.
//  Using DFD = 0 causes an error in CUMF when it calls BETA_INC.
//  The lower bound was set to the more reasonable value of 1.
//  JVB, 14 April 2007.
//
//
  else if ( 4 == *which )
  {

    T12 = 1.0;
    T13 = inf;
    T14 = atol;
    T15 = tol;
    dstinv ( &T12, &T13, &K4, &K4, &K5, &T14, &T15 );

    *status = 0;
    *dfd = 5.0;
    fx = 0.0;
    dinvr ( status, dfd, &fx, &qleft, &qhi );

    while ( *status == 1 )
    {
      cumf ( f, dfn, dfd, &cum, &ccum );

      if ( *p <= *q ) 
      {
        fx = cum - *p;
      }
      else
      {
        fx = ccum - *q;
      }
      dinvr ( status, dfd, &fx, &qleft, &qhi );
    }

    if ( *status == -1 )
    {
      if ( qleft )
      {
        *status = 1;
        *bound = 1.0;
      }
      else
      {
        *status = 2;
        *bound = inf;
      }
    }
  }

  return;
# undef tol
# undef atol
# undef zero
# undef inf
}
//****************************************************************************80

void cdffnc ( int *which, double *p, double *q, double *f, double *dfn,
  double *dfd, double *phonc, int *status, double *bound )

//****************************************************************************80
//
//  Purpose:
// 
//    CDFFNC evaluates the CDF of the Noncentral F distribution.
//
//  Discussion:
//
//    This routine originally used 1.0E+300 as the upper bound for the 
//    interval in which many of the missing parameters are to be sought.
//    Since the underlying rootfinder routine needs to evaluate the 
//    function at this point, it is no surprise that the program was
//    experiencing overflows.  A less extravagant upper bound
//    is being tried for now!
//
//
//    This routine calculates any one parameter of the Noncentral F distribution 
//    given the others.
//
//    The value P of the cumulative distribution function is calculated 
//    directly.
//
//    Computation of the other parameters involves a seach for a value that
//    produces the desired value of P.  The search relies on the
//    monotonicity of P with respect to the other parameters.
//
//    The computation time required for this routine is proportional
//    to the noncentrality parameter PNONC.  Very large values of
//    this parameter can consume immense computer resources.  This is
//    why the search range is bounded by 10,000.
//
//    The value of the cumulative noncentral F distribution is not
//    necessarily monotone in either degree of freedom.  There thus
//    may be two values that provide a given CDF value.  This routine
//    assumes monotonicity and will find an arbitrary one of the two
//    values.
//
//    The CDF of the noncentral F distribution can be evaluated 
//    within Mathematica by commands such as:
//
//      Needs["Statistics`ContinuousDistributions`"]
//      CDF [ NoncentralFRatioDistribution [ DFN, DFD, PNONC ], X ]
//
//  Modified:
//
//    15 June 2004
//
//  Reference:
//
//    Milton Abramowitz and Irene Stegun,
//    Handbook of Mathematical Functions 
//    1966, Formula 26.6.20.
//
//    Stephen Wolfram,
//    The Mathematica Book,
//    Fourth Edition,
//    Wolfram Media / Cambridge University Press, 1999.
//
//  Parameters:
//
//    Input, int *WHICH, indicates which argument is to be calculated
//    from the others.
//    1: Calculate P and Q from F, DFN, DFD and PNONC;
//    2: Calculate F from P, Q, DFN, DFD and PNONC;
//    3: Calculate DFN from P, Q, F, DFD and PNONC;
//    4: Calculate DFD from P, Q, F, DFN and PNONC;
//    5: Calculate PNONC from P, Q, F, DFN and DFD.
//
//    Input/output, double *P, the integral from 0 to F of 
//    the noncentral F-density.  If P is an input value it should
//    lie in the range [0,1) (Not including 1!).
//
//    Dummy, double *Q, is not used by this subroutine,
//    and is only included for similarity with the other routines.
//    Its input value is not checked.  If P is to be computed, the
//    Q is set to 1 - P.
//
//    Input/output, double *F, the upper limit of integration 
//    of the noncentral F-density.  If this is an input value, it should
//    lie in the range: [0, +infinity).  If it is an output value, it
//    will be searched for in the range: [0,1.0D+30].
//
//    Input/output, double *DFN, the number of degrees of freedom 
//    of the numerator sum of squares.  If this is an input value, it should
//    lie in the range: (0, +infinity).  If it is an output value, it will
//    be searched for in the range: [ 1.0, 1.0D+30].
//
//    Input/output, double *DFD, the number of degrees of freedom 
//    of the denominator sum of squares.  If this is an input value, it should
//    be in range: (0, +infinity).  If it is an output value, it will be
//    searched for in the range [1.0, 1.0D+30].
//
//    Input/output, double *PNONC, the noncentrality parameter
//    If this is an input value, it should be nonnegative.
//    If it is an output value, it will be searched for in the range: [0,1.0D+4].
//
//    Output, int *STATUS, reports the status of the computation.
//     0, if the calculation completed correctly;
//    -I, if the input parameter number I is out of range;
//    +1, if the answer appears to be lower than lowest search bound;
//    +2, if the answer appears to be higher than greatest search bound;
//    +3, if P + Q /= 1.
//
//    Output, double *BOUND, is only defined if STATUS is nonzero.
//    If STATUS is negative, then this is the value exceeded by parameter I.
//    if STATUS is 1 or 2, this is the search bound that was exceeded.
//
{
# define tent4 1.0e4
# define tol (1.0e-8)
# define atol (1.0e-50)
# define zero (1.0e-300)
# define one (1.0e0-1.0e-16)
# define inf 1.0e300

  static double K1 = 0.0e0;
  static double K3 = 0.5e0;
  static double K4 = 5.0e0;
  static double fx,cum,ccum;
  static unsigned long qhi,qleft;
  static double T2,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15,T16,T17;

  *status = 0;
  *bound = 0.0;
//
//     Check arguments
//
    if(!(*which < 1 || *which > 5)) goto S30;
    if(!(*which < 1)) goto S10;
    *bound = 1.0e0;
    goto S20;
S10:
    *bound = 5.0e0;
S20:
    *status = -1;
    return;
S30:
    if(*which == 1) goto S70;
//
//     P
//
    if(!(*p < 0.0e0 || *p > one)) goto S60;
    if(!(*p < 0.0e0)) goto S40;
    *bound = 0.0e0;
    goto S50;
S40:
    *bound = one;
S50:
    *status = -2;
    return;
S70:
S60:
    if(*which == 2) goto S90;
//
//     F
//
    if(!(*f < 0.0e0)) goto S80;
    *bound = 0.0e0;
    *status = -4;
    return;
S90:
S80:
    if(*which == 3) goto S110;
//
//     DFN
//
    if(!(*dfn <= 0.0e0)) goto S100;
    *bound = 0.0e0;
    *status = -5;
    return;
S110:
S100:
    if(*which == 4) goto S130;
//
//     DFD
//
    if(!(*dfd <= 0.0e0)) goto S120;
    *bound = 0.0e0;
    *status = -6;
    return;
S130:
S120:
    if(*which == 5) goto S150;
//
//     PHONC
//
    if(!(*phonc < 0.0e0)) goto S140;
    *bound = 0.0e0;
    *status = -7;
    return;
S150:
S140:
//
//     Calculate ANSWERS
//
    if(1 == *which) {
//
//     Calculating P
//
        cumfnc(f,dfn,dfd,phonc,p,q);
        *status = 0;
    }
    else if(2 == *which) {
//
//     Calculating F
//
        *f = 5.0e0;
        T2 = inf;
        T5 = atol;
        T6 = tol;
        dstinv(&K1,&T2,&K3,&K3,&K4,&T5,&T6);
        *status = 0;
        dinvr(status,f,&fx,&qleft,&qhi);
S160:
        if(!(*status == 1)) goto S170;
        cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
        fx = cum-*p;
        dinvr(status,f,&fx,&qleft,&qhi);
        goto S160;
S170:
        if(!(*status == -1)) goto S200;
        if(!qleft) goto S180;
        *status = 1;
        *bound = 0.0e0;
        goto S190;
S180:
        *status = 2;
        *bound = inf;
S200:
S190:
        ;
    }
    else if(3 == *which) {
//
//     Calculating DFN
//
        *dfn = 5.0e0;
        T7 = zero;
        T8 = inf;
        T9 = atol;
        T10 = tol;
        dstinv(&T7,&T8,&K3,&K3,&K4,&T9,&T10);
        *status = 0;
        dinvr(status,dfn,&fx,&qleft,&qhi);
S210:
        if(!(*status == 1)) goto S220;
        cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
        fx = cum-*p;
        dinvr(status,dfn,&fx,&qleft,&qhi);
        goto S210;
S220:
        if(!(*status == -1)) goto S250;
        if(!qleft) goto S230;
        *status = 1;
        *bound = zero;
        goto S240;
S230:
        *status = 2;
        *bound = inf;
S250:
S240:
        ;
    }
    else if(4 == *which) {
//
//     Calculating DFD
//
        *dfd = 5.0e0;
        T11 = zero;
        T12 = inf;
        T13 = atol;
        T14 = tol;
        dstinv(&T11,&T12,&K3,&K3,&K4,&T13,&T14);
        *status = 0;
        dinvr(status,dfd,&fx,&qleft,&qhi);
S260:
        if(!(*status == 1)) goto S270;
        cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
        fx = cum-*p;
        dinvr(status,dfd,&fx,&qleft,&qhi);
        goto S260;
S270:
        if(!(*status == -1)) goto S300;
        if(!qleft) goto S280;
        *status = 1;
        *bound = zero;
        goto S290;
S280:
        *status = 2;
        *bound = inf;
S300:
S290:
        ;
    }
    else if(5 == *which) {
//
//     Calculating PHONC
//
        *phonc = 5.0e0;
        T15 = tent4;
        T16 = atol;
        T17 = tol;
        dstinv(&K1,&T15,&K3,&K3,&K4,&T16,&T17);
        *status = 0;
        dinvr(status,phonc,&fx,&qleft,&qhi);
S310:
        if(!(*status == 1)) goto S320;
        cumfnc(f,dfn,dfd,phonc,&cum,&ccum);
        fx = cum-*p;
        dinvr(status,phonc,&fx,&qleft,&qhi);
        goto S310;
S320:
        if(!(*status == -1)) goto S350;
        if(!qleft) goto S330;
        *status = 1;
        *bound = 0.0e0;
        goto S340;
S330:
        *status = 2;
        *bound = tent4;
S340:
        ;
    }
S350:
    return;
# undef tent4
# undef tol
# undef atol
# undef zero
# undef one
# undef inf
}
//****************************************************************************80

void cdfgam ( int *which, double *p, double *q, double *x, double *shape,
  double *scale, int *status, double *bound )

//****************************************************************************80
//
//  Purpose:
// 
//    CDFGAM evaluates the CDF of the Gamma Distribution.
//
//  Discussion:
//
//    This routine calculates any one parameter of the Gamma distribution 
//    given the others.
//
//    The cumulative distribution function P is calculated directly.
//
//    Computation of the other parameters involves a seach for a value that
//    produces the desired value of P.  The search relies on the
//    monotonicity of P with respect to the other parameters.
//
//    The gamma density is proportional to T**(SHAPE - 1) * EXP(- SCALE * T)
//
//  Reference:
//
//    Armido DiDinato and Alfred Morris, 
//    Computation of the incomplete gamma function ratios and their inverse,
//    ACM Transactions on Mathematical Software,
//    Volume 12, 1986, pages 377-393.
//
//  Parameters:
//
//    Input, int *WHICH, indicates which argument is to be calculated
//    from the others.
//    1: Calculate P and Q from X, SHAPE and SCALE;
//    2: Calculate X from P, Q, SHAPE and SCALE;
//    3: Calculate SHAPE from P, Q, X and SCALE;
//    4: Calculate SCALE from P, Q, X and SHAPE.
//
//    Input/output, double *P, the integral from 0 to X of the 
//    Gamma density.  If this is an input value, it should lie in the 
//    range: [0,1].
//
//    Input/output, double *Q, equal to 1-P.  If Q is an input
//    value, it should lie in the range [0,1].  If Q is an output value,
//    it will lie in the range [0,1].
//
//    Input/output, double *X, the upper limit of integration of 
//    the Gamma density.  If this is an input value, it should lie in the
//    range: [0, +infinity).  If it is an output value, it will lie in
//    the range: [0,1E300].
//
//    Input/output, double *SHAPE, the shape parameter of the 
//    Gamma density.  If this is an input value, it should lie in the range: 
//    (0, +infinity).  If it is an output value, it will be searched for
//    in the range: [1.0D-300,1.0D+300].
//
//    Input/output, double *SCALE, the scale parameter of the 
//    Gamma density.  If this is an input value, it should lie in the range
//    (0, +infinity).  If it is an output value, it will be searched for
//    in the range: (1.0D-300,1.0D+300].
//
//    Output, int *STATUS, reports the status of the computation.
//     0, if the calculation completed correctly;
//    -I, if the input parameter number I is out of range;
//    +1, if the answer appears to be lower than lowest search bound;
//    +2, if the answer appears to be higher than greatest search bound;
//    +3, if P + Q /= 1;
//    +10, if the Gamma or inverse Gamma routine cannot compute the answer.  
//    This usually happens only for X and SHAPE very large (more than 1.0D+10.
//
//    Output, double *BOUND, is only defined if STATUS is nonzero.
//    If STATUS is negative, then this is the value exceeded by parameter I.
//    if STATUS is 1 or 2, this is the search bound that was exceeded.
//
{
# define tol (1.0e-8)
# define atol (1.0e-50)
# define zero (1.0e-300)
# define inf 1.0e300

  static int K1 = 1;
  static double K5 = 0.5e0;
  static double K6 = 5.0e0;
  static double xx,fx,xscale,cum,ccum,pq,porq;
  static int ierr;
  static unsigned long qhi,qleft,qporq;
  static double T2,T3,T4,T7,T8,T9;

  *status = 0;
  *bound = 0.0;
//
//     Check arguments
//
    if(!(*which < 1 || *which > 4)) goto S30;
    if(!(*which < 1)) goto S10;
    *bound = 1.0e0;
    goto S20;
S10:
    *bound = 4.0e0;
S20:
    *status = -1;
    return;
S30:
    if(*which == 1) goto S70;
//
//     P
//
    if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
    if(!(*p < 0.0e0)) goto S40;
    *bound = 0.0e0;
    goto S50;
S40:
    *bound = 1.0e0;
S50:
    *status = -2;
    return;
S70:
S60:
    if(*which == 1) goto S110;
//
//     Q
//
    if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
    if(!(*q <= 0.0e0)) goto S80;
    *bound = 0.0e0;
    goto S90;
S80:
    *bound = 1.0e0;
S90:
    *status = -3;
    return;
S110:
S100:
    if(*which == 2) goto S130;
//
//     X
//
    if(!(*x < 0.0e0)) goto S120;
    *bound = 0.0e0;
    *status = -4;
    return;
S130:
S120:
    if(*which == 3) goto S150;
//
//     SHAPE
//
    if(!(*shape <= 0.0e0)) goto S140;
    *bound = 0.0e0;
    *status = -5;
    return;
S150:
S140:
    if(*which == 4) goto S170;
//
//     SCALE
//
    if(!(*scale <= 0.0e0)) goto S160;
    *bound = 0.0e0;
    *status = -6;
    return;
S170:
S160:
    if(*which == 1) goto S210;
//
//     P + Q
//
    pq = *p+*q;
    if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*dpmpar(&K1))) goto S200;
    if(!(pq < 0.0e0)) goto S180;
    *bound = 0.0e0;
    goto S190;
S180:
    *bound = 1.0e0;
S190:
    *status = 3;
    return;
S210:
S200:
    if(*which == 1) goto S240;
//
//     Select the minimum of P or Q
//
    qporq = *p <= *q;
    if(!qporq) goto S220;
    porq = *p;
    goto S230;
S220:
    porq = *q;
S240:
S230:
//
//     Calculate ANSWERS
//
    if(1 == *which) {
//
//     Calculating P
//
        *status = 0;
        xscale = *x**scale;
        cumgam(&xscale,shape,p,q);
        if(porq > 1.5e0) *status = 10;
    }
    else if(2 == *which) {
//
//     Computing X
//
        T2 = -1.0e0;
        gamma_inc_inv ( shape, &xx, &T2, p, q, &ierr );
        if(ierr < 0.0e0) {
            *status = 10;
            return;
        }
        else  {
            *x = xx/ *scale;
            *status = 0;
        }
    }
    else if(3 == *which) {
//
//     Computing SHAPE
//
        *shape = 5.0e0;
        xscale = *x**scale;
        T3 = zero;
        T4 = inf;
        T7 = atol;
        T8 = tol;
        dstinv(&T3,&T4,&K5,&K5,&K6,&T7,&T8);
        *status = 0;
        dinvr(status,shape,&fx,&qleft,&qhi);
S250:
        if(!(*status == 1)) goto S290;
        cumgam(&xscale,shape,&cum,&ccum);
        if(!qporq) goto S260;
        fx = cum-*p;
        goto S270;
S260:
        fx = ccum-*q;
S270:
        if(!(qporq && cum > 1.5e0 || !qporq && ccum > 1.5e0)) goto S280;
        *status = 10;
        return;
S280:
        dinvr(status,shape,&fx,&qleft,&qhi);
        goto S250;
S290:
        if(!(*status == -1)) goto S320;
        if(!qleft) goto S300;
        *status = 1;
        *bound = zero;
        goto S310;
S300:
        *status = 2;
        *bound = inf;
S320:
S310:
        ;
    }
    else if(4 == *which) {
//
//     Computing SCALE
//
        T9 = -1.0e0;
        gamma_inc_inv ( shape, &xx, &T9, p, q, &ierr );
        if(ierr < 0.0e0) {
            *status = 10;
            return;
        }
        else  {
            *scale = xx/ *x;
            *status = 0;
        }
    }
    return;
# undef tol
# undef atol
# undef zero
# undef inf
}
//****************************************************************************80

void cdfnbn ( int *which, double *p, double *q, double *s, double *xn,
  double *pr, double *ompr, int *status, double *bound )

//****************************************************************************80
//
//  Purpose:
// 
//    CDFNBN evaluates the CDF of the Negative Binomial distribution
//
//  Discussion:
//
//    This routine calculates any one parameter of the negative binomial
//    distribution given values for the others.
//
//    The cumulative negative binomial distribution returns the
//    probability that there will be F or fewer failures before the
//    S-th success in binomial trials each of which has probability of
//    success PR.
//
//    The individual term of the negative binomial is the probability of
//    F failures before S successes and is
//    Choose( F, S+F-1 ) * PR^(S) * (1-PR)^F
//
//    Computation of other parameters involve a seach for a value that
//    produces the desired value of P.  The search relies on the
//    monotonicity of P with respect to the other parameters.
//
//  Reference:
//
//    Milton Abramowitz and Irene Stegun,
//    Handbook of Mathematical Functions 
//    1966, Formula 26.5.26.
//
//  Parameters:
//
//    Input, int WHICH, indicates which argument is to be calculated
//    from the others.
//    1: Calculate P and Q from F, S, PR and OMPR;
//    2: Calculate F from P, Q, S, PR and OMPR;
//    3: Calculate S from P, Q, F, PR and OMPR;
//    4: Calculate PR and OMPR from P, Q, F and S.
//
//    Input/output, double P, the cumulation from 0 to F of 
//    the negative binomial distribution.  If P is an input value, it
//    should lie in the range [0,1].
//
//    Input/output, double Q, equal to 1-P.  If Q is an input
//    value, it should lie in the range [0,1].  If Q is an output value,
//    it will lie in the range [0,1].
//
//    Input/output, double F, the upper limit of cumulation of 
//    the binomial distribution.  There are F or fewer failures before 
//    the S-th success.  If this is an input value, it may lie in the
//    range [0,+infinity), and if it is an output value, it will be searched
//    for in the range [0,1.0D+300].
//
//    Input/output, double S, the number of successes.
//    If this is an input value, it should lie in the range: [0, +infinity).
//    If it is an output value, it will be searched for in the range: 
//    [0, 1.0D+300].
//
//    Input/output, double PR, the probability of success in each 
//    binomial trial.  Whether an input or output value, it should lie in the
//    range [0,1].
//
//    Input/output, double OMPR, the value of (1-PR).  Whether an
//    input or output value, it should lie in the range [0,1].  
//
//    Output, int STATUS, reports the status of the computation.
//     0, if the calculation completed correctly;
//    -I, if the input parameter number I is out of range;
//    +1, if the answer appears to be lower than lowest search bound;
//    +2, if the answer appears to be higher than greatest search bound;
//    +3, if P + Q /= 1;
//    +4, if PR + OMPR /= 1.
//
//    Output, double BOUND, is only defined if STATUS is nonzero.
//    If STATUS is negative, then this is the value exceeded by parameter I.
//    if STATUS is 1 or 2, this is the search bound that was exceeded.
//
{
# define tol (1.0e-8)
# define atol (1.0e-50)
# define inf 1.0e300
# define one 1.0e0

  static int K1 = 1;
  static double K2 = 0.0e0;
  static double K4 = 0.5e0;
  static double K5 = 5.0e0;
  static double K11 = 1.0e0;
  static double fx,xhi,xlo,pq,prompr,cum,ccum;
  static unsigned long qhi,qleft,qporq;
  static double T3,T6,T7,T8,T9,T10,T12,T13;

  *status = 0;
  *bound = 0.0;
//
//     Check arguments
//
    if(!(*which < 1 || *which > 4)) goto S30;
    if(!(*which < 1)) goto S10;
    *bound = 1.0e0;
    goto S20;
S10:
    *bound = 4.0e0;
S20:
    *status = -1;
    return;
S30:
    if(*which == 1) goto S70;
//
//     P
//
    if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
    if(!(*p < 0.0e0)) goto S40;
    *bound = 0.0e0;
    goto S50;
S40:
    *bound = 1.0e0;
S50:
    *status = -2;
    return;
S70:
S60:
    if(*which == 1) goto S110;
//
//     Q
//
    if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
    if(!(*q <= 0.0e0)) goto S80;
    *bound = 0.0e0;
    goto S90;
S80:
    *bound = 1.0e0;
S90:
    *status = -3;
    return;
S110:
S100:
    if(*which == 2) goto S130;
//
//     S
//
    if(!(*s < 0.0e0)) goto S120;
    *bound = 0.0e0;
    *status = -4;
    return;
S130:
S120:
    if(*which == 3) goto S150;
//
//     XN
//
    if(!(*xn < 0.0e0)) goto S140;
    *bound = 0.0e0;
    *status = -5;
    return;
S150:
S140:
    if(*which == 4) goto S190;
//
//     PR
//
    if(!(*pr < 0.0e0 || *pr > 1.0e0)) goto S180;
    if(!(*pr < 0.0e0)) goto S160;
    *bound = 0.0e0;
    goto S170;
S160:
    *bound = 1.0e0;
S170:
    *status = -6;
    return;
S190:
S180:
    if(*which == 4) goto S230;
//
//     OMPR
//
    if(!(*ompr < 0.0e0 || *ompr > 1.0e0)) goto S220;
    if(!(*ompr < 0.0e0)) goto S200;
    *bound = 0.0e0;
    goto S210;
S200:
    *bound = 1.0e0;
S210:
    *status = -7;
    return;
S230:
S220:
    if(*which == 1) goto S270;
//
//     P + Q
//
    pq = *p+*q;
    if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*dpmpar(&K1))) goto S260;
    if(!(pq < 0.0e0)) goto S240;
    *bound = 0.0e0;
    goto S250;
S240:
    *bound = 1.0e0;
S250:
    *status = 3;
    return;
S270:
S260:
    if(*which == 4) goto S310;
//
//     PR + OMPR
//
    prompr = *pr+*ompr;
    if(!(fabs(prompr-0.5e0-0.5e0) > 3.0e0*dpmpar(&K1))) goto S300;
    if(!(prompr < 0.0e0)) goto S280;
    *bound = 0.0e0;
    goto S290;
S280:
    *bound = 1.0e0;
S290:
    *status = 4;
    return;
S310:
S300:
    if(!(*which == 1)) qporq = *p <= *q;
//
//     Select the minimum of P or Q
//     Calculate ANSWERS
//
    if(1 == *which) {
//
//     Calculating P
//
        cumnbn(s,xn,pr,ompr,p,q);
        *status = 0;
    }
    else if(2 == *which) {
//
//     Calculating S
//
        *s = 5.0e0;
        T3 = inf;
        T6 = atol;
        T7 = tol;
        dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
        *status = 0;
        dinvr(status,s,&fx,&qleft,&qhi);
S320:
        if(!(*status == 1)) goto S350;
        cumnbn(s,xn,pr,ompr,&cum,&ccum);
        if(!qporq) goto S330;
        fx = cum-*p;
        goto S340;
S330:
        fx = ccum-*q;
S340:
        dinvr(status,s,&fx,&qleft,&qhi);
        goto S320;
S350:
        if(!(*status == -1)) goto S380;
        if(!qleft) goto S360;
        *status = 1;
        *bound = 0.0e0;
        goto S370;
S360:
        *status = 2;
        *bound = inf;
S380:
S370:
        ;
    }
    else if(3 == *which) {
//
//     Calculating XN
//
        *xn = 5.0e0;
        T8 = inf;
        T9 = atol;
        T10 = tol;
        dstinv(&K2,&T8,&K4,&K4,&K5,&T9,&T10);
        *status = 0;
        dinvr(status,xn,&fx,&qleft,&qhi);
S390:
        if(!(*status == 1)) goto S420;
        cumnbn(s,xn,pr,ompr,&cum,&ccum);
        if(!qporq) goto S400;
        fx = cum-*p;
        goto S410;
S400:
        fx = ccum-*q;
S410:
        dinvr(status,xn,&fx,&qleft,&qhi);
        goto S390;
S420:
        if(!(*status == -1)) goto S450;
        if(!qleft) goto S430;
        *status = 1;
        *bound = 0.0e0;
        goto S440;
S430:
        *status = 2;
        *bound = inf;
S450:
S440:
        ;
    }
    else if(4 == *which) {
//
//     Calculating PR and OMPR
//
        T12 = atol;
        T13 = tol;
        dstzr(&K2,&K11,&T12,&T13);
        if(!qporq) goto S480;
        *status = 0;
        dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
        *ompr = one-*pr;
S460:
        if(!(*status == 1)) goto S470;
        cumnbn(s,xn,pr,ompr,&cum,&ccum);
        fx = cum-*p;
        dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
        *ompr = one-*pr;
        goto S460;
S470:
        goto S510;
S480:
        *status = 0;
        dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
        *pr = one-*ompr;
S490:
        if(!(*status == 1)) goto S500;
        cumnbn(s,xn,pr,ompr,&cum,&ccum);
        fx = ccum-*q;
        dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
        *pr = one-*ompr;
        goto S490;
S510:
S500:
        if(!(*status == -1)) goto S540;
        if(!qleft) goto S520;
        *status = 1;
        *bound = 0.0e0;
        goto S530;
S520:
        *status = 2;
        *bound = 1.0e0;
S530:
        ;
    }
S540:
    return;
# undef tol
# undef atol
# undef inf
# undef one
}
//****************************************************************************80

void cdfnor ( int *which, double *p, double *q, double *x, double *mean,
  double *sd, int *status, double *bound )

//****************************************************************************80
//
//  Purpose:
// 
//    CDFNOR evaluates the CDF of the Normal distribution.
//
//  Discussion:
//
//    A slightly modified version of ANORM from SPECFUN
//    is used to calculate the cumulative standard normal distribution.
//
//    The rational functions from pages 90-95 of Kennedy and Gentle
//    are used as starting values to Newton's Iterations which 
//    compute the inverse standard normal.  Therefore no searches are
//    necessary for any parameter.
//
//    For X < -15, the asymptotic expansion for the normal is used  as
//    the starting value in finding the inverse standard normal.
//
//    The normal density is proportional to
//    exp( - 0.5D+00 * (( X - MEAN)/SD)**2)
//
//  Reference:
//
//    Milton Abramowitz and Irene Stegun,
//    Handbook of Mathematical Functions 
//    1966, Formula 26.2.12.
//
//    William Cody,
//    Algorithm 715: SPECFUN - A Portable FORTRAN Package of
//      Special Function Routines and Test Drivers,
//    ACM Transactions on Mathematical Software,
//    Volume 19, pages 22-32, 1993.
//
//    Kennedy and Gentle,
//    Statistical Computing,
//    Marcel Dekker, NY, 1980,
//    QA276.4  K46
//
//  Parameters:
//
//    Input, int *WHICH, indicates which argument is to be calculated
//    from the others.
//    1: Calculate P and Q from X, MEAN and SD;
//    2: Calculate X from P, Q, MEAN and SD;
//    3: Calculate MEAN from P, Q, X and SD;
//    4: Calculate SD from P, Q, X and MEAN.
//
//    Input/output, double *P, the integral from -infinity to X 
//    of the Normal density.  If this is an input or output value, it will
//    lie in the range [0,1].
//
//    Input/output, double *Q, equal to 1-P.  If Q is an input
//    value, it should lie in the range [0,1].  If Q is an output value,
//    it will lie in the range [0,1].
//
//    Input/output, double *X, the upper limit of integration of 
//    the Normal density.
//
//    Input/output, double *MEAN, the mean of the Normal density.
//
//    Input/output, double *SD, the standard deviation of the 
//    Normal density.  If this is an input value, it should lie in the
//    range (0,+infinity).
//
//    Output, int *STATUS, the status of the calculation.
//    0, if calculation completed correctly;
//    -I, if input parameter number I is out of range;
//    1, if answer appears to be lower than lowest search bound;
//    2, if answer appears to be higher than greatest search bound;
//    3, if P + Q /= 1.
//
//    Output, double *BOUND, is only defined if STATUS is nonzero.
//    If STATUS is negative, then this is the value exceeded by parameter I.
//    if STATUS is 1 or 2, this is the search bound that was exceeded.
//
{
  static int K1 = 1;
  static double z,pq;

  *status = 0;
  *bound = 0.0;
//
//     Check arguments
//
    *status = 0;
    if(!(*which < 1 || *which > 4)) goto S30;
    if(!(*which < 1)) goto S10;
    *bound = 1.0e0;
    goto S20;
S10:
    *bound = 4.0e0;
S20:
    *status = -1;
    return;
S30:
    if(*which == 1) goto S70;
//
//     P
//
    if(!(*p <= 0.0e0 || *p > 1.0e0)) goto S60;
    if(!(*p <= 0.0e0)) goto S40;
    *bound = 0.0e0;
    goto S50;
S40:
    *bound = 1.0e0;
S50:
    *status = -2;
    return;
S70:
S60:
    if(*which == 1) goto S110;
//
//     Q
//
    if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
    if(!(*q <= 0.0e0)) goto S80;
    *bound = 0.0e0;
    goto S90;
S80:
    *bound = 1.0e0;
S90:
    *status = -3;
    return;
S110:
S100:
    if(*which == 1) goto S150;
//
//     P + Q
//
    pq = *p+*q;
    if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*dpmpar(&K1))) goto S140;
    if(!(pq < 0.0e0)) goto S120;
    *bound = 0.0e0;
    goto S130;
S120:
    *bound = 1.0e0;
S130:
    *status = 3;
    return;
S150:
S140:
    if(*which == 4) goto S170;
//
//     SD
//
    if(!(*sd <= 0.0e0)) goto S160;
    *bound = 0.0e0;
    *status = -6;
    return;
S170:
S160:
//
//     Calculate ANSWERS
//
    if(1 == *which) {
//
//     Computing P
//
        z = (*x-*mean)/ *sd;
        cumnor(&z,p,q);
    }
    else if(2 == *which) {
//
//     Computing X
//
        z = dinvnr(p,q);
        *x = *sd*z+*mean;
    }
    else if(3 == *which) {
//
//     Computing the MEAN
//
        z = dinvnr(p,q);
        *mean = *x-*sd*z;
    }
    else if(4 == *which) {
//
//     Computing SD
//
        z = dinvnr(p,q);
        *sd = (*x-*mean)/z;
    }
    return;
}
//****************************************************************************80

void cdfpoi ( int *which, double *p, double *q, double *s, double *xlam,
  int *status, double *bound )

//****************************************************************************80
//
//  Purpose:
// 
//    CDFPOI evaluates the CDF of the Poisson distribution.
//
//  Discussion:
//
//    This routine calculates any one parameter of the Poisson distribution 
//    given the others.
//
//    The value P of the cumulative distribution function is calculated 
//    directly.
//
//    Computation of other parameters involve a seach for a value that
//    produces the desired value of P.  The search relies on the
//    monotonicity of P with respect to the other parameters.
//
//  Reference:
//
//    Milton Abramowitz and Irene Stegun,
//    Handbook of Mathematical Functions 
//    1966, Formula 26.4.21.
//
//  Parameters:
//
//    Input, int *WHICH, indicates which argument is to be calculated
//    from the others.
//    1: Calculate P and Q from S and XLAM;
//    2: Calculate A from P, Q and XLAM;
//    3: Calculate XLAM from P, Q and S.
//
//    Input/output, double *P, the cumulation from 0 to S of the 
//    Poisson density.  Whether this is an input or output value, it will
//    lie in the range [0,1].
//
//    Input/output, double *Q, equal to 1-P.  If Q is an input
//    value, it should lie in the range [0,1].  If Q is an output value,
//    it will lie in the range [0,1].
//
//    Input/output, double *S, the upper limit of cumulation of
//    the Poisson CDF.  If this is an input value, it should lie in
//    the range: [0, +infinity).  If it is an output value, it will be
//    searched for in the range: [0,1.0D+300].
//
//    Input/output, double *XLAM, the mean of the Poisson 
//    distribution.  If this is an input value, it should lie in the range
//    [0, +infinity).  If it is an output value, it will be searched for
//    in the range: [0,1E300].
//
//    Output, int *STATUS, reports the status of the computation.
//     0, if the calculation completed correctly;
//    -I, if the input parameter number I is out of range;
//    +1, if the answer appears to be lower than lowest search bound;
//    +2, if the answer appears to be higher than greatest search bound;
//    +3, if P + Q /= 1.
//
//    Output, double *BOUND, is only defined if STATUS is nonzero.
//    If STATUS is negative, then this is the value exceeded by parameter I.
//    if STATUS is 1 or 2, this is the search bound that was exceeded.
//
{
# define tol (1.0e-8)
# define atol (1.0e-50)
# define inf 1.0e300

  static int K1 = 1;
  static double K2 = 0.0e0;
  static double K4 = 0.5e0;
  static double K5 = 5.0e0;
  static double fx,cum,ccum,pq;
  static unsigned long qhi,qleft,qporq;
  static double T3,T6,T7,T8,T9,T10;

  *status = 0;
  *bound = 0.0;
//
//     Check arguments
//
    if(!(*which < 1 || *which > 3)) goto S30;
    if(!(*which < 1)) goto S10;
    *bound = 1.0e0;
    goto S20;
S10:
    *bound = 3.0e0;
S20:
    *status = -1;
    return;
S30:
    if(*which == 1) goto S70;
//
//     P
//
    if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
    if(!(*p < 0.0e0)) goto S40;
    *bound = 0.0e0;
    goto S50;
S40:
    *bound = 1.0e0;
S50:
    *status = -2;
    return;
S70:
S60:
    if(*which == 1) goto S110;
//
//     Q
//
    if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
    if(!(*q <= 0.0e0)) goto S80;
    *bound = 0.0e0;
    goto S90;
S80:
    *bound = 1.0e0;
S90:
    *status = -3;
    return;
S110:
S100:
    if(*which == 2) goto S130;
//
//     S
//
    if(!(*s < 0.0e0)) goto S120;
    *bound = 0.0e0;
    *status = -4;
    return;
S130:
S120:
    if(*which == 3) goto S150;
//
//     XLAM
//
    if(!(*xlam < 0.0e0)) goto S140;
    *bound = 0.0e0;
    *status = -5;
    return;
S150:
S140:
    if(*which == 1) goto S190;
//
//     P + Q
//
    pq = *p+*q;
    if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*dpmpar(&K1))) goto S180;
    if(!(pq < 0.0e0)) goto S160;
    *bound = 0.0e0;
    goto S170;
S160:
    *bound = 1.0e0;
S170:
    *status = 3;
    return;
S190:
S180:
    if(!(*which == 1)) qporq = *p <= *q;
//
//     Select the minimum of P or Q
//     Calculate ANSWERS
//
    if(1 == *which) {
//
//     Calculating P
//
        cumpoi(s,xlam,p,q);
        *status = 0;
    }
    else if(2 == *which) {
//
//     Calculating S
//
        *s = 5.0e0;
        T3 = inf;
        T6 = atol;
        T7 = tol;
        dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
        *status = 0;
        dinvr(status,s,&fx,&qleft,&qhi);
S200:
        if(!(*status == 1)) goto S230;
        cumpoi(s,xlam,&cum,&ccum);
        if(!qporq) goto S210;
        fx = cum-*p;
        goto S220;
S210:
        fx = ccum-*q;
S220:
        dinvr(status,s,&fx,&qleft,&qhi);
        goto S200;
S230:
        if(!(*status == -1)) goto S260;
        if(!qleft) goto S240;
        *status = 1;
        *bound = 0.0e0;
        goto S250;
S240:
        *status = 2;
        *bound = inf;
S260:
S250:
        ;
    }
    else if(3 == *which) {
//
//     Calculating XLAM
//
        *xlam = 5.0e0;
        T8 = inf;
        T9 = atol;
        T10 = tol;
        dstinv(&K2,&T8,&K4,&K4,&K5,&T9,&T10);
        *status = 0;
        dinvr(status,xlam,&fx,&qleft,&qhi);
S270:
        if(!(*status == 1)) goto S300;
        cumpoi(s,xlam,&cum,&ccum);
        if(!qporq) goto S280;
        fx = cum-*p;
        goto S290;
S280:
        fx = ccum-*q;
S290:
        dinvr(status,xlam,&fx,&qleft,&qhi);
        goto S270;
S300:
        if(!(*status == -1)) goto S330;
        if(!qleft) goto S310;
        *status = 1;
        *bound = 0.0e0;
        goto S320;
S310:
        *status = 2;
        *bound = inf;
S320:
        ;
    }
S330:
    return;
# undef tol
# undef atol
# undef inf
}
//****************************************************************************80

void cdft ( int *which, double *p, double *q, double *t, double *df,
  int *status, double *bound )

//****************************************************************************80
//
//  Purpose:
// 
//    CDFT evaluates the CDF of the T distribution.
//
//  Discussion:
//
//    This routine calculates any one parameter of the T distribution 
//    given the others.
//
//    The value P of the cumulative distribution function is calculated 
//    directly.
//
//    Computation of other parameters involve a seach for a value that
//    produces the desired value of P.   The search relies on the
//    monotonicity of P with respect to the other parameters.
//
//    The original version of this routine allowed the search interval
//    to extend from -1.0E+300 to +1.0E+300, which is fine until you 
//    try to evaluate a function at such a point!
//
//  Reference:
//
//    Milton Abramowitz and Irene Stegun,
//    Handbook of Mathematical Functions 
//    1966, Formula 26.5.27.
//
//  Parameters:
//
//    Input, int *WHICH, indicates which argument is to be calculated
//    from the others.
//    1 : Calculate P and Q from T and DF;
//    2 : Calculate T from P, Q and DF;
//    3 : Calculate DF from P, Q and T.
//
//    Input/output, double *P, the integral from -infinity to T of 
//    the T-density.  Whether an input or output value, this will lie in the
//    range [0,1].
//
//    Input/output, double *Q, equal to 1-P.  If Q is an input
//    value, it should lie in the range [0,1].  If Q is an output value,
//    it will lie in the range [0,1].
//
//    Input/output, double *T, the upper limit of integration of 
//    the T-density.  If this is an input value, it may have any value.
//    It it is an output value, it will be searched for in the range
//    [ -1.0D+30, 1.0D+30 ].
//
//    Input/output, double *DF, the number of degrees of freedom
//    of the T distribution.  If this is an input value, it should lie
//    in the range: (0 , +infinity).  If it is an output value, it will be
//    searched for in the range: [1, 1.0D+10].
//
//    Output, int *STATUS, reports the status of the computation.
//     0, if the calculation completed correctly;
//    -I, if the input parameter number I is out of range;
//    +1, if the answer appears to be lower than lowest search bound;
//    +2, if the answer appears to be higher than greatest search bound;
//    +3, if P + Q /= 1.
//
//    Output, double *BOUND, is only defined if STATUS is nonzero.
//    If STATUS is negative, then this is the value exceeded by parameter I.
//    if STATUS is 1 or 2, this is the search bound that was exceeded.
//
{
# define tol (1.0e-8)
# define atol (1.0e-50)
# define zero (1.0e-300)
# define inf 1.0e30
# define maxdf 1.0e10

  static int K1 = 1;
  static double K4 = 0.5e0;
  static double K5 = 5.0e0;
  static double fx,cum,ccum,pq;
  static unsigned long qhi,qleft,qporq;
  static double T2,T3,T6,T7,T8,T9,T10,T11;

  *status = 0;
  *bound = 0.0;
//
//     Check arguments
//
    if(!(*which < 1 || *which > 3)) goto S30;
    if(!(*which < 1)) goto S10;
    *bound = 1.0e0;
    goto S20;
S10:
    *bound = 3.0e0;
S20:
    *status = -1;
    return;
S30:
    if(*which == 1) goto S70;
//
//     P
//
    if(!(*p <= 0.0e0 || *p > 1.0e0)) goto S60;
    if(!(*p <= 0.0e0)) goto S40;
    *bound = 0.0e0;
    goto S50;
S40:
    *bound = 1.0e0;
S50:
    *status = -2;
    return;
S70:
S60:
    if(*which == 1) goto S110;
//
//     Q
//
    if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
    if(!(*q <= 0.0e0)) goto S80;
    *bound = 0.0e0;
    goto S90;
S80:
    *bound = 1.0e0;
S90:
    *status = -3;
    return;
S110:
S100:
    if(*which == 3) goto S130;
//
//     DF
//
    if(!(*df <= 0.0e0)) goto S120;
    *bound = 0.0e0;
    *status = -5;
    return;
S130:
S120:
    if(*which == 1) goto S170;
//
//     P + Q
//
    pq = *p+*q;
    if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*dpmpar(&K1))) goto S160;
    if(!(pq < 0.0e0)) goto S140;
    *bound = 0.0e0;
    goto S150;
S140:
    *bound = 1.0e0;
S150:
    *status = 3;
    return;
S170:
S160:
    if(!(*which == 1)) qporq = *p <= *q;
//
//     Select the minimum of P or Q
//     Calculate ANSWERS
//
    if(1 == *which) {
//
//     Computing P and Q
//
        cumt(t,df,p,q);
        *status = 0;
    }
    else if(2 == *which) {
//
//     Computing T
//     .. Get initial approximation for T
//
        *t = dt1(p,q,df);
        T2 = -inf;
        T3 = inf;
        T6 = atol;
        T7 = tol;
        dstinv(&T2,&T3,&K4,&K4,&K5,&T6,&T7);
        *status = 0;
        dinvr(status,t,&fx,&qleft,&qhi);
S180:
        if(!(*status == 1)) goto S210;
        cumt(t,df,&cum,&ccum);
        if(!qporq) goto S190;
        fx = cum-*p;
        goto S200;
S190:
        fx = ccum-*q;
S200:
        dinvr(status,t,&fx,&qleft,&qhi);
        goto S180;
S210:
        if(!(*status == -1)) goto S240;
        if(!qleft) goto S220;
        *status = 1;
        *bound = -inf;
        goto S230;
S220:
        *status = 2;
        *bound = inf;
S240:
S230:
        ;
    }
    else if(3 == *which) {
//
//     Computing DF
//
        *df = 5.0e0;
        T8 = zero;
        T9 = maxdf;
        T10 = atol;
        T11 = tol;
        dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
        *status = 0;
        dinvr(status,df,&fx,&qleft,&qhi);
S250:
        if(!(*status == 1)) goto S280;
        cumt(t,df,&cum,&ccum);
        if(!qporq) goto S260;
        fx = cum-*p;
        goto S270;
S260:
        fx = ccum-*q;
S270:
        dinvr(status,df,&fx,&qleft,&qhi);
        goto S250;
S280:
        if(!(*status == -1)) goto S310;
        if(!qleft) goto S290;
        *status = 1;
        *bound = zero;
        goto S300;
S290:
        *status = 2;
        *bound = maxdf;
S300:
        ;
    }
S310:
    return;
# undef tol
# undef atol
# undef zero
# undef inf
# undef maxdf
}
//****************************************************************************80**

void chi_noncentral_cdf_values ( int *n_data, double *x, double *lambda, 
  int *df, double *cdf )

//****************************************************************************80**
//
//  Purpose:
//
//    CHI_NONCENTRAL_CDF_VALUES returns values of the noncentral chi CDF.
//
//  Discussion:
//
//    The CDF of the noncentral chi square distribution can be evaluated
//    within Mathematica by commands such as:
//
//      Needs["Statistics`ContinuousDistributions`"]
//      CDF [ NoncentralChiSquareDistribution [ DF, LAMBDA ], X ]
//
//  Modified:
//
//    12 June 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Stephen Wolfram,
//    The Mathematica Book,
//    Fourth Edition,
//    Wolfram Media / Cambridge University Press, 1999.
//
//  Parameters:
//
//    Input/output, int *N_DATA.  The user sets N_DATA to 0 before the
//    first call.  On each call, the routine increments N_DATA by 1, and
//    returns the corresponding data; when there is no more data, the
//    output value of N_DATA will be 0 again.
//
//    Output, double *X, the argument of the function.
//
//    Output, double *LAMBDA, the noncentrality parameter.
//
//    Output, int *DF, the number of degrees of freedom.
//
//    Output, double *CDF, the noncentral chi CDF.
//
{
# define N_MAX 27

  double cdf_vec[N_MAX] = {
    0.839944E+00, 0.695906E+00, 0.535088E+00, 
    0.764784E+00, 0.620644E+00, 0.469167E+00,
    0.307088E+00, 0.220382E+00, 0.150025E+00, 
    0.307116E-02, 0.176398E-02, 0.981679E-03, 
    0.165175E-01, 0.202342E-03, 0.498448E-06, 
    0.151325E-01, 0.209041E-02, 0.246502E-03, 
    0.263684E-01, 0.185798E-01, 0.130574E-01, 
    0.583804E-01, 0.424978E-01, 0.308214E-01, 
    0.105788E+00, 0.794084E-01, 0.593201E-01 };
  int df_vec[N_MAX] = {
      1,   2,   3, 
      1,   2,   3, 
      1,   2,   3,
      1,   2,   3, 
     60,  80, 100, 
      1,   2,   3, 
     10,  10,  10, 
     10,  10,  10, 
     10,  10,  10 };
  double lambda_vec[N_MAX] = { 
     0.5E+00,  0.5E+00,  0.5E+00, 
     1.0E+00,  1.0E+00,  1.0E+00, 
     5.0E+00,  5.0E+00,  5.0E+00, 
    20.0E+00, 20.0E+00, 20.0E+00, 
    30.0E+00, 30.0E+00, 30.0E+00, 
     5.0E+00,  5.0E+00,  5.0E+00, 
     2.0E+00,  3.0E+00,  4.0E+00, 
     2.0E+00,  3.0E+00,  4.0E+00, 
     2.0E+00,  3.0E+00,  4.0E+00 };
  double x_vec[N_MAX] = {
     3.000E+00,  3.000E+00,  3.000E+00, 
     3.000E+00,  3.000E+00,  3.000E+00, 
     3.000E+00,  3.000E+00,  3.000E+00, 
     3.000E+00,  3.000E+00,  3.000E+00, 
    60.000E+00, 60.000E+00, 60.000E+00, 
     0.050E+00,  0.050E+00,  0.050E+00, 
     4.000E+00,  4.000E+00,  4.000E+00, 
     5.000E+00,  5.000E+00,  5.000E+00, 
     6.000E+00,  6.000E+00,  6.000E+00 };

  if ( *n_data < 0 )
  {
    *n_data = 0;
  }

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *x = 0.0E+00;
    *lambda = 0.0E+00;
    *df = 0;
    *cdf = 0.0E+00;
  }
  else
  {
    *x = x_vec[*n_data-1];
    *lambda = lambda_vec[*n_data-1];
    *df = df_vec[*n_data-1];
    *cdf = cdf_vec[*n_data-1];
  }

  return;
# undef N_MAX
}
//****************************************************************************80***

void chi_square_cdf_values ( int *n_data, int *a, double *x, double *fx )

//****************************************************************************80***
//
//  Purpose: 
//
//    CHI_SQUARE_CDF_VALUES returns some values of the Chi-Square CDF.
//
//  Discussion:
//
//    The value of CHI_CDF ( DF, X ) can be evaluated in Mathematica by 
//    commands like:
//
//      Needs["Statistics`ContinuousDistributions`"]
//      CDF[ChiSquareDistribution[DF], X ]
//
//  Modified:
//
//    11 June 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz and Irene Stegun,
//    Handbook of Mathematical Functions,
//    US Department of Commerce, 1964.
//
//    Stephen Wolfram,
//    The Mathematica Book,
//    Fourth Edition,
//    Wolfram Media / Cambridge University Press, 1999.
//
//  Parameters:
//
//    Input/output, int *N_DATA.  The user sets N_DATA to 0 before the
//    first call.  On each call, the routine increments N_DATA by 1, and
//    returns the corresponding data; when there is no more data, the
//    output value of N_DATA will be 0 again.
//
//    Output, int *A, the parameter of the function.
//
//    Output, double *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 21

  int a_vec[N_MAX] = { 
     1,  2,  1,  2, 
     1,  2,  3,  4, 
     1,  2,  3,  4, 
     5,  3,  3,  3, 
     3,  3, 10, 10,
    10 };
  double fx_vec[N_MAX] = { 
    0.0796557E+00, 0.00498752E+00, 0.112463E+00,    0.00995017E+00, 
    0.472911E+00,  0.181269E+00,   0.0597575E+00,   0.0175231E+00, 
    0.682689E+00,  0.393469E+00,   0.198748E+00,    0.090204E+00, 
    0.0374342E+00, 0.427593E+00,   0.608375E+00,    0.738536E+00, 
    0.828203E+00,  0.88839E+00,    0.000172116E+00, 0.00365985E+00, 
    0.0185759E+00 };
  double x_vec[N_MAX] = { 
    0.01E+00, 0.01E+00, 0.02E+00, 0.02E+00, 
    0.40E+00, 0.40E+00, 0.40E+00, 0.40E+00, 
    1.00E+00, 1.00E+00, 1.00E+00, 1.00E+00, 
    1.00E+00, 2.00E+00, 3.00E+00, 4.00E+00, 
    5.00E+00, 6.00E+00, 1.00E+00, 2.00E+00, 
    3.00E+00 };

  if ( *n_data < 0 )
  {
    *n_data = 0;
  }

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *a = 0;
    *x = 0.0E+00;
    *fx = 0.0E+00;
  }
  else
  {
    *a = a_vec[*n_data-1];
    *x = x_vec[*n_data-1];
    *fx = fx_vec[*n_data-1];
  }
  return;
# undef N_MAX
}
//****************************************************************************80

void cumbet ( double *x, double *y, double *a, double *b, double *cum,
  double *ccum )

//****************************************************************************80
// 
//  Purpose:
// 
//    CUMBET evaluates the cumulative incomplete beta distribution.
//
//  Discussion:
//
//    This routine calculates the CDF to X of the incomplete beta distribution
//    with parameters A and B.  This is the integral from 0 to x
//    of (1/B(a,b))*f(t)) where f(t) = t**(a-1) * (1-t)**(b-1)
//
//  Modified:
//
//    14 March 2006
//
//  Reference:
//
//    A R Didonato and Alfred Morris,
//    Algorithm 708:
//    Significant Digit Computation of the Incomplete Beta Function Ratios. 
//    ACM Transactions on Mathematical Software, 
//    Volume 18, Number 3, September 1992, pages 360-373.
//
//  Parameters:
//
//    Input, double *X, the upper limit of integration.
//
//    Input, double *Y, the value of 1-X.
//
//    Input, double *A, *B, the parameters of the distribution.
//
//    Output, double *CUM, *CCUM, the values of the cumulative
//    density function and complementary cumulative density function.
//
{
  static int ierr;

  if ( *x <= 0.0 )
  {
    *cum = 0.0;
    *ccum = 1.0;
  }
  else if ( *y <= 0.0 ) 
  {
    *cum = 1.0;
    *ccum = 0.0;
  }
  else
  {
    beta_inc ( a, b, x, y, cum, ccum, &ierr );
  }
  return;
}
//****************************************************************************80

void cumbin ( double *s, double *xn, double *pr, double *ompr,
  double *cum, double *ccum )

//****************************************************************************80
// 
//  Purpose:
// 
//    CUMBIN evaluates the cumulative binomial distribution.
//
//  Discussion:
//
//    This routine returns the probability of 0 to S successes in XN binomial
//    trials, each of which has a probability of success, PR.
//
//  Modified:
//
//    14 March 2006
//
//  Reference:
//
//    Milton Abramowitz and Irene Stegun,
//    Handbook of Mathematical Functions 
//    1966, Formula 26.5.24.
//
//  Parameters:
//
//    Input, double *S, the upper limit of summation.
//
//    Input, double *XN, the number of trials.
//
//    Input, double *PR, the probability of success in one trial.
//
//    Input, double *OMPR, equals ( 1 - PR ).
//
//    Output, double *CUM, the cumulative binomial distribution.
//
//    Output, double *CCUM, the complement of the cumulative 
//    binomial distribution.
//
{
  static double T1,T2;

  if ( *s < *xn )
  {
    T1 = *s + 1.0;
    T2 = *xn - *s;
    cumbet ( pr, ompr, &T1, &T2, ccum, cum );
  }
  else
  {
    *cum = 1.0;
    *ccum = 0.0;
  }
  return;
}
//****************************************************************************80

void cumchi ( double *x, double *df, double *cum, double *ccum )

//****************************************************************************80
//
//  Purpose:
// 
//    CUMCHI evaluates the cumulative chi-square distribution.
//
//  Parameters:
//
//    Input, double *X, the upper limit of integration.
//
//    Input, double *DF, the degrees of freedom of the
//    chi-square distribution.
//
//    Output, double *CUM, the cumulative chi-square distribution.
//
//    Output, double *CCUM, the complement of the cumulative 
//    chi-square distribution.
//
{
  static double a;
  static double xx;

  a = *df * 0.5;
  xx = *x * 0.5;
  cumgam ( &xx, &a, cum, ccum );
  return;
}
//****************************************************************************80

void cumchn ( double *x, double *df, double *pnonc, double *cum,
  double *ccum )

//****************************************************************************80
// 
//  Purpose:
// 
//    CUMCHN evaluates the cumulative noncentral chi-square distribution.
//
//  Discussion:
//
//    Calculates the cumulative noncentral chi-square
//    distribution, i.e., the probability that a random variable
//    which follows the noncentral chi-square distribution, with
//    noncentrality parameter PNONC and continuous degrees of
//    freedom DF, is less than or equal to X.
//
//  Reference:
//
//    Milton Abramowitz and Irene Stegun,
//    Handbook of Mathematical Functions 
//    1966, Formula 26.4.25.
//
//  Parameters:
//
//    Input, double *X, the upper limit of integration.
//
//    Input, double *DF, the number of degrees of freedom.
//
//    Input, double *PNONC, the noncentrality parameter of 
//    the noncentral chi-square distribution.
//
//    Output, double *CUM, *CCUM, the CDF and complementary
//    CDF of the noncentral chi-square distribution.
//
//  Local Parameters:
//
//    Local, double EPS, the convergence criterion.  The sum 
//    stops when a term is less than EPS*SUM.
//
//    Local, int NTIRED, the maximum number of terms to be evaluated
//    in each sum.
//
//    Local, bool QCONV, is TRUE if convergence was achieved, that is,
//    the program did not stop on NTIRED criterion.
//
{
# define dg(i) (*df+2.0e0*(double)(i))
# define qsmall(xx) (int)(sum < 1.0e-20 || (xx) < eps*sum)
# define qtired(i) (int)((i) > ntired)

  static double eps = 1.0e-5;
  static int ntired = 1000;
  static double adj,centaj,centwt,chid2,dfd2,lcntaj,lcntwt,lfact,pcent,pterm,sum,
    sumadj,term,wt,xnonc;
  static int i,icent,iterb,iterf;
  static double T1,T2,T3;

    if(!(*x <= 0.0e0)) goto S10;
    *cum = 0.0e0;
    *ccum = 1.0e0;
    return;
S10:
    if(!(*pnonc <= 1.0e-10)) goto S20;
//
//     When non-centrality parameter is (essentially) zero,
//     use cumulative chi-square distribution
//
    cumchi(x,df,cum,ccum);
    return;
S20:
    xnonc = *pnonc/2.0e0;
//
//     The following code calculates the weight, chi-square, and
//     adjustment term for the central term in the infinite series.
//     The central term is the one in which the poisson weight is
//     greatest.  The adjustment term is the amount that must
//     be subtracted from the chi-square to move up two degrees
//     of freedom.
//
    icent = fifidint(xnonc);
    if(icent == 0) icent = 1;
    chid2 = *x/2.0e0;
//
//     Calculate central weight term
//
    T1 = (double)(icent+1);
    lfact = gamma_log ( &T1 );
    lcntwt = -xnonc+(double)icent*log(xnonc)-lfact;
    centwt = exp(lcntwt);
//
//     Calculate central chi-square
//
    T2 = dg(icent);
    cumchi(x,&T2,&pcent,ccum);
//
//     Calculate central adjustment term
//
    dfd2 = dg(icent)/2.0e0;
    T3 = 1.0e0+dfd2;
    lfact = gamma_log ( &T3 );
    lcntaj = dfd2*log(chid2)-chid2-lfact;
    centaj = exp(lcntaj);
    sum = centwt*pcent;
//
//     Sum backwards from the central term towards zero.
//     Quit whenever either
//     (1) the zero term is reached, or
//     (2) the term gets small relative to the sum, or
//     (3) More than NTIRED terms are totaled.
//
    iterb = 0;
    sumadj = 0.0e0;
    adj = centaj;
    wt = centwt;
    i = icent;
    goto S40;
S30:
    if(qtired(iterb) || qsmall(term) || i == 0) goto S50;
S40:
    dfd2 = dg(i)/2.0e0;
//
//     Adjust chi-square for two fewer degrees of freedom.
//     The adjusted value ends up in PTERM.
//
    adj = adj*dfd2/chid2;
    sumadj = sumadj + adj;
    pterm = pcent+sumadj;
//
//     Adjust poisson weight for J decreased by one
//
    wt *= ((double)i/xnonc);
    term = wt*pterm;
    sum = sum + term;
    i -= 1;
    iterb = iterb + 1;
    goto S30;
S50:
    iterf = 0;
//
//     Now sum forward from the central term towards infinity.
//     Quit when either
//     (1) the term gets sma