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).
!
!  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 ( kind = 4 ) a
  real ( kind = 4 ) algdiv
  real ( kind = 4 ) alnrel
  real ( kind = 4 ) b
  real ( kind = 4 ) c
  real ( kind = 4 ), parameter :: c0 =  0.833333333333333E-01
  real ( kind = 4 ), parameter :: c1 = -0.277777777760991E-02
  real ( kind = 4 ), parameter :: c2 =  0.793650666825390E-03
  real ( kind = 4 ), parameter :: c3 = -0.595202931351870E-03
  real ( kind = 4 ), parameter :: c4 =  0.837308034031215E-03
  real ( kind = 4 ), parameter :: c5 = -0.165322962780713E-02
  real ( kind = 4 ) d
  real ( kind = 4 ) h
  real ( kind = 4 ) s11
  real ( kind = 4 ) s3
  real ( kind = 4 ) s5
  real ( kind = 4 ) s7
  real ( kind = 4 ) s9
  real ( kind = 4 ) t
  real ( kind = 4 ) u
  real ( kind = 4 ) v
  real ( kind = 4 ) w
  real ( kind = 4 ) x
  real ( kind = 4 ) 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).
!
!  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 ( kind = 4 ) a
  real ( kind = 4 ) alnrel
  real ( kind = 4 ), parameter :: p1 = -0.129418923021993E+01
  real ( kind = 4 ), parameter :: p2 =  0.405303492862024E+00
  real ( kind = 4 ), parameter :: p3 = -0.178874546012214E-01
  real ( kind = 4 ), parameter :: q1 = -0.162752256355323E+01
  real ( kind = 4 ), parameter :: q2 =  0.747811014037616E+00
  real ( kind = 4 ), parameter :: q3 = -0.845104217945565E-01
  real ( kind = 4 ) t
  real ( kind = 4 ) t2
  real ( kind = 4 ) w
  real ( kind = 4 ) 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 = 1.0E+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) for
!     a <= min(eps,eps*b), b*x <= 1, and x <= 0.5. used when
!     a is very small. use only if above inequalities are satisfied.
!
!  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 ( kind = 4 ) a
  real ( kind = 4 ) aj
  real ( kind = 4 ) apser
  real ( kind = 4 ) b
  real ( kind = 4 ) bx
  real ( kind = 4 ) c
  real ( kind = 4 ) eps
  real ( kind = 4 ), parameter :: g = 0.577215664901533E+00
  real ( kind = 4 ) j
  real ( kind = 4 ) psi
  real ( kind = 4 ) s
  real ( kind = 4 ) t
  real ( kind = 4 ) tol
  real ( kind = 4 ) 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.
!
!  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 ( kind = 4 ) a
  real ( kind = 4 ) a0(21)
  real ( kind = 4 ) b
  real ( kind = 4 ) b0(21)
  real ( kind = 4 ) bcorr
  real ( kind = 4 ) basym
  real ( kind = 4 ) bsum
  real ( kind = 4 ) c(21)
  real ( kind = 4 ) d(21)
  real ( kind = 4 ) dsum
  real ( kind = 4 ), parameter :: e0 = 1.12837916709551E+00
  real ( kind = 4 ), parameter :: e1 = 0.353553390593274E+00
  real ( kind = 4 ) eps
  real ( kind = 4 ) erfc1
  real ( kind = 4 ) f
  real ( kind = 4 ) h
  real ( kind = 4 ) h2
  real ( kind = 4 ) hn
  integer ( kind = 4 ) i
  integer ( kind = 4 ) im1
  integer ( kind = 4 ) imj
  integer ( kind = 4 ) j
  real ( kind = 4 ) j0
  real ( kind = 4 ) j1
  real ( kind = 4 ) lambda
  integer ( kind = 4 ) m
  integer ( kind = 4 ) mm1
  integer ( kind = 4 ) mmj
  integer ( kind = 4 ) n
  integer ( kind = 4 ) np1
  integer ( kind = 4 ), parameter :: num = 20
  real ( kind = 4 ) r
  real ( kind = 4 ) r0
  real ( kind = 4 ) r1
  real ( kind = 4 ) rlog1
  real ( kind = 4 ) s
  real ( kind = 4 ) sum2
  real ( kind = 4 ) t
  real ( kind = 4 ) t0
  real ( kind = 4 ) t1
  real ( kind = 4 ) u
  real ( kind = 4 ) w
  real ( kind = 4 ) w0
  real ( kind = 4 ) z
  real ( kind = 4 ) z0
  real ( kind = 4 ) z2
  real ( kind = 4 ) zn
  real ( kind = 4 ) 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, kind = 4 )
    np1 = n + 1
    s = s + hn
    a0(np1) = 2.0 * r1 * s / real ( n + 3, kind = 4 )

    do i = n, np1

      r = -0.5 * real ( i + 1, kind = 4 )
      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)  where
