function algdiv ( a, b )

!*****************************************************************************80
!
!! 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).
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  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, real ( kind = 8 ) A, B, define the arguments.
!
!    Output, real ( kind = 8 ) ALGDIV, the value of ln(Gamma(B)/Gamma(A+B)).
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) algdiv
  real ( kind = 8 ) alnrel
  real ( kind = 8 ) b
  real ( kind = 8 ) c
  real ( kind = 8 ), parameter :: c0 =  0.833333333333333D-01
  real ( kind = 8 ), parameter :: c1 = -0.277777777760991D-02
  real ( kind = 8 ), parameter :: c2 =  0.793650666825390D-03
  real ( kind = 8 ), parameter :: c3 = -0.595202931351870D-03
  real ( kind = 8 ), parameter :: c4 =  0.837308034031215D-03
  real ( kind = 8 ), parameter :: c5 = -0.165322962780713D-02
  real ( kind = 8 ) d
  real ( kind = 8 ) h
  real ( kind = 8 ) s11
  real ( kind = 8 ) s3
  real ( kind = 8 ) s5
  real ( kind = 8 ) s7
  real ( kind = 8 ) s9
  real ( kind = 8 ) t
  real ( kind = 8 ) u
  real ( kind = 8 ) v
  real ( kind = 8 ) w
  real ( kind = 8 ) x
  real ( kind = 8 ) x2

  if ( b < a ) then
    h = b / a
    c = 1.0D+00 / ( 1.0D+00 + h )
    x = h / ( 1.0D+00 + h )
    d = a + ( b - 0.5D+00 )
  else
    h = a / b
    c = h / ( 1.0D+00 + h )
    x = 1.0D+00 / ( 1.0D+00 + h )
    d = b + ( a - 0.5D+00 )
  end if
!
!  Set SN = (1 - X**N)/(1 - X).
!
  x2 = x * x
  s3 = 1.0D+00 + ( x + x2 )
  s5 = 1.0D+00 + ( x + x2 * s3 )
  s7 = 1.0D+00 + ( x + x2 * s5 )
  s9 = 1.0D+00 + ( x + x2 * s7 )
  s11 = 1.0D+00 + ( x + x2 * s9 )
!
!  Set W = DEL(B) - DEL(A + B).
!
  t = ( 1.0D+00 / b )**2
  w = (((( &
          c5 * s11  * t &
        + c4 * s9 ) * t &
        + c3 * s7 ) * t &
        + c2 * s5 ) * t &
        + c1 * s3 ) * t &
        + c0

  w = w * ( c / b )
!
!  Combine the results.
!
  u = d * alnrel ( a / b )
  v = a * ( log ( b ) - 1.0D+00 )

  if ( v < u ) then
    algdiv = ( w - v ) - u
  else
    algdiv = ( w - u ) - v
  end if

  return
end
function alnrel ( a )

!*****************************************************************************80
!
!! ALNREL evaluates the function ln ( 1 + A ).
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  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, real ( kind = 8 ) A, the argument.
!
!    Output, real ( kind = 8 ) ALNREL, the value of ln ( 1 + A ).
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) alnrel
  real ( kind = 8 ), parameter :: p1 = -0.129418923021993D+01
  real ( kind = 8 ), parameter :: p2 =  0.405303492862024D+00
  real ( kind = 8 ), parameter :: p3 = -0.178874546012214D-01
  real ( kind = 8 ), parameter :: q1 = -0.162752256355323D+01
  real ( kind = 8 ), parameter :: q2 =  0.747811014037616D+00
  real ( kind = 8 ), parameter :: q3 = -0.845104217945565D-01
  real ( kind = 8 ) t
  real ( kind = 8 ) t2
  real ( kind = 8 ) w
  real ( kind = 8 ) x

  if ( abs ( a ) <= 0.375D+00 ) then

    t = a / ( a +  2.0D+00  )
    t2 = t * t

    w = ((( p3 * t2 + p2 ) * t2 + p1 ) * t2 + 1.0D+00 ) &
      / ((( q3 * t2 + q2 ) * t2 + q1 ) * t2 + 1.0D+00 )

    alnrel =  2.0D+00  * t * w

  else

    x = 1.0D+00 + real ( a, kind = 8 )
    alnrel = log ( x )

  end if

  return
end
function apser ( a, b, x, eps )

!*****************************************************************************80
!
!! 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.
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  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, real ( kind = 8 ) A, B, X, the parameters of the
!    incomplete beta ratio.
!
!    Input, real ( kind = 8 ) EPS, a tolerance.
!
!    Output, real ( kind = 8 ) APSER, the computed value of the
!    incomplete beta ratio.
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) aj
  real ( kind = 8 ) apser
  real ( kind = 8 ) b
  real ( kind = 8 ) bx
  real ( kind = 8 ) c
  real ( kind = 8 ) eps
  real ( kind = 8 ), parameter :: g = 0.577215664901533D+00
  real ( kind = 8 ) j
  real ( kind = 8 ) psi
  real ( kind = 8 ) s
  real ( kind = 8 ) t
  real ( kind = 8 ) tol
  real ( kind = 8 ) x

  bx = b * x
  t = x - bx

  if ( b * eps <= 0.02D+00 ) then
    c = log ( x ) + psi ( b ) + g + t
  else
    c = log ( bx ) + g + t
  end if

  tol = 5.0D+00 * eps * abs ( c )
  j = 1.0D+00
  s = 0.0D+00

  do

    j = j + 1.0D+00
    t = t * ( x - bx / j )
    aj = t / j
    s = s + aj

    if ( abs ( aj ) <= tol ) then
      exit
    end if

  end do

  apser = -a * ( c + s )

  return
end
function bcorr ( a0, b0 )

!*****************************************************************************80
!
!! 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 ).
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  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, real ( kind = 8 ) A0, B0, the arguments.
!    It is assumed that 8 <= A0 and 8 <= B0.
!
!    Output, real ( kind = 8 ) BCORR, the value of the function.
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) a0
  real ( kind = 8 ) b
  real ( kind = 8 ) b0
  real ( kind = 8 ) bcorr
  real ( kind = 8 ) c
  real ( kind = 8 ), parameter :: c0 =  0.833333333333333D-01
  real ( kind = 8 ), parameter :: c1 = -0.277777777760991D-02
  real ( kind = 8 ), parameter :: c2 =  0.793650666825390D-03
  real ( kind = 8 ), parameter :: c3 = -0.595202931351870D-03
  real ( kind = 8 ), parameter :: c4 =  0.837308034031215D-03
  real ( kind = 8 ), parameter :: c5 = -0.165322962780713D-02
  real ( kind = 8 ) h
  real ( kind = 8 ) s11
  real ( kind = 8 ) s3
  real ( kind = 8 ) s5
  real ( kind = 8 ) s7
  real ( kind = 8 ) s9
  real ( kind = 8 ) t
  real ( kind = 8 ) w
  real ( kind = 8 ) x
  real ( kind = 8 ) x2

  a = min ( a0, b0 )
  b = max ( a0, b0 )

  h = a / b
  c = h / ( 1.0D+00 + h )
  x = 1.0D+00 / ( 1.0D+00 + h )
  x2 = x * x
!
!  Set SN = (1 - X**N)/(1 - X)
!
  s3 = 1.0D+00 + ( x + x2 )
  s5 = 1.0D+00 + ( x + x2 * s3 )
  s7 = 1.0D+00 + ( x + x2 * s5 )
  s9 = 1.0D+00 + ( x + x2 * s7 )
  s11 = 1.0D+00 + ( x + x2 * s9 )
!
!  Set W = DEL(B) - DEL(A + B)
!
  t = ( 1.0D+00 / b )**2

  w = (((( &
       c5 * s11  * t &
     + c4 * s9 ) * t &
     + c3 * s7 ) * t &
     + c2 * s5 ) * t &
     + c1 * s3 ) * t &
     + c0

  w = w * ( c / b )
!
!  Compute  DEL(A) + W.
!
  t = ( 1.0D+00 / a )**2

  bcorr = ((((( &
         c5   * t &
       + c4 ) * t &
       + c3 ) * t &
       + c2 ) * t &
       + c1 ) * t &
       + c0 ) / a + w

  return
end
function beta ( a, b )

!*****************************************************************************80
!
!! BETA evaluates the beta function.
!
!  Modified:
!
!    03 December 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) A, B, the arguments of the beta function.
!
!    Output, real ( kind = 8 ) BETA, the value of the beta function.
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  real ( kind = 8 ) beta
  real ( kind = 8 ) beta_log

  beta = exp ( beta_log ( a, b ) )

  return
end
function beta_asym ( a, b, lambda, eps )

!*****************************************************************************80
!
!! BETA_ASYM computes an asymptotic expansion for IX(A,B), for large A and B.
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  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, real ( kind = 8 ) 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, real ( kind = 8 ) LAMBDA, the value of ( A + B ) * Y - B.
!    It is assumed that 0 <= LAMBDA.
!
!    Input, real ( kind = 8 ) EPS, the tolerance.
!
  implicit none

  integer, parameter :: num = 20

  real ( kind = 8 ) a
  real ( kind = 8 ) a0(num+1)
  real ( kind = 8 ) b
  real ( kind = 8 ) b0(num+1)
  real ( kind = 8 ) bcorr
  real ( kind = 8 ) beta_asym
  real ( kind = 8 ) bsum
  real ( kind = 8 ) c(num+1)
  real ( kind = 8 ) d(num+1)
  real ( kind = 8 ) dsum
  real ( kind = 8 ), parameter :: e0 = 1.12837916709551D+00
  real ( kind = 8 ), parameter :: e1 = 0.353553390593274D+00
  real ( kind = 8 ) eps
  real ( kind = 8 ) error_fc
  real ( kind = 8 ) f
  real ( kind = 8 ) h
  real ( kind = 8 ) h2
  real ( kind = 8 ) hn
  integer i
  integer j
  real ( kind = 8 ) j0
  real ( kind = 8 ) j1
  real ( kind = 8 ) lambda
  integer m
  integer mm1
  integer mmj
  integer n
  integer np1
  real ( kind = 8 ) r
  real ( kind = 8 ) r0
  real ( kind = 8 ) r1
  real ( kind = 8 ) rlog1
  real ( kind = 8 ) s
  real ( kind = 8 ) sum1
  real ( kind = 8 ) t
  real ( kind = 8 ) t0
  real ( kind = 8 ) t1
  real ( kind = 8 ) u
  real ( kind = 8 ) w
  real ( kind = 8 ) w0
  real ( kind = 8 ) z
  real ( kind = 8 ) z0
  real ( kind = 8 ) z2
  real ( kind = 8 ) zn
  real ( kind = 8 ) znm1

  beta_asym = 0.0D+00

  if ( a < b ) then
    h = a / b
    r0 = 1.0D+00 / ( 1.0D+00 + h )
    r1 = ( b - a ) / b
    w0 = 1.0D+00 / sqrt ( a * ( 1.0D+00 + h ))
  else
    h = b / a
    r0 = 1.0D+00 / ( 1.0D+00 + h )
    r1 = ( b - a ) / a
    w0 = 1.0D+00 / sqrt ( b * ( 1.0D+00 + h ))
  end if

  f = a * rlog1 ( - lambda / a ) + b * rlog1 ( lambda / b )
  t = exp ( - f )
  if ( t == 0.0D+00 ) then
    return
  end if

  z0 = sqrt ( f )
  z = 0.5D+00 * ( z0 / e1 )
  z2 = f + f

  a0(1) = ( 2.0D+00 / 3.0D+00 ) * r1
  c(1) = -0.5D+00 * a0(1)
  d(1) = -c(1)
  j0 = ( 0.5D+00 / e0 ) * error_fc ( 1, z0 )
  j1 = e1
  sum1 = j0 + d(1) * w0 * j1

  s = 1.0D+00
  h2 = h * h
  hn = 1.0D+00
  w = w0
  znm1 = z
  zn = z2

  do n = 2, num, 2

    hn = h2 * hn
    a0(n) = 2.0D+00 * r0 * ( 1.0D+00 + h * hn ) &
      / ( n + 2.0D+00 )
    np1 = n + 1
    s = s + hn
    a0(np1) = 2.0D+00 * r1 * s / ( n + 3.0D+00 )

    do i = n, np1

      r = -0.5D+00 * ( i + 1.0D+00 )
      b0(1) = r * a0(1)
      do m = 2, i
        bsum = 0.0D+00
        mm1 = m - 1
        do j = 1, mm1
          mmj = m - j
          bsum = bsum + ( j * r - mmj ) * a0(j) * b0(mmj)
        end do
        b0(m) = r * a0(m) + bsum / m
      end do

      c(i) = b0(i) / ( i + 1.0D+00 )

      dsum = 0.0
      do j = 1, i-1
        dsum = dsum + d(i-j) * c(j)
      end do
      d(i) = - ( dsum + c(i) )

    end do

    j0 = e1 * znm1 + ( n - 1.0D+00 ) * j0
    j1 = e1 * zn + n * j1
    znm1 = z2 * znm1
    zn = z2 * zn
    w = w0 * w
    t0 = d(n) * w * j0
    w = w0 * w
    t1 = d(np1) * w * j1
    sum1 = sum1 + ( t0 + t1 )

    if ( ( abs ( t0 ) + abs ( t1 )) <= eps * sum1 ) then
      u = exp ( - bcorr ( a, b ) )
      beta_asym = e0 * t * u * sum1
      return
    end if

  end do

  u = exp ( - bcorr ( a, b ) )
  beta_asym = e0 * t * u * sum1

  return
