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:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

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

  if ( b < a ) then
    h = b / a
    c = 1.0E+00 / ( 1.0E+00 + h )
    x = h / ( 1.0E+00 + h )
    d = a + ( b - 0.5E+00 )
  else
    h = a / b
    c = h / ( 1.0E+00 + h )
    x = 1.0E+00 / ( 1.0E+00 + h )
    d = b + ( a - 0.5E+00 )
  end if
!
!  Set sn = (1 - x**n) / ( 1 - x )
!
  x2 = x * x
  s3 = 1.0E+00 + ( x + x2)
  s5 = 1.0E+00 + ( x + x2 * s3 )
  s7 = 1.0E+00 + ( x + x2 * s5 )
  s9 = 1.0E+00 + ( x + x2 * s7 )
  s11 = 1.0E+00 + ( x + x2 * s9 )
!
!  Set w = del(b) - del(a + b)
!
  t = ( 1.0E+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.0E+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:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

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

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

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

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

    alnrel = 2.0E+00 * t * w

  else

    x = real ( 1.0D+00 + dble ( a ) )
    alnrel = log ( x )

  end if

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

!*****************************************************************************80
!
!! apser() yields the incomplete beta ratio i(sub(1-x))(b,a) using a series.
!
!  Discussion:
!
!    The algorithm is used for a <= min(eps,eps*b), b*x <= 1, and x <= 0.5. 
!    It is used when a is very small. 
!    Use only if above inequalities are satisfied.
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  real a
  real aj
  real apser
  real b
  real bx
  real c
  real eps
  real, parameter :: g = 0.577215664901533E+00
  real j
  real psi
  real s
  real t
  real tol
  real x

  bx = b * x
  t = x - bx

  if ( b * eps <= 2.0E-02 ) then
    c = log ( x) + psi ( b ) + g + t
  else
    c = log ( bx ) + g + t
  end if

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

  do

    j = j + 1.0E+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 basym ( a, b, lambda, eps )

!*****************************************************************************80
!
!! basym() uses an asymptotic expansion for Ix(A,B) for large A and B.
!
!  Discussion:
!
!     lambda = (a + b)*y - b  and eps is the tolerance used.
!     it is assumed that lambda is nonnegative and that
!     a and b are greater than or equal to 15.
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  real a
  real a0(21)
  real b
  real b0(21)
  real bcorr
  real basym
  real bsum
  real c(21)
  real d(21)
  real dsum
  real, parameter :: e0 = 1.12837916709551E+00
  real, parameter :: e1 = 0.353553390593274E+00
  real eps
  real erfc1
  real f
  real h
  real h2
  real hn
  integer i
  integer im1
  integer imj
  integer j
  real j0
  real j1
  real lambda
  integer m
  integer mm1
  integer mmj
  integer n
  integer np1
  integer, parameter :: num = 20
  real r
  real r0
  real r1
  real rlog1
  real s
  real sum2
  real t
  real t0
  real t1
  real u
  real w
  real w0
  real z
  real z0
  real z2
  real zn
  real znm1
!
!  num is the maximum value that n can take in the do loop
!  ending at statement 50. it is required that num be even.
!  the arrays a0, b0, c, d have dimension num + 1.
!
!     e0 = 2/sqrt(pi)
!     e1 = 2**(-3/2)

  basym = 0.0E+00

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

  f = a * rlog1 ( - lambda / a ) + b * rlog1 ( lambda / b )
  t = exp ( - f )

  if ( t == 0.0E+00 ) then
    return
  end if

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

  a0(1) = ( 2.0E+00 / 3.0E+00 ) * r1
  c(1) = - 0.5E+00 * a0(1)
  d(1) = - c(1)
  j0 = ( 0.5E+00 / e0 ) * erfc1 ( 1, z0 )
  j1 = e1
  sum2 = j0 + d(1) * w0 * j1

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

  do n = 2, num, 2

    hn = h2 * hn
    a0(n) = 2.0 * r0 * ( 1.0 + h * hn ) / real ( n + 2 )
    np1 = n + 1
    s = s + hn
    a0(np1) = 2.0 * r1 * s / real ( n + 3 )

    do i = n, np1

      r = -0.5 * real ( i + 1 )
      b0(1) = r * a0(1)

      do m = 2, i

        bsum = 0.0
        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.0 )

      dsum = 0.0
      im1 = i - 1

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

      d(i) = -(dsum + c(i))

    end do

    j0 = e1*znm1 + ( n - 1.0 ) * 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
    sum2 = sum2 + (t0 + t1)

    if ((abs(t0) + abs(t1)) <= eps * sum2 ) then
      exit
    end if

  end do

  u = exp ( - bcorr ( a, b ) )
  basym = e0 * t * u * sum2

  return
end
function bcorr ( a0, b0 )

!*****************************************************************************80
!
!! bcorr() evaluates  del(a0) + del(b0) - del(a0 + b0).
!
!  Discussion:
!
!    Here ln(gamma(a)) = (a - 0.5)*ln(a) - a + 0.5*ln(2*pi) + del(a).
!    It is assumed that a0 .ge. 8 and b0 .ge. 8.
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

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

  a = min ( a0, b0 )
  b = max ( a0, b0 )
  h = a / b
  c = h/(1.0 + h)
  x = 1.0/(1.0 + h)
  x2 = x*x
!
!  Set sn = (1 - x**n)/(1 - x)
!
  s3 = 1.0 + (x + x2)
  s5 = 1.0 + (x + x2*s3)
  s7 = 1.0 + (x + x2*s5)
  s9 = 1.0 + (x + x2*s7)
  s11 = 1.0 + (x + x2*s9)
!
!  Set w = del(b) - del(a + b)
!
  t = (1.0/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.0/a)**2
  bcorr = (((((c5*t + c4)*t + c3)*t + c2)*t + c1)*t + c0)/a + w

  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