!     ln(gamma(a)) = (a - 0.5)*ln(a) - a + 0.5*ln(2*pi) + del(a).
!     it is assumed that a0 .ge. 8 and b0 .ge. 8.
!
!  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 ( kind = 4 ) a
  real ( kind = 4 ) a0
  real ( kind = 4 ) b
  real ( kind = 4 ) b0
  real ( kind = 4 ) bcorr
  real ( kind = 4 ) c
  real ( kind = 4 ), parameter :: c0 =  0.833333333333333e-01
  real ( kind = 4 ), parameter :: c1 = -0.277777777760991e-02
  real ( kind = 4 ), parameter :: c2 =  0.793650666825390e-03
  real ( kind = 4 ), parameter :: c3 = -0.595202931351870e-03
  real ( kind = 4 ), parameter :: c4 =  0.837308034031215e-03
  real ( kind = 4 ), parameter :: c5 = -0.165322962780713e-02
  real ( kind = 4 ) h
  real ( kind = 4 ) s11
  real ( kind = 4 ) s3
  real ( kind = 4 ) s5
  real ( kind = 4 ) s7
  real ( kind = 4 ) s9
  real ( kind = 4 ) t
  real ( kind = 4 ) w
  real ( kind = 4 ) x
  real ( kind = 4 ) 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 ]
!
!  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 ( kind = 4 ) 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 = 4 ) A, B, the parameters of the function.
!
!    Output, real ( kind = 4 ) X, the argument of the function.
!
!    Output, real ( kind = 4 ) FX, the value of the function.
!
  implicit none

  integer ( kind = 4 ), parameter :: n_max = 42

  real ( kind = 4 ) a
  real ( kind = 4 ), 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 ( kind = 4 ) b
  real ( kind = 4 ), 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 ( kind = 4 ) fx
  real ( kind = 4 ), 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 ( kind = 4 ) n_data
  real ( kind = 4 ) x
  real ( kind = 4 ), 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.
!
!  Modified:
!
!    15 February 2004
!
!  Reference:
!
!    Milton Abramowitz and Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964.
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer ( kind = 4 ) 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 ( kind = 4 ) X, Y, the arguments of the function.
!
!    Output, real ( kind = 4 ) FXY, the value of the function.
!
  implicit none

  integer ( kind = 4 ), parameter :: nmax = 17

  real ( kind = 4 ), 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 ( kind = 4 ) fxy
  integer ( kind = 4 ) n
  real ( kind = 4 ) x
  real ( kind = 4 ), 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 ( kind = 4 ) y
  real ( kind = 4 ), 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.
!
!  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 Parameters:
!
!    Local, real ( kind = 4 ) E, the value of Log ( 2 * PI ) / 2.
!
  implicit none

  real ( kind = 4 ) a
  real ( kind = 4 ) a0
  real ( kind = 4 ) algdiv
  real ( kind = 4 ) alnrel
  real ( kind = 4 ) b
  real ( kind = 4 ) b0
  real ( kind = 4 ) bcorr
  real ( kind = 4 ) betaln
  real ( kind = 4 ) c
  real ( kind = 4 ), parameter :: e = 0.918938533204673E+00
  real ( kind = 4 ) gamln
  real ( kind = 4 ) gsumln
  real ( kind = 4 ) h
  integer ( kind = 4 ) i
  integer ( kind = 4 ) n
  real ( kind = 4 ) u
  real ( kind = 4 ) v
  real ( kind = 4 ) w
  real ( kind = 4 ) 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 = 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 = 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 = 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.