end
function beta_frac ( a, b, x, y, lambda, eps )

!*****************************************************************************80
!
!! BETA_FRAC evaluates a continued fraction expansion for IX(A,B).
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  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, real ( kind = 8 ) 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, real ( kind = 8 ) X, Y.  X is the argument of the
!    function, and should satisy 0 <= X <= 1.  Y should equal 1 - X.
!
!    Input, real ( kind = 8 ) LAMBDA, the value of ( A + B ) * Y - B.
!
!    Input, real ( kind = 8 ) EPS, a tolerance.
!
!    Output, real ( kind = 8 ) BETA_FRAC, the value of the continued
!    fraction approximation for IX(A,B).
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) alpha
  real ( kind = 8 ) an
  real ( kind = 8 ) anp1
  real ( kind = 8 ) b
  real ( kind = 8 ) beta
  real ( kind = 8 ) beta_frac
  real ( kind = 8 ) beta_rcomp
  real ( kind = 8 ) bn
  real ( kind = 8 ) bnp1
  real ( kind = 8 ) c
  real ( kind = 8 ) c0
  real ( kind = 8 ) c1
  real ( kind = 8 ) e
  real ( kind = 8 ) eps
  real ( kind = 8 ) lambda
  real ( kind = 8 ) n
  real ( kind = 8 ) p
  real ( kind = 8 ) r
  real ( kind = 8 ) r0
  real ( kind = 8 ) s
  real ( kind = 8 ) t
  real ( kind = 8 ) w
  real ( kind = 8 ) x
  real ( kind = 8 ) y
  real ( kind = 8 ) yp1

  beta_frac = beta_rcomp ( a, b, x, y )

  if ( beta_frac == 0.0D+00 ) then
    return
  end if

  c = 1.0D+00 + lambda
  c0 = b / a
  c1 = 1.0D+00 + 1.0D+00 / a
  yp1 = y + 1.0D+00

  n = 0.0D+00
  p = 1.0D+00
  s = a + 1.0D+00
  an = 0.0D+00
  bn = 1.0D+00
  anp1 = 1.0D+00
  bnp1 = c / c1
  r = c1 / c
!
!  Continued fraction calculation.
!
  do

    n = n + 1.0D+00
    t = n / a
    w = n * ( b - n ) * x
    e = a / s
    alpha = ( p * ( p + c0 ) * e * e ) * ( w * x )
    e = ( 1.0D+00 + t ) / ( c1 + t + t )
    beta = n + w / s + e * ( c + n * yp1 )
    p = 1.0D+00 + t
    s = s +  2.0D+00 
!
!  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 ( abs ( r - r0 ) <= eps * r ) then
      beta_frac = beta_frac * r
      exit
    end if
!
!  Rescale AN, BN, ANP1, and BNP1.
!
    an = an / bnp1
    bn = bn / bnp1
    anp1 = r
    bnp1 = 1.0D+00

  end do

  return
end
subroutine beta_grat ( a, b, x, y, w, eps, ierr )

!*****************************************************************************80
!
!! BETA_GRAT evaluates an asymptotic expansion for IX(A,B).
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  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, real ( kind = 8 ) 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, real ( kind = 8 ) X, Y.  X is the argument of the
!    function, and should satisy 0 <= X <= 1.  Y should equal 1 - X.
!
!    Input/output, real ( kind = 8 ) W, a quantity to which the
!    result of the computation is to be added on output.
!
!    Input, real ( kind = 8 ) EPS, a tolerance.
!
!    Output, integer IERR, an error flag, which is 0 if no error
!    was detected.
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) algdiv
  real ( kind = 8 ) alnrel
  real ( kind = 8 ) b
  real ( kind = 8 ) bm1
  real ( kind = 8 ) bp2n
  real ( kind = 8 ) c(30)
  real ( kind = 8 ) cn
  real ( kind = 8 ) coef
  real ( kind = 8 ) d(30)
  real ( kind = 8 ) dj
  real ( kind = 8 ) eps
  real ( kind = 8 ) gam1
  integer i
  integer ierr
  real ( kind = 8 ) j
  real ( kind = 8 ) l
  real ( kind = 8 ) lnx
  integer n
  real ( kind = 8 ) n2
  real ( kind = 8 ) nu
  real ( kind = 8 ) p
  real ( kind = 8 ) q
  real ( kind = 8 ) r
  real ( kind = 8 ) s
  real ( kind = 8 ) sum1
  real ( kind = 8 ) t
  real ( kind = 8 ) t2
  real ( kind = 8 ) u
  real ( kind = 8 ) v
  real ( kind = 8 ) w
  real ( kind = 8 ) x
  real ( kind = 8 ) y
  real ( kind = 8 ) z

  bm1 = ( b - 0.5D+00 ) - 0.5D+00
  nu = a + 0.5D+00 * bm1

  if ( y <= 0.375D+00 ) then
    lnx = alnrel ( - y )
  else
    lnx = log ( x )
  end if

  z = -nu * lnx

  if ( b * z == 0.0D+00 ) then
    ierr = 1
    return
  end if
!
!  Computation of the expansion.
!
!  Set R = EXP(-Z)*Z**B/GAMMA(B)
!
  r = b * ( 1.0D+00 + gam1 ( b ) ) * exp ( b * log ( z ))
  r = r * exp ( a * lnx ) * exp ( 0.5D+00 * bm1 * lnx )
  u = algdiv ( b, a ) + b * log ( nu )
  u = r * exp ( - u )

  if ( u == 0.0D+00 ) then
    ierr = 1
    return
  end if

  call gamma_rat1 ( b, z, r, p, q, eps )

  v = 0.25D+00 * ( 1.0D+00 / nu )**2
  t2 = 0.25D+00 * lnx * lnx
  l = w / u
  j = q / r
  sum1 = j
  t = 1.0D+00
  cn = 1.0D+00
  n2 = 0.0D+00

  do n = 1, 30

    bp2n = b + n2
    j = ( bp2n * ( bp2n + 1.0D+00 ) * j &
      + ( z + bp2n + 1.0D+00 ) * t ) * v
    n2 = n2 +  2.0D+00 
    t = t * t2
    cn = cn / ( n2 * ( n2 + 1.0D+00 ))
    c(n) = cn
    s = 0.0D+00

    coef = b - n
    do i = 1, n-1
      s = s + coef * c(i) * d(n-i)
      coef = coef + b
    end do

    d(n) = bm1 * cn + s / n
    dj = d(n) * j
    sum1 = sum1 + dj

    if ( sum1 <= 0.0D+00 ) then
      ierr = 1
      return
    end if

    if ( abs ( dj ) <= eps * ( sum1 + l ) ) then
      ierr = 0
      w = w + u * sum1
      return
    end if

  end do

  ierr = 0
  w = w + u * sum1

  return
end
subroutine beta_inc ( a, b, x, y, w, w1, ierr )

!*****************************************************************************80
!
!! BETA_INC evaluates the incomplete beta function IX(A,B).
!
!  Author:
!
!    Alfred Morris,
!    Naval Surface Weapons Center,
!    Dahlgren, Virginia.
!
!  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, real ( kind = 8 ) A, B, the parameters of the function.
!    A and B should be nonnegative.
!
!    Input, real ( kind = 8 ) X, Y.  X is the argument of the
!    function, and should satisy 0 <= X <= 1.  Y should equal 1 - X.
!
!    Output, real ( kind = 8 ) W, W1, the values of IX(A,B) and
!    1-IX(A,B).
!
!    Output, integer 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.
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) a0
  real ( kind = 8 ) apser
  real ( kind = 8 ) b
  real ( kind = 8 ) b0
  real ( kind = 8 ) beta_asym
  real ( kind = 8 ) beta_frac
  real ( kind = 8 ) beta_pser
  real ( kind = 8 ) beta_up
  real ( kind = 8 ) eps
  real ( kind = 8 ) fpser
  integer ierr
  integer ierr1
  integer ind
  real ( kind = 8 ) lambda
  integer n
  real ( kind = 8 ) t
  real ( kind = 8 ) w
  real ( kind = 8 ) w1
  real ( kind = 8 ) x
  real ( kind = 8 ) x0
  real ( kind = 8 ) y
  real ( kind = 8 ) y0
  real ( kind = 8 ) z

  eps = epsilon ( eps )
  w = 0.0D+00
  w1 = 0.0D+00

  if ( a < 0.0D+00 .or. b < 0.0D+00 ) then
    ierr = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'BETA_INC - Fatal error!'
    write ( *, '(a,i8)' ) '  IERR = ', ierr
    return
  end if

  if ( a == 0.0D+00 .and. b == 0.0D+00 ) then
    ierr = 2
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'BETA_INC - Fatal error!'
    write ( *, '(a,i8)' ) '  IERR = ', ierr
    return
  end if

  if ( x < 0.0D+00 .or. 1.0D+00 < x ) then
    ierr = 3
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'BETA_INC - Fatal error!'
    write ( *, '(a,i8)' ) '  IERR = ', ierr
    return
  end if

  if ( y < 0.0D+00 .or. 1.0D+00 < y ) then
    ierr = 4
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'BETA_INC - Fatal error!'
    write ( *, '(a,i8)' ) '  IERR = ', ierr
    return
  end if

  z = ( ( x + y ) - 0.5D+00 ) - 0.5D+00

  if ( 3.0D+00 * eps < abs ( z ) ) then
    ierr = 5
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'BETA_INC - Fatal error!'
    write ( *, '(a,i8)' ) '  IERR = ', ierr
    return
  end if

  ierr = 0

  if ( x == 0.0D+00 ) then
    w = 0.0D+00
    w1 = 1.0D+00
    if ( a == 0.0D+00 ) then
      ierr = 6
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'BETA_INC - Fatal error!'
      write ( *, '(a,i8)' ) '  IERR = ', ierr
    end if
    return
  end if

  if ( y == 0.0D+00 ) then
    if ( b == 0.0D+00 ) then
      ierr = 7
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'BETA_INC - Fatal error!'
      write ( *, '(a,i8)' ) '  IERR = ', ierr
      return
    end if
    w = 1.0D+00
    w1 = 0.0D+00
    return
  end if

  if ( a == 0.0D+00 ) then
    w = 1.0D+00
    w1 = 0.0D+00
    return
  end if

  if ( b == 0.0D+00 ) then
    w = 0.0D+00
    w1 = 1.0D+00
    return
  end if

  eps = max ( eps, 1.0D-15 )

  if ( max ( a, b ) < 0.001D+00 * eps ) then
    go to 260
  end if

  ind = 0
  a0 = a
  b0 = b
  x0 = x
  y0 = y

  if ( 1.0D+00 < min ( a0, b0 ) ) then
    go to 40
  end if
!
!  Procedure for A0 <= 1 or B0 <= 1
!
  if ( 0.5D+00 < x ) then
    ind = 1
    a0 = b
    b0 = a
    x0 = y
    y0 = x
  end if

  if ( b0 < min ( eps, eps * a0 ) ) then
    go to 90
  end if

  if ( a0 < min ( eps, eps * b0 ) .and. b0 * x0 <= 1.0D+00 ) then
    go to 100
  end if

  if ( 1.0D+00 < max ( a0, b0 ) ) then
    go to 20
  end if

  if ( min ( 0.2D+00, b0 ) <= a0 ) then
    go to 110
  end if

  if ( x0**a0 <= 0.9D+00 ) then
    go to 110
  end if

  if ( 0.3D+00 <= x0 ) then
    go to 120
  end if

  n = 20
  go to 140

20 continue

  if ( b0 <= 1.0D+00 ) then
    go to 110
  end if

  if ( 0.3D+00 <= x0 ) then
    go to 120
  end if

  if ( 0.1D+00 <= x0 ) then
    go to 30
  end if

  if ( ( x0 * b0 )**a0 <= 0.7D+00 ) then
    go to 110
  end if

30 continue

  if ( 15.0D+00 < b0 ) then
    go to 150
  end if

  n = 20
  go to 140
!
!  PROCEDURE for 1 < A0 and 1 < B0.
!
40 continue

  if ( a <= b ) then
    lambda = a - ( a + b ) * x
  else
    lambda = ( a + b ) * y - b
  end if

  if ( lambda < 0.0D+00 ) then
    ind = 1
    a0 = b
    b0 = a
    x0 = y
    y0 = x
    lambda = abs ( lambda )
  end if

70 continue

  if ( b0 < 40.0D+00 .and. b0 * x0 <= 0.7D+00 ) then
    go to 110
  end if

  if ( b0 < 40.0D+00 ) then
    go to 160
  end if

  if ( b0 < a0 ) then
    go to 80
  end if

  if ( a0 <= 100.0D+00 ) then
    go to 130
  end if

  if ( 0.03D+00 * a0 < lambda ) then
    go to 130
  end if

  go to 200