!
!    The incomplete Beta function is also sometimes called the
!    "modified" Beta function, or the "normalized" Beta function
!    or the Beta CDF (cumulative density function.
!
!    In Mathematica, the function can be evaluated by:
!
!      BETA[X,A,B] / BETA[A,B]
!
!    The function can also be evaluated by using the Statistics package:
!
!      Needs["Statistics`ContinuousDistributions`"]
!      dist = BetaDistribution [ a, b ]
!      CDF [ dist, x ]
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    04 January 2006
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz and Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964.
!
!    Karl Pearson,
!    Tables of the Incomplete Beta Function,
!    Cambridge University Press, 1968.
!
!    Stephen Wolfram,
!    The Mathematica Book,
!    Fourth Edition,
!    Wolfram Media / Cambridge University Press, 1999.
!
!  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 A, B, the parameters of the function.
!
!    Output, real X, the argument of the function.
!
!    Output, real FX, the value of the function.
!
  implicit none

  integer, parameter :: n_max = 42

  real a
  real, save, dimension ( n_max ) :: a_vec = (/ &
     0.5E+00, &
     0.5E+00, &
     0.5E+00, &
     1.0E+00, &
     1.0E+00, &
     1.0E+00, &
     1.0E+00, &
     1.0E+00, &
     2.0E+00, &
     2.0E+00, &
     2.0E+00, &
     2.0E+00, &
     2.0E+00, &
     2.0E+00, &
     2.0E+00, &
     2.0E+00, &
     2.0E+00, &
     5.5E+00, &
    10.0E+00, &
    10.0E+00, &
    10.0E+00, &
    10.0E+00, &
    20.0E+00, &
    20.0E+00, &
    20.0E+00, &
    20.0E+00, &
    20.0E+00, &
    30.0E+00, &
    30.0E+00, &
    40.0E+00, &
     0.1E+01, &
     0.1E+01, &
     0.1E+01, &
     0.1E+01, &
     0.1E+01, &
     0.1E+01, &
     0.1E+01, &
     0.1E+01, &
     0.2E+01, &
     0.3E+01, &
     0.4E+01, &
     0.5E+01 /)
  real b
  real, save, dimension ( n_max ) :: b_vec = (/ &
     0.5E+00, &
     0.5E+00, &
     0.5E+00, &
     0.5E+00, &
     0.5E+00, &
     0.5E+00, &
     0.5E+00, &
     1.0E+00, &
     2.0E+00, &
     2.0E+00, &
     2.0E+00, &
     2.0E+00, &
     2.0E+00, &
     2.0E+00, &
     2.0E+00, &
     2.0E+00, &
     2.0E+00, &
     5.0E+00, &
     0.5E+00, &
     5.0E+00, &
     5.0E+00, &
    10.0E+00, &
     5.0E+00, &
    10.0E+00, &
    10.0E+00, &
    20.0E+00, &
    20.0E+00, &
    10.0E+00, &
    10.0E+00, &
    20.0E+00, &
     0.5E+00, &
     0.5E+00, &
     0.5E+00, &
     0.5E+00, &
     0.2E+01, &
     0.3E+01, &
     0.4E+01, &
     0.5E+01, &
     0.2E+01, &
     0.2E+01, &
     0.2E+01, &
     0.2E+01 /)
  real fx
  real, save, dimension ( n_max ) :: fx_vec = (/ &
    0.6376856085851985E-01, &
    0.2048327646991335E+00, &
    0.1000000000000000E+01, &
    0.0000000000000000E+00, &
    0.5012562893380045E-02, &
    0.5131670194948620E-01, &
    0.2928932188134525E+00, &
    0.5000000000000000E+00, &
    0.2800000000000000E-01, &
    0.1040000000000000E+00, &
    0.2160000000000000E+00, &
    0.3520000000000000E+00, &
    0.5000000000000000E+00, &
    0.6480000000000000E+00, &
    0.7840000000000000E+00, &
    0.8960000000000000E+00, &
    0.9720000000000000E+00, &
    0.4361908850559777E+00, &
    0.1516409096347099E+00, &
    0.8978271484375000E-01, &
    0.1000000000000000E+01, &
    0.5000000000000000E+00, &
    0.4598773297575791E+00, &
    0.2146816102371739E+00, &
    0.9507364826957875E+00, &
    0.5000000000000000E+00, &
    0.8979413687105918E+00, &
    0.2241297491808366E+00, &
    0.7586405487192086E+00, &
    0.7001783247477069E+00, &
    0.5131670194948620E-01, &
    0.1055728090000841E+00, &
    0.1633399734659245E+00, &
    0.2254033307585166E+00, &
    0.3600000000000000E+00, &
    0.4880000000000000E+00, &
    0.5904000000000000E+00, &
    0.6723200000000000E+00, &
    0.2160000000000000E+00, &
    0.8370000000000000E-01, &
    0.3078000000000000E-01, &
    0.1093500000000000E-01 /)
  integer n_data
  real x
  real, save, dimension ( n_max ) :: x_vec = (/ &
    0.01E+00, &
    0.10E+00, &
    1.00E+00, &
    0.00E+00, &
    0.01E+00, &
    0.10E+00, &
    0.50E+00, &
    0.50E+00, &
    0.10E+00, &
    0.20E+00, &
    0.30E+00, &
    0.40E+00, &
    0.50E+00, &
    0.60E+00, &
    0.70E+00, &
    0.80E+00, &
    0.90E+00, &
    0.50E+00, &
    0.90E+00, &
    0.50E+00, &
    1.00E+00, &
    0.50E+00, &
    0.80E+00, &
    0.60E+00, &
    0.80E+00, &
    0.50E+00, &
    0.60E+00, &
    0.70E+00, &
    0.80E+00, &
    0.70E+00, &
    0.10E+00, &
    0.20E+00, &
    0.30E+00, &
    0.40E+00, &
    0.20E+00, &
    0.20E+00, &
    0.20E+00, &
    0.20E+00, &
    0.30E+00, &
    0.30E+00, &
    0.30E+00, &
    0.30E+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.0E+00
    b = 0.0E+00
    x = 0.0E+00
    fx = 0.0E+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 beta_log_values ( n, x, y, fxy )

!*****************************************************************************80
!
!! beta_log_values() returns some values of the Beta function for testing.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    15 February 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz and Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964.
!
!  Parameters:
!
!    Input/output, integer N.
!    On input, if N is 0, the first test data is returned, and N is set
!    to the index of the test data.  On each subsequent call, N is
!    incremented and that test data is returned.  When there is no more
!    test data, N is set to 0.
!
!    Output, real X, Y, the arguments of the function.
!
!    Output, real FXY, the value of the function.
!
  implicit none

  integer, parameter :: nmax = 17

  real, save, dimension ( nmax ) :: fxvec = (/ &
    1.609437912E+00,  0.9162907319E+00,  0.5108256238E+00,  0.2231435513E+00, &
    1.609437912E+00,  0.9162907319E+00,  0.0000000000E+00, -1.791759469E+00, &
   -3.401197382E+00, -4.941642423E+00,  -6.445719819E+00,  -3.737669618E+00, &
   -5.123963979E+00, -6.222576268E+00,  -7.138867000E+00,  -7.927324360E+00, &
   -9.393661429E+00 /)
  real fxy
  integer n
  real x
  real, save, dimension ( nmax ) :: xvec = (/ &
    0.2E+00, 0.4E+00, 0.6E+00, 0.8E+00, &
    1.0E+00, 1.0E+00, 1.0E+00, 2.0E+00, &
    3.0E+00, 4.0E+00, 5.0E+00, 6.0E+00, &
    6.0E+00, 6.0E+00, 6.0E+00, 6.0E+00, &
    7.0E+00 /)
  real y
  real, save, dimension ( nmax ) :: yvec = (/ &
    1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, &
    0.2E+00, 0.4E+00, 1.0E+00, 2.0E+00, &
    3.0E+00, 4.0E+00, 5.0E+00, 2.0E+00, &
    3.0E+00, 4.0E+00, 5.0E+00, 6.0E+00, &
    7.0E+00 /)

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

  n = n + 1

  if ( nmax < n ) then
    n = 0
    x = 0.0E+00
    y = 0.0E+00
    fxy = 0.0E+00
    return
  end if

  x = xvec(n)
  y = yvec(n)
  fxy = fxvec(n)

  return
end
function betaln ( a0, b0 )

!*****************************************************************************80
!
!! betaln() evaluates the logarithm of the Beta function.
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
!  Local:
!
!    real E, the value of Log ( 2 * PI ) / 2.
!
  implicit none

  real a
  real a0
  real algdiv
  real alnrel
  real b
  real b0
  real bcorr
  real betaln
  real c
  real, parameter :: e = 0.918938533204673E+00
  real gamln
  real gsumln
  real h
  integer i
  integer n
  real u
  real v
  real w
  real z

  a = min ( a0, b0 )
  b = max ( a0, b0 )
  if (a .ge. 8.0) go to 60
  if (a .ge. 1.0) go to 20
!
!  a < 1
!
  if ( b < 8.0 ) then
    betaln = gamln ( a ) + ( gamln ( b ) - gamln ( a + b ) )
    return
  else
    betaln = gamln ( a ) + algdiv ( a, b )
    return
  end if
!
!  procedure when 1 <= a < 8
!
   20 if (a .gt. 2.0) go to 30
      if (b .gt. 2.0) go to 21
         betaln = gamln(a) + gamln(b) - gsumln(a,b)
         return
   21 w = 0.0
      if (b < 8.0) go to 40
         betaln = gamln(a) + algdiv(a,b)
         return
!
! reduction of a when b <= 1000
!
   30 if (b .gt. 1000.0) go to 50
      n = int ( a - 1.0 )
      w = 1.0
      do i = 1,n
         a = a - 1.0
         h = a/b
         w = w * (h/(1.0 + h))
      end do
      w = log ( w )
      if (b < 8.0) go to 40
      betaln = w + gamln(a) + algdiv(a,b)
      return
!
!  Reduction of b when b < 8
!
   40 n = int ( b - 1.0 )
      z = 1.0
      do i = 1,n
         b = b - 1.0
         z = z * (b/(a + b))
      end do
      betaln = w + log ( z ) + (gamln(a) + (gamln(b) - gsumln(a,b)))
      return
!
!  reduction of a when b .gt. 1000
!
   50 n = int ( a - 1.0 )
      w = 1.0
      do i = 1,n
         a = a - 1.0
         w = w * (a/(1.0 + a/b))
      end do
      betaln = ( log ( w ) - n* log ( b ) ) + (gamln(a) + algdiv(a,b))
      return
!
!  procedure when a .ge. 8
!
   60 w = bcorr(a,b)
      h = a/b
      c = h/(1.0 + h)
      u = -(a - 0.5) * log ( c )
      v = b*alnrel(h)
      if (u <= v) go to 61
         betaln = (((-0.5* log ( b ) + e) + w) - v) - u
         return
   61 betaln = (((-0.5* log ( b ) + e) + w) - u) - v

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

!*****************************************************************************80
!
!! bfrac() uses a continued fraction expansion for ix(a,b) when a,b .gt. 1.
!
!  Discussion:
!
!    It is assumed that  lambda = (a + b)*y - b.
!
!  Modified:
!
!    28 August 2004
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  real a
  real alpha
  real an
  real anp1
  real b
  real beta
  real bfrac
  real bn
  real bnp1
  real brcomp
  real c
  real c0
  real c1
  real e
  real eps
  real lambda
  real n
  real p
  real r
  real r0
  real s
  real t
  real w
  real x
  real y
  real yp1

  bfrac = brcomp ( a, b, x, y )

  if ( bfrac == 0.0 ) then
    return
  end if

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

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

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

  end do
!
!  Termination.
!
  bfrac = bfrac * r

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

!*****************************************************************************80
!
!! bgrat() uses an asymptotic expansion for Ix(a,b) when a is larger than b.
!
!  Discussion:
!
!    The result of the expansion is added to w. it is assumed
!    that a .ge. 15 and b <= 1.  eps is the tolerance used.
!    ierr is a variable that reports the status of the results.
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

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

  bm1 = (b - 0.5) - 0.5
  nu = a + 0.5*bm1
  if (y .gt. 0.375) go to 10
  lnx = alnrel(-y)
  go to 11
10 continue
  lnx = log ( x )
11 continue
  z = -nu*lnx
  if (b*z == 0.0) go to 100
!
!  computation of the expansion
!  set r = exp(-z)*z**b/gamma(b)
!
      r = b*(1.0 + gam1(b))*exp ( b * log ( z ) )
      r = r*exp(a*lnx)*exp(0.5*bm1*lnx)
      u = algdiv(b,a) + b* log ( nu )
      u = r*exp(-u)
      if (u == 0.0) go to 100
      call grat1(b,z,r,p,q,eps)

      v = 0.25*(1.0/nu)**2
      t2 = 0.25*lnx*lnx
      l = w/u
      j = q/r
      sum1 = j
      t = 1.0
      cn = 1.0
      n2 = 0.0
      do n = 1,30
         bp2n = b + n2
         j = (bp2n*(bp2n + 1.0)*j + (z + bp2n + 1.0)*t)*v
         n2 = n2 + 2.0
         t = t*t2
         cn = cn/(n2*(n2 + 1.0))
         c(n) = cn
         s = 0.0
         if (n == 1) go to 21
            nm1 = n - 1
            coef = b - n
            do i = 1,nm1
               s = s + coef*c(i)*d(n-i)
               coef = coef + b
            end do
   21    d(n) = bm1*cn + s/n
         dj = d(n)*j
         sum1 = sum1 + dj

         if ( sum1 <= 0.0 ) then
           go to 100
         end if

         if (abs(dj) <= eps*(sum1 + l)) then
           go to 30
         end if

      end do
!
!  Add the results to w
!
   30 ierr = 0
      w = w + u * sum1
      return
!
!  the expansion cannot be computed
!
  100 ierr = 1

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

!*****************************************************************************80
!
!! bpser() uses the power series expansion for evaluating ix(a,b).
!
!  Discussion:
!
!    It is assumed that b <= 1 or b*x <= 0.7.  eps is the tolerance used.
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  real a
  real a0
  real algdiv
  real apb
  real b
  real b0
  real betaln
  real bpser
  real c
  real eps
  real gam1
  real gamln1
  integer i
  integer m
  real n
  real sum1
  real t
  real tol
  real u
  real w
  real x
  real z

  bpser = 0.0
  if ( x == 0.0 ) then
    return
  end if
!
!  compute the factor x**a/(a*beta(a,b))
!
      a0 = min ( a, b )
      if (a0 < 1.0) go to 10
         z = a* log ( x ) - betaln(a,b)
         bpser = exp(z)/a
         go to 70
   10 b0 = max ( a, b )
      if (b0 .ge. 8.0) go to 60
      if (b0 .gt. 1.0) go to 40
!
!  procedure for a0 < 1 and b0 <= 1
!
      bpser = x**a
      if (bpser == 0.0) return

      apb = a + b
      if (apb .gt. 1.0) go to 20
         z = 1.0 + gam1(apb)
         go to 30
   20 u = real ( dble(a) + dble(b) - 1.d0 )
      z = (1.0 + gam1(u))/apb

   30 c = (1.0 + gam1(a))*(1.0 + gam1(b))/z
      bpser = bpser*c*(b/apb)
      go to 70
!
!  procedure for a0 < 1 and 1 < b0 < 8
!
   40 u = gamln1(a0)
      m = int ( b0 - 1.0 )
      if (m < 1) go to 50
      c = 1.0
      do i = 1,m
         b0 = b0 - 1.0
         c = c*(b0/(a0 + b0))
      end do
      u = log ( c ) + u

   50 z = a* log ( x ) - u
      b0 = b0 - 1.0
      apb = a0 + b0
      if (apb .gt. 1.0) go to 51
         t = 1.0 + gam1(apb)
         go to 52
   51 u = real ( dble(a0) + dble(b0) - 1.d0 )
      t = (1.0 + gam1(u))/apb
   52 bpser = exp(z)*(a0/a)*(1.0 + gam1(b0))/t
      go to 70
!
!  procedure for a0 < 1 and b0 .ge. 8
!
   60 u = gamln1(a0) + algdiv(a0,b0)
      z = a* log ( x ) - u
      bpser = (a0/a)*exp(z)
   70 if (bpser == 0.0 .or. a <= 0.1*eps) return
!
!  Compute the series
!
  sum1 = 0.0
  n = 0.0
  c = 1.0
  tol = eps/a

  do
    n = n + 1.0
    c = c*(0.5 + (0.5 - b/n))*x
    w = c/(a + n)
    sum1 = sum1 + w
    if ( abs ( w ) <= tol ) then
      exit
    end if
  end do

  bpser = bpser*(1.0 + a*sum1)

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

!*****************************************************************************80
!
!! bratio() evaluates the incomplete beta function Ix(A,B).
!
!  Discussion:
!
!    It is assumed that X <= 1
!    and Y = 1 - X.  BRATIO assigns W and W1 the values
!
!                      W  = ix(a,b)
!                      W1 = 1 - ix(a,b)
!
!     ierr is a variable that reports the status of the results.
!     if no input errors are detected then ierr is set to 0 and
!     w and w1 are computed. otherwise, if an error is detected,
!     then w and w1 are assigned the value 0 and ierr is set to
!     one of the following values ...
!
!        ierr = 1  if a or b is negative
!        ierr = 2  if a = b = 0
!        ierr = 3  if x < 0 or x .gt. 1
!        ierr = 4  if y < 0 or y .gt. 1
!        ierr = 5  if x + y /= 1
!        ierr = 6  if x = a = 0
!        ierr = 7  if y = b = 0
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
!  Parameters:
!
!    Input, real A, B, the parameters of the function.  A and B
!    should be nonnegative.
!
  implicit none

  real a
  real a0
  real apser
  real b
  real b0
  real basym
  real bfrac
  real bpser
  real bup
  real eps
  real fpser
  integer ierr
  integer ierr1
  integer ind
  real lambda
  integer n
  real t
  real w
  real w1
  real x
  real x0
  real y
  real y0
  real z

  eps = epsilon ( eps )

  w = 0.0
  w1 = 0.0
  if (a < 0.0 .or. b < 0.0) go to 300
  if (a == 0.0 .and. b == 0.0) go to 310
  if (x < 0.0 .or. x .gt. 1.0) go to 320
      if (y < 0.0 .or. y .gt. 1.0) go to 330

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

  if ( abs(z) .gt. 3.0*eps) then
    ierr = 5
    return
  end if

  ierr = 0

  if (x == 0.0) go to 200
  if (y == 0.0) go to 210
  if (a == 0.0) go to 211
  if (b == 0.0) go to 201

      eps = max ( eps, 1.e-15 )
      if ( max ( a, b ) < 1.e-3*eps) go to 230

      ind = 0
      a0 = a
      b0 = b
      x0 = x
      y0 = y
      if ( min ( a0, b0 ) .gt. 1.0) go to 30
!
!  procedure for a0 <= 1 or b0 <= 1
!
      if (x <= 0.5) go to 10
      ind = 1
      a0 = b
      b0 = a
      x0 = y
      y0 = x

   10 if (b0 < min ( eps, eps * a0 ) ) go to 80
      if (a0 < min ( eps, eps * b0 ) .and. b0*x0 <= 1.0) go to 90
      if ( max ( a0, b0 ) .gt. 1.0) go to 20
      if (a0 .ge. min ( 0.2, b0 ) ) go to 100
      if (x0**a0 <= 0.9) go to 100
      if (x0 .ge. 0.3) go to 110
      n = 20
      go to 130

   20 if (b0 <= 1.0) go to 100
      if (x0 .ge. 0.3) go to 110
      if (x0 .ge. 0.1) go to 21
      if ((x0*b0)**a0 <= 0.7) go to 100
   21 if (b0 .gt. 15.0) go to 131
      n = 20
      go to 130
!
!  procedure for a0 .gt. 1 and b0 .gt. 1
!
   30 if (a .gt. b) go to 31
         lambda = a - (a + b)*x
         go to 32
   31 lambda = (a + b)*y - b

   32 if (lambda .ge. 0.0) go to 40
      ind = 1
      a0 = b
      b0 = a
      x0 = y
      y0 = x
      lambda = abs(lambda)

   40 if (b0 < 40.0 .and. b0*x0 <= 0.7) go to 100
      if (b0 < 40.0) go to 140
      if (a0 .gt. b0) go to 50
         if (a0 <= 100.0) go to 120
         if (lambda .gt. 0.03*a0) go to 120
         go to 180
   50 if (b0 <= 100.0) go to 120
      if (lambda .gt. 0.03*b0) go to 120
      go to 180
!
!  evaluation of the appropriate algorithm
!
   80 continue
      w = fpser(a0, b0, x0, eps)
      w1 = 0.5 + (0.5 - w)
      go to 220

   90 continue
      w1 = apser(a0, b0, x0, eps)
      w = 0.5 + (0.5 - w1)
      go to 220

  100 continue
      w = bpser(a0, b0, x0, eps)
      w1 = 0.5 + (0.5 - w)
      go to 220

  110 continue
      w1 = bpser(b0, a0, y0, eps)
      w = 0.5 + (0.5 - w1)
      go to 220

  120 continue
      w = bfrac(a0, b0, x0, y0, lambda, 15.0*eps)
      w1 = 0.5 + (0.5 - w)
      go to 220

  130 continue
      w1 = bup(b0, a0, y0, x0, n, eps)
      b0 = b0 + n
  131 call bgrat(b0, a0, y0, x0, w1, 15.0*eps, ierr1)
      w = 0.5 + (0.5 - w1)
      go to 220

  140 n = int ( b0 )
      b0 = b0 - n

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

      w = bup(b0, a0, y0, x0, n, eps)
      if (x0 .gt. 0.7) go to 150
      w = w + bpser(a0, b0, x0, eps)
      w1 = 0.5 + (0.5 - w)
      go to 220

  150 if ( 15.0 < a0 ) go to 151
         n = 20
         w = w + bup(a0, b0, x0, y0, n, eps)
         a0 = a0 + n
  151 call bgrat(a0, b0, x0, y0, w, 15.0*eps, ierr1)
      w1 = 0.5 + (0.5 - w)
      go to 220

  180 w = basym(a0, b0, lambda, 100.0*eps)
      w1 = 0.5 + (0.5 - w)
      go to 220
!
!  termination of the procedure
!
  200 if (a == 0.0) go to 350
  201 w = 0.0
      w1 = 1.0
      return

  210 if (b == 0.0) go to 360
  211 w = 1.0
      w1 = 0.0
      return

  220 if (ind == 0) return
      t = w
      w = w1
      w1 = t
      return
!
!  procedure for a and b < 1.e-3*eps
!
  230 w = b/(a + b)
      w1 = a/(a + b)
      return
!
!  Error return
!
  300 ierr = 1
      return
  310 ierr = 2
      return
  320 ierr = 3
      return
  330 ierr = 4
      return

  350 ierr = 6
      return
  360 ierr = 7

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

!*****************************************************************************80
!
!! brcmp1() evaluates exp(mu) * (x**a*y**b/beta(a,b)).
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  real a
  real a0
  real algdiv
  real alnrel
  real apb
  real b
  real b0
  real bcorr
  real betaln
  real brcmp1
  real c
  real, parameter :: const = 0.398942280401433
  real e
  real esum
  real gam1
  real gamln1
  real h
  integer i
  real lambda
  real lnx
  real lny
  integer mu
  integer n
  real rlog1
  real t
  real u
  real v
  real x
  real x0
  real y
  real y0
  real z

  a0 = min ( a, b )
      if (a0 .ge. 8.0) go to 100

      if (x .gt. 0.375) go to 10
         lnx = log ( x )
         lny = alnrel(-x)
         go to 20
   10 if (y .gt. 0.375) go to 11
         lnx = alnrel(-y)
         lny = log ( y )
         go to 20
   11 lnx = log ( x )
      lny = log ( y )

   20 z = a*lnx + b*lny
      if (a0 < 1.0) go to 30
      z = z - betaln(a,b)
      brcmp1 = esum(mu,z)
      return
!-----------------------------------------------------------------------
!  procedure for a < 1 or b < 1
!-----------------------------------------------------------------------
   30 b0 = max ( a, b )
      if (b0 .ge. 8.0) go to 80
      if (b0 .gt. 1.0) go to 60
!
!  algorithm for b0 <= 1
!
      brcmp1 = esum(mu,z)
      if (brcmp1 == 0.0) return

      apb = a + b
      if ( 1.0 < apb ) go to 40
         z = 1.0 + gam1(apb)
         go to 50
   40 u = real ( dble(a) + dble(b) - 1.d0 )
      z = (1.0 + gam1(u))/apb

   50 c = (1.0 + gam1(a))*(1.0 + gam1(b))/z
      brcmp1 = brcmp1*(a0*c)/(1.0 + a0/b0)
      return
!
!  algorithm for 1 < b0 < 8
!
   60 u = gamln1(a0)
      n = int ( b0 - 1.0 )
      if (n < 1) go to 70
      c = 1.0
      do 61 i = 1,n
         b0 = b0 - 1.0
         c = c*(b0/(a0 + b0))
   61 continue
      u = log ( c ) + u

   70 z = z - u
      b0 = b0 - 1.0
      apb = a0 + b0
      if (apb .gt. 1.0) go to 71
         t = 1.0 + gam1(apb)
         go to 72
   71 u = real ( dble(a0) + dble(b0) - 1.d0 )
      t = (1.0 + gam1(u))/apb
   72 brcmp1 = a0*esum(mu,z)*(1.0 + gam1(b0))/t
      return
!
!  algorithm for b0 .ge. 8
!
   80 u = gamln1(a0) + algdiv(a0,b0)
      brcmp1 = a0*esum(mu,z - u)
      return
!-----------------------------------------------------------------------
!  procedure for a .ge. 8 and b .ge. 8
!-----------------------------------------------------------------------
  100 if (a .gt. b) go to 101
         h = a/b
         x0 = h/(1.0 + h)
         y0 = 1.0/(1.0 + h)
         lambda = a - (a + b)*x
         go to 110
  101 h = b/a
      x0 = 1.0/(1.0 + h)
      y0 = h/(1.0 + h)
      lambda = (a + b)*y - b

  110 e = -lambda/a
      if (abs(e) .gt. 0.6) go to 111
         u = rlog1(e)
         go to 120
  111 u = e - log ( x / x0 )

  120 e = lambda/b
      if (abs(e) .gt. 0.6) go to 121
         v = rlog1(e)
         go to 130
  121 v = e - log ( y / y0 )

  130 z = esum(mu,-(a*u + b*v))
      brcmp1 = const*sqrt(b*x0)*z*exp(-bcorr(a,b))

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

!*****************************************************************************80
!
!! brcomp() evaluates X**a * y**b / beta(a,b).
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  real a
  real a0
  real algdiv
  real alnrel
  real apb
  real b
  real b0
  real bcorr
  real betaln
  real brcomp
  real c
  real, parameter :: const = 0.398942280401433
  real e
  real gam1
  real gamln1
  real h
  integer i
  real lambda
  real lnx
  real lny
  integer n
  real rlog1
  real t
  real u
  real v
  real x
  real x0
  real y
  real y0
  real z

  brcomp = 0.0

  if ( x == 0.0 .or. y == 0.0 ) then
    return
  end if

  a0 = min ( a, b )

      if ( 8.0 <= a0 ) go to 100

      if (x .gt. 0.375) go to 10
         lnx = log ( x )
         lny = alnrel(-x)
         go to 20
   10 if (y .gt. 0.375) go to 11
         lnx = alnrel(-y)
         lny = log ( y )
         go to 20
   11 lnx = log ( x )
      lny = log ( y )

   20 z = a*lnx + b*lny

      if ( 1.0 <= a ) then
        z = z - betaln(a,b)
        brcomp = exp(z)
        return
      end if
!-----------------------------------------------------------------------
!  procedure for a < 1 or b < 1
!-----------------------------------------------------------------------
      b0 = max ( a, b )
      if (b0 .ge. 8.0) go to 80
      if ( 1.0 < b0 ) go to 60
!
!  algorithm for b0 <= 1
!
      brcomp = exp(z)
      if (brcomp == 0.0) return

      apb = a + b
      if (apb .gt. 1.0) go to 40
         z = 1.0 + gam1(apb)
         go to 50
   40 u = real ( dble(a) + dble(b) - 1.d0 )
      z = (1.0 + gam1(u))/apb

   50 c = (1.0 + gam1(a))*(1.0 + gam1(b))/z
      brcomp = brcomp*(a0*c)/(1.0 + a0/b0)
      return
!
!  algorithm for 1 < b0 < 8
!
   60 u = gamln1(a0)
      n = int ( b0 - 1.0 )
      if (n < 1) go to 70
      c = 1.0
      do i = 1,n
         b0 = b0 - 1.0
         c = c*(b0/(a0 + b0))
      end do
      u = log ( c ) + u

   70 z = z - u
      b0 = b0 - 1.0
      apb = a0 + b0
      if (apb .gt. 1.0) go to 71
         t = 1.0 + gam1(apb)
         go to 72
   71 u = real ( dble(a0) + dble(b0) - 1.d0 )
      t = (1.0 + gam1(u))/apb
   72 brcomp = a0*exp(z)*(1.0 + gam1(b0))/t
      return
!
!  algorithm for b0 .ge. 8
!
   80 u = gamln1(a0) + algdiv(a0,b0)
      brcomp = a0*exp(z - u)
      return
!-----------------------------------------------------------------------
!  procedure for a .ge. 8 and b .ge. 8
!-----------------------------------------------------------------------
  100 if (a .gt. b) go to 101
         h = a/b
         x0 = h/(1.0 + h)
         y0 = 1.0/(1.0 + h)
         lambda = a - (a + b)*x
         go to 110
  101 h = b/a
      x0 = 1.0/(1.0 + h)
      y0 = h/(1.0 + h)
      lambda = (a + b)*y - b

  110 e = -lambda/a
      if (abs(e) .gt. 0.6) go to 111
         u = rlog1(e)
         go to 120
  111 u = e - log ( x / x0 )

  120 e = lambda/b
      if (abs(e) .gt. 0.6) go to 121
         v = rlog1(e)
         go to 130
  121 v = e - log ( y / y0 )

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

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

!*****************************************************************************80
!
!! bup() evaluates ix(a,b) - ix(a+n,b) where n is a positive integer.
!
!  Discussion:
!
!    eps is the tolerance used.
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  real a
  real ap1
  real apb
  real b
  real brcmp1
  real bup
  real d
  real eps
  real exparg
  integer i
  integer k
  integer kp1
  real l
  integer mu
  integer n
  integer nm1
  real r
  real t
  real w
  real x
  real y
!
!  obtain the scaling factor exp(-mu) and
!  exp(mu)*(x**a*y**b/beta(a,b))/a
!
  apb = a + b
  ap1 = a + 1.0
  mu = 0
  d = 1.0
  if (n == 1 .or. a < 1.0) go to 10
      if (apb < 1.1*ap1) go to 10
         mu = int ( abs(exparg(1)) )
         k = int ( exparg(0) )
         if (k < mu) mu = k
         t = mu
         d = exp(-t)

10 continue

  bup = brcmp1(mu,a,b,x,y)/a
      if (n == 1 .or. bup == 0.0) return
      nm1 = n - 1
      w = d
!
!  let k be the index of the maximum term
!
      k = 0
      if (b <= 1.0) go to 40
      if ( 1.e-4 < y ) go to 20
         k = nm1
         go to 30
   20 r = (b - 1.0)*x/y - a
      if (r < 1.0) go to 40
      k = nm1
      t = nm1
      if (r < t) k = int ( r )
!
!  add the increasing terms of the series
!
   30 do i = 1,k
         l = i - 1
         d = ((apb + l)/(ap1 + l))*x*d
         w = w + d
      end do
      if (k == nm1) go to 50
!
!  add the remaining terms of the series
!
   40 kp1 = k + 1
      do i = kp1,nm1
         l = i - 1
         d = ((apb + l)/(ap1 + l))*x*d
         w = w + d
         if (d <= eps*w) go to 50
      end do
!
!  terminate the procedure
!
   50 bup = bup*w

  return
end
function erfc1 ( ind, x )

!*****************************************************************************80
!
!! erfc1() evaluates the complementary error function
!
!  Discussion:
!
!    erfc1(ind,x) = erfc(x)            if ind = 0
!    erfc1(ind,x) = exp(x*x)*erfc(x)   otherwise
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  real a(5)
  real ax
  real b(3)
  real bot
  real, parameter :: c = 0.564189583547756
  real e
  real erfc1
  real exparg
  integer ind
  real p(8)
  real q(8)
  real r(5)
  real s(4)
  real t
  real top
  double precision w
  real x

      data a(1) /.771058495001320e-04/, a(2)/-.133733772997339e-02/, &
           a(3) /.323076579225834e-01/, a(4) /.479137145607681e-01/, &
           a(5) /.128379167095513e+00/
      data b(1) /.301048631703895e-02/, b(2) /.538971687740286e-01/, &
           b(3) /.375795757275549e+00/

      data p(1)/-1.36864857382717e-07/, p(2) /5.64195517478974e-01/, &
           p(3) /7.21175825088309e+00/, p(4) /4.31622272220567e+01/, &
           p(5) /1.52989285046940e+02/, p(6) /3.39320816734344e+02/, &
           p(7) /4.51918953711873e+02/, p(8) /3.00459261020162e+02/
      data q(1) /1.00000000000000e+00/, q(2) /1.27827273196294e+01/, &
           q(3) /7.70001529352295e+01/, q(4) /2.77585444743988e+02/, &
           q(5) /6.38980264465631e+02/, q(6) /9.31354094850610e+02/, &
           q(7) /7.90950925327898e+02/, q(8) /3.00459260956983e+02/
      data r(1) /2.10144126479064e+00/, r(2) /2.62370141675169e+01/, &
           r(3) /2.13688200555087e+01/, r(4) /4.65807828718470e+00/, &
           r(5) /2.82094791773523e-01/
      data s(1) /9.41537750555460e+01/, s(2) /1.87114811799590e+02/, &
           s(3) /9.90191814623914e+01/, s(4) /1.80124575948747e+01/

!
!  abs(x) <= 0.5
!
      ax = abs(x)
      if (ax .gt. 0.5) go to 10
      t = x*x
      top = ((((a(1)*t + a(2))*t + a(3))*t + a(4))*t + a(5)) + 1.0
      bot = ((b(1)*t + b(2))*t + b(3))*t + 1.0
      erfc1 = 0.5 + (0.5 - x*(top/bot))
      if (ind /= 0) erfc1 = exp(t) * erfc1
      return
!
!  0.5 < abs(x) <= 4
!
   10 if (ax .gt. 4.0) go to 20
      top = ((((((p(1)*ax + p(2))*ax + p(3))*ax + p(4))*ax + p(5))*ax &
                          + p(6))*ax + p(7))*ax + p(8)
      bot = ((((((q(1)*ax + q(2))*ax + q(3))*ax + q(4))*ax + q(5))*ax &
                          + q(6))*ax + q(7))*ax + q(8)
      erfc1 = top/bot
      go to 40
!
!  abs(x) .gt. 4
!
   20 if (x <= -5.6) go to 50
      if (ind /= 0) go to 30
      if ( 100.0 < x ) go to 60
      if ( -exparg(1) < x * x ) go to 60

   30 t = (1.0/x)**2
      top = (((r(1)*t + r(2))*t + r(3))*t + r(4))*t + r(5)
      bot = (((s(1)*t + s(2))*t + s(3))*t + s(4))*t + 1.0
      erfc1 = (c - t*top/bot)/ax
!
!  final assembly
!
   40 if (ind == 0) go to 41
         if (x < 0.0) erfc1 = 2.0*exp(x*x) - erfc1
         return
   41 w = dble(x)*dble(x)
      t = real ( w )
      e = real ( w - dble(t) )
      erfc1 = ((0.5 + (0.5 - e)) * exp(-t)) * erfc1
      if (x < 0.0) erfc1 = 2.0 - erfc1
      return
!
!  limit value for large negative x
!
   50 erfc1 = 2.0
      if (ind /= 0) erfc1 = 2.0*exp(x*x)
      return
!
!  limit value for large positive x when ind = 0
!
   60 erfc1 = 0.0

  return
end
function r4_erf ( x )

!*****************************************************************************80
!
!! r4_erf() evaluates the real error function.
!
!  Discussion:
!
!    This function was originally named "erf(x)", but this conflicts with
!    the corresponding intrinsic.
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  real a(5)
  real ax
  real b(3)
  real bot
  real, parameter :: c = 0.564189583547756
  real r4_erf
  real p(8)
  real q(8)
  real r(5)
  real s(4)
  real t
  real top
  real x
  real x2

      data a(1) /.771058495001320e-04/, a(2)/-.133733772997339e-02/, &
           a(3) /.323076579225834e-01/, a(4) /.479137145607681e-01/, &
           a(5) /.128379167095513e+00/
      data b(1) /.301048631703895e-02/, b(2) /.538971687740286e-01/, &
           b(3) /.375795757275549e+00/
      data p(1)/-1.36864857382717e-07/, p(2) /5.64195517478974e-01/, &
           p(3) /7.21175825088309e+00/, p(4) /4.31622272220567e+01/, &
           p(5) /1.52989285046940e+02/, p(6) /3.39320816734344e+02/, &
           p(7) /4.51918953711873e+02/, p(8) /3.00459261020162e+02/
      data q(1) /1.00000000000000e+00/, q(2) /1.27827273196294e+01/, &
           q(3) /7.70001529352295e+01/, q(4) /2.77585444743988e+02/, &
           q(5) /6.38980264465631e+02/, q(6) /9.31354094850610e+02/, &
           q(7) /7.90950925327898e+02/, q(8) /3.00459260956983e+02/
      data r(1) /2.10144126479064e+00/, r(2) /2.62370141675169e+01/, &
           r(3) /2.13688200555087e+01/, r(4) /4.65807828718470e+00/, &
           r(5) /2.82094791773523e-01/
      data s(1) /9.41537750555460e+01/, s(2) /1.87114811799590e+02/, &
           s(3) /9.90191814623914e+01/, s(4) /1.80124575948747e+01/

  ax = abs ( x )

  if ( ax <= 0.5 ) then

    t = x*x
    top = ((((a(1)*t + a(2))*t + a(3))*t + a(4))*t + a(5)) + 1.0
    bot = ((b(1)*t + b(2))*t + b(3))*t + 1.0
    r4_erf = x*(top/bot)
  
  else if ( ax <= 4.0 ) then

    top = ((((((p(1)*ax + p(2))*ax + p(3))*ax + p(4))*ax + p(5))*ax &
                        + p(6))*ax + p(7))*ax + p(8)
    bot = ((((((q(1)*ax + q(2))*ax + q(3))*ax + q(4))*ax + q(5))*ax &
                        + q(6))*ax + q(7))*ax + q(8)
    r4_erf = 0.5 + (0.5 - exp(-x*x)*top/bot)
    if ( x < 0.0 ) then
      r4_erf = -r4_erf
    end if

  else if ( ax < 5.8 ) then

    x2 = x*x
    t = 1.0/x2
    top = (((r(1)*t + r(2))*t + r(3))*t + r(4))*t + r(5)
    bot = (((s(1)*t + s(2))*t + s(3))*t + s(4))*t + 1.0
    r4_erf = (c - top/(x2*bot)) / ax
    r4_erf = 0.5 + (0.5 - exp(-x2)*r4_erf)

    if ( x < 0.0 ) then
      r4_erf = -r4_erf
    end if

  else

   r4_erf = sign ( 1.0, x )

  end if

  return
end
subroutine erf_values ( n_data, x, fx )

!*****************************************************************************80
!
!! erf_values() returns some values of the ERF or "error" function for testing.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    17 April 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer N_DATA.
!    On input, if N_DATA is 0, the first test data is returned, and
!    N_DATA is set to the index of the test data.  On each subsequent
!    call, N_DATA is incremented and that test data is returned.  When
!    there is no more test data, N_DATA is set to 0.
!
!    Output, real X, the argument of the function.
!
!    Output, real FX, the value of the function.
!
  implicit none

  integer, parameter :: nmax = 21

  real, save, dimension ( nmax ) :: bvec = (/ &
    0.0000000000E+00, 0.1124629160E+00, 0.2227025892E+00, 0.3286267595E+00, &
    0.4283923550E+00, 0.5204998778E+00, 0.6038560908E+00, 0.6778011938E+00, &
    0.7421009647E+00, 0.7969082124E+00, 0.8427007929E+00, 0.8802050696E+00, &
    0.9103139782E+00, 0.9340079449E+00, 0.9522851198E+00, 0.9661051465E+00, &
    0.9763483833E+00, 0.9837904586E+00, 0.9890905016E+00, 0.9927904292E+00, &
    0.9953222650E+00 /)
  real fx
  integer n_data
  real x
  real, save, dimension ( nmax ) :: xvec = (/ &
    0.0E+00, 0.1E+00, 0.2E+00, 0.3E+00, &
    0.4E+00, 0.5E+00, 0.6E+00, 0.7E+00, &
    0.8E+00, 0.9E+00, 1.0E+00, 1.1E+00, &
    1.2E+00, 1.3E+00, 1.4E+00, 1.5E+00, &
    1.6E+00, 1.7E+00, 1.8E+00, 1.9E+00, &
    2.0E+00 /)

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

  n_data = n_data + 1

  if ( nmax < n_data ) then
    n_data = 0
    x = 0.0E+00
    fx = 0.0E+00
  else
    x = xvec(n_data)
    fx = bvec(n_data)
  end if

  return
end
function esum ( mu, x )

!*****************************************************************************80
!
!! esum() evaluates exp(mu + x).
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  real esum
  integer mu
  real w
  real x

  if (x .gt. 0.0) go to 10

      if (mu < 0) go to 20
         w = mu + x
         if (w .gt. 0.0) go to 20
         esum = exp(w)
         return

   10 if (mu .gt. 0) go to 20
         w = mu + x
         if (w < 0.0) go to 20
         esum = exp(w)
         return

   20 w = mu

  esum = exp(w)*exp(x)

  return
end
function exparg ( l )

!*****************************************************************************80
!
!! exparg() reports the largest safe arguments for EXP(X).
!
!  Discussion:
!
!    if l = 0 then  exparg(l) = the largest positive w for which
!    exp(w) can be computed.
!
!    if l is nonzero then  exparg(l) = the largest negative w for
!    which the computed value of exp(w) is nonzero.
!
!    Only an approximate value for exparg(l) is needed.
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  integer b
  real exparg
  integer ipmpar
  integer l
  integer m
  real lnb

  b = ipmpar(4)
  if (b /= 2) go to 10
         lnb = .69314718055995
         go to 50
   10 if (b /= 8) go to 20
         lnb = 2.0794415416798
         go to 50
   20 if (b /= 16) go to 30
         lnb = 2.7725887222398
         go to 50
   30 lnb = log ( float ( b ) )

   50 if (l == 0) go to 60
         m = ipmpar(6) - 1
         exparg = 0.99999 * (m * lnb)
         return
   60 m = ipmpar(7)
      exparg = 0.99999 * (m * lnb)

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

!*****************************************************************************80
!
!! fpser() evaluates Ix(A,B) for small B and X <= 0.5.
!
!  Discussion:
!
!    for b < min(eps,eps*a) and x <= 0.5.
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  real a
  real an
  real b
  real c
  real eps
  real exparg
  real fpser
  real s
  real t
  real tol
  real x
!
!  Set FPSER = X**A.
!
  fpser = 1.0

  if ( 1.0E-03 * eps < a ) then
    fpser = 0.0
    t = a * log ( x )
    if ( t < exparg ( 1 ) ) then
      return
    end if
    fpser = exp ( t )
  end if
!
!  Note that 1/b(a,b) = b
!
  fpser = ( b / a ) * fpser
  tol = eps/a
  an = a + 1.0
  t = x
  s = t/an

  do

    an = an + 1.0
    t = x * t
    c = t / an
    s = s + c

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

  end do

  fpser = fpser * ( 1.0E+00 + a * s )

  return
end
function gam1 ( a )

!*****************************************************************************80
!
!! gam1() computes 1/gamma(a+1) - 1  for -0.5 <= a <= 1.5
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  real a
  real bot
  real d
  real gam1
  real p(7)
  real q(5)
  real r(9)
  real s1
  real s2
  real t
  real top
  real w

      data p(1)/ .577215664901533e+00/, p(2)/-.409078193005776e+00/, &
           p(3)/-.230975380857675e+00/, p(4)/ .597275330452234e-01/, &
           p(5)/ .766968181649490e-02/, p(6)/-.514889771323592e-02/, &
           p(7)/ .589597428611429e-03/

      data q(1)/ .100000000000000e+01/, q(2)/ .427569613095214e+00/, &
           q(3)/ .158451672430138e+00/, q(4)/ .261132021441447e-01/, &
           q(5)/ .423244297896961e-02/

      data r(1)/-.422784335098468e+00/, r(2)/-.771330383816272e+00/, &
           r(3)/-.244757765222226e+00/, r(4)/ .118378989872749e+00/, &
           r(5)/ .930357293360349e-03/, r(6)/-.118290993445146e-01/, &
           r(7)/ .223047661158249e-02/, r(8)/ .266505979058923e-03/, &
           r(9)/-.132674909766242e-03/

      data s1  / .273076135303957e+00/, s2  / .559398236957378e-01/

      t = a
      d = a - 0.5
      if (d .gt. 0.0) then
        t = d - 0.5
      end if

      if ( t == 0.0 ) then

        gam1 = 0.0

      else if ( 0.0 < t ) then

        top = (((((p(7)*t + p(6))*t + p(5))*t + p(4))*t + p(3))*t &
                        + p(2))*t + p(1)
        bot = (((q(5)*t + q(4))*t + q(3))*t + q(2))*t + 1.0
        w = top/bot
        if (d <= 0.0) then
          gam1 = a*w
        else
          gam1 = (t/a)*((w - 0.5) - 0.5)
        end if

      else

        top = (((((((r(9)*t + r(8))*t + r(7))*t + r(6))*t + r(5))*t &
                          + r(4))*t + r(3))*t + r(2))*t + r(1)
        bot = (s2*t + s1)*t + 1.0
        w = top/bot
        if ( d <= 0.0 ) then
          gam1 = a*((w + 0.5) + 0.5)
        else
          gam1 = t*w/a
        end if

      end if

  return
end
function gamln1 ( a )

!*****************************************************************************80
!
!! gamln1() evaluates ln(gamma(1 + a)) for -0.2 <= A <= 1.25
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  real a
  real gamln1
  real p0
  real p1
  real p2
  real p3
  real p4
  real p5
  real p6
  real q1
  real q2
  real q3
  real q4
  real q5
  real q6
  real r0
  real r1
  real r2
  real r3
  real r4
  real r5
  real s1
  real s2
  real s3
  real s4
  real s5
  real w
  real x

      data p0/ .577215664901533e+00/, p1/ .844203922187225e+00/, &
           p2/-.168860593646662e+00/, p3/-.780427615533591e+00/, &
           p4/-.402055799310489e+00/, p5/-.673562214325671e-01/, &
           p6/-.271935708322958e-02/
      data q1/ .288743195473681e+01/, q2/ .312755088914843e+01/, &
           q3/ .156875193295039e+01/, q4/ .361951990101499e+00/, &
           q5/ .325038868253937e-01/, q6/ .667465618796164e-03/

      data r0/.422784335098467e+00/,  r1/.848044614534529e+00/, &
           r2/.565221050691933e+00/,  r3/.156513060486551e+00/, &
           r4/.170502484022650e-01/,  r5/.497958207639485e-03/
      data s1/.124313399877507e+01/,  s2/.548042109832463e+00/, &
           s3/.101552187439830e+00/,  s4/.713309612391000e-02/, &
           s5/.116165475989616e-03/

  if ( a < 0.6 ) then
      w = ((((((p6*a + p5)*a + p4)*a + p3)*a + p2)*a + p1)*a + p0)/ &
          ((((((q6*a + q5)*a + q4)*a + q3)*a + q2)*a + q1)*a + 1.0)
      gamln1 = -a*w
  else
    x = (a - 0.5) - 0.5
      w = (((((r5*x + r4)*x + r3)*x + r2)*x + r1)*x + r0)/ &
          (((((s5*x + s4)*x + s3)*x + s2)*x + s1)*x + 1.0)
      gamln1 = x*w
  end if

  return
end
function gamln ( a )

!*****************************************************************************80
!
!! gamln() evaluates ln(gamma(a)) for positive A.
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none
!
!     d = 0.5*(ln(2*pi) - 1)
!
  real a
  real c0
  real c1
  real c2
  real c3
  real c4
  real c5
  real, parameter :: d = 0.418938533204673
  real gamln
  real gamln1
  integer i
  integer n
  real t
  real w


      data c0/.833333333333333e-01/, c1/-.277777777760991e-02/, &
           c2/.793650666825390e-03/, c3/-.595202931351870e-03/, &
           c4/.837308034031215e-03/, c5/-.165322962780713e-02/

  if ( a <= 0.8 ) then

    gamln = gamln1 ( a ) - log ( a )

  else if ( a <= 2.25 ) then

    t = (a - 0.5) - 0.5
    gamln = gamln1 ( t )

  else if ( a < 10.0 ) then

    n = int ( a - 1.25 )
    t = a
    w = 1.0
    do i = 1, n
      t = t - 1.0
      w = t*w
    end do
    gamln = gamln1(t - 1.0) + log ( w )

  else

    t = (1.0/a)**2
    w = (((((c5*t + c4)*t + c3)*t + c2)*t + c1)*t + c0)/a
    gamln = (d + w) + (a - 0.5)*( log ( a ) - 1.0)

  end if

  return
end
subroutine gamma_inc_values ( n_data, a, x, fx )

!*****************************************************************************80
!
!! gamma_inc_values() returns some values of the incomplete Gamma function.
!
!  Discussion:
!
!    The (normalized) incomplete Gamma function P(A,X) is defined as:
!
!      PN(A,X) = 1/Gamma(A) * Integral ( 0 <= T <= X ) T**(A-1) * exp(-T) dT.
!
!    With this definition, for all A and X,
!
!      0 <= PN(A,X) <= 1
!
!    and
!
!      PN(A,INFINITY) = 1.0
!
!    In Mathematica, the function can be evaluated by:
!
!      1 - GammaRegularized[A,X]
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    28 August 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz and Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964.
!
!    Stephen Wolfram,
!    The Mathematica Book,
!    Fourth Edition,
!    Wolfram Media / Cambridge University Press, 1999.
!
!  Parameters:
!
!    Input/output, 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 A, the parameter of the function.
!
!    Output, real X, the argument of the function.
!
!    Output, real FX, the value of the function.
!
  implicit none

  integer, parameter :: n_max = 20

  real a
  real, save, dimension ( n_max ) :: a_vec = (/ &
     0.10E+00, &
     0.10E+00, &
     0.10E+00, &
     0.50E+00, &
     0.50E+00, &
     0.50E+00, &
     0.10E+01, &
     0.10E+01, &
     0.10E+01, &
     0.11E+01, &
     0.11E+01, &
     0.11E+01, &
     0.20E+01, &
     0.20E+01, &
     0.20E+01, &
     0.60E+01, &
     0.60E+01, &
     0.11E+02, &
     0.26E+02, &
     0.41E+02 /)
  real fx
  real, save, dimension ( n_max ) :: fx_vec = (/ &
     0.7382350532339351E+00, &
     0.9083579897300343E+00, &
     0.9886559833621947E+00, &
     0.3014646416966613E+00, &
     0.7793286380801532E+00, &
     0.9918490284064973E+00, &
     0.9516258196404043E-01, &
     0.6321205588285577E+00, &
     0.9932620530009145E+00, &
     0.7205974576054322E-01, &
     0.5891809618706485E+00, &
     0.9915368159845525E+00, &
     0.1392920235749422E+00, &
     0.7768698398515702E+00, &
     0.9990881180344455E+00, &
     0.4202103819530612E-01, &
     0.9796589705830716E+00, &
     0.9226039842296429E+00, &
     0.4470785799755852E+00, &
     0.7444549220718699E+00 /)
  integer n_data
  real x
  real, save, dimension ( n_max ) :: x_vec = (/ &
     0.30E-01, &
     0.30E+00, &
     0.15E+01, &
     0.75E-01, &
     0.75E+00, &
     0.35E+01, &
     0.10E+00, &
     0.10E+01, &
     0.50E+01, &
     0.10E+00, &
     0.10E+01, &
     0.50E+01, &
     0.15E+00, &
     0.15E+01, &
     0.70E+01, &
     0.25E+01, &
     0.12E+02, &
     0.16E+02, &
     0.25E+02, &
     0.45E+02 /)

  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.0E+00
    x = 0.0E+00
    fx = 0.0E+00
  else
    a = a_vec(n_data)
    x = x_vec(n_data)
    fx = fx_vec(n_data)
  end if

  return
end
subroutine gamma_log_values ( n_data, x, fx )

!*****************************************************************************80
!
!! gamma_log_values() returns some values of the Log Gamma function for testing.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    17 April 2001
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz and Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964.
!
!  Parameters:
!
!    Input/output, integer N_DATA.
!    On input, if N_DATA is 0, the first test data is returned, and
!    N_DATA is set to the index of the test data.  On each subsequent
!    call, N_DATA is incremented and that test data is returned.  When
!    there is no more test data, N_DATA is set to 0.
!
!    Output, real X, the argument of the function.
!
!    Output, real FX, the value of the function.
!
  implicit none

  integer, parameter :: nmax = 18

  real, save, dimension ( nmax ) :: bvec = (/ &
     1.524064183E+00,    0.7966780066E+00,   0.3982337117E+00,  &
     0.1520599127E+00,   0.000000000E+00,   -0.04987246543E+00, &
    -0.08537410945E+00, -0.1081747934E+00,  -0.1196128950E+00,  & 
    -0.1207822040E+00,  -0.1125917658E+00,  -0.09580771625E+00, &
    -0.07108385116E+00, -0.03898428380E+00,  0.000000000E+00,   &
    12.80182743E+00,    39.33988571E+00,    71.25704193E+00 /)
  real fx
  integer n_data
  real x
  real, save, dimension ( nmax ) :: xvec = (/ &
    0.2E+00,  0.4E+00,  0.6E+00,  0.8E+00, &
    1.0E+00,  1.1E+00,  1.2E+00,  1.3E+00, &
    1.4E+00,  1.5E+00,  1.6E+00,  1.7E+00, &
    1.8E+00,  1.9E+00,  2.0E+00, 10.0E+00, &
   20.0E+00, 30.0E+00 /)

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

  n_data = n_data + 1

  if ( nmax < n_data ) then
    n_data = 0
    x = 0.0E+00
    fx = 0.0E+00
  else
    x = xvec(n_data)
    fx = bvec(n_data)
  end if

  return
end
subroutine grat1 ( a, x, r, p, q, eps )

!*****************************************************************************80
!
!! grat1() evaluates the incomplete Gamma ratio functions P(A,X) and Q(A,X).
!
!  Discussion:
!
!    It is assumed that a <= 1.  eps is the tolerance to be used.
!    the input argument r has the value e**(-x)*x**a/gamma(a).
!
!  Modified:
!
!    28 August 2004
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  real a
  real a2n
  real a2nm1
  real am0
  real an
  real an0
  real b2n
  real b2nm1
  real c
  real cma
  real eps
  real erf
  real erfc1
  real g
  real gam1
  real h
  real j
  real l
  real p
  real q
  real r
  real rexp
  real sum2
  real t
  real tol
  real w
  real x
  real z

  if (a*x == 0.0) go to 130
  if (a == 0.5) go to 120
  if (x < 1.1) go to 10
  go to 50
!
!  Taylor series for p(a,x)/x**a
!
10 continue

      an = 3.0
      c = x
      sum2 = x/(a + 3.0)
      tol = 0.1 * eps / ( a + 1.0 )

      do
        an = an + 1.0
        c = -c*(x/an)
        t = c/(a + an)
        sum2 = sum2 + t
        if ( abs ( t ) <= tol ) then
          exit
        end if
      end do

      j = a*x*((sum2/6.0 - 0.5/(a + 2.0))*x + 1.0/(a + 1.0))

      z = a * log ( x )
      h = gam1(a)
      g = 1.0 + h
      if (x < 0.25) go to 20
         if (a < x/2.59) go to 40
         go to 30
   20 if (z .gt. -0.13394) go to 40

   30 w = exp(z)
      p = w*g*(0.5 + (0.5 - j))
      q = 0.5 + (0.5 - p)
      return

   40 l = rexp(z)
      w = 0.5 + (0.5 + l)
      q = (w*j - l)*g - h
      if (q < 0.0) go to 110
      p = 0.5 + (0.5 - q)
      return
!
!  continued fraction expansion
!
   50 a2nm1 = 1.0
      a2n = 1.0
      b2nm1 = x
      b2n = x + (1.0 - a)
      c = 1.0
   51    a2nm1 = x*a2n + c*a2nm1
         b2nm1 = x*b2n + c*b2nm1
         am0 = a2nm1/b2nm1
         c = c + 1.0
         cma = c - a
         a2n = a2nm1 + cma*a2n
         b2n = b2nm1 + cma*b2n
         an0 = a2n/b2n
         if (abs(an0 - am0) .ge. eps*an0) go to 51
      q = r*an0
      p = 0.5 + (0.5 - q)
      return
!
!  special cases
!
  100 p = 0.0
      q = 1.0
      return

  110 p = 1.0
      q = 0.0
      return

  120 if (x .ge. 0.25) go to 121
      p = erf(sqrt(x))
      q = 0.5 + (0.5 - p)
      return
  121 q = erfc1(0,sqrt(x))
      p = 0.5 + (0.5 - q)
      return

  130 if (x <= a) go to 100
      go to 110

end
function gsumln ( a, b )

!*****************************************************************************80
!
!! gsumln() evaluates the function Log ( Gamma ( A + B ) ) in a special range.
!
!  Discussion:
!
!    The range is for 1 <= a <= 2  and  1 <= b <= 2
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  real a
  real alnrel
  real b
  real gamln1
  real gsumln
  real x

  x = real ( dble ( a ) + dble ( b ) - 2.0d0 )

  if ( x <= 0.25 ) then
    gsumln = gamln1 ( 1.0 + x )
  else if (x <= 1.25 ) then
    gsumln = gamln1 ( x ) + alnrel ( x )
  else
    gsumln = gamln1 ( x - 1.0 ) + log ( x * ( 1.0 + x ) )
  end if

  return
end
function ipmpar ( i )

!*****************************************************************************80
!
!! ipmpar() provides the integer machine constants.
!
!  Discussion:
!
!    it is assumed that the argument i is an integer
!    having one of the values 1-10. ipmpar(i) has the value ...
!
!    integers.
!
!     assume integers are represented in the n-digit, base-a form
!
!               sign ( x(n-1)*a**(n-1) + ... + x(1)*a + x(0) )
!
!               where 0 <= x(i) < a for i=0,...,n-1.
!
!     ipmpar(1) = a, the base.
!
!     ipmpar(2) = n, the number of base-a digits.
!
!     ipmpar(3) = a**n - 1, the largest magnitude.
!
!    floating-point numbers.
!
!     it is assumed that the single and double precision floating
!     point arithmetics have the same base, say b, and that the
!     nonzero numbers are represented in the form
!
!               sign (b**e) * (x(1)/b + ... + x(m)/b**m)
!
!               where x(i) = 0,1,...,b-1 for i=1,...,m,
!               x(1) .ge. 1, and emin <= e <= emax.
!
!     ipmpar(4) = b, the base.
!
!    single-precision
!
!     ipmpar(5) = m, the number of base-b digits.
!
!     ipmpar(6) = emin, the smallest exponent e.
!
!     ipmpar(7) = emax, the largest exponent e.
!
!    double-precision
!
!     ipmpar(8) = m, the number of base-b digits.
!
!     ipmpar(9) = emin, the smallest exponent e.
!
!     ipmpar(10) = emax, the largest exponent e.
!
!     to define this function for the computer being used, activate
!     the data statments for the computer by removing the c from
!     column 1. (all the other data statements should have c in
!     column 1.)
!
!     ipmpar is an adaptation of the function i1mach, written by
!     p.a. fox, a.d. hall, and n.l. schryer (bell laboratories).
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
  integer imach(10)
  integer ipmpar
!
!     machine constants for amdahl machines.
!
!     data imach( 1) /   2 /
!     data imach( 2) /  31 /
!     data imach( 3) / 2147483647 /
!     data imach( 4) /  16 /
!     data imach( 5) /   6 /
!     data imach( 6) / -64 /
!     data imach( 7) /  63 /
!     data imach( 8) /  14 /
!     data imach( 9) / -64 /
!     data imach(10) /  63 /
!
!     machine constants for the at&t 3b series, at&t
!     pc 7300, and at&t 6300.
!
!     data imach( 1) /     2 /
!     data imach( 2) /    31 /
!     data imach( 3) / 2147483647 /
!     data imach( 4) /     2 /
!     data imach( 5) /    24 /
!     data imach( 6) /  -125 /
!     data imach( 7) /   128 /
!     data imach( 8) /    53 /
!     data imach( 9) / -1021 /
!     data imach(10) /  1024 /
!
!     machine constants for the burroughs 1700 system.
!
!     data imach( 1) /    2 /
!     data imach( 2) /   33 /
!     data imach( 3) / 8589934591 /
!     data imach( 4) /    2 /
!     data imach( 5) /   24 /
!     data imach( 6) / -256 /
!     data imach( 7) /  255 /
!     data imach( 8) /   60 /
!     data imach( 9) / -256 /
!     data imach(10) /  255 /
!
!     machine constants for the burroughs 5700 system.
!
!     data imach( 1) /    2 /
!     data imach( 2) /   39 /
!     data imach( 3) / 549755813887 /
!     data imach( 4) /    8 /
!     data imach( 5) /   13 /
!     data imach( 6) /  -50 /
!     data imach( 7) /   76 /
!     data imach( 8) /   26 /
!     data imach( 9) /  -50 /
!     data imach(10) /   76 /
!
!     machine constants for the burroughs 6700/7700 systems.
!
!     data imach( 1) /      2 /
!     data imach( 2) /     39 /
!     data imach( 3) / 549755813887 /
!     data imach( 4) /      8 /
!     data imach( 5) /     13 /
!     data imach( 6) /    -50 /
!     data imach( 7) /     76 /
!     data imach( 8) /     26 /
!     data imach( 9) / -32754 /
!     data imach(10) /  32780 /
!
!     machine constants for the cdc 6000/7000 series
!     60 bit arithmetic, and the cdc cyber 995 64 bit
!     arithmetic (nos operating system).
!
!     data imach( 1) /    2 /
!     data imach( 2) /   48 /
!     data imach( 3) / 281474976710655 /
!     data imach( 4) /    2 /
!     data imach( 5) /   48 /
!     data imach( 6) / -974 /
!     data imach( 7) / 1070 /
!     data imach( 8) /   95 /
!     data imach( 9) / -926 /
!     data imach(10) / 1070 /
!
!     machine constants for the cdc cyber 995 64 bit
!     arithmetic (nos/ve operating system).
!
!     data imach( 1) /     2 /
!     data imach( 2) /    63 /
!     data imach( 3) / 9223372036854775807 /
!     data imach( 4) /     2 /
!     data imach( 5) /    48 /
!     data imach( 6) / -4096 /
!     data imach( 7) /  4095 /
!     data imach( 8) /    96 /
!     data imach( 9) / -4096 /
!     data imach(10) /  4095 /
!
!     machine constants for the cray 1, xmp, 2, and 3.
!
!     data imach( 1) /     2 /
!     data imach( 2) /    63 /
!     data imach( 3) / 9223372036854775807 /
!     data imach( 4) /     2 /
!     data imach( 5) /    47 /
!     data imach( 6) / -8189 /
!     data imach( 7) /  8190 /
!     data imach( 8) /    94 /
!     data imach( 9) / -8099 /
!     data imach(10) /  8190 /
!
!     machine constants for the data general eclipse s/200.
!
!     data imach( 1) /    2 /
!     data imach( 2) /   15 /
!     data imach( 3) / 32767 /
!     data imach( 4) /   16 /
!     data imach( 5) /    6 /
!     data imach( 6) /  -64 /
!     data imach( 7) /   63 /
!     data imach( 8) /   14 /
!     data imach( 9) /  -64 /
!     data imach(10) /   63 /
!
!     machine constants for the harris 220.
!
!     data imach( 1) /    2 /
!     data imach( 2) /   23 /
!     data imach( 3) / 8388607 /
!     data imach( 4) /    2 /
!     data imach( 5) /   23 /
!     data imach( 6) / -127 /
!     data imach( 7) /  127 /
!     data imach( 8) /   38 /
!     data imach( 9) / -127 /
!     data imach(10) /  127 /
!
!     machine constants for the honeywell 600/6000
!     and dps 8/70 series.
!
!     data imach( 1) /    2 /
!     data imach( 2) /   35 /
!     data imach( 3) / 34359738367 /
!     data imach( 4) /    2 /
!     data imach( 5) /   27 /
!     data imach( 6) / -127 /
!     data imach( 7) /  127 /
!     data imach( 8) /   63 /
!     data imach( 9) / -127 /
!     data imach(10) /  127 /
!
!     machine constants for the hp 2100
!     3 word double precision option with ftn4
!
!     data imach( 1) /    2 /
!     data imach( 2) /   15 /
!     data imach( 3) / 32767 /
!     data imach( 4) /    2 /
!     data imach( 5) /   23 /
!     data imach( 6) / -128 /
!     data imach( 7) /  127 /
!     data imach( 8) /   39 /
!     data imach( 9) / -128 /
!     data imach(10) /  127 /
!
!     machine constants for the hp 2100
!     4 word double precision option with ftn4
!
!     data imach( 1) /    2 /
!     data imach( 2) /   15 /
!     data imach( 3) / 32767 /
!     data imach( 4) /    2 /
!     data imach( 5) /   23 /
!     data imach( 6) / -128 /
!     data imach( 7) /  127 /
!     data imach( 8) /   55 /
!     data imach( 9) / -128 /
!     data imach(10) /  127 /
!
!     machine constants for the hp 9000.
!
!     data imach( 1) /     2 /
!     data imach( 2) /    31 /
!     data imach( 3) / 2147483647 /
!     data imach( 4) /     2 /
!     data imach( 5) /    24 /
!     data imach( 6) /  -126 /
!     data imach( 7) /   128 /
!     data imach( 8) /    53 /
!     data imach( 9) / -1021 /
!     data imach(10) /  1024 /
!
!     machine constants for the ibm 360/370 series,
!     the icl 2900, the itel as/6, the xerox sigma
!     5/7/9 and the sel systems 85/86.
!
!     data imach( 1) /    2 /
!     data imach( 2) /   31 /
!     data imach( 3) / 2147483647 /
!     data imach( 4) /   16 /
!     data imach( 5) /    6 /
!     data imach( 6) /  -64 /
!     data imach( 7) /   63 /
!     data imach( 8) /   14 /
!     data imach( 9) /  -64 /
!     data imach(10) /   63 /
!
!     machine constants for the ibm pc.
!
      data imach( 1) /     2 /
      data imach( 2) /    31 /
      data imach( 3) / 2147483647 /
      data imach( 4) /     2 /
      data imach( 5) /    24 /
      data imach( 6) /  -125 /
      data imach( 7) /   128 /
      data imach( 8) /    53 /
      data imach( 9) / -1021 /
      data imach(10) /  1024 /
!
!     machine constants for the macintosh ii - absoft
!     macfortran ii.
!
!     data imach( 1) /     2 /
!     data imach( 2) /    31 /
!     data imach( 3) / 2147483647 /
!     data imach( 4) /     2 /
!     data imach( 5) /    24 /
!     data imach( 6) /  -125 /
!     data imach( 7) /   128 /
!     data imach( 8) /    53 /
!     data imach( 9) / -1021 /
!     data imach(10) /  1024 /
!
!     machine constants for the microvax - vms fortran.
!
!     data imach( 1) /    2 /
!     data imach( 2) /   31 /
!     data imach( 3) / 2147483647 /
!     data imach( 4) /    2 /
!     data imach( 5) /   24 /
!     data imach( 6) / -127 /
!     data imach( 7) /  127 /
!     data imach( 8) /   56 /
!     data imach( 9) / -127 /
!     data imach(10) /  127 /
!
!     machine constants for the pdp-10 (ka processor).
!
!     data imach( 1) /    2 /
!     data imach( 2) /   35 /
!     data imach( 3) / 34359738367 /
!     data imach( 4) /    2 /
!     data imach( 5) /   27 /
!     data imach( 6) / -128 /
!     data imach( 7) /  127 /
!     data imach( 8) /   54 /
!     data imach( 9) / -101 /
!     data imach(10) /  127 /
!
!     machine constants for the pdp-10 (ki processor).
!
!     data imach( 1) /    2 /
!     data imach( 2) /   35 /
!     data imach( 3) / 34359738367 /
!     data imach( 4) /    2 /
!     data imach( 5) /   27 /
!     data imach( 6) / -128 /
!     data imach( 7) /  127 /
!     data imach( 8) /   62 /
!     data imach( 9) / -128 /
!     data imach(10) /  127 /
!
!     machine constants for the pdp-11 fortran supporting
!     32-bit integer arithmetic.
!
!     data imach( 1) /    2 /
!     data imach( 2) /   31 /
!     data imach( 3) / 2147483647 /
!     data imach( 4) /    2 /
!     data imach( 5) /   24 /
!     data imach( 6) / -127 /
!     data imach( 7) /  127 /
!     data imach( 8) /   56 /
!     data imach( 9) / -127 /
!     data imach(10) /  127 /
!
!     machine constants for the sequent balance 8000.
!
!     data imach( 1) /     2 /
!     data imach( 2) /    31 /
!     data imach( 3) / 2147483647 /
!     data imach( 4) /     2 /
!     data imach( 5) /    24 /
!     data imach( 6) /  -125 /
!     data imach( 7) /   128 /
!     data imach( 8) /    53 /
!     data imach( 9) / -1021 /
!     data imach(10) /  1024 /
!
!     machine constants for the silicon graphics iris-4d
!     series (mips r3000 processor).
!
!     data imach( 1) /     2 /
!     data imach( 2) /    31 /
!     data imach( 3) / 2147483647 /
!     data imach( 4) /     2 /
!     data imach( 5) /    24 /
!     data imach( 6) /  -125 /
!     data imach( 7) /   128 /
!     data imach( 8) /    53 /
!     data imach( 9) / -1021 /
!     data imach(10) /  1024 /
!
!     machine constants for the sun 3.
!
!     data imach( 1) /     2 /
!     data imach( 2) /    31 /
!     data imach( 3) / 2147483647 /
!     data imach( 4) /     2 /
!     data imach( 5) /    24 /
!     data imach( 6) /  -125 /
!     data imach( 7) /   128 /
!     data imach( 8) /    53 /
!     data imach( 9) / -1021 /
!     data imach(10) /  1024 /
!
!     machine constants for the univac 1100 series.
!
!     data imach( 1) /    2 /
!     data imach( 2) /   35 /
!     data imach( 3) / 34359738367 /
!     data imach( 4) /    2 /
!     data imach( 5) /   27 /
!     data imach( 6) / -128 /
!     data imach( 7) /  127 /
!     data imach( 8) /   60 /
!     data imach( 9) /-1024 /
!     data imach(10) / 1023 /
!
!     machine constants for the vax 11/780.
!
!     data imach( 1) /    2 /
!     data imach( 2) /   31 /
!     data imach( 3) / 2147483647 /
!     data imach( 4) /    2 /
!     data imach( 5) /   24 /
!     data imach( 6) / -127 /
!     data imach( 7) /  127 /
!     data imach( 8) /   56 /
!     data imach( 9) / -127 /
!     data imach(10) /  127 /
!
      ipmpar = imach(i)

  return
end
function psi ( xx )

!*****************************************************************************80
!
!! psi() evaluates the Psi or Digamma function.
!
!  Discussion:
!
!    psi(xx) is assigned the value 0 when the digamma function cannot
!    be computed.
!
!    the main computation involves evaluation of rational chebyshev
!    approximations published in math. comp. 27, 123-127(1973) by
!    cody, strecok and thacher.
!
!    psi was written at argonne national laboratory for the funpack
!    package of special function subroutines. psi was modified by
!    a.h. morris (nswc).
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
!  Local:
!
!    real ( kind = rk ) DX0, zero of PSI.
!
!    real PIOV4 = pi/4.
!
!    real XMAX1, the smallest positive floating point 
!    constant with entirely integer representation.  Also used as negative of 
!    lower bound on acceptable negative arguments and as the positive 
!    argument beyond which PSI may be represented as log(x).
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real aug
  real den
  real ( kind = rk ), parameter :: dx0 = 1.461632144968362341262659542325721325d0
  integer i
  integer ipmpar
  integer m
  integer n
  integer nq
  real p1(7)
  real p2(4)
  real, parameter :: piov4 = 0.785398163397448e0
  real psi
  real q1(6)
  real q2(4)
  real sgn
  real upper
  real w
  real x
  real xmax1
  real xmx0
  real xsmall
  real xx
  real z
!
!  coefficients for rational approximation of
!  psi(x) / (x - x0),  0.5 <= x <= 3.0
!
      data p1(1)/.895385022981970e-02/,  p1(2)/.477762828042627e+01/, &
           p1(3)/.142441585084029e+03/,  p1(4)/.118645200713425e+04/, &
           p1(5)/.363351846806499e+04/,  p1(6)/.413810161269013e+04/, &
           p1(7)/.130560269827897e+04/
      data q1(1)/.448452573429826e+02/,  q1(2)/.520752771467162e+03/, &
           q1(3)/.221000799247830e+04/,  q1(4)/.364127349079381e+04/, &
           q1(5)/.190831076596300e+04/,  q1(6)/.691091682714533e-05/
!
!  coefficients for rational approximation of
!  psi(x) - ln(x) + 1 / (2*x),  x .gt. 3.0
!
      data p2(1)/-.212940445131011e+01/, p2(2)/-.701677227766759e+01/, &
           p2(3)/-.448616543918019e+01/, p2(4)/-.648157123766197e+00/
      data q2(1)/ .322703493791143e+02/, q2(2)/ .892920700481861e+02/, &
           q2(3)/ .546117738103215e+02/, q2(4)/ .777788548522962e+01/
!
!  machine dependent constants ...
!
!  xsmall = absolute argument below which pi*cotan(pi*x)
!  may be represented by 1/x.
!
      xmax1 = ipmpar(3)
      xmax1 =  min ( xmax1, 1.0 / epsilon ( xmax1 ) )
      xsmall = 1.e-9

      x = xx
      aug = 0.0e0
      if (x .ge. 0.5e0) go to 200
!
!  x < 0.5,  use reflection formula
!  psi(1-x) = psi(x) + pi * cotan(pi*x)
!
      if ( xsmall < abs(x) ) go to 100
      if (x == 0.0e0) go to 400
!
!  0 < abs(x) <= xsmall.  use 1/x as a substitute
!  for  pi*cotan(pi*x)
!
      aug = -1.0e0 / x
      go to 150
!
!  reduction of argument for cotan
!
  100 w = - x
      sgn = piov4
      if ( 0.0e0 < w ) go to 120
      w = - w
      sgn = -sgn
!
!  Error exit if x <= -xmax1
!
  120 if (w .ge. xmax1) go to 400
      nq = int(w)
      w = w - float(nq)
      nq = int(w*4.0e0)
      w = 4.0e0 * (w - float(nq) * .25e0)
!
!  w is now related to the fractional part of  4.0 * x.
!  adjust argument to correspond to values in first
!  quadrant and determine sign
!
      n = nq / 2
      if ((n+n) /= nq) then
        w = 1.0e0 - w
      end if
      z = piov4 * w
      m = n / 2
      if ((m+m) /= n) then 
        sgn = - sgn
      end if
!
!  Determine final value for  -pi*cotan(pi*x)
!
      n = (nq + 1) / 2
      m = n / 2
      m = m + m
      if (m /= n) go to 140
!
!  Check for singularity
!
      if (z == 0.0e0) go to 400
!
!  use cos/sin as a substitute for cotan, and
!  sin/cos as a substitute for tan
!
      aug = sgn * ((cos(z) / sin(z)) * 4.0e0)
      go to 150
  140 aug = sgn * ((sin(z) / cos(z)) * 4.0e0)
  150 x = 1.0e0 - x
  200 if (x .gt. 3.0e0) go to 300
!
!  0.5 <= x <= 3.0
!
      den = x
      upper = p1(1) * x

      do i = 1, 5
         den = (den + q1(i)) * x
         upper = (upper + p1(i+1)) * x
      end do

      den = (upper + p1(7)) / (den + q1(6))
      xmx0 = real ( dble(x) - dx0 )
      psi = den * xmx0 + aug
      return
!
!  if x .ge. xmax1, psi = ln(x)
!
  300 if (x .ge. xmax1) go to 350
!
!  3.0 < x < xmax1
!
      w = 1.0e0 / (x * x)
      den = w
      upper = p2(1) * w

      do i = 1, 3
         den = (den + q2(i)) * w
         upper = (upper + p2(i+1)) * w
      end do

      aug = upper / (den + q2(4)) - 0.5e0 / x + aug
  350 psi = aug + log ( x )
      return
!
!  error return
!
  400 psi = 0.0e0

  return
end
subroutine psi_values ( n, x, fx )

!*****************************************************************************80
!
!! psi_values() returns some values of the Psi or Digamma function for testing.
!
!  Discussion:
!
!    PSI(X) = d LN ( GAMMA ( X ) ) / d X = GAMMA'(X) / GAMMA(X)
!
!    PSI(1) = - Euler's constant.
!
!    PSI(X+1) = PSI(X) + 1 / X.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    17 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz and Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964.
!
!  Parameters:
!
!    Input/output, integer N.
!    On input, if N is 0, the first test data is returned, and N is set
!    to the index of the test data.  On each subsequent call, N is
!    incremented and that test data is returned.  When there is no more
!    test data, N is set to 0.
!
!    Output, real X, the argument of the function.
!
!    Output, real FX, the value of the function.
!
  implicit none

  integer, parameter :: nmax = 11

  real fx
  real, save, dimension ( nmax ) :: fxvec = (/ &
    -0.5772156649E+00, -0.4237549404E+00, -0.2890398966E+00, &
    -0.1691908889E+00, -0.0613845446E+00, -0.0364899740E+00, &
     0.1260474528E+00,  0.2085478749E+00,  0.2849914333E+00, &  
     0.3561841612E+00,  0.4227843351E+00 /)
  integer n
  real x
  real, save, dimension ( nmax ) :: xvec = (/ &
    1.0E+00,  1.1E+00,  1.2E+00,  &
    1.3E+00,  1.4E+00,  1.5E+00,  &
    1.6E+00,  1.7E+00,  1.8E+00,  &
    1.9E+00,  2.0E+00 /)

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

  n = n + 1

  if ( nmax < n ) then
    n = 0
    x = 0.0E+00
    fx = 0.0E+00
    return
  end if

  x = xvec(n)
  fx = fxvec(n)

  return
end
function rexp ( x )

!*****************************************************************************80
!
!! rexp() evaluates the function Exp(X) - 1.
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  real, parameter :: p1 =  0.914041914819518e-09
  real, parameter :: p2 =  0.238082361044469e-01
  real, parameter :: q1 = -0.499999999085958e+00
  real, parameter :: q2 =  0.107141568980644e+00
  real, parameter :: q3 = -0.119041179760821e-01
  real, parameter :: q4 =  0.595130811860248e-03
  real rexp
  real w
  real x

  if ( abs ( x ) <= 0.15 ) then

    rexp = x*(((p2*x + p1)*x + 1.0) / ((((q4*x + q3)*x + q2)*x &
      + q1)*x + 1.0))

  else if ( x < 0.0 ) then

    w = exp ( x )
    rexp = ( w - 0.5 ) - 0.5

  else

    w = exp ( x )
    rexp = w * ( 0.5 + ( 0.5 - 1.0 / w ) )

  end if

  return
end
function rlog1 ( x )

!*****************************************************************************80
!
!! rlog1() evaluates the function X - Log ( 1 + X ).
!
!  Author:
!
!    Original Fortran77 version by Alfred Morris.
!
!  Reference:
!
!    Armido Didonato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, Number 3, 1992, pages 360-373.
!
  implicit none

  real, parameter :: a = 0.566749439387324e-01
  real, parameter :: b = 0.456512608815524e-01
  real h
  real, parameter :: p0 =  0.333333333333333e+00
  real, parameter :: p1 = -0.224696413112536e+00
  real, parameter :: p2 =  0.620886815375787e-02
  real, parameter :: q1 = -0.127408923933623e+01
  real, parameter :: q2 =  0.354508718369557e+00
  real r
  real rlog1
  real t
  real w
  real w1
  real x

  if ( x < -0.39 ) then
    w = ( x + 0.5 ) + 0.5
    rlog1 = x - log ( w )
  else if ( x < -0.18 ) then
    h = real ( ( dble ( x ) + 0.3d0 ) / 0.7 )
    w1 = a - h*0.3
    r = h/(h + 2.0)
    t = r * r
    w = ( ( p2 * t + p1 ) * t + p0 ) / ( ( q2 * t + q1 ) * t + 1.0 )
    rlog1 = 2.0*t*(1.0/(1.0 - r) - r*w) + w1
  else if ( x <= 0.18 ) then
    h = x
    w1 = 0.0
    r = h/(h + 2.0)
    t = r*r
    w = ((p2*t + p1)*t + p0)/((q2*t + q1)*t + 1.0 )
    rlog1 = 2.0*t*(1.0/(1.0 - r) - r*w) + w1
  else if ( x <= 0.57 ) then
    h = real ( 0.75d0 * dble ( x ) - 0.25d0 )
    w1 = b + h/3.0
    r = h/(h + 2.0)
    t = r*r
    w = ((p2*t + p1)*t + p0)/((q2*t + q1)*t + 1.0)
    rlog1 = 2.0 * t * ( 1.0 / ( 1.0 - r ) - r * w ) + w1
  else
    w = ( x + 0.5 ) + 0.5
    rlog1 = x - log ( w )
  end if

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! timestamp() prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    31 May 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    18 May 2013
!
!  Author:
!
!    John Burkardt
!
  implicit none

  character ( len = 8 ) ampm
  integer d
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  integer values(8)
  integer y

  call date_and_time ( values = values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end