!
!     it is assumed that  lambda = (a + b)*y - b.
!
!  Modified:
!
!    28 August 2004
!
!  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 ( kind = 4 ) a
  real ( kind = 4 ) alpha
  real ( kind = 4 ) an
  real ( kind = 4 ) anp1
  real ( kind = 4 ) b
  real ( kind = 4 ) beta
  real ( kind = 4 ) bfrac
  real ( kind = 4 ) bn
  real ( kind = 4 ) bnp1
  real ( kind = 4 ) brcomp
  real ( kind = 4 ) c
  real ( kind = 4 ) c0
  real ( kind = 4 ) c1
  real ( kind = 4 ) e
  real ( kind = 4 ) eps
  real ( kind = 4 ) lambda
  real ( kind = 4 ) n
  real ( kind = 4 ) p
  real ( kind = 4 ) r
  real ( kind = 4 ) r0
  real ( kind = 4 ) s
  real ( kind = 4 ) t
  real ( kind = 4 ) w
  real ( kind = 4 ) x
  real ( kind = 4 ) y
  real ( kind = 4 ) 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.
!
!  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 ( kind = 4 ) a
  real ( kind = 4 ) algdiv
  real ( kind = 4 ) alnrel
  real ( kind = 4 ) b
  real ( kind = 4 ) bm1
  real ( kind = 4 ) bp2n
  real ( kind = 4 ) c(30)
  real ( kind = 4 ) cn
  real ( kind = 4 ) coef
  real ( kind = 4 ) d(30)
  real ( kind = 4 ) dj
  real ( kind = 4 ) eps
  real ( kind = 4 ) gam1
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ierr
  real ( kind = 4 ) j
  real ( kind = 4 ) l
  real ( kind = 4 ) lnx
  integer ( kind = 4 ) n
  real ( kind = 4 ) n2
  integer ( kind = 4 ) nm1
  real ( kind = 4 ) nu
  real ( kind = 4 ) p
  real ( kind = 4 ) q
  real ( kind = 4 ) r
  real ( kind = 4 ) s
  real ( kind = 4 ) sum1
  real ( kind = 4 ) t
  real ( kind = 4 ) t2
  real ( kind = 4 ) u
  real ( kind = 4 ) v
  real ( kind = 4 ) w
  real ( kind = 4 ) x
  real ( kind = 4 ) y
  real ( kind = 4 ) 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) when b <= 1
!     or b*x <= 0.7.  eps is the tolerance used.
!
!  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 ( kind = 4 ) a
  real ( kind = 4 ) a0
  real ( kind = 4 ) algdiv
  real ( kind = 4 ) apb
  real ( kind = 4 ) b
  real ( kind = 4 ) b0
  real ( kind = 4 ) betaln
  real ( kind = 4 ) bpser
  real ( kind = 4 ) c
  real ( kind = 4 ) eps
  real ( kind = 4 ) gam1
  real ( kind = 4 ) gamln1
  integer ( kind = 4 ) i
  integer ( kind = 4 ) m
  real ( kind = 4 ) n
  real ( kind = 4 ) sum1
  real ( kind = 4 ) t
  real ( kind = 4 ) tol
  real ( kind = 4 ) u
  real ( kind = 4 ) w
  real ( kind = 4 ) x
  real ( kind = 4 ) 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 = 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 = 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 = 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:
!
!    Alfred Morris
!    Naval Surface Warfare Center
!    Dahlgren, Virginia
!
!  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 ( kind = 4 ) A, B, the parameters of the function.  A and B
!    should be nonnegative.
!
  implicit none

  real ( kind = 4 ) a
  real ( kind = 4 ) a0
  real ( kind = 4 ) apser
  real ( kind = 4 ) b
  real ( kind = 4 ) b0
  real ( kind = 4 ) basym
  real ( kind = 4 ) bfrac
  real ( kind = 4 ) bpser
  real ( kind = 4 ) bup
  real ( kind = 4 ) eps
  real ( kind = 4 ) fpser
  integer ( kind = 4 ) ierr
  integer ( kind = 4 ) ierr1
  integer ( kind = 4 ) ind
  real ( kind = 4 ) lambda
  integer ( kind = 4 ) n
  real ( kind = 4 ) t
  real ( kind = 4 ) w
  real ( kind = 4 ) w1
  real ( kind = 4 ) x
  real ( kind = 4 ) x0
  real ( kind = 4 ) y
  real ( kind = 4 ) y0
  real ( kind = 4 ) 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 = b0
      b0 = b0 - n

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

  141 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)).