80 continue

  if ( b0 <= 100.0D+00 ) then
    go to 130
  end if

  if ( 0.03D+00 * b0 < lambda ) then
    go to 130
  end if

  go to 200
!
!  Evaluation of the appropriate algorithm.
!
90 continue

  w = fpser ( a0, b0, x0, eps )
  w1 = 0.5D+00 + ( 0.5D+00 - w )
  go to 250

100 continue

  w1 = apser ( a0, b0, x0, eps )
  w = 0.5D+00 + ( 0.5D+00 - w1 )
  go to 250

110 continue

  w = beta_pser ( a0, b0, x0, eps )
  w1 = 0.5D+00 + ( 0.5D+00 - w )
  go to 250

120 continue

  w1 = beta_pser ( b0, a0, y0, eps )
  w = 0.5D+00 + ( 0.5D+00 - w1 )
  go to 250

130 continue

  w = beta_frac ( a0, b0, x0, y0, lambda, 15.0D+00 * eps )
  w1 = 0.5D+00 + ( 0.5D+00 - w )
  go to 250

140 continue

  w1 = beta_up ( b0, a0, y0, x0, n, eps )
  b0 = b0 + n

150 continue

  call beta_grat ( b0, a0, y0, x0, w1, 15.0D+00 * eps, ierr1 )
  w = 0.5D+00 + ( 0.5D+00 - w1 )
  go to 250

160 continue

  n = b0
  b0 = b0 - n

  if ( b0 == 0.0D+00 ) then
    n = n - 1
    b0 = 1.0D+00
  end if

170 continue

  w = beta_up ( b0, a0, y0, x0, n, eps )

  if ( x0 <= 0.7D+00 ) then
    w = w + beta_pser ( a0, b0, x0, eps )
    w1 = 0.5D+00 + ( 0.5D+00 - w )
    go to 250
  end if

  if ( a0 <= 15.0D+00 ) then
    n = 20
    w = w + beta_up ( a0, b0, x0, y0, n, eps )
    a0 = a0 + n
  end if

190 continue

  call beta_grat ( a0, b0, x0, y0, w, 15.0D+00 * eps, ierr1 )
  w1 = 0.5D+00 + ( 0.5D+00 - w )
  go to 250

200 continue

  w = beta_asym ( a0, b0, lambda, 100.0D+00 * eps )
  w1 = 0.5D+00 + ( 0.5D+00 - w )
  go to 250
!
!  Termination of the procedure.
!
250 continue

  if ( ind /= 0 ) then
    t = w
    w = w1
    w1 = t
  end if

  return
!
!  Procedure for A and B < 0.001 * EPS
!
260 continue

  w = b / ( a + b )
  w1 = a / ( a + b )

  return
end
subroutine beta_inc_values ( n_data, a, b, x, fx )

!*****************************************************************************80
!
!! 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:
!
!    17 February 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz, 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, integer 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, real ( kind = 8 ) A, B, X, the arguments of the function.
!
!    Output, real ( kind = 8 ) FX, the value of the function.
!
  implicit none

  integer, parameter :: n_max = 30

  real ( kind = 8 ) a
  real ( kind = 8 ), save, dimension ( n_max ) :: a_vec = (/ &
     0.5D+00,  0.5D+00,  0.5D+00,  1.0D+00, &
     1.0D+00,  1.0D+00,  1.0D+00,  1.0D+00, &
     2.0D+00,  2.0D+00,  2.0D+00,  2.0D+00, &
     2.0D+00,  2.0D+00,  2.0D+00,  2.0D+00, &
     2.0D+00,  5.5D+00, 10.0D+00, 10.0D+00, &
    10.0D+00, 10.0D+00, 20.0D+00, 20.0D+00, &
    20.0D+00, 20.0D+00, 20.0D+00, 30.0D+00, &
    30.0D+00, 40.0D+00 /)
  real ( kind = 8 ) b
  real ( kind = 8 ), save, dimension ( n_max ) :: b_vec = (/ &
     0.5D+00,  0.5D+00,  0.5D+00,  0.5D+00, &
     0.5D+00,  0.5D+00,  0.5D+00,  1.0D+00, &
     2.0D+00,  2.0D+00,  2.0D+00,  2.0D+00, &
     2.0D+00,  2.0D+00,  2.0D+00,  2.0D+00, &
     2.0D+00,  5.0D+00,  0.5D+00,  5.0D+00, &
     5.0D+00, 10.0D+00,  5.0D+00, 10.0D+00, &
    10.0D+00, 20.0D+00, 20.0D+00, 10.0D+00, &
    10.0D+00, 20.0D+00 /)
  real ( kind = 8 ) fx
  real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ &
    0.0637686D+00, 0.2048328D+00, 1.0000000D+00, 0.0D+00,       &
    0.0050126D+00, 0.0513167D+00, 0.2928932D+00, 0.5000000D+00, &
    0.028D+00,     0.104D+00,     0.216D+00,     0.352D+00,     &
    0.500D+00,     0.648D+00,     0.784D+00,     0.896D+00,     &
    0.972D+00,     0.4361909D+00, 0.1516409D+00, 0.0897827D+00, &
    1.0000000D+00, 0.5000000D+00, 0.4598773D+00, 0.2146816D+00, &
    0.9507365D+00, 0.5000000D+00, 0.8979414D+00, 0.2241297D+00, &
    0.7586405D+00, 0.7001783D+00 /)
  integer n_data
  real ( kind = 8 ) x
  real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ &
    0.01D+00, 0.10D+00, 1.00D+00, 0.0D+00,  &
    0.01D+00, 0.10D+00, 0.50D+00, 0.50D+00, &
    0.1D+00,  0.2D+00,  0.3D+00,  0.4D+00,  &
    0.5D+00,  0.6D+00,  0.7D+00,  0.8D+00,  &
    0.9D+00,  0.50D+00, 0.90D+00, 0.50D+00, &
    1.00D+00, 0.50D+00, 0.80D+00, 0.60D+00, &
    0.80D+00, 0.50D+00, 0.60D+00, 0.70D+00, &
    0.80D+00, 0.70D+00 /)

  if ( n_data < 0 ) then
    n_data = 0
  end if

  n_data = n_data + 1

  if ( n_max < n_data ) then
    n_data = 0
    a = 0.0D+00
    b = 0.0D+00
    x = 0.0D+00
    fx = 0.0D+00
  else
    a = a_vec(n_data)
    b = b_vec(n_data)
    x = x_vec(n_data)
    fx = fx_vec(n_data)
  end if

  return
end
function beta_log ( a0, b0 )

!*****************************************************************************80
!
!! BETA_LOG evaluates the logarithm of the beta function.
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  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, real ( kind = 8 ) A0, B0, the parameters of the function.
!    A0 and B0 should be nonnegative.
!
!    Output, real ( kind = 8 ) BETA_LOG, the value of the logarithm
!    of the Beta function.
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) a0
  real ( kind = 8 ) algdiv
  real ( kind = 8 ) alnrel
  real ( kind = 8 ) b
  real ( kind = 8 ) b0
  real ( kind = 8 ) bcorr
  real ( kind = 8 ) beta_log
  real ( kind = 8 ) c
  real ( kind = 8 ), parameter :: e = 0.918938533204673D+00
  real ( kind = 8 ) gamma_log
  real ( kind = 8 ) gsumln
  real ( kind = 8 ) h
  integer i
  integer n
  real ( kind = 8 ) u
  real ( kind = 8 ) v
  real ( kind = 8 ) w
  real ( kind = 8 ) z

  a = min ( a0, b0 )
  b = max ( a0, b0 )
!
!  8 < A.
!
  if ( 8.0D+00 <= a ) then

    w = bcorr ( a, b )
    h = a / b
    c = h / ( 1.0D+00 + h )
    u = - ( a - 0.5D+00 ) * log ( c )
    v = b * alnrel ( h )

    if ( v < u ) then
      beta_log = ((( -0.5D+00 * log ( b ) + e ) + w ) - v ) - u
    else
      beta_log = ((( -0.5D+00 * log ( b ) + e ) + w ) - u ) - v
    end if

    return
  end if
!
!  Procedure when A < 1
!
  if ( a < 1.0D+00 ) then

    if ( b < 8.0D+00 ) then
      beta_log = gamma_log ( a ) + ( gamma_log ( b ) - gamma_log ( a + b ) )
    else
      beta_log = gamma_log ( a ) + algdiv ( a, b )
    end if

    return

  end if
!
!  Procedure when 1 <= A < 8
!
  if ( 2.0D+00 < a ) then
    go to 40
  end if

  if ( b <= 2.0D+00 ) then
    beta_log = gamma_log ( a ) + gamma_log ( b ) - gsumln ( a, b )
    return
  end if

  w = 0.0D+00

  if ( b < 8.0D+00 ) then
    go to 60
  end if

  beta_log = gamma_log ( a ) + algdiv ( a, b )
  return

40 continue
!
!  Reduction of A when 1000 < B.
!
  if ( 1000.0D+00 < b ) then

    n = a - 1.0D+00
    w = 1.0D+00
    do i = 1, n
      a = a - 1.0D+00
      w = w * ( a / ( 1.0D+00 + a / b ))
    end do

    beta_log = ( log ( w ) - n * log ( b ) ) &
      + ( gamma_log ( a ) + algdiv ( a, b ) )

    return
  end if

  n = a - 1.0D+00
  w = 1.0D+00
  do i = 1, n
    a = a - 1.0D+00
    h = a / b
    w = w * ( h / ( 1.0D+00 + h ) )
  end do
  w = log ( w )

  if ( 8.0D+00 <= b ) then
    beta_log = w + gamma_log ( a ) + algdiv ( a, b )
    return
  end if
!
!  Reduction of B when B < 8.
!
60 continue

  n = b - 1.0D+00
  z = 1.0D+00
  do i = 1, n
    b = b - 1.0D+00
    z = z * ( b / ( a + b ))
  end do

  beta_log = w + log ( z ) + ( gamma_log ( a ) + ( gamma_log ( b ) &
    - gsumln ( a, b ) ) )

  return
end
function beta_pser ( a, b, x, eps )

!*****************************************************************************80
!
!! 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.
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  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, real ( kind = 8 ) A, B, the parameters.
!
!    Input, real ( kind = 8 ) X, the point where the function
!    is to be evaluated.
!
!    Input, real ( kind = 8 ) EPS, the tolerance.
!
!    Output, real ( kind = 8 ) BETA_PSER, the approximate value of IX(A,B)(X).
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) a0
  real ( kind = 8 ) algdiv
  real ( kind = 8 ) apb
  real ( kind = 8 ) b
  real ( kind = 8 ) b0
  real ( kind = 8 ) beta_log
  real ( kind = 8 ) beta_pser
  real ( kind = 8 ) c
  real ( kind = 8 ) eps
  real ( kind = 8 ) gam1
  real ( kind = 8 ) gamma_ln1
  integer i
  integer m
  real ( kind = 8 ) n
  real ( kind = 8 ) sum1
  real ( kind = 8 ) t
  real ( kind = 8 ) tol
  real ( kind = 8 ) u
  real ( kind = 8 ) w
  real ( kind = 8 ) x
  real ( kind = 8 ) z

  beta_pser = 0.0D+00

  if ( x == 0.0D+00 ) then
    return
  end if
!
!  Compute the factor X**A/(A*BETA(A,B))
!
  a0 = min ( a, b )

  if ( 1.0D+00 <= a0 ) then

    z = a * log ( x ) - beta_log ( a, b )
    beta_pser = exp ( z ) / a
  
  else

    b0 = max ( a, b )

    if ( b0 <= 1.0D+00 ) then

      beta_pser = x**a
      if ( beta_pser == 0.0D+00 ) then
        return
      end if

      apb = a + b

      if ( apb <= 1.0D+00 ) then
        z = 1.0D+00 + gam1 ( apb )
      else
        u = a + b - 1.0D+00
        z = ( 1.0D+00 + gam1 ( u ) ) / apb
      end if

      c = ( 1.0D+00 + gam1 ( a ) ) &
        * ( 1.0D+00 + gam1 ( b ) ) / z
      beta_pser = beta_pser * c * ( b / apb )

    else if ( b0 < 8.0D+00 ) then

      u = gamma_ln1 ( a0 )
      m = b0 - 1.0D+00

      c = 1.0D+00
      do i = 1, m
        b0 = b0 - 1.0D+00
        c = c * ( b0 / ( a0 + b0 ))
      end do

      u = log ( c ) + u
      z = a * log ( x ) - u
      b0 = b0 - 1.0D+00
      apb = a0 + b0

      if ( apb <= 1.0D+00 ) then
        t = 1.0D+00 + gam1 ( apb )
      else
        u = a0 + b0 - 1.0D+00
        t = ( 1.0D+00 + gam1 ( u ) ) / apb
      end if

      beta_pser = exp ( z ) * ( a0 / a ) &
        * ( 1.0D+00 + gam1 ( b0 )) / t

    else if ( 8.0D+00 <= b0 ) then

      u = gamma_ln1 ( a0 ) + algdiv ( a0, b0 )
      z = a * log ( x ) - u
      beta_pser = ( a0 / a ) * exp ( z )

    end if

  end if

  if ( beta_pser == 0.0D+00 .or. a <= 0.1D+00 * eps ) then
    return
  end if
!
!  Compute the series.
!
  sum1 = 0.0D+00
  n = 0.0D+00
  c = 1.0D+00
  tol = eps / a

  do

    n = n + 1.0D+00
    c = c * ( 0.5D+00 + ( 0.5D+00 - b / n ) ) * x
    w = c / ( a + n )
    sum1 = sum1 + w

    if ( abs ( w ) <= tol ) then
      exit
    end if

  end do

  beta_pser = beta_pser * ( 1.0D+00 + a * sum1 )

  return
end
function beta_rcomp ( a, b, x, y )

!*****************************************************************************80
!
!! BETA_RCOMP evaluates X**A * Y**B / Beta(A,B).
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  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, real ( kind = 8 ) A, B, the parameters of the Beta function.
!    A and B should be nonnegative.
!
!    Input, real ( kind = 8 ) X, Y, define the numerator of the fraction.
!
!    Output, real ( kind = 8 ) BETA_RCOMP, the value of X**A * Y**B / Beta(A,B).
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) a0
  real ( kind = 8 ) algdiv
  real ( kind = 8 ) alnrel
  real ( kind = 8 ) apb
  real ( kind = 8 ) b
  real ( kind = 8 ) b0
  real ( kind = 8 ) bcorr
  real ( kind = 8 ) beta_log
  real ( kind = 8 ) beta_rcomp
  real ( kind = 8 ) c
  real ( kind = 8 ), parameter :: const = 0.398942280401433D+00
  real ( kind = 8 ) e
  real ( kind = 8 ) gam1
  real ( kind = 8 ) gamma_ln1
  real ( kind = 8 ) h
  integer i
  real ( kind = 8 ) lambda
  real ( kind = 8 ) lnx
  real ( kind = 8 ) lny
  integer n
  real ( kind = 8 ) rlog1
  real ( kind = 8 ) t
  real ( kind = 8 ) u
  real ( kind = 8 ) v
  real ( kind = 8 ) x
  real ( kind = 8 ) x0
  real ( kind = 8 ) y
  real ( kind = 8 ) y0
  real ( kind = 8 ) z

  beta_rcomp = 0.0D+00
  if ( x == 0.0D+00 .or. y == 0.0D+00 ) then
    return
  end if

  a0 = min ( a, b )

  if ( a0 < 8.0D+00 ) then

    if ( x <= 0.375D+00 ) then
      lnx = log ( x )
      lny = alnrel ( - x )      
    else if ( y <= 0.375D+00 ) then
      lnx = alnrel ( - y )
      lny = log ( y )
    else
      lnx = log ( x )
      lny = log ( y )
    end if

    z = a * lnx + b * lny

    if ( 1.0D+00 <= a0 ) then
      z = z - beta_log ( a, b )
      beta_rcomp = exp ( z )
      return
    end if
!
!  Procedure for A < 1 or B < 1
!
    b0 = max ( a, b )

    if ( b0 <= 1.0D+00 ) then

      beta_rcomp = exp ( z )
      if ( beta_rcomp == 0.0D+00 ) then
        return
      end if

      apb = a + b

      if ( apb <= 1.0D+00 ) then
        z = 1.0D+00 + gam1 ( apb )
      else
        u = a + b - 1.0D+00
        z = ( 1.0D+00 + gam1 ( u ) ) / apb
      end if

      c = ( 1.0D+00 + gam1 ( a ) ) &
        * ( 1.0D+00 + gam1 ( b ) ) / z
      beta_rcomp = beta_rcomp * ( a0 * c ) &
        / ( 1.0D+00 + a0 / b0 )

    else if ( b0 < 8.0D+00 ) then

      u = gamma_ln1 ( a0 )
      n = b0 - 1.0D+00

      c = 1.0D+00
      do i = 1, n
        b0 = b0 - 1.0D+00
        c = c * ( b0 / ( a0 + b0 ))
      end do
      u = log ( c ) + u

      z = z - u
      b0 = b0 - 1.0D+00
      apb = a0 + b0

      if ( apb <= 1.0D+00 ) then
        t = 1.0D+00 + gam1 ( apb )
      else
        u = a0 + b0 - 1.0D+00
        t = ( 1.0D+00 + gam1 ( u ) ) / apb
      end if

      beta_rcomp = a0 * exp ( z ) * ( 1.0D+00 + gam1 ( b0 ) ) / t

    else if ( 8.0D+00 <= b0 ) then

      u = gamma_ln1 ( a0 ) + algdiv ( a0, b0 )
      beta_rcomp = a0 * exp ( z - u ) 

    end if

  else

    if ( a <= b ) then
      h = a / b
      x0 = h / ( 1.0D+00 + h )
      y0 = 1.0D+00 / (  1.0D+00 + h )
      lambda = a - ( a + b ) * x
    else
      h = b / a
      x0 = 1.0D+00 / ( 1.0D+00 + h )
      y0 = h / ( 1.0D+00 + h )
      lambda = ( a + b ) * y - b
    end if

    e = -lambda / a

    if ( abs ( e ) <= 0.6D+00 ) then
      u = rlog1 ( e )
    else
      u = e - log ( x / x0 )
    end if

    e = lambda / b
 
    if ( abs ( e ) <= 0.6D+00 ) then
      v = rlog1 ( e )
    else
      v = e - log ( y / y0 )
    end if

    z = exp ( - ( a * u + b * v ) )
    beta_rcomp = const * sqrt ( b * x0 ) * z * exp ( - bcorr ( a, b ))

  end if

  return
end
function beta_rcomp1 ( mu, a, b, x, y )

!*****************************************************************************80
!
!! BETA_RCOMP1 evaluates exp(MU) * X**A * Y**B / Beta(A,B).
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  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, integer MU, ?
!
!    Input, real ( kind = 8 ) A, B, the parameters of the Beta function.
!    A and B should be nonnegative.
!
!    Input, real ( kind = 8 ) X, Y, quantities whose powers form part of
!    the expression.
!
!    Output, real ( kind = 8 ) BETA_RCOMP1, the value of
!    exp(MU) * X**A * Y**B / Beta(A,B).
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) a0
  real ( kind = 8 ) algdiv
  real ( kind = 8 ) alnrel
  real ( kind = 8 ) apb
  real ( kind = 8 ) b
  real ( kind = 8 ) b0
  real ( kind = 8 ) bcorr
  real ( kind = 8 ) beta_log
  real ( kind = 8 ) beta_rcomp1
  real ( kind = 8 ) c
  real ( kind = 8 ), parameter :: const = 0.398942280401433D+00
  real ( kind = 8 ) e
  real ( kind = 8 ) esum
  real ( kind = 8 ) gam1
  real ( kind = 8 ) gamma_ln1
  real ( kind = 8 ) h
  integer i
  real ( kind = 8 ) lambda
  real ( kind = 8 ) lnx
  real ( kind = 8 ) lny
  integer mu
  integer n
  real ( kind = 8 ) rlog1
  real ( kind = 8 ) t
  real ( kind = 8 ) u
  real ( kind = 8 ) v
  real ( kind = 8 ) x
  real ( kind = 8 ) x0
  real ( kind = 8 ) y
  real ( kind = 8 ) y0
  real ( kind = 8 ) z

  a0 = min ( a, b )
!
!  Procedure for 8 <= A and 8 <= B.
!
  if ( 8.0D+00 <= a0 ) then

    if ( a <= b ) then
      h = a / b
      x0 = h / ( 1.0D+00 + h )
      y0 = 1.0D+00 / ( 1.0D+00 + h )
      lambda = a - ( a + b ) * x
    else
      h = b / a
      x0 = 1.0D+00 / ( 1.0D+00 + h )
      y0 = h / ( 1.0D+00 + h )
      lambda = ( a + b ) * y - b
    end if

    e = -lambda / a

    if ( abs ( e ) <= 0.6D+00 ) then
      u = rlog1 ( e )
    else
      u = e - log ( x / x0 )
    end if

    e = lambda / b

    if ( abs ( e ) <= 0.6D+00 ) then
      v = rlog1 ( e )
    else
      v = e - log ( y / y0 )
    end if

    z = esum ( mu, - ( a * u + b * v ))
    beta_rcomp1 = const * sqrt ( b * x0 ) * z * exp ( - bcorr ( a, b ) )
!
!  Procedure for A < 8 or B < 8.
!
  else

    if ( x <= 0.375D+00 ) then
      lnx = log ( x )
      lny = alnrel ( - x )
    else if ( y <= 0.375D+00 ) then
      lnx = alnrel ( - y )
      lny = log ( y )
    else
      lnx = log ( x )
      lny = log ( y )
    end if
  
    z = a * lnx + b * lny

    if ( 1.0D+00 <= a0 ) then
      z = z - beta_log ( a, b )
      beta_rcomp1 = esum ( mu, z )
      return
    end if
!
!  Procedure for A < 1 or B < 1.
!
    b0 = max ( a, b )

    if ( 8.0D+00 <= b0 ) then
      u = gamma_ln1 ( a0 ) + algdiv ( a0, b0 )
      beta_rcomp1 = a0 * esum ( mu, z-u )
      return
    end if

    if ( 1.0D+00 < b0 ) then
!
!  Algorithm for 1 < B0 < 8
!
      u = gamma_ln1 ( a0 )
      n = b0 - 1.0D+00

      c = 1.0D+00
      do i = 1, n
        b0 = b0 - 1.0D+00
        c = c * ( b0 / ( a0 + b0 ) )
      end do
      u = log ( c ) + u

      z = z - u
      b0 = b0 - 1.0D+00
      apb = a0 + b0

      if ( apb <= 1.0D+00 ) then
        t = 1.0D+00 + gam1 ( apb )
      else
        u = a0 + b0 - 1.0D+00
        t = ( 1.0D+00 + gam1 ( u ) ) / apb
      end if

      beta_rcomp1 = a0 * esum ( mu, z ) &
        * ( 1.0D+00 + gam1 ( b0 ) ) / t
!
!  Algorithm for B0 <= 1
!
    else

      beta_rcomp1 = esum ( mu, z )
      if ( beta_rcomp1 == 0.0D+00 ) then
        return
      end if

      apb = a + b

      if ( apb <= 1.0D+00 ) then
        z = 1.0D+00 + gam1 ( apb )
      else
        u = real ( a, kind = 8 ) + real ( b, kind = 8 ) - 1.0D+00
        z = ( 1.0D+00 + gam1 ( u )) / apb
      end if

      c = ( 1.0D+00 + gam1 ( a ) ) &
        * ( 1.0D+00 + gam1 ( b ) ) / z
      beta_rcomp1 = beta_rcomp1 * ( a0 * c ) / ( 1.0D+00 + a0 / b0 )

    end if

  end if

  return
end
function beta_up ( a, b, x, y, n, eps )

!*****************************************************************************80
!
!! BETA_UP evaluates IX(A,B) - IX(A+N,B) where N is a positive integer.
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  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, real ( kind = 8 ) A, B, the parameters of the function.
!    A and B should be nonnegative.
!
!    Input, real ( kind = 8 ) X, Y, ?
!
!    Input, integer N, the increment to the first argument of IX.
!
!    Input, real ( kind = 8 ) EPS, the tolerance.
!
!    Output, real ( kind = 8 ) BETA_UP, the value of IX(A,B) - IX(A+N,B).
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) ap1
  real ( kind = 8 ) apb
  real ( kind = 8 ) b
  real ( kind = 8 ) beta_rcomp1
  real ( kind = 8 ) beta_up
  real ( kind = 8 ) d
  real ( kind = 8 ) eps
  real ( kind = 8 ) exparg
  integer i
  integer k
  real ( kind = 8 ) l
  integer mu
  integer n
  real ( kind = 8 ) r
  real ( kind = 8 ) t
  real ( kind = 8 ) w
  real ( kind = 8 ) x
  real ( kind = 8 ) y
!
!  Obtain the scaling factor EXP(-MU) AND
!  EXP(MU)*(X**A*Y**B/BETA(A,B))/A
!
  apb = a + b
  ap1 = a + 1.0D+00
  mu = 0
  d = 1.0D+00

  if ( n /= 1 ) then

    if ( 1.0D+00 <= a ) then

      if ( 1.1D+00 * ap1 <= apb ) then
        mu = abs ( exparg ( 1 ) )
        k = exparg ( 0 )
        if ( k < mu ) then
          mu = k
        end if
        t = mu
        d = exp ( - t )
      end if

    end if

  end if

  beta_up = beta_rcomp1 ( mu, a, b, x, y ) / a

  if ( n == 1 .or. beta_up == 0.0D+00 ) then
    return
  end if

  w = d
!
!  Let K be the index of the maximum term.
!
  k = 0

  if ( 1.0D+00 < b ) then

    if ( y <= 0.0001D+00 ) then

      k = n - 1

    else

      r = ( b - 1.0D+00 ) * x / y - a

      if ( 1.0D+00 <= r ) then
        k = n - 1
        t = n - 1
        if ( r < t ) then
          k = r
        end if
      end if

    end if