!
!  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 ( kind = 4 ) a
  real ( kind = 4 ) a0
  real ( kind = 4 ) algdiv
  real ( kind = 4 ) alnrel
  real ( kind = 4 ) apb
  real ( kind = 4 ) b
  real ( kind = 4 ) b0
  real ( kind = 4 ) bcorr
  real ( kind = 4 ) betaln
  real ( kind = 4 ) brcmp1
  real ( kind = 4 ) c
  real ( kind = 4 ), parameter :: const = 0.398942280401433
  real ( kind = 4 ) e
  real ( kind = 4 ) esum
  real ( kind = 4 ) gam1
  real ( kind = 4 ) gamln1
  real ( kind = 4 ) h
  integer ( kind = 4 ) i
  real ( kind = 4 ) lambda
  real ( kind = 4 ) lnx
  real ( kind = 4 ) lny
  integer ( kind = 4 ) mu
  integer ( kind = 4 ) n
  real ( kind = 4 ) rlog1
  real ( kind = 4 ) t
  real ( kind = 4 ) u
  real ( kind = 4 ) v
  real ( kind = 4 ) x
  real ( kind = 4 ) x0
  real ( kind = 4 ) y
  real ( kind = 4 ) y0
  real ( kind = 4 ) 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 = 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 = 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 = 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).
!
!  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 ( kind = 4 ) a
  real ( kind = 4 ) a0
  real ( kind = 4 ) algdiv
  real ( kind = 4 ) alnrel
  real ( kind = 4 ) apb
  real ( kind = 4 ) b
  real ( kind = 4 ) b0
  real ( kind = 4 ) bcorr
  real ( kind = 4 ) betaln
  real ( kind = 4 ) brcomp
  real ( kind = 4 ) c
  real ( kind = 4 ), parameter :: const = 0.398942280401433
  real ( kind = 4 ) e
  real ( kind = 4 ) gam1
  real ( kind = 4 ) gamln1
  real ( kind = 4 ) h
  integer ( kind = 4 ) i
  real ( kind = 4 ) lambda
  real ( kind = 4 ) lnx
  real ( kind = 4 ) lny
  integer ( kind = 4 ) n
  real ( kind = 4 ) rlog1
  real ( kind = 4 ) t
  real ( kind = 4 ) u
  real ( kind = 4 ) v
  real ( kind = 4 ) x
  real ( kind = 4 ) x0
  real ( kind = 4 ) y
  real ( kind = 4 ) y0
  real ( kind = 4 ) 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
!-----------------------------------------------------------------------
   30 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 = 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 = 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 = 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.
!
!     eps is the tolerance used.
!
!  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 ( kind = 4 ) a
  real ( kind = 4 ) ap1
  real ( kind = 4 ) apb
  real ( kind = 4 ) b
  real ( kind = 4 ) brcmp1
  real ( kind = 4 ) bup
  real ( kind = 4 ) d
  real ( kind = 4 ) eps
  real ( kind = 4 ) exparg
  integer ( kind = 4 ) i
  integer ( kind = 4 ) k
  integer ( kind = 4 ) kp1
  real ( kind = 4 ) l
  integer ( kind = 4 ) mu
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nm1
  real ( kind = 4 ) r
  real ( kind = 4 ) t
  real ( kind = 4 ) w
  real ( kind = 4 ) x
  real ( kind = 4 ) 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 = abs(exparg(1))
         k = 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 = 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