!
!  Add the increasing terms of the series.
!
    do i = 1, k
      l = i - 1
      d = ( ( apb + l ) / ( ap1 + l ) ) * x * d
      w = w + d
    end do

  end if
!
!  Add the remaining terms of the series.
!
  do i = k+1, n-1
    l = i - 1
    d = ( ( apb + l ) / ( ap1 + l ) ) * x * d
    w = w + d
    if ( d <= eps * w ) then
      beta_up = beta_up * w
      return
    end if
  end do

  beta_up = beta_up * w

  return
end
subroutine binomial_cdf_values ( n_data, a, b, x, fx )

!*****************************************************************************80
!
!! 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:
!
!    27 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz, 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, integer 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, integer A, real ( kind = 8 ) B, integer X, the arguments 
!    of the function.
!
!    Output, real ( kind = 8 ) FX, the value of the function.
!
  implicit none

  integer, parameter :: n_max = 17

  integer a
  integer, save, dimension ( n_max ) :: a_vec = (/ &
     2,  2,  2,  2, &
     2,  4,  4,  4, &
     4, 10, 10, 10, &
    10, 10, 10, 10, &
    10 /)
  real ( kind = 8 ) b
  real ( kind = 8 ), save, dimension ( n_max ) :: b_vec = (/ &
    0.05D+00, 0.05D+00, 0.05D+00, 0.50D+00, &
    0.50D+00, 0.25D+00, 0.25D+00, 0.25D+00, &
    0.25D+00, 0.05D+00, 0.10D+00, 0.15D+00, &
    0.20D+00, 0.25D+00, 0.30D+00, 0.40D+00, &
    0.50D+00 /)
  real ( kind = 8 ) fx
  real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ &
    0.9025D+00, 0.9975D+00, 1.0000D+00, 0.2500D+00, &
    0.7500D+00, 0.3164D+00, 0.7383D+00, 0.9492D+00, &
    0.9961D+00, 0.9999D+00, 0.9984D+00, 0.9901D+00, &
    0.9672D+00, 0.9219D+00, 0.8497D+00, 0.6331D+00, &
    0.3770D+00 /)
  integer n_data
  integer x
  integer, save, dimension ( n_max ) :: x_vec = (/ &
     0, 1, 2, 0, &
     1, 0, 1, 2, &
     3, 4, 4, 4, &
     4, 4, 4, 4, &
     4 /)

  if ( n_data < 0 ) then
    n_data = 0
  end if

  n_data = n_data + 1

  if ( n_max < n_data ) then
    n_data = 0
    a = 0
    b = 0.0D+00
    x = 0
    fx = 0.0D+00
  else
    a = a_vec(n_data)
    b = b_vec(n_data)
    x = x_vec(n_data)
    fx = fx_vec(n_data)
  end if

  return
end
subroutine cdfbet ( which, p, q, x, y, a, b, status, bound )

!*****************************************************************************80
!
!! 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:
!
!    14 April 2007
!
!  Author:
!
!    Barry Brown, James Lovato, Kathy Russell
!
!  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, integer 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, real ( kind = 8 ) P, the integral from 0 to X of the
!    chi-square distribution.  Input range: [0, 1].
!
!    Input/output, real ( kind = 8 ) Q, equals 1-P.  Input range: [0, 1].
!
!    Input/output, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, integer 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, real ( kind = 8 ) 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.
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ), parameter :: atol = 1.0D-10
  real ( kind = 8 ) b
  real ( kind = 8 ) bound
  real ( kind = 8 ) ccum
  real ( kind = 8 ) cum
  real ( kind = 8 ) fx
  real ( kind = 8 ), parameter :: inf = 1.0D+300
  real ( kind = 8 ) p
  real ( kind = 8 ) q
  logical qhi
  logical qleft
  integer status
  real ( kind = 8 ), parameter :: tol = 1.0D-08
  integer which
  real ( kind = 8 ) x
  real ( kind = 8 ) xhi
  real ( kind = 8 ) xlo
  real ( kind = 8 ) y

  status = 0
  bound = 0.0D+00

  if ( which < 1 ) then
    bound = 1.0D+00
    status = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CDFBET - Fatal error!'
    write ( *, '(a)' ) '  The input parameter WHICH is out of range.'
    write ( *, '(a)' ) '  Legal values are between 1 and 4.'
    return
  end if

  if ( 4 < which ) then
    bound = 4.0D+00
    status = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CDFBET - Fatal error!'
    write ( *, '(a)' ) '  The input parameter WHICH is out of range.'
    write ( *, '(a)' ) '  Legal values are between 1 and 4.'
    return
  end if
!
!  Unless P is to be computed, make sure it is legal.
!
  if ( which /= 1 ) then
    if ( p < 0.0D+00 ) then
      bound = 0.0D+00
      status = -2
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter P is out of range.'
      return
    else if ( 1.0D+00 < p ) then
      bound = 1.0D+00
      status = -2
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter P is out of range.'
      return
    end if
  end if
!
!  Unless Q is to be computed, make sure it is legal.
!
  if ( which /= 1 ) then
    if ( q < 0.0D+00 ) then
      bound = 0.0D+00
      status = -3
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter Q is out of range.'
      return
    else if ( 1.0D+00 < q ) then
      bound = 1.0D+00
      status = -3
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter Q is out of range.'
      return
    end if
  end if
!
!  Unless X is to be computed, make sure it is legal.
!
  if ( which /= 2 ) then
    if ( x < 0.0D+00 ) then
      bound = 0.0D+00
      status = -4
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter X is out of range.'
      return
    else if ( 1.0D+00 < x ) then
      bound = 1.0D+00
      status = -4
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter X is out of range.'
      return
    end if
  end if
!
!  Unless Y is to be computed, make sure it is legal.
!
  if ( which /= 2 ) then
    if ( y < 0.0D+00 ) then
      bound = 0.0D+00
      status = -5
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter Y is out of range.'
      return
    else if ( 1.0D+00 < y ) then
      bound = 1.0D+00
      status = -5
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter Y is out of range.'
      return
    end if
  end if
!
!  Unless A is to be computed, make sure it is legal.
!
  if ( which /= 3 ) then
    if ( a <= 0.0D+00 ) then
      bound = 0.0D+00
      status = -6
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter A is out of range.'
      return
    end if
  end if
!
!  Unless B is to be computed, make sure it is legal.
!
  if ( which /= 4 ) then
    if ( b <= 0.0D+00 ) then
      bound = 0.0D+00
      status = -7
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter B is out of range.'
      return
    end if
  end if
!
!  Check that P + Q = 1.
!
  if ( which /= 1 ) then
    if ( 3.0D+00 * epsilon ( p ) < abs ( ( p + q ) - 1.0D+00 ) ) then
      status = 3
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  P + Q /= 1.'
      return
    end if
  end if
!
!  Check that X + Y = 1.
!
  if ( which /= 2 ) then
    if ( 3.0D+00 * epsilon ( x ) < abs ( ( x + y ) - 1.0D+00 ) ) then
      status = 4
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  X + Y /= 1.'
      return
    end if
  end if
!
!  Compute P and Q.
!
  if ( which == 1 ) then

    call cumbet ( x, y, a, b, p, q )
    status = 0
!
!  Compute X and Y.
!
  else if ( which == 2 ) then

    call dstzr ( 0.0D+00, 1.0D+00, atol, tol )

    if ( p <= q ) then

      status = 0
      fx = 0.0D+00
      call dzror ( status, x, fx, xlo, xhi, qleft, qhi )
      y = 1.0D+00 - x

      do while ( status == 1 )
 
        call cumbet ( x, y, a, b, cum, ccum )
        fx = cum - p
        call dzror ( status, x, fx, xlo, xhi, qleft, qhi )
        y = 1.0D+00 - x
 
      end do

    else

      status = 0
      fx = 0.0D+00
      call dzror ( status, y, fx, xlo, xhi, qleft, qhi )
      x = 1.0D+00 - y

      do while ( status == 1 )

        call cumbet ( x, y, a, b, cum, ccum )
        fx = ccum - q
        call dzror ( status, y, fx, xlo, xhi, qleft, qhi )
        x = 1.0D+00 - y

      end do

    end if

    if ( status == -1 ) then
      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFBET - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      else
        status = 2
        bound = 1.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFBET - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      end if
    end if
!
!  Compute A.
!
  else if ( which == 3 ) then

    call dstinv ( 0.0D+00, inf, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    a = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, a, fx, qleft, qhi )

    do while ( status == 1 )

      call cumbet ( x, y, a, b, cum, ccum )

      if ( p <= q ) then
        fx = cum - p
      else
        fx = ccum - q
      end if

      call dinvr ( status, a, fx, qleft, qhi )

    end do

    if ( status == -1 ) then

      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFBET - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      else
        status = 2
        bound = inf
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFBET - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      end if

    end if
!
!  Compute B.
!
  else if ( which == 4 ) then

    call dstinv ( 0.0D+00, inf, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    b = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, b, fx, qleft, qhi )

    do while ( status == 1 )

      call cumbet ( x, y, a, b, cum, ccum )

      if ( p <= q ) then
        fx = cum - p
      else
        fx = ccum - q
      end if

      call dinvr ( status, b, fx, qleft, qhi )

    end do

    if ( status == -1 ) then
      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFBET - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      else
        status = 2
        bound = inf
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFBET - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      end if
    end if

  end if

  return
end
subroutine cdfbin ( which, p, q, s, xn, pr, ompr, status, bound )

!*****************************************************************************80
!
!! 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:
!
!    14 April 2007
!
!  Author:
!
!    Barry Brown, James Lovato, Kathy Russell
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,  
!    Handbook of Mathematical Functions 
!    1966, Formula 26.5.24.
!
!  Parameters:
!
!    Input, integer 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) S, the number of successes observed.
!    Whether this is an input or output value, it should lie in the
!    range [0,XN].
!
!    Input/output, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, integer 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, real ( kind = 8 ) 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.
!
  implicit none

  real ( kind = 8 ), parameter :: atol = 1.0D-10
  real ( kind = 8 ) bound
  real ( kind = 8 ) ccum
  real ( kind = 8 ) cum
  real ( kind = 8 ) fx
  real ( kind = 8 ), parameter :: inf = 1.0D+300
  real ( kind = 8 ) ompr
  real ( kind = 8 ) p
  real ( kind = 8 ) pr
  real ( kind = 8 ) q
  logical qhi
  logical qleft
  real ( kind = 8 ) s
  integer status
  real ( kind = 8 ), parameter :: tol = 1.0D-08
  integer which
  real ( kind = 8 ) xhi
  real ( kind = 8 ) xlo
  real ( kind = 8 ) xn

  status = 0
  bound = 0.0D+00
!
!  Check the arguments.
!
  if ( which < 1 ) then
    bound = 1.0D+00
    status = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CDFBET - Fatal error!'
    write ( *, '(a)' ) '  The input parameter WHICH is out of range.'
    write ( *, '(a)' ) '  Legal values are between 1 and 4.'
    return
  end if

  if ( 4 < which ) then
    bound = 4.0D+00
    status = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CDFBET - Fatal error!'
    write ( *, '(a)' ) '  The input parameter WHICH is out of range.'
    write ( *, '(a)' ) '  Legal values are between 1 and 4.'
    return
  end if
!
!  Unless P is to be computed, make sure it is legal.
!
  if ( which /= 1 ) then
    if ( p < 0.0D+00 ) then
      status = -2
      bound = 0.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter P is out of range.'
      return
    else if ( 1.0D+00 < p ) then
      status = -2
      bound = 1.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter P is out of range.'
      return
    end if
  end if
!
!  Unless Q is to be computed, make sure it is legal.
!
  if ( which /= 1 ) then
    if ( q < 0.0D+00 ) then
      status = -3
      bound = 0.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter Q is out of range.'
      return
    else if ( 1.0D+00 < q ) then
      status = -3
      bound = 1.0
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter Q is out of range.'
      return
    end if
  end if
!
!  Unless XN is to be computed, make sure it is legal.
!
  if ( which /= 3 ) then
    if ( xn <= 0.0D+00 ) then
      bound = 0.0D+00
      status = -5
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter XN is out of range.'
      return
    end if
  end if
!
!  Unless S is to be computed, make sure it is legal.
!
  if ( which /= 2 ) then
    if ( s < 0.0D+00 ) then
      bound = 0.0D+00
      status = -4
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter S is out of range.'
      return
    else if ( which /= 3 .and. xn < s ) then
      bound = xn
      status = -4
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter S is out of range.'
      return
    end if
  end if
!
!  Unless PR is to be computed, make sure it is legal.
!
  if ( which /= 4 ) then
    if ( pr < 0.0D+00 ) then
      bound = 0.0D+00
      status = -6
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter PR is out of range.'
      return
    else if ( 1.0D+00 < pr ) then
      bound = 1.0D+00
      status = -6
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter PR is out of range.'
      return
    end if
  end if
!
!  Unless OMPR is to be computed, make sure it is legal.
!
  if ( which /= 4 ) then
    if ( ompr < 0.0D+00 ) then
      bound = 0.0D+00
      status = -7
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter OMPR is out of range.'
      return
    else if ( 1.0D+00 < ompr ) then
      bound = 1.0D+00
      status = -7
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  Input parameter OMPR is out of range.'
      return
    end if
  end if
!
!  Check that P + Q = 1.
!
  if ( which /= 1 ) then
    if ( 3.0D+00 * epsilon ( p ) < abs ( ( p + q ) - 1.0D+00 ) ) then
      status = 3
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  P + Q /= 1.'
      return
    end if
  end if
!
!  Check that PR + OMPR = 1.
!
  if ( which /= 4 ) then
    if ( 3.0D+00 * epsilon ( 1.0D+00 ) &
      < abs ( ( pr + ompr ) - 1.0D+00 ) ) then
      status = 4
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFBET - Fatal error!'
      write ( *, '(a)' ) '  PR + OMPR /= 1.'
      return
    end if
  end if
!
!  Calculate P and Q.
!
  if ( which == 1 ) then

    call cumbin ( s, xn, pr, ompr, p, q )
    status = 0
!
!  Calculate S.
!
  else if ( which == 2 ) then

    call dstinv ( 0.0D+00, xn, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    s = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, s, fx, qleft, qhi )

    do while ( status == 1 )

      call cumbin ( s, xn, pr, ompr, cum, ccum )

      if ( p <= q ) then
        fx = cum - p
      else
        fx = ccum - q   
      end if

      call dinvr ( status, s, fx, qleft, qhi )

    end do

    if ( status == -1 ) then
      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFBIN - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      else
        status = 2
        bound = xn
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFBIN - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      end if
    end if
!
!  Calculate XN.
!
  else if ( which == 3 ) then

    call dstinv ( 0.0D+00, inf, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    xn = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, xn, fx, qleft, qhi )

    do while ( status == 1 )

      call cumbin ( s, xn, pr, ompr, cum, ccum )

      if ( p <= q ) then
        fx = cum - p
      else
        fx = ccum - q
      end if

      call dinvr ( status, xn, fx, qleft, qhi )

    end do

    if ( status == -1 ) then
      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFBIN - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
        return
      else
        status = 2
        bound = inf
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFBIN - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
        return
      end if
    end if
!
!  Calculate PR and OMPR.
!
  else if ( which == 4 ) then

    call dstzr ( 0.0D+00, 1.0D+00, atol, tol )

    if ( p <= q ) then

      status = 0
      call dzror ( status, pr, fx, xlo, xhi, qleft, qhi )
      ompr = 1.0D+00 - pr

      do while ( status == 1 )

        call cumbin ( s, xn, pr, ompr, cum, ccum )
        fx = cum - p
        call dzror ( status, pr, fx, xlo, xhi, qleft, qhi )
        ompr = 1.0D+00 - pr

      end do

    else

      status = 0
      call dzror ( status, ompr, fx, xlo, xhi, qleft, qhi )
      pr = 1.0D+00 - ompr

      do while ( status == 1 )

        call cumbin ( s, xn, pr, ompr, cum, ccum )
        fx = ccum - q
        call dzror ( status, ompr, fx, xlo, xhi, qleft, qhi )
        pr = 1.0D+00 - ompr

      end do

    end if

    if ( status == -1 ) then
      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFBIN - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      else
        status = 2
        bound = 1.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFBIN - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      end if
    end if

  end if

  return
end
subroutine cdfchi ( which, p, q, x, df, status, bound )

!*****************************************************************************80
!
!! 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 ]
!
!  Modified:
!
!    14 April 2007
!
!  Author:
!
!    Barry Brown, James Lovato, Kathy Russell
!
!  Reference:
!
!    Milton Abramowitz, 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, integer 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, integer 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, real ( kind = 8 ) 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.
!
  implicit none

  real ( kind = 8 ), parameter :: atol = 1.0D-10
  real ( kind = 8 ) bound
  real ( kind = 8 ) ccum
  real ( kind = 8 ) cum
  real ( kind = 8 ) df
  real ( kind = 8 ) fx
  real ( kind = 8 ), parameter :: inf = 1.0D+300
  real ( kind = 8 ) p
  real ( kind = 8 ) porq
  real ( kind = 8 ) q
  logical qhi
  logical qleft
  integer status
  real ( kind = 8 ), parameter :: tol = 1.0D-08
  integer which
  real ( kind = 8 ) x

  status = 0
  bound = 0.0D+00
!
!  Check the arguments.
!
  if ( which < 1 ) then
    bound = 1.0
    status = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CDFCHI - Fatal error!'
    write ( *, '(a)' ) '  The input parameter WHICH is out of range.'
    write ( *, '(a)' ) '  Legal values are between 1 and 3.'
    return
  end if

  if ( 3 < which ) then
    bound = 3.0
    status = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CDFCHI - Fatal error!'
    write ( *, '(a)' ) '  The input parameter WHICH is out of range.'
    write ( *, '(a)' ) '  Legal values are between 1 and 3.'
    return
  end if
!
!  Unless P is to be computed, make sure it is legal.
!
  if ( which /= 1 ) then
    if ( p < 0.0D+00 ) then
      status = -2
      bound = 0.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFCHI - Fatal error!'
      write ( *, '(a)' ) '  Input parameter P is out of range.'
      return
    else if ( 1.0D+00 < p ) then
      status = -2
      bound = 1.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFCHI - Fatal error!'
      write ( *, '(a)' ) '  Input parameter P is out of range.'
      return
    end if
  end if
!
!  Unless Q is to be computed, make sure it is legal.
!
  if ( which /= 1 ) then
    if ( q < 0.0D+00 ) then
      status = -3
      bound = 0.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFCHI - Fatal error!'
      write ( *, '(a)' ) '  Input parameter Q is out of range.'
      return
    else if ( 1.0D+00 < q ) then
      status = -3
      bound = 1.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFCHI - Fatal error!'
      write ( *, '(a)' ) '  Input parameter Q is out of range.'
      return
    end if
  end if
!
!  Unless X is to be computed, make sure it is legal.
!
  if ( which /= 2 ) then
    if ( x < 0.0D+00 ) then
      bound = 0.0D+00
      status = -4
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFCHI - Fatal error!'
      write ( *, '(a)' ) '  Input parameter X is out of range.'
      return
    end if
  end if
!
!  Unless DF is to be computed, make sure it is legal.
!
  if ( which /= 3 ) then
    if ( df <= 0.0D+00 ) then
      bound = 0.0D+00
      status = -5
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFCHI - Fatal error!'
      write ( *, '(a)' ) '  Input parameter DF is out of range.'
      return
    end if
  end if
!
!  Check that P + Q = 1.
!
  if ( which /= 1 ) then
    if ( 3.0D+00 * epsilon ( p ) < abs ( ( p + q ) - 1.0D+00 ) ) then
      status = 3
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFCHI - Fatal error!'
      write ( *, '(a)' ) '  P + Q /= 1.'
      return
    end if
  end if
!
!  Select the minimum of P and Q.
!
  if ( which /= 1 ) then
    porq = min ( p, q )
  end if
!
!  Calculate P and Q.
!
  if ( which == 1 ) then

    status = 0
    call cumchi ( x, df, p, q )

    if ( 1.5D+00 < porq ) then
      status = 10
      return
    end if
!
!  Calculate X.
!
  else if ( which == 2 ) then

    call dstinv ( 0.0D+00, inf, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    x = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, x, fx, qleft, qhi )

    do while ( status == 1 )

      call cumchi ( x, df, cum, ccum )

      if ( p <= q ) then
        fx = cum - p
      else
        fx = ccum - q
      end if

      if ( 1.5D+00 < fx + porq ) then
        status = 10
        return
      end if

      call dinvr ( status, x, fx, qleft, qhi )

    end do

    if ( status == -1 ) then
      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFCHI - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      else
        status = 2
        bound = inf
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFCHI - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      end if

    end if
!
!  Calculate DF.
!
  else if ( which == 3 ) then

    call dstinv ( 0.0D+00, inf, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    df = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, df, fx, qleft, qhi )

    do while ( status == 1 )

      call cumchi ( x, df, cum, ccum )

      if ( p <= q ) then
        fx = cum - p
      else
        fx = ccum - q
      end if

      if ( 1.5D+00 < fx + porq ) then
        status = 10
        return
      end if

      call dinvr ( status, df, fx, qleft, qhi )

    end do

    if ( status == -1 ) then
      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFCHI - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      else
        status = 2
        bound = inf
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFCHI - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      end if
    end if

  end if

  return
end
subroutine cdfchn ( which, p, q, x, df, pnonc, status, bound )

!*****************************************************************************80
!
!! 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 ]
!
!  Modified:
!
!    14 April 2007
!
!  Author:
!
!    Barry Brown, James Lovato, Kathy Russell
!
!  Reference:
!
!    Milton Abramowitz, 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, integer 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, integer 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, real ( kind = 8 ) 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.
!
  implicit none

  real ( kind = 8 ), parameter :: atol=1.0D-50
  real ( kind = 8 ) bound
  real ( kind = 8 ) ccum
  real ( kind = 8 ) cum
  real ( kind = 8 ) df
  real ( kind = 8 ) fx
  real ( kind = 8 ), parameter :: inf=1.0D+300
  real ( kind = 8 ) p
  real ( kind = 8 ) pnonc
  real ( kind = 8 ) q
  logical qhi
  logical qleft
  integer status
  real ( kind = 8 ), parameter :: tent4=1.0D+04
  real ( kind = 8 ), parameter :: tol=1.0D-08
  integer which
  real ( kind = 8 ) x

  status = 0
  bound = 0.0D+00
!
!  Check the arguments.
!
  if ( which < 1 ) then
    bound = 1.0D+00
    status = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CDFCHN - Fatal error!'
    write ( *, '(a)' ) '  The input parameter WHICH is out of range.'
    write ( *, '(a)' ) '  Legal values are between 1 and 4.'
    return
  end if

  if ( 4 < which ) then
    bound = 4.0D+00
    status = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CDFCHN - Fatal error!'
    write ( *, '(a)' ) '  The input parameter WHICH is out of range.'
    write ( *, '(a)' ) '  Legal values are between 1 and 4.'
    return
  end if
!
!  Unless P is to be computed, make sure it is legal.
!
  if ( which /= 1 ) then
    if ( p < 0.0D+00 ) then
      status = -2
      bound = 0.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFCHN - Fatal error!'
      write ( *, '(a)' ) '  Input parameter P is out of range.'
      return
    else if ( 1.0D+00 < p ) then
      status = -2
      bound = 1.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFCHN - Fatal error!'
      write ( *, '(a)' ) '  Input parameter P is out of range.'
      return
    end if
  end if
!
!  Unless X is to be computed, make sure it is legal.
!
  if ( which /= 2 ) then
    if ( x < 0.0D+00 ) then
      bound = 0.0D+00
      status = -4
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFCHN - Fatal error!'
      write ( *, '(a)' ) '  Input parameter X is out of range.'
      return
    end if
  end if
!
!  Unless DF is to be computed, make sure it is legal.
!
  if ( which /= 3 ) then
    if ( df <= 0.0D+00 ) then
      bound = 0.0D+00
      status = -5
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFCHN - Fatal error!'
      write ( *, '(a)' ) '  Input parameter DF is out of range.'
      return
    end if
  end if
!
!  Unless PNONC is to be computed, make sure it is legal.
!
  if ( which /= 4 ) then
    if ( pnonc < 0.0D+00 ) then
      bound = 0.0D+00
      status = -6
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFCHN - Fatal error!'
      write ( *, '(a)' ) '  Input parameter PNONC is out of range.'
      return
    end if
  end if
!
!  Calculate P and Q.
!
  if ( which == 1 ) then

    call cumchn ( x, df, pnonc, p, q )
    status = 0
!
!  Calculate X.
!
  else if ( which == 2 ) then

    call dstinv ( 0.0D+00, inf, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    x = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, x, fx, qleft, qhi )

    do while ( status == 1 )

      call cumchn ( x, df, pnonc, cum, ccum )
      fx = cum - p
      call dinvr ( status, x, fx, qleft, qhi )

    end do

    if ( status == -1 ) then
      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFCHN - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      else
        status = 2
        bound = inf
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFCHN - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      end if
    end if
!
!  Calculate DF.
!
  else if ( which == 3 ) then

    call dstinv ( 0.0D+00, inf, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    df = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, df, fx, qleft, qhi )

    do while ( status == 1 )

      call cumchn ( x, df, pnonc, cum, ccum )
      fx = cum - p
      call dinvr ( status, df, fx, qleft, qhi )

    end do

    if ( status == -1 )then
      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFCHN - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      else
        status = 2
        bound = inf
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFCHN - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      end if
    end if
!
!  Calculate PNONC.
!
  else if ( which == 4 ) then

    call dstinv ( 0.0D+00, tent4, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    pnonc = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, pnonc, fx, qleft, qhi )

    do while ( status == 1 )

      call cumchn ( x, df, pnonc, cum, ccum )
      fx = cum - p
      call dinvr ( status, pnonc, fx, qleft, qhi )

    end do

    if ( status == -1 ) then
      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFCHN - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      else
        status = 2
        bound = tent4
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFCHN - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      end if
    end if

  end if

  return