!
!          erfc1(ind,x) = erfc(x)            if ind = 0
!          erfc1(ind,x) = exp(x*x)*erfc(x)   otherwise
!
!  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 ( kind = 4 ) a(5)
  real ( kind = 4 ) ax
  real ( kind = 4 ) b(3)
  real ( kind = 4 ) bot
  real ( kind = 4 ), parameter :: c = 0.564189583547756
  real ( kind = 4 ) e
  real ( kind = 4 ) erfc1
  real ( kind = 4 ) exparg
  integer ( kind = 4 ) ind
  real ( kind = 4 ) p(8)
  real ( kind = 4 ) q(8)
  real ( kind = 4 ) r(5)
  real ( kind = 4 ) s(4)
  real ( kind = 4 ) t
  real ( kind = 4 ) top
  double precision w
  real ( kind = 4 ) 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 = w
      e = 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 erf ( x )

!*****************************************************************************80
!
!! ERF evaluates the real error function.
!
!  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 ( kind = 4 ) a(5)
  real ( kind = 4 ) ax
  real ( kind = 4 ) b(3)
  real ( kind = 4 ) bot
  real ( kind = 4 ), parameter :: c = 0.564189583547756
  real ( kind = 4 ) erf
  real ( kind = 4 ) p(8)
  real ( kind = 4 ) q(8)
  real ( kind = 4 ) r(5)
  real ( kind = 4 ) s(4)
  real ( kind = 4 ) t
  real ( kind = 4 ) top
  real ( kind = 4 ) x
  real ( kind = 4 ) 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
    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)
    erf = 0.5 + (0.5 - exp(-x*x)*top/bot)
    if ( x < 0.0 ) then
      erf = -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
    erf = (c - top/(x2*bot)) / ax
    erf = 0.5 + (0.5 - exp(-x2)*erf)

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

  else

   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.
!
!  Modified:
!
!    17 April 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer ( kind = 4 ) 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 ( kind = 4 ) X, the argument of the function.
!
!    Output, real ( kind = 4 ) FX, the value of the function.
!
  implicit none

  integer ( kind = 4 ), parameter :: nmax = 21

  real ( kind = 4 ), 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 ( kind = 4 ) fx
  integer ( kind = 4 ) n_data
  real ( kind = 4 ) x
  real ( kind = 4 ), 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).
!
!  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 ( kind = 4 ) esum
  integer ( kind = 4 ) mu
  real ( kind = 4 ) w
  real ( kind = 4 ) 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).
!
!    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.
!
!    note... only an approximate value for exparg(l) is needed.
!
!  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 ( kind = 4 ) b
  real ( kind = 4 ) exparg
  integer ( kind = 4 ) ipmpar
  integer ( kind = 4 ) l
  integer ( kind = 4 ) m
  real ( kind = 4 ) 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.
!
!          for b < min(eps,eps*a) and x <= 0.5.
!
!  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 ( kind = 4 ) a
  real ( kind = 4 ) an
  real ( kind = 4 ) b
  real ( kind = 4 ) c
  real ( kind = 4 ) eps
  real ( kind = 4 ) exparg
  real ( kind = 4 ) fpser
  real ( kind = 4 ) s
  real ( kind = 4 ) t
  real ( kind = 4 ) tol
  real ( kind = 4 ) 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
!
!  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 ( kind = 4 ) a
  real ( kind = 4 ) bot
  real ( kind = 4 ) d
  real ( kind = 4 ) gam1
  real ( kind = 4 ) p(7)
  real ( kind = 4 ) q(5)
  real ( kind = 4 ) r(9)
  real ( kind = 4 ) s1
  real ( kind = 4 ) s2
  real ( kind = 4 ) t
  real ( kind = 4 ) top
  real ( kind = 4 ) 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) t = d - 0.5
      if (t) 30,10,20

   10 gam1 = 0.0
      return

   20 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 .gt. 0.0) go to 21
         gam1 = a*w
         return
   21 gam1 = (t/a)*((w - 0.5) - 0.5)
      return

   30 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 ( 0.0 < d ) go to 31
         gam1 = a*((w + 0.5) + 0.5)
         return
   31 gam1 = t*w/a

  return
end
function gamln1 ( a )