end
subroutine cdff ( which, p, q, f, dfn, dfd, status, bound )

!*****************************************************************************80
!
!! 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 degrees 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
!
!  Author:
!
!    Barry Brown, James Lovato, Kathy Russell
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,  
!    Handbook of Mathematical Functions 
!    1966, Formula 26.6.2.
!
!  Parameters:
!
!    Input, integer 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, integer 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, real ( kind = 8 ) 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.
!
  implicit none

  real ( kind = 8 ), parameter :: atol = 1.0D-10
  real ( kind = 8 ) bound
  real ( kind = 8 ) bound_hi
  real ( kind = 8 ) bound_lo
  real ( kind = 8 ) ccum
  real ( kind = 8 ) cum
  real ( kind = 8 ) dfd
  real ( kind = 8 ) dfn
  real ( kind = 8 ) f
  real ( kind = 8 ) fx
  real ( kind = 8 ), parameter :: inf = 1.0D+300
  real ( kind = 8 ) p
  real ( kind = 8 ) q
  logical qhi
  logical qleft
  integer status
  real ( kind = 8 ), parameter :: tol = 1.0D-08
  integer which

  status = 0
  bound = 0.0D+00
!
!  Check the arguments.
!
  if ( which < 1 ) then
    bound = 1.0D+00
    status = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CDFF - Fatal error!'
    write ( *, '(a)' ) '  The input parameter WHICH is out of range.'
    write ( *, '(a)' ) '  Legal values are between 1 and 4.'
    return
  end if

  if ( 4 < which ) then
    bound = 4.0D+00
    status = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CDFF - Fatal error!'
    write ( *, '(a)' ) '  The input parameter WHICH is out of range.'
    write ( *, '(a)' ) '  Legal values are between 1 and 4.'
    return
  end if
!
!  Unless P is to be computed, make sure it is legal.
!
  if ( which /= 1 ) then
    if ( p < 0.0D+00 ) then
      status = -2
      bound = 0.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFF - Fatal error!'
      write ( *, '(a)' ) '  Input parameter P is out of range.'
      return
    else if ( 1.0D+00 < p ) then
      status = -2
      bound = 1.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFF - Fatal error!'
      write ( *, '(a)' ) '  Input parameter P is out of range.'
      return
    end if
  end if
!
!  Unless Q is to be computed, make sure it is legal.
!
  if ( which /= 1 ) then
    if ( q < 0.0D+00 ) then
      status = -3
      bound = 0.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFF - Fatal error!'
      write ( *, '(a)' ) '  Input parameter Q is out of range.'
      return
    else if ( 1.0D+00 < q ) then
      status = -3
      bound = 1.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFF - Fatal error!'
      write ( *, '(a)' ) '  Input parameter Q is out of range.'
      return
    end if
  end if
!
!  Unless F is to be computed, make sure it is legal.
!
  if ( which /= 2 ) then
    if ( f < 0.0D+00 ) then
      bound = 0.0D+00
      status = -4
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFF - Fatal error!'
      write ( *, '(a)' ) '  Input parameter F is out of range.'
      return
    end if
  end if
!
!  Unless DFN is to be computed, make sure it is legal.
!
  if ( which /= 3 ) then
    if ( dfn <= 0.0D+00 ) then
      bound = 0.0D+00
      status = -5
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFF - Fatal error!'
      write ( *, '(a)' ) '  Input parameter DFN is out of range.'
      return
    end if
  end if
!
!  Unless DFD is to be computed, make sure it is legal.
!
  if ( which /= 4 ) then
    if ( dfd <= 0.0D+00 ) then
      bound = 0.0D+00
      status = -6
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFF - Fatal error!'
      write ( *, '(a)' ) '  Input parameter DFD is out of range.'
      return
    end if
  end if
!
!  Check that P + Q = 1.
!
  if ( which /= 1 ) then
    if ( 3.0D+00 * epsilon ( p ) < abs ( ( p + q ) - 1.0D+00 ) ) then
      status = 3
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFF - Fatal error!'
      write ( *, '(a)' ) '  P + Q /= 1.'
      return
    end if
  end if
!
!  Calculate P and Q.
!
  if ( which == 1 ) then

    call cumf ( f, dfn, dfd, p, q )
    status = 0
!
!  Calculate F.
!
  else if ( which == 2 ) then

    call dstinv ( 0.0D+00, inf, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    f = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, f, fx, qleft, qhi )

    do while ( status == 1 )

      call cumf ( f, dfn, dfd, cum, ccum )

      if ( p <= q ) then
        fx = cum - p
      else
        fx = ccum - q
      end if

      call dinvr ( status, f, fx, qleft, qhi )

    end do

    if ( status == -1 ) then

      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFF - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      else
        status = 2
        bound = inf
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFF - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      end if

    end if
!
!  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 ( which == 3 ) then

    bound_lo = 1.0D+00
    bound_hi = inf

    call dstinv ( bound_lo, bound_hi, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    dfn = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, dfn, fx, qleft, qhi )

    do while ( status == 1 )

      call cumf ( f, dfn, dfd, cum, ccum )

      if ( p <= q ) then
        fx = cum - p
      else
        fx = ccum - q
      end if
      
      call dinvr ( status, dfn, fx, qleft, qhi )

    end do

    if ( status == -1 ) then

      if ( qleft ) then
        status = 1
        bound = bound_lo
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFF - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound_lo
        return
      else
        status = 2
        bound = bound_hi
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFF - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound_hi
        return
      end if

    end if
!
!  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 ( which == 4 ) then

    bound_lo = 1.0D+00
    bound_hi = inf

    call dstinv ( bound_lo, bound_hi, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    dfd = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, dfd, fx, qleft, qhi )

    do while ( status == 1 )

      call cumf ( f, dfn, dfd, cum, ccum )

      if ( p <= q ) then
        fx = cum - p
      else
        fx = ccum - q
      end if

      call dinvr ( status, dfd, fx, qleft, qhi )

    end do

    if ( status == -1 ) then
      if ( qleft ) then
        status = 1
        bound = bound_lo
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFF - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound_lo
      else
        status = 2
        bound = bound_hi
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFF - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound_hi
      end if
    end if

  end if

  return
end
subroutine cdffnc ( which, p, q, f, dfn, dfd, pnonc, status, bound )

!*****************************************************************************80
!
!! CDFFNC evaluates the CDF of the Noncentral F distribution.
!
!  Discussion:
!
!    This routine originally used 1.0D+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:
!
!    14 April 2007
!
!  Author:
!
!    Barry Brown, James Lovato, Kathy Russell
!
!  Reference:
!
!    Milton Abramowitz, 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, integer 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, integer 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, real ( kind = 8 ) 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.
!
  implicit none

  real ( kind = 8 ), parameter :: atol = 1.0D-10
  real ( kind = 8 ) bound
  real ( kind = 8 ) ccum
  real ( kind = 8 ) cum
  real ( kind = 8 ) dfd
  real ( kind = 8 ) dfn
  real ( kind = 8 ) f
  real ( kind = 8 ) fx
  real ( kind = 8 ), parameter :: inf = 1.0D+30
  real ( kind = 8 ) p
  real ( kind = 8 ) pnonc
  real ( kind = 8 ) q
  logical qhi
  logical qleft
  integer status
  real ( kind = 8 ), parameter :: tent4 = 1.0D+04
  real ( kind = 8 ), parameter :: tol = 1.0D-08
  integer which

  status = 0
  bound = 0.0D+00
!
!  Check the arguments.
!
  if ( which < 1 ) then
    bound = 1.0D+00
    status = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CDFFNC - Fatal error!'
    write ( *, '(a)' ) '  The input parameter WHICH is out of range.'
    write ( *, '(a)' ) '  Legal values are between 1 and 5.'
    return
  end if

  if ( 5 < which ) then
    bound = 5.0D+00
    status = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CDFFNC - Fatal error!'
    write ( *, '(a)' ) '  The input parameter WHICH is out of range.'
    write ( *, '(a)' ) '  Legal values are between 1 and 5.'
    return
  end if
!
!  Unless P is to be computed, make sure it is legal.
!
  if ( which /= 1 ) then
    if ( p < 0.0D+00 ) then
      status = -2
      bound = 0.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFFNC - Fatal error!'
      write ( *, '(a)' ) '  Input parameter P is out of range.'
      return
    else if ( 1.0D+00 < p ) then
      status = -2
      bound = 1.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFFNC - Fatal error!'
      write ( *, '(a)' ) '  Input parameter P is out of range.'
      return
    end if
  end if
!
!  Unless F is to be computed, make sure it is legal.
!
  if ( which /= 2 ) then
    if ( f < 0.0D+00 ) then
      bound = 0.0D+00
      status = -4
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFFNC - Fatal error!'
      write ( *, '(a)' ) '  Input parameter F is out of range.'
      return
    end if
  end if
!
!  Unless DFN is to be computed, make sure it is legal.
!
  if ( which /= 3 ) then
    if ( dfn <= 0.0D+00 ) then
      bound = 0.0D+00
      status = -5
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFFNC - Fatal error!'
      write ( *, '(a)' ) '  Input parameter DFN is out of range.'
      return
    end if
  end if
!
!  Unless DFD is to be computed, make sure it is legal.
!
  if ( which /= 4 ) then
    if ( dfd <= 0.0D+00 ) then
      bound = 0.0D+00
      status = -6
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFFNC - Fatal error!'
      write ( *, '(a)' ) '  Input parameter DFD is out of range.'
      return
    end if
  end if
!
!  Unless PNONC is to be computed, make sure it is legal.
!
  if ( which /= 5 ) then
    if ( pnonc < 0.0D+00 ) then
      bound = 0.0D+00
      status = -7
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFFNC - Fatal error!'
      write ( *, '(a)' ) '  Input parameter PNONC is out of range.'
      return
    end if
  end if
!
! Calculate P and Q.
!
  if ( which == 1 ) then

    call cumfnc ( f, dfn, dfd, pnonc, p, q )
    status = 0
!
!  Calculate F.
!
  else if ( which == 2 ) then

    call dstinv ( 0.0D+00, inf, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0  
    f = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, f, fx, qleft, qhi )

    do while ( status == 1 )
 
      call cumfnc ( f, dfn, dfd, pnonc, cum, ccum )
      fx = cum - p
      call dinvr ( status, f, fx, qleft, qhi )

    end do

    if ( status == -1 ) then
      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFFNC - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
        return
      else
        status = 2
        bound = inf
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFFNC - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
        return
      end if
    end if
!
!  Calculate DFN.
!
  else if ( which == 3 ) then

    call dstinv ( 1.0D+00, inf, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    dfn = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, dfn, fx, qleft, qhi )

    do while ( status == 1 )

      call cumfnc ( f, dfn, dfd, pnonc, cum, ccum )
      fx = cum - p

      call dinvr ( status, dfn, fx, qleft, qhi )

    end do

    if ( status == -1 ) then
      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFFNC - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      else
        status = 2
        bound = inf
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFFNC - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      end if
    end if
!
!  Calculate DFD.
!
  else if ( which == 4 ) then

    call dstinv ( 1.0D+00, inf, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    dfd = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, dfd, fx, qleft, qhi )

    do while ( status == 1 )

      call cumfnc ( f, dfn, dfd, pnonc, cum, ccum )
      fx = cum - p
      call dinvr ( status, dfd, fx, qleft, qhi )

    end do

    if ( status == -1 ) then
      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFFNC - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      else
        status = 2
        bound = inf
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFFNC - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      end if
    end if
!
!  Calculate PNONC.
!
  else if ( which == 5 ) then

    call dstinv ( 0.0D+00, tent4, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    pnonc = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, pnonc, fx, qleft, qhi )

    do while ( status == 1 )

      call cumfnc ( f, dfn, dfd, pnonc, cum, ccum )
      fx = cum - p

      call dinvr ( status, pnonc, fx, qleft, qhi )

    end do

    if ( status == -1 ) then
      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFFNC - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      else
        status = 2
        bound = tent4
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFFNC - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      end if
    end if

  end if

  return
end
subroutine cdfgam ( which, p, q, x, shape, scale, status, bound )

!*****************************************************************************80
!
!! 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)
!
!  Modified:
!
!    14 April 2007
!
!  Author:
!
!    Barry Brown, James Lovato, Kathy Russell
!
!  Reference:
!
!    Armido DiDinato, Alfred Morris, 
!    Algorithm 654:
!    Computation of the incomplete gamma function ratios and their inverse,
!    ACM Transactions on Mathematical Software,
!    Volume 12, 1986, pages 377-393.
!
!  Parameters:
!
!    Input, integer 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, integer 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, real ( kind = 8 ) 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.
!
  implicit none

  real ( kind = 8 ), parameter :: atol = 1.0D-10
  real ( kind = 8 ) bound
  real ( kind = 8 ) ccum
  real ( kind = 8 ) cum
  real ( kind = 8 ) fx
  integer ierr
  real ( kind = 8 ), parameter :: inf=1.0D+300
  real ( kind = 8 ) p
  real ( kind = 8 ) porq
  real ( kind = 8 ) q
  logical qhi
  logical qleft
  real ( kind = 8 ) scale
  real ( kind = 8 ) shape
  real ( kind = 8 ), parameter :: tol = 1.0D-08
  integer status,which
  real ( kind = 8 ) x
  real ( kind = 8 ) xscale
  real ( kind = 8 ) xx

  status = 0
  bound = 0.0D+00