!*****************************************************************************80
!
!! GAMLN1 evaluates ln(gamma(1 + a)) for -0.2 <= A <= 1.25
!
!  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 ( kind = 4 ) a
  real ( kind = 4 ) gamln1
  real ( kind = 4 ) p0
  real ( kind = 4 ) p1
  real ( kind = 4 ) p2
  real ( kind = 4 ) p3
  real ( kind = 4 ) p4
  real ( kind = 4 ) p5
  real ( kind = 4 ) p6
  real ( kind = 4 ) q1
  real ( kind = 4 ) q2
  real ( kind = 4 ) q3
  real ( kind = 4 ) q4
  real ( kind = 4 ) q5
  real ( kind = 4 ) q6
  real ( kind = 4 ) r0
  real ( kind = 4 ) r1
  real ( kind = 4 ) r2
  real ( kind = 4 ) r3
  real ( kind = 4 ) r4
  real ( kind = 4 ) r5
  real ( kind = 4 ) s1
  real ( kind = 4 ) s2
  real ( kind = 4 ) s3
  real ( kind = 4 ) s4
  real ( kind = 4 ) s5
  real ( kind = 4 ) w
  real ( kind = 4 ) 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:
!
!    Alfred Morris
!    Naval Surface Warfare Center
!    Dahlgren, Virginia
!
!  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 ( kind = 4 ) a
  real ( kind = 4 ) c0
  real ( kind = 4 ) c1
  real ( kind = 4 ) c2
  real ( kind = 4 ) c3
  real ( kind = 4 ) c4
  real ( kind = 4 ) c5
  real ( kind = 4 ), parameter :: d = 0.418938533204673
  real ( kind = 4 ) gamln
  real ( kind = 4 ) gamln1
  integer ( kind = 4 ) i
  integer ( kind = 4 ) n
  real ( kind = 4 ) t
  real ( kind = 4 ) 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 = 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]
!
!  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 ( kind = 4 ) 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 = 4 ) A, the parameter of the function.
!
!    Output, real ( kind = 4 ) X, the argument of the function.
!
!    Output, real ( kind = 4 ) FX, the value of the function.
!
  implicit none

  integer ( kind = 4 ), parameter :: n_max = 20

  real ( kind = 4 ) a
  real ( kind = 4 ), 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 ( kind = 4 ) fx
  real ( kind = 4 ), 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 ( kind = 4 ) n_data
  real ( kind = 4 ) x
  real ( kind = 4 ), 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.
!
!  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 ( kind = 4 ) 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 ( kind = 4 ) X, the argument of the function.
!
!    Output, real ( kind = 4 ) FX, the value of the function.
!
  implicit none

  integer ( kind = 4 ), parameter :: nmax = 18

  real ( kind = 4 ), 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 ( kind = 4 ) fx
  integer ( kind = 4 ) n_data
  real ( kind = 4 ) x
  real ( kind = 4 ), 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
!
!  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 ( kind = 4 ) a
  real ( kind = 4 ) a2n
  real ( kind = 4 ) a2nm1
  real ( kind = 4 ) am0
  real ( kind = 4 ) an
  real ( kind = 4 ) an0
  real ( kind = 4 ) b2n
  real ( kind = 4 ) b2nm1
  real ( kind = 4 ) c
  real ( kind = 4 ) cma
  real ( kind = 4 ) eps
  real ( kind = 4 ) erf
  real ( kind = 4 ) erfc1
  real ( kind = 4 ) g
  real ( kind = 4 ) gam1
  real ( kind = 4 ) h
  real ( kind = 4 ) j
  real ( kind = 4 ) l
  real ( kind = 4 ) p
  real ( kind = 4 ) q
  real ( kind = 4 ) r
  real ( kind = 4 ) rexp
  real ( kind = 4 ) sum2
  real ( kind = 4 ) t
  real ( kind = 4 ) tol
  real ( kind = 4 ) w
  real ( kind = 4 ) x
  real ( kind = 4 ) 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.
!
!          for 1 <= a <= 2  and  1 <= b <= 2
!
!  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 ( kind = 4 ) a
  real ( kind = 4 ) alnrel
  real ( kind = 4 ) b
  real ( kind = 4 ) gamln1
  real ( kind = 4 ) gsumln
  real ( kind = 4 ) x

  x = 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 ( kind = 4 ) machine constants for the computer