!
!  Check the arguments.
!
  if ( which < 1 ) then
    bound = 1.0D+00
    status = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CDFGAM - Fatal error!'
    write ( *, '(a)' ) '  The input parameter WHICH is out of range.'
    write ( *, '(a)' ) '  Legal values are between 1 and 4.'
    return
  end if

  if ( 4 < which ) then
    bound = 4.0D+00
    status = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CDFGAM - Fatal error!'
    write ( *, '(a)' ) '  The input parameter WHICH is out of range.'
    write ( *, '(a)' ) '  Legal values are between 1 and 4.'
    return
  end if
!
!  Unless P is to be computed, make sure it is legal.
!
  if ( which /= 1 ) then
    if ( p < 0.0D+00 ) then
      status = -2
      bound = 0.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFGAM - Fatal error!'
      write ( *, '(a)' ) '  Input parameter P is out of range.'
      return
    else if ( 1.0D+00 < p ) then
      status = -2
      bound = 1.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFGAM - Fatal error!'
      write ( *, '(a)' ) '  Input parameter P is out of range.'
      return
    end if
  end if
!
!  Unless Q is to be computed, make sure it is legal.
!
  if ( which /= 1 ) then
    if ( q < 0.0D+00 ) then
      status = -3
      bound = 0.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFGAM - Fatal error!'
      write ( *, '(a)' ) '  Input parameter Q is out of range.'
      return
    else if ( 1.0D+00 < q ) then
      status = -3
      bound = 1.0
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFGAM - Fatal error!'
      write ( *, '(a)' ) '  Input parameter Q is out of range.'
      return
    end if
  end if
!
!  Unless X is to be computed, make sure it is legal.
!
  if ( which /= 2 ) then
    if ( x < 0.0D+00 ) then
      bound = 0.0D+00
      status = -4
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFGAM - Fatal error!'
      write ( *, '(a)' ) '  Input parameter X is out of range.'
      return
    end if
  end if
!
!  Unless SHAPE is to be computed, make sure it is legal.
!
  if ( which /= 3 ) then
    if ( shape <= 0.0D+00 ) then
      bound = 0.0D+00
      status = -5
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFGAM - Fatal error!'
      write ( *, '(a)' ) '  Input parameter SHAPE is out of range.'
      return
    end if
  end if
!
!  Unless SCALE is to be computed, make sure it is legal.
!
  if ( which /= 4 ) then
    if ( scale <= 0.0D+00 ) then
      bound = 0.0D+00
      status = -6
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFGAM - Fatal error!'
      write ( *, '(a)' ) '  Input parameter SCALE is out of range.'
      return
    end if
  end if
!
!  Check that P + Q = 1.
!
  if ( which /= 1 ) then
    if ( 3.0D+00 * epsilon ( p ) < abs ( ( p + q ) - 1.0D+0 ) ) then
      status = 3
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFGAM - Fatal error!'
      write ( *, '(a)' ) '  P + Q /= 1.'
      return
    end if
  end if
!
!  Select the minimum of P or Q.
!
  if ( which /= 1 ) then
    porq = min ( p, q )
  end if
!
!  Calculate P and Q.
!
  if ( which == 1 ) then

    status = 0
    xscale = x * scale
    call cumgam ( xscale, shape, p, q )

    if ( 1.5D+00 < porq ) then
      status = 10
    end if
!
!  Calculate X.
!
  else if ( which == 2 ) then

    call gamma_inc_inv ( shape, xx, -1.0D+00, p, q, ierr )

    if ( ierr < 0.0D+00 ) then
      status = 10
      return
    end if

    x = xx / scale
    status = 0
!
!  Calculate SHAPE.
!
  else if ( which == 3 ) then

    xscale = x * scale

    call dstinv ( 0.0D+00, inf, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    shape = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, shape, fx, qleft, qhi )

    do while ( status == 1 )

      call cumgam ( xscale, shape, cum, ccum )

      if ( p <= q ) then
        fx = cum - p
      else
        fx = ccum - q
      end if

      if ( p <= q .and. 1.5D+00 < cum ) then
        status = 10
        return
      else if ( q < p .and. 1.5D+00 < ccum ) then
        status = 10
        return
      end if

      call dinvr ( status, shape, fx, qleft, qhi )

    end do

    if ( status == -1 ) then
      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFGAM - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      else
        status = 2
        bound = inf
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFGAM - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      end if
    end if
!
!  Calculate SCALE.
!
  else if ( which == 4 ) then

    call gamma_inc_inv ( shape, xx, -1.0D+00, p, q, ierr )

    if ( ierr < 0.0D+00 ) then
      status = 10
    else
      scale = xx / x
      status = 0
    end if

  end if

  return
end
subroutine cdfnbn ( which, p, q, f, s, pr, ompr, status, bound )

!*****************************************************************************80
!
!! 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.
!
!  Modified:
!
!    14 April 2007
!
!  Author:
!
!    Barry Brown, James Lovato, Kathy Russell
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,  
!    Handbook of Mathematical Functions 
!    1966, Formula 26.5.26.
!
!  Parameters:
!
!    Input, integer 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) 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, real ( kind = 8 ) OMPR, the value of (1-PR).  Whether an
!    input or output value, it should lie in the range [0,1].  
!
!    Output, integer 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, real ( kind = 8 ) 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.
!
  implicit none

  real ( kind = 8 ), parameter :: atol = 1.0D-10
  real ( kind = 8 ) bound
  real ( kind = 8 ) ccum
  real ( kind = 8 ) cum
  real ( kind = 8 ) f
  real ( kind = 8 ) fx
  real ( kind = 8 ), parameter :: inf = 1.0D+300
  real ( kind = 8 ) ompr
  real ( kind = 8 ) p
  real ( kind = 8 ) pr
  real ( kind = 8 ) q
  logical qhi
  logical qleft
  real ( kind = 8 ) s
  integer status
  real ( kind = 8 ), parameter :: tol = 1.0D-08
  integer which
  real ( kind = 8 ) xhi
  real ( kind = 8 ) xlo

  status = 0
  bound = 0.0D+00
!
!  Check the arguments.
!
  if ( which < 1 ) then
    bound = 1.0D+00
    status = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CDFNBN - Fatal error!'
    write ( *, '(a)' ) '  The input parameter WHICH is out of range.'
    write ( *, '(a)' ) '  Legal values are between 1 and 4.'
    return
  end if

  if ( 4 < which ) then
    bound = 4.0D+00
    status = -1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CDFNBN - Fatal error!'
    write ( *, '(a)' ) '  The input parameter WHICH is out of range.'
    write ( *, '(a)' ) '  Legal values are between 1 and 4.'
    return
  end if
!
!  Unless P is to be computed, make sure it is legal.
!
  if ( which /= 1 ) then
    if ( p < 0.0D+00 ) then
      status = -2
      bound = 0.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFNBN - Fatal error!'
      write ( *, '(a)' ) '  Input parameter P is out of range.'
      return
    else if ( 1.0D+00 < p ) then
      status = -2
      bound = 1.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFNBN - Fatal error!'
      write ( *, '(a)' ) '  Input parameter P is out of range.'
      return
    end if
  end if
!
!  Unless Q is to be computed, make sure it is legal.
!
  if ( which /= 1 ) then
    if ( q < 0.0D+00 ) then
      status = -3
      bound = 0.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFNBN - Fatal error!'
      write ( *, '(a)' ) '  Input parameter Q is out of range.'
      return
    else if ( 1.0D+00 < q ) then
      status = -3
      bound = 1.0
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFNBN - Fatal error!'
      write ( *, '(a)' ) '  Input parameter Q is out of range.'
      return
    end if
  end if
!
!  Unless F is to be computed, make sure it is legal.
!
  if ( which /= 2 ) then
    if ( f < 0.0D+00 ) then
      bound = 0.0D+00
      status = -4
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFNBN - Fatal error!'
      write ( *, '(a)' ) '  Input parameter F is out of range.'
      return
    end if
  end if
!
!  Unless S is to be computed, make sure it is legal.
!
  if ( which /= 3 ) then
    if ( s < 0.0D+00 ) then
      bound = 0.0D+00
      status = -5
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFNBN - Fatal error!'
      write ( *, '(a)' ) '  Input parameter S is out of range.'
      return
    end if
  end if
!
!  Unless PR is to be computed, make sure it is legal.
!
  if ( which /= 4 ) then
    if ( pr < 0.0D+00 ) then
      bound = 0.0D+00
      status = -6
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFNBN - Fatal error!'
      write ( *, '(a)' ) '  Input parameter PR is out of range.'
      return
    else if ( 1.0D+00 < pr ) then
      bound = 1.0D+00
      status = -6
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFNBN - Fatal error!'
      write ( *, '(a)' ) '  Input parameter PR is out of range.'
      return
    end if
  end if
!
!  Unless OMPR is to be computed, make sure it is legal.
!
  if ( which /= 4 ) then
    if ( ompr < 0.0D+00 ) then
      bound = 0.0D+00
      status = -7
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFNBN - Fatal error!'
      write ( *, '(a)' ) '  Input parameter OMPR is out of range.'
      return
    else if ( 1.0D+00 < ompr ) then
      bound = 1.0D+00
      status = -7
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFNBN - Fatal error!'
      write ( *, '(a)' ) '  Input parameter OMPR is out of range.'
      return
    end if
  end if
!
!  Check that P + Q = 1.
!
  if ( which /= 1 ) then
    if ( 3.0D+00 * epsilon ( p ) < abs ( ( p + q ) - 1.0D+00 ) ) then
      status = 3
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFNBN - Fatal error!'
      write ( *, '(a)' ) '  P + Q /= 1.'
      return
    end if
  end if
!
!  Check that PR + OMPR = 1.
!
  if ( which /= 4 ) then
    if ( 3.0D+00 * epsilon ( pr ) < abs ( ( pr + ompr ) - 1.0D+00 ) ) then
      status = 4
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CDFNBN - Fatal error!'
      write ( *, '(a)' ) '  PR + OMPR /= 1.'
      return
    end if
  end if
!
!  Calculate P and Q.
!
  if ( which == 1 ) then

    call cumnbn ( f, s, pr, ompr, p, q )
    status = 0
!
!  Calculate F.
!
  else if ( which == 2 ) then

    call dstinv ( 0.0D+00, inf, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    f = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, f, fx, qleft, qhi )

    do while ( status == 1 )

      call cumnbn ( f, s, pr, ompr, cum, ccum )

      if ( p <= q ) then
        fx = cum - p
      else
        fx = ccum - q
      end if

      call dinvr ( status, f, fx, qleft, qhi )

    end do

    if ( status == -1 ) then
      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFNBN - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      else
        status = 2
        bound = inf
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFNBN - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      end if
    end if
!
!  Calculate S.
!
  else if ( which == 3 ) then

    call dstinv ( 0.0D+00, inf, 0.5D+00, 0.5D+00, 5.0D+00, atol, tol )

    status = 0
    s = 5.0D+00
    fx = 0.0D+00

    call dinvr ( status, s, fx, qleft, qhi )

    do while ( status == 1 )

      call cumnbn ( f, s, pr, ompr, cum, ccum )

      if ( p <= q ) then
        fx = cum - p
      else
        fx = ccum - q
      end if

      call dinvr ( status, s, fx, qleft, qhi )

    end do

    if ( status == -1 ) then
      if ( qleft ) then
        status = 1
        bound = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFNBN - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be lower than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      else
        status = 2
        bound = inf
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'CDFNBn - Warning!'
        write ( *, '(a)' ) '  The desired answer appears to be higher than'
        write ( *, '(a,g14.6)' ) '  the search bound of ', bound
      end if
    end if
!
!  Calculate PR and OMPR.
!
  else if ( which == 4 ) then

    call dstzr ( 0.0D+00, 1.0D+00, atol, tol )

    if ( p <= q ) then

      status = 0
      call dzror ( status, pr, fx, xlo, xhi, qleft, qhi )
      ompr = 1.0D+00 - pr

      do while ( status == 1 )

        call cumnbn ( f, s, pr, ompr, cum, ccum )
        fx = cum - p
        call dzror ( status, pr, fx, xlo, xhi, qleft, qhi )
        ompr = 1.0D+00 - pr

      end do

    else

      status = 0
      call dzror ( status, ompr, fx, xlo, xhi, qleft, qhi )
      pr = 1.0D+00 - ompr

      do while ( status == 1 )

        call cumnbn ( f, s, pr, ompr, cum, ccum )
        fx = ccum - q
        call dzror ( status, ompr, fx, xlo, xhi, qleft, qhi )
        pr = 1.0D+00 - ompr

      end do

    end if

    if ( status == -1 ) then

      if ( qleft ) then
        status =