!     that is used. it is assumed that the argument i is an integer ( kind = 4 )
!     having one of the values 1-10. ipmpar(i) has the value ...
!
!  integer ( kind = 4 )s.
!
!     assume integer ( kind = 4 )s 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).
!     ipmpar was formed by a.h. morris (nswc). the constants are
!     from bell laboratories, nswc, and other sources.
!
      integer ( kind = 4 ) imach(10)
  integer ( kind = 4 ) 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 ( kind = 4 ) 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.
!
!     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).
!
!  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 parameters:
!
!    Local, real ( kind = 8 ) DX0, zero of PSI.
!
!    Local, real ( kind = 4 ) PIOV4 = pi/4.
!
!    Local, real ( kind = 4 ) 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

  real ( kind = 4 ) aug
  real ( kind = 4 ) den
  real ( kind = 8 ), parameter :: dx0 = 1.461632144968362341262659542325721325d0
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ipmpar
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nq
  real ( kind = 4 ) p1(7)
  real ( kind = 4 ) p2(4)
  real ( kind = 4 ), parameter :: piov4 = 0.785398163397448e0
  real ( kind = 4 ) psi
  real ( kind = 4 ) q1(6)
  real ( kind = 4 ) q2(4)
  real ( kind = 4 ) sgn
  real ( kind = 4 ) upper
  real ( kind = 4 ) w
  real ( kind = 4 ) x
  real ( kind = 4 ) xmax1
  real ( kind = 4 ) xmx0
  real ( kind = 4 ) xsmall
  real ( kind = 4 ) xx
  real ( kind = 4 ) 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 = 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.
!
!  Modified:
!
!    17 May 2001
!
!  Reference:
!
!    Milton Abramowitz and Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964.
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer ( kind = 4 ) 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 ( kind = 4 ) X, the argument of the function.
!
!    Output, real ( kind = 4 ) FX, the value of the function.
!
  implicit none

  integer ( kind = 4 ), parameter :: nmax = 11

  real ( kind = 4 ) fx
  real ( kind = 4 ), 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 ( kind = 4 ) n
  real ( kind = 4 ) x
  real ( kind = 4 ), 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.
!
!  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 ( kind = 4 ), parameter :: p1 =  0.914041914819518e-09
  real ( kind = 4 ), parameter :: p2 =  0.238082361044469e-01
  real ( kind = 4 ), parameter :: q1 = -0.499999999085958e+00
  real ( kind = 4 ), parameter :: q2 =  0.107141568980644e+00
  real ( kind = 4 ), parameter :: q3 = -0.119041179760821e-01
  real ( kind = 4 ), parameter :: q4 =  0.595130811860248e-03
  real ( kind = 4 ) rexp
  real ( kind = 4 ) w
  real ( kind = 4 ) 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 ).
!
!  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 ( kind = 4 ), parameter :: a = 0.566749439387324e-01
  real ( kind = 4 ), parameter :: b = 0.456512608815524e-01
  real ( kind = 4 ) h
  real ( kind = 4 ), parameter :: p0 =  0.333333333333333e+00
  real ( kind = 4 ), parameter :: p1 = -0.224696413112536e+00
  real ( kind = 4 ), parameter :: p2 =  0.620886815375787e-02
  real ( kind = 4 ), parameter :: q1 = -0.127408923933623e+01
  real ( kind = 4 ), parameter :: q2 =  0.354508718369557e+00
  real ( kind = 4 ) r
  real ( kind = 4 ) rlog1
  real ( kind = 4 ) t
  real ( kind = 4 ) w
  real ( kind = 4 ) w1
  real ( kind = 4 ) x

  if ( x < -0.39 ) then
    w = ( x + 0.5 ) + 0.5
    rlog1 = x - log ( w )
  else if ( x < -0.18 ) then
    h = ( 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 = 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 GNU LGPL license.
!
!  Modified:
!
!    18 May 2013
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer ( kind = 4 ) d
  integer ( kind = 4 ) h
  integer ( kind = 4 ) m
  integer ( kind = 4 ) mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer ( kind = 4 ) n
  integer ( kind = 4 ) s
  integer ( kind = 4 ) values(8)
  integer ( kind = 4 ) 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