/* ---------------------------------------------------------------
!     ALGORITHM AS 27  APPL. STATIST. VOL.19, NO.1
!
!     Calculate the upper tail area under Student's t-distribution
!
!     Translated from Algol by Alan Miller
! --------------------------------------------------------------*/
      float V, X, TT;
      double A1, A2, A3, A4, A5, B1, B2, C1, C2, C3, C4, C5, D1, D2;
      double E1, E2, E3, E4, E5, F1, F2, G1, G2, G3, G4, G5, H1, H2;
      double I1, I2, I3, I4, I5, J1, J2;
      bool POS;

      A1=0.09979441; A2=-0.581821; A3=1.390993; A4=-1.222452; A5=2.151185;
      B1=5.537409; B2=11.42343;
      C1=0.04431742; C2=-0.2206018; C3=-0.03317253; C4=5.679969; C5=-12.96519;
      D1=5.166733; D2=13.49862;
      E1=0.009694901; E2=-0.1408854; E3=1.88993; E4=-12.75532; E5=25.77532;
      F1=4.233736; F2=14.3963;
      G1=-9.187228E-5; G2=0.03789901; G3=-1.280346; G4=9.249528; G5=-19.08115;
      H1=2.777816; H2=16.46132;
      I1=5.79602E-4; I2=-0.02763334; I3=0.4517029; I4=-2.657697; I5=5.127212;
      J1=0.5657187; J2=21.83269;


!*****************************************************************************80
!
!! GAMMA_SAMPLE samples the Gamma PDF.
!
!  Modified:
!
!    13 February 1999
!
!  Reference:
!
!    Joachim Ahrens, Ulrich Dieter,
!    Computer Methods for Sampling from Gamma, Beta, Poisson and
!    Binomial Distributions,
!    Computing,
!    Volume 12, Number 3, September 1974, pages 223-246.
!
!    Joachim Ahrens, Ulrich Dieter,
!    Generating Gamma Variates by a Modified Rejection Technique,
!    Communications of the ACM,
!    Volume 25, Number 1, January 1982, pages 47-54.
!
!  Parameters:
!
!    Input, real A, B, the parameters of the PDF.
!    0.0 < A,
!    0.0 < B.
!
!    Output, real X, a sample of the PDF.
!
  implicit none

  real, parameter :: a1 =  0.3333333E+00
  real, parameter :: a2 = -0.2500030E+00
  real, parameter :: a3 =  0.2000062E+00
  real, parameter :: a4 = -0.1662921E+00
  real, parameter :: a5 =  0.1423657E+00
  real, parameter :: a6 = -0.1367177E+00
  real, parameter :: a7 =  0.1233795E+00
  real, parameter :: e1 = 1.0E+00
  real, parameter :: e2 = 0.4999897E+00
  real, parameter :: e3 = 0.1668290E+00
  real, parameter :: e4 = 0.0407753E+00
  real, parameter :: e5 = 0.0102930E+00
  real, parameter :: q1 =  0.04166669E+00
  real, parameter :: q2 =  0.02083148E+00
  real, parameter :: q3 =  0.00801191E+00
  real, parameter :: q4 =  0.00144121E+00
  real, parameter :: q5 = -0.00007388E+00
  real, parameter :: q6 =  0.00024511E+00
  real, parameter :: q7 =  0.00024240E+00

!
!  Allow A = 0.
!
  if ( a == 0.0E+00 ) then
    x = 0.0E+00
    return
  end if
!
!  A < 1.
!
  if ( a < 1.0E+00 ) then

    do

      call r4_random ( 0.0E+00, 1.0E+00, p )
      p = ( 1.0E+00 + 0.3678794E+00 * a ) * p

      call exponential_01_sample ( e )

      if ( 1.0E+00 <= p ) then

        x = - log ( ( 1.0E+00 + 0.3678794E+00 * a - p ) / a )

        if ( ( 1.0E+00 - a ) * log ( x ) <= e ) then
          x = x / b
          return
        end if

      else

        x = exp ( log ( p ) / a )

        if ( x <= e ) then
          x = x / b
          return
        end if

      end if

    end do
!
!  1 <= A.
!
  else

    s2 = a - 0.5E+00
    s = sqrt ( a - 0.5E+00 )
    d = sqrt ( 32.0E+00 ) - 12.0E+00 * sqrt ( a - 0.5E+00 )

    call normal_01_sample ( t )
    x = ( sqrt ( a - 0.5E+00 ) + 0.5E+00 * t )**2

    if ( 0.0E+00 <= t ) then
      x = x / b
      return
    end if

    call r4_random ( 0.0E+00, 1.0E+00, u )

    if ( d * u <= t**3 ) then
      x = x / b
      return
    end if

    r = 1.0E+00 / a
    q0 = ( ( ( ( ( ( &
           q7   * r &
         + q6 ) * r &
         + q5 ) * r &
         + q4 ) * r &
         + q3 ) * r &
         + q2 ) * r &
         + q1 ) * r 

    if ( a <= 3.686E+00 ) then
      bcoef = 0.463E+00 + s - 0.178E+00 * s2
      si = 1.235E+00
      c = 0.195E+00 / s - 0.079E+00 + 0.016E+00 * s
    else if ( a <= 13.022E+00 ) then
      bcoef = 1.654E+00 + 0.0076E+00 * s2
      si = 1.68E+00 / s + 0.275E+00
      c = 0.062E+00 / s + 0.024E+00
    else
      bcoef = 1.77E+00
      si = 0.75E+00
      c = 0.1515E+00 / s
    end if

    if ( 0.0E+00 < sqrt ( a - 0.5E+00 ) + 0.5E+00 * t ) then

      v = 0.5E+00 * t / s

      if ( 0.25E+00 < abs ( v ) ) then
        q = q0 - s * t + 0.25E+00 * t**2 + 2.0E+00 * s2 * log ( 1.0E+00 + v )
      else
        q = q0 + 0.5E+00 * t**2 * ( ( ( ( ( ( &
               a7   * v &
             + a6 ) * v &
             + a5 ) * v &
             + a4 ) * v &
             + a3 ) * v &
             + a2 ) * v &
             + a1 ) * v
      end if

      if ( log ( 1.0E+00 - u ) <= q ) then
        x = x / b
        return
      end if

    end if

    do

      call exponential_01_sample ( e )

      call r4_random ( 0.0E+00, 1.0E+00, u )
      u = 2.0E+00 * u - 1.0E+00
      t = bcoef + sign ( si * e, u )

      if ( - 0.7187449E+00 <= t ) then

        v = 0.5E+00 * t / s

        if ( 0.25E+00 < abs ( v ) ) then
          q = q0 - s * t + 0.25E+00 * t**2 + 2.0E+00 * s2 * log ( 1.0E+00 + v )
        else
          q = q0 + 0.5E+00 * t**2 * ( ( ( ( ( ( &
                 a7   * v &
               + a6 ) * v &
               + a5 ) * v &
               + a4 ) * v &
               + a3 ) * v &
               + a2 ) * v &
               + a1 ) * v
        end if

        if ( 0.0E+00 < q ) then

          if ( 0.5E+00 < q ) then
            w = exp ( q ) - 1.0E+00
          else
            w = ( ( ( ( &
                   e5   * q &
                 + e4 ) * q &
                 + e3 ) * q &
                 + e2 ) * q &
                 + e1 ) * q
          end if

          if ( c * abs ( u ) <= w * exp ( e - 0.5E+00 * t**2 ) ) then
            x = ( s + 0.5E+00 * t )**2 / b
            return
          end if

        end if

      end if

    end do

  end if

end
function gammad ( x, p, ifault )


!*****************************************************************************80
!
!! PPCHI2 evaluates the percentage points of the chi-squared PDF.
!
!  Modified:
!
!    30 March 1999
!
!  Reference:
!
!    Donald Best, DE Roberts,
!    Algorithm AS 91:
!    The Percentage Points of the Chi-Squared Distribution,
!    Applied Statistics,
!    Volume 24, Number 3, 1975, pages 385-390.
!
!  Parameters:
!
!    Input, real P, a value of the chi-squared cumulative probability
!    density function.
!    0.000002 <= P <= 0.999998.
!
!    Input, real V, the parameter of the chi-squared probability density
!    function.  0 < V.
!
!    Output, integer IFAULT, error flag.
!    0, no error detected.
!    1, P < PMIN or PMAX < P.
!    2, V <= 0.0.
!    3, an error occurred in the incomplete Gamma function routine.
!    4, the maximum number of iterations were taken without convergence.
!    5, an error occurred in the log Gamma routine
!
  implicit none

  real, parameter :: aa = 0.6931471806E+00
  real, parameter :: c1 = 0.01E+00
  real, parameter :: c2 = 0.222222E+00
  real, parameter :: c3 = 0.32E+00
  real, parameter :: c4 = 0.4E+00
  real, parameter :: c5 = 1.24E+00
  real, parameter :: c6 = 2.2E+00
  real, parameter :: c7 = 4.67E+00
  real, parameter :: c8 = 6.66E+00
  real, parameter :: c9 = 6.73E+00
  real, parameter :: c10 = 13.32E+00
  real, parameter :: c11 = 60.0E+00
  real, parameter :: c12 = 70.0E+00
  real, parameter :: c13 = 84.0E+00
  real, parameter :: c14 = 105.0E+00
  real, parameter :: c15 = 120.0E+00
  real, parameter :: c16 = 127.0E+00
  real, parameter :: c17 = 140.0E+00
  real, parameter :: c18 = 175.0E+00
  real, parameter :: c19 = 210.0E+00
  real, parameter :: c20 = 252.0E+00
  real, parameter :: c21 = 264.0E+00
  real, parameter :: c22 = 294.0E+00
  real, parameter :: c23 = 346.0E+00
  real, parameter :: c24 = 420.0E+00
  real, parameter :: c25 = 462.0E+00
  real, parameter :: c26 = 606.0E+00
  real, parameter :: c27 = 672.0E+00
  real, parameter :: c28 = 707.0E+00
  real, parameter :: c29 = 735.0E+00
  real, parameter :: c30 = 889.0E+00
  real, parameter :: c31 = 932.0E+00
  real, parameter :: c32 = 966.0E+00
  real, parameter :: c33 = 1141.0E+00
  real, parameter :: c34 = 1182.0E+00
  real, parameter :: c35 = 1278.0E+00
  real, parameter :: c36 = 1740.0E+00
  real, parameter :: c37 = 2520.0E+00
  real, parameter :: c38 = 5040.0E+00
  real, parameter :: e = 0.0000005E+00
  integer, parameter :: maxit = 20
  real, parameter :: pmax = 0.999998E+00
  real, parameter :: pmin = 0.000002E+00

  ifault2 = 0
!
!  Check the input.
!
  if ( p < pmin ) then
    ifault = 1
    ppchi2 = - 1.0E+00
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PPCHI2 - Fatal error!'
    write ( *, '(a)' ) '  P < PMIN.'
    stop
  end if

  if ( pmax < p ) then
    ifault = 1
    ppchi2 = - 1.0E+00
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PPCHI2 - Fatal error!'
    write ( *, '(a)' ) '  PMAX < P.'
    stop
  end if

  if ( v <= 0.0E+00 ) then
    ifault = 2
    ppchi2 = - 1.0E+00
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PPCHI2 - Fatal error!'
    write ( *, '(a)' ) '  V <= 0.0.'
    stop
  end if

  ifault = 0
  xx = 0.5E+00 * v
  c = xx - 1.0E+00
!
!  Compute Log ( Gamma ( V/2 ) ).
!
  g = alngam ( v / 2.0E+00, ifault )

  if ( ifault /= 0 ) then
    ifault = 5
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PPCHI2 - Fatal error!'
    write ( *, '(a)' ) '  ALNGAM returns error code.'
    stop
  end if
!
!  Starting approximation for small chi-squared.
!
  if ( v < - c5 * log ( p ) ) then

    ch = ( p * xx * exp ( g + xx * aa ) )**( 1.0E+00 / xx )

    if ( ch < e ) then
      ifault = 0
      ppchi2 = ch
      return
    end if
!
!  Starting approximation for V less than or equal to 0.32.
!
  else if ( v <= c3 ) then

    ch = c4
    a = log ( 1.0E+00 - p )

    do

      q = ch
      p1 = 1.0E+00 + ch * ( c7 + ch )
      p2 = ch * ( c9 + ch * ( c8 + ch ) )

      t = - 0.5E+00 + ( c7 + 2.0E+00 * ch ) / p1 - ( c9 + ch * ( c10 + 3.0E+00 * ch ) ) / p2

      ch = ch - ( 1.0E+00 - exp ( a + g + 0.5E+00 * ch + c * aa ) * p2 / p1 ) / t

      if ( abs ( q / ch - 1.0E+00 ) <= c1 ) then
        exit
      end if

    end do
!
!  Call to algorithm AS 111.
!  Note that P has been tested above.
!  AS 241 could be used as an alternative.
!
  else

    x = ppnd ( p, ifault2 )
!
!  Starting approximation using Wilson and Hilferty estimate.
!
    p1 = c2 / v
    ch = v * ( x * sqrt ( p1 ) + 1.0E+00 - p1 )**3
!
!  Starting approximation for P tending to 1.
!
    if ( c6 * v + 6.0E+00 < ch ) then
      ch = - 2.0E+00 * ( log ( 1.0E+00 - p ) - c * log ( 0.5E+00 * ch ) + g )
    end if

  end if
!
!  Call to algorithm AS 239 and calculation of seven term Taylor series.
!
  do i = 1, maxit

    q = ch
    p1 = 0.5E+00 * ch
    p2 = p - gammad ( p1, xx, ifault2 )

    if ( ifault2 /= 0 ) then
      ppchi2 = - 1.0E+00
      ifault = 3
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'PPCHI2 - Fatal error!'
      write ( *, '(a)' ) '  GAMMAD returns error code.'
      stop
    end if

    t = p2 * exp ( xx * aa + g + p1 - c * log ( ch ) )
    b = t / ch
    a = 0.5E+00 * t - b * c

    s1 = &
           ( c19 + a &
         * ( c17 + a &
         * ( c14 + a &
         * ( c13 + a &
         * ( c12 + a &
         *   c11 ) ) ) ) ) / c24

    s2 = &
           ( c24 + a &
         * ( c29 + a &
         * ( c32 + a &
         * ( c33 + a &
         *   c35 ) ) ) ) / c37

    s3 = &
           ( c19 + a &
         * ( c25 + a &
         * ( c28 + a * &
             c31 ) ) ) / c37

    s4 = &
           ( c20 + a &
         * ( c27 + a &
         *   c34 ) + c &
         * ( c22 + a &
         * ( c30 + a &
         *   c36 ) ) ) / c38

    s5 = ( c13 + c21 * a + c * ( c18 + c26 * a ) ) / c37

    s6 = ( c15 + c * ( c23 + c16 * c ) ) / c38

    ch = ch + t * ( 1.0E+00 + 0.5E+00 * t * s1 - b * c &
         * ( s1 - b &
         * ( s2 - b &
         * ( s3 - b &
         * ( s4 - b &
         * ( s5 - b &
         *   s6 ) ) ) ) ) )

    if ( e < abs ( q / ch - 1.0E+00 ) ) then
      ifault = 0
      ppchi2 = ch
      return
    end if

  end do

  ifault = 4
  ppchi2 = ch
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'PPCHI2 - Warning!'
  write ( *, '(a)' ) '  Convergence not reached.'

  return
end
function ppnd ( p, ifault )

!*****************************************************************************80
!
!! PPND produces the normal deviate value corresponding to lower tail area = P.
!
!  Modified:
!
!    28 March 1999
!
!  Reference:
!
!    JD Beasley, SG Springer,
!    Algorithm AS 111:
!    The Percentage Points of the Normal Distribution,
!    Applied Statistics,
!    Volume 26, Number 1, 1977, pages 118-121.
!
!  Parameters:
!
!    Input, real P, the value of the cumulative probability densitity function.
!    0 < P < 1.
!
!    Output, integer IFAULT, error flag.
!    0, no error.
!    1, P <= 0 or 1 <= P.  PPND is returned as 0.
!
!    Output, real PPND, the normal deviate value with the property that
!    the probability of a standard normal deviate being less than or
!    equal to PPND is P.
!
  implicit none

  real, parameter :: a0 = 2.50662823884E+00
  real, parameter :: a1 = -18.61500062529E+00
  real, parameter :: a2 = 41.39119773534E+00
  real, parameter :: a3 = -25.44106049637E+00
  real, parameter :: b1 = -8.47351093090E+00
  real, parameter :: b2 = 23.08336743743E+00
  real, parameter :: b3 = -21.06224101826E+00
  real, parameter :: b4 = 3.13082909833E+00
  real, parameter :: c0 = -2.78718931138E+00
  real, parameter :: c1 = -2.29796479134E+00
  real, parameter :: c2 = 4.85014127135E+00
  real, parameter :: c3 = 2.32121276858E+00
  real, parameter :: d1 = 3.54388924762E+00
  real, parameter :: d2 = 1.63706781897E+00
  real, parameter :: split = 0.42E+00

  integer ifault
  real p
  real ppnd
  real r

  ifault = 0
!
!  0.08 < P < 0.92
!
  if ( abs ( p - 0.5E+00 ) <= split ) then

    r = ( p - 0.5E+00 )**2

    ppnd = ( p - 0.5E+00 ) * ( ( ( &
           a3   * r &
         + a2 ) * r &
         + a1 ) * r &
         + a0 ) / ( ( ( ( &
           b4   * r &
         + b3 ) * r &
         + b2 ) * r &
         + b1 ) * r &
         + 1.0E+00 )
!
!  P < 0.08 or 0.92 < P, 
!  R = min ( P, 1-P )
!
  else if ( 0.0E+00 < p .and. p < 1.0E+00 ) then

    if ( 0.5E+00 < p ) then
      r = sqrt ( - log ( 1.0E+00 - p ) )
    else
      r = sqrt ( - log ( p ) )
    end if

    ppnd = ( ( ( &
           c3   * r &
         + c2 ) * r &
         + c1 ) * r &
         + c0 ) / ( ( &
           d2   * r &
         + d1 ) * r &
         + 1.0E+00 )

    if ( p < 0.5E+00 ) then
      ppnd = - ppnd
    end if
!
!  P <= 0.0 or 1.0 <= P.
!
  else

    ifault = 1
    ppnd = 0.0E+00

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PPND - Warning!'
    write ( *, '(a)' ) '  P <= 0 or 1 <= P.'
    write ( *, '(a)' ) '  PPND value would be infinite.'
    
  end if

  return
end
function ppnd16 ( p, ifault )

!*****************************************************************************80
!
!! PPND16 produces the normal deviate value corresponding to lower tail area = P.
!
!  Discussion:
!
!    The result is accurate to about 1 part in 10**16.
!
!  Modified:
!
!    28 March 1999
!
!  Reference:
!
!    Michael Wichura,
!    Algorithm AS 241:
!    The Percentage Points of the Normal Distribution,
!    Applied Statistics,
!    Volume 37, Number 3, 1988, pages 477-484.
!
!  Parameters:
!
!    Input, double precision P, the value of the cumulative probability 
!    densitity function.  0 < P < 1.
!
!    Output, integer IFAULT, error flag.
!    0, no error.
!    1, P <= 0 or 1 <= P.
!
!    Output, double precision PPND16, the normal deviate value with the 
!    property that the probability of a standard normal deviate being 
!    less than or equal to PPND16 is P.
!
  implicit none

  double precision, parameter :: a0 = 3.3871328727963666080d+00
  double precision, parameter :: a1 = 1.3314166789178437745d+02
  double precision, parameter :: a2 = 1.9715909503065514427d+03
  double precision, parameter :: a3 = 1.3731693765509461125d+04
  double precision, parameter :: a4 = 4.5921953931549871457d+04
  double precision, parameter :: a5 = 6.7265770927008700853d+04
  double precision, parameter :: a6 = 3.3430575583588128105d+04
  double precision, parameter :: a7 = 2.5090809287301226727d+03
  double precision, parameter :: b1 = 4.2313330701600911252d+01
  double precision, parameter :: b2 = 6.8718700749205790830d+02
  double precision, parameter :: b3 = 5.3941960214247511077d+03
  double precision, parameter :: b4 = 2.1213794301586595867d+04
  double precision, parameter :: b5 = 3.9307895800092710610d+04
  double precision, parameter :: b6 = 2.8729085735721942674d+04
  double precision, parameter :: b7 = 5.2264952788528545610d+03
  double precision, parameter :: c0 = 1.42343711074968357734d+00
  double precision, parameter :: c1 = 4.63033784615654529590d+00
  double precision, parameter :: c2 = 5.76949722146069140550d+00
  double precision, parameter :: c3 = 3.64784832476320460504d+00
  double precision, parameter :: c4 = 1.27045825245236838258d+00
  double precision, parameter :: c5 = 2.41780725177450611770d-01
  double precision, parameter :: c6 = 2.27238449892691845833d-02
  double precision, parameter :: c7 = 7.74545014278341407640d-04
  double precision, parameter :: const1 = 0.180625d+00
  double precision, parameter :: const2 = 1.6d+00
  double precision, parameter :: d1 = 2.05319162663775882187d+00
  double precision, parameter :: d2 = 1.67638483018380384940d+00
  double precision, parameter :: d3 = 6.89767334985100004550d-01
  double precision, parameter :: d4 = 1.48103976427480074590d-01
  double precision, parameter :: d5 = 1.51986665636164571966d-02
  double precision, parameter :: d6 = 5.47593808499534494600d-04
  double precision, parameter :: d7 = 1.05075007164441684324d-09
  double precision, parameter :: e0 = 6.65790464350110377720d+00
  double precision, parameter :: e1 = 5.46378491116411436990d+00
  double precision, parameter :: e2 = 1.78482653991729133580d+00
  double precision, parameter :: e3 = 2.96560571828504891230d-01
  double precision, parameter :: e4 = 2.65321895265761230930d-02
  double precision, parameter :: e5 = 1.24266094738807843860d-03
  double precision, parameter :: e6 = 2.71155556874348757815d-05
  double precision, parameter :: e7 = 2.01033439929228813265d-07
  double precision, parameter :: f1 = 5.99832206555887937690d-01
  double precision, parameter :: f2 = 1.36929880922735805310d-01
  double precision, parameter :: f3 = 1.48753612908506148525d-02
  double precision, parameter :: f4 = 7.86869131145613259100d-04
  double precision, parameter :: f5 = 1.84631831751005468180d-05
  double precision, parameter :: f6 = 1.42151175831644588870d-07
  double precision, parameter :: f7 = 2.04426310338993978564d-15
  double precision, parameter :: split1 = 0.425d+00
  double precision, parameter :: split2 = 5.d+00

  integer ifault
  double precision p
  double precision ppnd16
  double precision q
  double precision r

  ifault = 0
  q = p - 0.5D+00

  if ( abs ( q ) <= split1 ) then

    r = const1 - q**2

    ppnd16 = q * ((((((( &
           a7   * r &
         + a6 ) * r &
         + a5 ) * r &
         + a4 ) * r &
         + a3 ) * r &
         + a2 ) * r &
         + a1 ) * r &
         + a0 ) / ((((((( &
           b7   * r &
         + b6 ) * r &
         + b5 ) * r &
         + b4 ) * r &
         + b3 ) * r &
         + b2 ) * r &
         + b1 ) * r &
         + 1.0D+00 )

  else

    if ( q < 0.0D+00 ) then
      r = p
    else
      r = 1.0D+00 - p
    end if

    if ( r <= 0.0D+00 ) then
      ifault = 1
      ppnd16 = 0.0D+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'PPND16 - Warning!'
      write ( *, '(a)' ) '  P <= 0 or 1 <= P.'
      write ( *, '(a)' ) '  PPND16 value would be infinite.'
      return
    end if

    r = sqrt ( - log ( r ) )

    if ( r <= split2 ) then

      r = r - const2

      ppnd16 = ((((((( &
             c7   * r &
           + c6 ) * r &
           + c5 ) * r &
           + c4 ) * r &
           + c3 ) * r &
           + c2 ) * r &
           + c1 ) * r &
           + c0 ) / ((((((( &
             d7   * r &
           + d6 ) * r &
           + d5 ) * r &
           + d4 ) * r &
           + d3 ) * r &
           + d2 ) * r &
           + d1 ) * r &
           + 1.0D+00 )

    else

      r = r - split2

      ppnd16 = ((((((( &
             e7   * r &
           + e6 ) * r &
           + e5 ) * r &
           + e4 ) * r &
           + e3 ) * r &
           + e2 ) * r &
           + e1 ) * r &
           + e0 ) / ((((((( &
             f7   * r &
           + f6 ) * r &
           + f5 ) * r &
           + f4 ) * r &
           + f3 ) * r &
           + f2 ) * r &
           + f1 ) * r &
           + 1.0D+00 )

    end if

    if ( q < 0.0D+00 ) then
      ppnd16 = - ppnd16
    end if

  end if

  return
end
function ppnd7 ( p, ifault )

!*****************************************************************************80
!
!! PPND7 produces the normal deviate value corresponding to lower tail area = P.
!
!  Discussion:
!
!    The result is accurate to about 1 part in 10**7.
!
!  Modified:
!
!    31 March 1999
!
!  Reference:
!
!    Michael Wichura,
!    Algorithm AS 241:
!    The Percentage Points of the Normal Distribution,
!    Applied Statistics,
!    Volume 37, Number 3, 1988, pages 477-484.
!
!  Parameters:
!
!    Input, real P, the value of the cumulative probability densitity function.
!    0 < P < 1.
!
!    Output, integer IFAULT, error flag.
!    0, no error.
!    1, P <= 0 or 1 <= P.
!
!    Output, real PPND7, the normal deviate value with the property that
!    the probability of a standard normal deviate being less than or
!    equal to PPND7 is P.
!
  implicit none

  real, parameter :: a0 = 3.3871327179E+00
  real, parameter :: a1 = 50.434271938E+00
  real, parameter :: a2 = 159.29113202E+00
  real, parameter :: a3 = 59.109374720E+00
  real, parameter :: b1 = 17.895169469E+00
  real, parameter :: b2 = 78.757757664E+00
  real, parameter :: b3 = 67.187563600E+00
  real, parameter :: c0 = 1.4234372777E+00
  real, parameter :: c1 = 2.7568153900E+00
  real, parameter :: c2 = 1.3067284816E+00
  real, parameter :: c3 = 0.17023821103E+00
  real, parameter :: const1 = 0.180625E+00
  real, parameter :: const2 = 1.6E+00
  real, parameter :: d1 = 0.73700164250E+00
  real, parameter :: d2 = 0.12021132975E+00
  real, parameter :: e0 = 6.6579051150E+00
  real, parameter :: e1 = 3.0812263860E+00
  real, parameter :: e2 = 0.42868294337E+00
  real, parameter :: e3 = 0.017337203997E+00
  real, parameter :: f1 = 0.24197894225E+00
  real, parameter :: f2 = 0.012258202635E+00
  real, parameter :: split1 = 0.425E+00
  real, parameter :: split2 = 5.0E+00

  integer ifault
  real p
  real ppnd7
  real q
  real r

  ifault = 0
  q = p - 0.5E+00

  if ( abs ( q ) <= split1 ) then

    r = const1 - q**2

    ppnd7 = q * ((( &
           a3   * r &
         + a2 ) * r &
         + a1 ) * r &
         + a0 ) / ((( &
           b3   * r &
         + b2 ) * r &
         + b1 ) * r &
         + 1.0E+00 )

  else

    if ( q < 0.0E+00 ) then
      r = p
    else
      r = 1.0E+00 - p
    end if

    if ( r <= 0.0E+00 ) then
      ifault = 1
      ppnd7 = 0.0E+00
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'PPND7 - Fatal error!'
      write ( *, '(a)' ) '  P <= 0 or 1 <= P.'
      write ( *, '(a)' ) '  PPND7 value would be infinite.'
      return
    end if

    r = sqrt ( - log ( r ) )

    if ( r <= split2 ) then

      r = r - const2

      ppnd7 = ((( &
           c3   * r &
         + c2 ) * r &
         + c1 ) * r &
         + c0 ) / (( &
           d2   * r &
         + d1 ) * r &
         + 1.0E+00 )

    else

      r = r - split2

      ppnd7 = ((( &
           e3   * r &
         + e2 ) * r &
         + e1 ) * r &
         + e0 ) / (( &
           f2   * r &
         + f1 ) * r &
         + 1.0E+00 )

    end if

    if ( q < 0.0E+00 ) then
      ppnd7 = - ppnd7
    end if

  end if

  return
end
function psi ( xx )

!*****************************************************************************80
!
!! PSI evaluates the psi or digamma function, d/dx ln(gamma(x)).
!
!  Discussion:
!
!    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:
!
!    William Cody, Anthony Strecok, Henry Thacher,
!    Chebyshev Approximations for the Psi Function,
!    Mathematics of Computation,
!    Volume 27, Number 121, January 1973, pages 123-127.
!
!  Parameters:
!
!    Input, real XX, the argument of the psi function.
!
!    Output, real PSI, the value of the psi function.  PSI is assigned
!    the value 0 when the psi function is undefined.
!
  implicit none

  real, parameter :: dx0 = 1.461632144968362341262659542325721325E+00
  real, parameter :: piov4 = 0.785398163397448E+00

  real aug
  real den
  real eps
  integer i
  integer m
  integer n
  integer nq
  real, parameter, dimension ( 7 ) :: p1 = (/&
    0.895385022981970E-02, &
    0.477762828042627E+01, &
    0.142441585084029E+03, &
    0.118645200713425E+04, &
    0.363351846806499E+04, &
    0.413810161269013E+04, &
    0.130560269827897E+04 /)
  real, parameter, dimension ( 4 ) :: p2 = (/ &
    -0.212940445131011E+01, &
    -0.701677227766759E+01, &
    -0.448616543918019E+01, &
    -0.648157123766197E+00 /)
  real psi
  real, parameter, dimension ( 6 ) :: q1 = (/&
    0.448452573429826E+02, &
    0.520752771467162E+03, &
    0.221000799247830E+04, &
    0.364127349079381E+04, &
    0.190831076596300E+04, &
    0.691091682714533E-05 /)
  real, parameter, dimension ( 4 ) :: q2 = (/ &
    0.322703493791143E+02, &
    0.892920700481861E+02, &
    0.546117738103215E+02, &
    0.777788548522962E+01 /)
  real sgn
  real temp
  real upper
  real w
  real x
  real xmax1
  real xmx0
  real xsmall
  real xx
  real z
!
!  XMAX1 is the largest positive floating point constant with entirely
!  integer representation.  It is 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).
!
  xmax1 = dble ( 2147483647 )

  eps = epsilon ( eps )

  xmax1 = min ( xmax1, 1.0E+00 / eps )
!
!  XSMALL is the absolute argument below which PI*COTAN(PI*X)
!  may be represented by 1/X.
!
  xsmall = 1.0E-09

  x = xx
  aug = 0.0E+00

  if ( x == 0.0E+00 ) then
    psi = 0.0E+00
    return
  end if
!
!  X < 0.5,  Use reflection formula PSI(1-X) = PSI(X) + PI * COTAN(PI*X)
!
  if ( x < 0.5E+00 ) then
!
!  0 < ABS(X) <= XSMALL.  Use 1/X as a substitute for PI*COTAN(PI*X)
!
    if ( abs ( x ) <= xsmall ) then
      aug = - 1.0E+00 / x
      go to 40
    end if
!
!  Reduction of argument for cotan
!
    w = - x
    sgn = piov4

    if ( w <= 0.0E+00 ) then
      w = - w
      sgn = -sgn
    end if
!
!  Make an error exit if X <= -XMAX1
!
    if ( xmax1 <= w ) then
      psi = 0.0E+00
      return
    end if

    nq = int ( w )
    w = w - real ( nq )
    nq = int ( w * 4.0E+00 )
    w = 4.0E+00 * ( w - real ( nq ) * 0.25E+00 )
!
!  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.0E+00 - 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 ) then

      if ( z == 0.0E+00 ) then
        psi = 0.0E+00
       return
    end  if

      aug = 4.0E+00 * sgn * ( cos(z) / sin(z) )

    else

      aug = 4.0E+00 * sgn * ( sin(z) / cos(z) )

    end if

   40   continue

    x = 1.0E+00 - x

  end if
!
!  0.5 <= X <= 3.0
!
  if ( x <= 3.0E+00 ) then

    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
!
!  3.0 < X < XMAX1
!
  else if ( x < xmax1 ) then

    w = 1.0E+00 / x**2
    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.5E+00 / x + aug
    psi = aug + log ( x )
!
!  XMAX1 <= X
!
  else

    psi = aug + log ( x )

  end if

  return
end
subroutine r4_random ( rlo, rhi, r )


function algdiv ( a, b )

!*****************************************************************************80
!
!! ALGDIV computes ln ( Gamma ( B ) / Gamma ( A + B ) ) when 8 <= B.
!
!  Discussion:
!
!    In this algorithm, DEL(X) is the function defined by
!
!      ln ( Gamma(X) ) = ( X - 0.5 ) * ln ( X ) - X + 0.5 * ln ( 2 * PI ) 
!                      + DEL(X).
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  Reference:
!
!    Armido DiDinato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, 1993, pages 360-373.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) A, B, define the arguments.
!
!    Output, real ( kind = 8 ) ALGDIV, the value of ln(Gamma(B)/Gamma(A+B)).
!
  implicit none

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

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

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

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

  return
end
function alnrel ( a )

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

function bcorr ( a0, b0 )

!*****************************************************************************80
!
!! BCORR evaluates DEL(A0) + DEL(B0) - DEL(A0 + B0).
!
!  Discussion:
!
!    The function DEL(A) is a remainder term that is used in the expression:
!
!      ln ( Gamma ( A ) ) = ( A - 0.5 ) * ln ( A ) 
!        - A + 0.5 * ln ( 2 * PI ) + DEL ( A ),
!
!    or, in other words, DEL ( A ) is defined as:
!
!      DEL ( A ) = ln ( Gamma ( A ) ) - ( A - 0.5 ) * ln ( A ) 
!        + A + 0.5 * ln ( 2 * PI ).
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  Reference:
!
!    Armido DiDinato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, 1993, pages 360-373.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) A0, B0, the arguments.
!    It is assumed that 8 <= A0 and 8 <= B0.
!
!    Output, real ( kind = 8 ) BCORR, the value of the function.
!
  implicit none

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

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

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

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

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

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

  return
end
function beta ( a, b )


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

  integer, parameter :: n_max = 30

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

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

  n_data = n_data + 1

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

  return
end
function beta_log ( a0, b0 )


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

  integer, parameter :: n_max = 17

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

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

  n_data = n_data + 1

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

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


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

  integer, parameter :: n_max = 27

  real ( kind = 8 ) cdf
  real, save, dimension ( n_max ) :: cdf_vec = (/ &
    0.839944D+00, 0.695906D+00, 0.535088D+00, &
    0.764784D+00, 0.620644D+00, 0.469167D+00, &
    0.307088D+00, 0.220382D+00, 0.150025D+00, &
    0.307116D-02, 0.176398D-02, 0.981679D-03, &
    0.165175D-01, 0.202342D-03, 0.498448D-06, &
    0.151325D-01, 0.209041D-02, 0.246502D-03, &
    0.263684D-01, 0.185798D-01, 0.130574D-01, &
    0.583804D-01, 0.424978D-01, 0.308214D-01, &
    0.105788D+00, 0.794084D-01, 0.593201D-01 /)
  integer df
  integer, save, dimension ( n_max ) :: df_vec = (/ &
      1,   2,   3, &
      1,   2,   3, &
      1,   2,   3, &
      1,   2,   3, &
     60,  80, 100, &
      1,   2,   3, &
     10,  10,  10, &
     10,  10,  10, &
     10,  10,  10 /)
  real ( kind = 8 ) lambda
  real, save, dimension ( n_max ) :: lambda_vec = (/ &
     0.5D+00,  0.5D+00,  0.5D+00, &
     1.0D+00,  1.0D+00,  1.0D+00, &
     5.0D+00,  5.0D+00,  5.0D+00, &
    20.0D+00, 20.0D+00, 20.0D+00, &
    30.0D+00, 30.0D+00, 30.0D+00, &
     5.0D+00,  5.0D+00,  5.0D+00, &
     2.0D+00,  3.0D+00,  4.0D+00, &
     2.0D+00,  3.0D+00,  4.0D+00, &
     2.0D+00,  3.0D+00,  4.0D+00 /)
  integer n_data
  real ( kind = 8 ) x
  real, save, dimension ( n_max ) :: x_vec = (/ &
     3.000D+00,  3.000D+00,  3.000D+00, &
     3.000D+00,  3.000D+00,  3.000D+00, &
     3.000D+00,  3.000D+00,  3.000D+00, &
     3.000D+00,  3.000D+00,  3.000D+00, &
    60.000D+00, 60.000D+00, 60.000D+00, &
     0.050D+00,  0.050D+00,  0.050D+00, &
     4.000D+00,  4.000D+00,  4.000D+00, &
     5.000D+00,  5.000D+00,  5.000D+00, &
     6.000D+00,  6.000D+00,  6.000D+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
    x = 0.0D+00
    lambda = 0.0D+00
    df = 0
    cdf = 0.0D+00
  else
    x = x_vec(n_data)
    lambda = lambda_vec(n_data)
    df = df_vec(n_data)
    cdf = cdf_vec(n_data)
  end if

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

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

  integer, parameter :: n_max = 21

  integer a
  integer, save, dimension ( n_max ) :: a_vec = (/ &
     1,  2,  1,  2, &
     1,  2,  3,  4, &
     1,  2,  3,  4, &
     5,  3,  3,  3, &
     3,  3, 10, 10, &
    10 /)
  real ( kind = 8 ) fx
  real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ &
    0.0796557D+00, 0.00498752D+00, 0.112463D+00,    0.00995017D+00, &
    0.472911D+00,  0.181269D+00,   0.0597575D+00,   0.0175231D+00, &
    0.682689D+00,  0.393469D+00,   0.198748D+00,    0.090204D+00, &
    0.0374342D+00, 0.427593D+00,   0.608375D+00,    0.738536D+00, &
    0.828203D+00,  0.88839D+00,    0.000172116D+00, 0.00365985D+00, &
    0.0185759D+00 /)
  integer n_data
  real ( kind = 8 ) x
  real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ &
    0.01D+00, 0.01D+00, 0.02D+00, 0.02D+00, &
    0.40D+00, 0.40D+00, 0.40D+00, 0.40D+00, &
    1.00D+00, 1.00D+00, 1.00D+00, 1.00D+00, &
    1.00D+00, 2.00D+00, 3.00D+00, 4.00D+00, &
    5.00D+00, 6.00D+00, 1.00D+00, 2.00D+00, &
    3.00D+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
    x = 0.0D+00
    fx = 0.0D+00
  else
    a = a_vec(n_data)
    x = x_vec(n_data)
    fx = fx_vec(n_data)
  end if

  return
end
subroutine cumbet ( x, y, a, b, cum, ccum )

subroutine cumnor ( arg, cum, ccum )

!*****************************************************************************80
!
!! CUMNOR computes the cumulative normal distribution.
!
!  Discussion:
!
!    This function evaluates the normal distribution function:
!
!                              / x
!                     1       |       -t*t/2
!          P(x) = ----------- |      e       dt
!                 sqrt(2 pi)  |
!                             /-oo
!
!    This transportable program uses rational functions that
!    theoretically approximate the normal distribution function to
!    at least 18 significant decimal digits.  The accuracy achieved
!    depends on the arithmetic system, the compiler, the intrinsic
!    functions, and proper selection of the machine dependent
!    constants.
!
!  Author: 
!
!    William Cody
!    Mathematics and Computer Science Division
!    Argonne National Laboratory
!    Argonne, IL 60439
!
!  Reference:
!
!    William Cody,
!    Rational Chebyshev approximations for the error function,
!    Mathematics of Computation, 
!    1969, pages 631-637.
!
!    William Cody, 
!    Algorithm 715: 
!    SPECFUN - A Portable FORTRAN Package of Special Function Routines 
!    and Test Drivers,
!    ACM Transactions on Mathematical Software,
!    Volume 19, Number 1, 1993, pages 22-32.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) ARG, the upper limit of integration.
!
!    Output, real ( kind = 8 ) CUM, CCUM, the Normal density CDF and
!    complementary CDF.
!
!  Local Parameters:
!
!    Local, real ( kind = 8 ) EPS, the argument below which anorm(x) 
!    may be represented by 0.5 and above which  x*x  will not underflow.
!    A conservative value is the largest machine number X
!    such that   1.0D+00 + X = 1.0D+00   to machine precision.
!
  implicit none

  real ( kind = 8 ), parameter, dimension ( 5 ) :: a = (/ &
    2.2352520354606839287D+00, &
    1.6102823106855587881D+02, &
    1.0676894854603709582D+03, &
    1.8154981253343561249D+04, &
    6.5682337918207449113D-02 /)
  real ( kind = 8 ) arg
  real ( kind = 8 ), parameter, dimension ( 4 ) :: b = (/ &
    4.7202581904688241870D+01, &
    9.7609855173777669322D+02, &
    1.0260932208618978205D+04, &
    4.5507789335026729956D+04 /)
  real ( kind = 8 ), parameter, dimension ( 9 ) :: c = (/ &
    3.9894151208813466764D-01, &
    8.8831497943883759412D+00, &
    9.3506656132177855979D+01, &
    5.9727027639480026226D+02, &
    2.4945375852903726711D+03, &
    6.8481904505362823326D+03, &
    1.1602651437647350124D+04, &
    9.8427148383839780218D+03, &
    1.0765576773720192317D-08 /)
  real ( kind = 8 ) ccum
  real ( kind = 8 ) cum
  real ( kind = 8 ), parameter, dimension ( 8 ) :: d = (/ &
    2.2266688044328115691D+01, &
    2.3538790178262499861D+02, &
    1.5193775994075548050D+03, &
    6.4855582982667607550D+03, &
    1.8615571640885098091D+04, &
    3.4900952721145977266D+04, &
    3.8912003286093271411D+04, &
    1.9685429676859990727D+04 /)
  real ( kind = 8 ) del
  real ( kind = 8 ) eps
  integer i
  real ( kind = 8 ), parameter, dimension ( 6 ) :: p = (/ &
    2.1589853405795699D-01, &
    1.274011611602473639D-01, &
    2.2235277870649807D-02, &
    1.421619193227893466D-03, &
    2.9112874951168792D-05, &
    2.307344176494017303D-02 /)
  real ( kind = 8 ), parameter, dimension ( 5 ) :: q = (/ &
    1.28426009614491121D+00, &
    4.68238212480865118D-01, &
    6.59881378689285515D-02, &
    3.78239633202758244D-03, &
    7.29751555083966205D-05 /)
  real ( kind = 8 ), parameter :: root32 = 5.656854248D+00
  real ( kind = 8 ), parameter :: sixten = 16.0D+00
  real ( kind = 8 ) temp
  real ( kind = 8 ), parameter :: sqrpi = 3.9894228040143267794D-01
  real ( kind = 8 ), parameter :: thrsh = 0.66291D+00
  real ( kind = 8 ) x
  real ( kind = 8 ) xden
  real ( kind = 8 ) xnum
  real ( kind = 8 ) y
  real ( kind = 8 ) xsq
!
!  Machine dependent constants
!
  eps = epsilon ( 1.0D+00 ) * 0.5D+00

  x = arg
  y = abs ( x )

  if ( y <= thrsh ) then
!
!  Evaluate  anorm  for  |X| <= 0.66291
!
    if ( eps < y ) then
      xsq = x * x
    else
      xsq = 0.0D+00
    end if

    xnum = a(5) * xsq
    xden = xsq
    do i = 1, 3
      xnum = ( xnum + a(i) ) * xsq
      xden = ( xden + b(i) ) * xsq
    end do
    cum = x * ( xnum + a(4) ) / ( xden + b(4) )
    temp = cum
    cum = 0.5D+00 + temp
    ccum = 0.5D+00 - temp
!
!  Evaluate ANORM for 0.66291 <= |X| <= sqrt(32)
!
  else if ( y <= root32 ) then

    xnum = c(9) * y
    xden = y
    do i = 1, 7
      xnum = ( xnum + c(i) ) * y
      xden = ( xden + d(i) ) * y
    end do
    cum = ( xnum + c(8) ) / ( xden + d(8) )
    xsq = aint ( y * sixten ) / sixten
    del = ( y - xsq ) * ( y + xsq )
    cum = exp ( - xsq * xsq * 0.5D+00 ) * exp ( -del * 0.5D+00 ) * cum
    ccum = 1.0D+00 - cum

    if ( 0.0D+00 < x ) then
      call r8_swap ( cum, ccum )
    end if
!
!  Evaluate ANORM for sqrt(32) < |X|.
!
  else

    cum = 0.0D+00
    xsq = 1.0D+00 / ( x * x )
    xnum = p(6) * xsq
    xden = xsq
    do i = 1, 4
      xnum = ( xnum + p(i) ) * xsq
      xden = ( xden + q(i) ) * xsq
    end do

    cum = xsq * ( xnum + p(5) ) / ( xden + q(5) )
    cum = ( sqrpi - cum ) / y
    xsq = aint ( x * sixten ) / sixten
    del = ( x - xsq ) * ( x + xsq )
    cum = exp ( - xsq * xsq * 0.5D+00 ) &
      * exp ( - del * 0.5D+00 ) * cum
    ccum = 1.0D+00 - cum

    if ( 0.0D+00 < x ) then
      call r8_swap ( cum, ccum )
    end if

  end if

  if ( cum < tiny ( cum ) ) then
    cum = 0.0D+00
  end if

  if ( ccum < tiny ( ccum ) ) then
    ccum = 0.0D+00
  end if

  return
end
subroutine cumpoi ( s, xlam, cum, ccum )

function dexpm1 ( x )

!*****************************************************************************80
!
!! DEXPM1 evaluates the function EXP(X) - 1.
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  Reference:
!
!    Armido DiDinato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software, 
!    Volume 18, 1993, pages 360-373.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, the value at which exp(X)-1 is desired.
!
!    Output, real ( kind = 8 ) DEXPM1, the value of exp(X)-1.
!
  implicit none

  real ( kind = 8 ) bot
  real ( kind = 8 ) dexpm1
  real ( kind = 8 ), parameter :: p1 =  0.914041914819518D-09
  real ( kind = 8 ), parameter :: p2 =  0.238082361044469D-01
  real ( kind = 8 ), parameter :: q1 = -0.499999999085958D+00
  real ( kind = 8 ), parameter :: q2 =  0.107141568980644D+00
  real ( kind = 8 ), parameter :: q3 = -0.119041179760821D-01
  real ( kind = 8 ), parameter :: q4 =  0.595130811860248D-03
  real ( kind = 8 ) top
  real ( kind = 8 ) w
  real ( kind = 8 ) x

  if ( abs ( x ) <= 0.15D+00 ) then

    top = ( p2 * x + p1 ) * x + 1.0D+00
    bot = ((( q4 * x + q3 ) * x + q2 ) * x + q1 ) * x + 1.0D+00
    dexpm1 = x * ( top / bot )

  else

    w = exp ( x )

    if ( x <= 0.0D+00 ) then
      dexpm1 = ( w - 0.5D+00 ) - 0.5D+00
    else
      dexpm1 = w * ( 0.5D+00 &
        + ( 0.5D+00 - 1.0D+00 / w ))
    end if

  end if

  return
end
function dinvnr ( p, q )

function dlanor ( x )

!*****************************************************************************80
!
!! DLANOR evaluates the logarithm of the asymptotic Normal CDF.
!
!  Discussion:
!
!    This routine computes the logarithm of the cumulative normal distribution
!    from abs ( x ) to infinity for  5 <= abs ( X ).
!
!    The relative error at X = 5 is about 0.5D-5.
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,  
!    Handbook of Mathematical Functions 
!    1966, Formula 26.2.12.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, the value at which the Normal CDF is to be
!    evaluated.  It is assumed that 5 <= abs ( X ).
!
!    Output, real ( kind = 8 ) DLANOR, the logarithm of the asymptotic
!    Normal CDF.
!
  implicit none

  real ( kind = 8 ) alnrel
  real ( kind = 8 ) approx
  real ( kind = 8 ), save, dimension ( 0:11 ) :: coef = (/ &
    -1.0D+00,  3.0D+00,  -15.0D+00,  105.0D+00,  -945.0D+00,  &
    10395.0D+00, -135135.0D+00,  2027025.0D+00,  -34459425.0D+00, &
    654729075.0D+00, -13749310575D+00,  316234143225.0D+00 /)
  real ( kind = 8 ) correc
  real ( kind = 8 ), parameter :: dlsqpi = 0.91893853320467274177D+00
  real ( kind = 8 ) eval_pol
  real ( kind = 8 ) dlanor
  real ( kind = 8 ) x
  real ( kind = 8 ) xx
  real ( kind = 8 ) xx2

  xx = abs ( x )

  if ( abs ( x ) < 5.0D+00 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DLANOR - Fatal error!'
    write ( *, '(a)' ) '  The argument X is too small.'
  end if

  approx = - dlsqpi - 0.5D+00 * x**2 - log ( abs ( x ) )

  xx2 = xx * xx
  correc = eval_pol ( coef, 11, 1.0D+00 / xx2 ) / xx2
  correc = alnrel ( correc )

  dlanor = approx + correc

  return
end
function dstrem ( z )

!*****************************************************************************80
!
!! DSTREM computes the Sterling remainder ln ( Gamma ( Z ) ) - Sterling ( Z ).
!
!  Discussion:
!
!    This routine returns 
!
!      ln ( Gamma ( Z ) ) - Sterling ( Z )  
!
!    where Sterling(Z) is Sterling's approximation to ln ( Gamma ( Z ) ).
!
!    Sterling(Z) = ln ( sqrt ( 2 * PI ) ) + ( Z - 0.5 ) * ln ( Z ) - Z
!
!    If 6 <= Z, the routine uses 9 terms of a series in Bernoulli numbers,
!    with values calculated using Maple.
!
!    Otherwise, the difference is computed explicitly.
!
!  Modified:
!
!    14 June 2004
!
!  Parameters:
!
!    Input, real ( kind = 8 ) Z, the value at which the Sterling 
!    remainder is to be calculated.  Z must be positive.
!
!    Output, real ( kind = 8 ) DSTREM, the Sterling remainder.
!
  implicit none

  integer, parameter :: ncoef = 9

  real ( kind = 8 ), parameter, dimension ( 0:ncoef ) :: coef = (/ &
    0.0D+00, &
    0.0833333333333333333333333333333D+00, &
    -0.00277777777777777777777777777778D+00, &
    0.000793650793650793650793650793651D+00, &
    -0.000595238095238095238095238095238D+00, &
    0.000841750841750841750841750841751D+00, &
    -0.00191752691752691752691752691753D+00, &
    0.00641025641025641025641025641026D+00, &
    -0.0295506535947712418300653594771D+00, &
    0.179644372368830573164938490016D+00 /)
  real ( kind = 8 ) dstrem
  real ( kind = 8 ) eval_pol
  real ( kind = 8 ) gamma_log
  real ( kind = 8 ), parameter :: hln2pi = 0.91893853320467274178D+00
  real ( kind = 8 ) sterl
  real ( kind = 8 ) z

  if ( z <= 0.0D+00 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'DSTREM - Fatal error!'
    write ( *, '(a)' ) '  Zero or negative argument Z.'
    stop
  end if

  if ( 6.0D+00 < z ) then
    dstrem = eval_pol ( coef, ncoef, 1.0D+00 / z**2 ) * z
  else
    sterl = hln2pi + ( z - 0.5D+00 ) * log ( z ) - z
    dstrem = gamma_log ( z ) - sterl
  end if

  return
end
function dt1 ( p, q, df )

!*****************************************************************************80
!
!! DT1 computes an approximate inverse of the cumulative T distribution.
!
!  Discussion:
!
!    This routine returns the inverse of the T distribution function, that is,
!    the integral from 0 to INVT of the T density is P.  This is an
!    initial approximation.
!
!    Thanks to Charles Katholi for pointing out that the RESHAPE
!    function should not use a range in the "SHAPE" field (0:4,4), 
!    but simply the number of rows and columns (5,4), JVB, 04 May 2006.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) P, Q, the value whose inverse from the 
!    T distribution CDF is desired, and the value (1-P).
!
!    Input, real ( kind = 8 ) DF, the number of degrees of freedom of the 
!    T distribution.
!
!    Output, real ( kind = 8 ) DT1, the approximate value of X for which
!    the T density CDF with DF degrees of freedom has value P.
!
  implicit none

  real ( kind = 8 ), dimension(0:4,4) :: coef = reshape ( (/ &
       1.0D+00,     1.0D+00,    0.0D+00,   0.0D+00,  0.0D+00, &
       3.0D+00,    16.0D+00,    5.0D+00,   0.0D+00,  0.0D+00, &
     -15.0D+00,    17.0D+00,   19.0D+00,   3.0D+00,  0.0D+00, &
    -945.0D+00, -1920.0D+00, 1482.0D+00, 776.0D+00, 79.0D+00/), (/ 5, 4 /) )
  real ( kind = 8 ), parameter, dimension ( 4 ) :: denom = (/ &
    4.0D+00, 96.0D+00, 384.0D+00, 92160.0D+00 /)
  real ( kind = 8 ) denpow
  real ( kind = 8 ) eval_pol
  real ( kind = 8 ) df
  real ( kind = 8 ) dinvnr
  real ( kind = 8 ) dt1
  integer i
  integer, parameter, dimension ( 4 ) :: ideg = (/ 1, 2, 3, 4 /)
  real ( kind = 8 ) p
  real ( kind = 8 ) q
  real ( kind = 8 ) sum1
  real ( kind = 8 ) term
  real ( kind = 8 ) x
  real ( kind = 8 ) xp
  real ( kind = 8 ) xx

  x = abs ( dinvnr ( p, q ) )
  xx = x * x

  sum1 = x
  denpow = 1.0D+00
  do i = 1, 4
    term = eval_pol ( coef(0,i), ideg(i), xx ) * x
    denpow = denpow * df
    sum1 = sum1 + term / ( denpow * denom(i) )
  end do

  if ( 0.5D+00 <= p ) then
    xp = sum1
  else
    xp = -sum1
  end if

  dt1 = xp

  return
end
subroutine dzror ( status, x, fx, xlo, xhi, qleft, qhi )


!*****************************************************************************80
!
!! ERF_VALUES returns some values of the ERF or "error" function.
!
!  Definition:
!
!    ERF(X) = ( 2 / sqrt ( PI ) * integral ( 0 <= T <= X ) exp ( - T^2 ) dT
!
!  Modified:
!
!    17 April 2001
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964.
!
!  Parameters:
!
!    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
!    first call.  On each call, the routine increments N_DATA by 1, and
!    returns the corresponding data; when there is no more data, the
!    output value of N_DATA will be 0 again.
!
!    Output, real ( kind = 8 ) X, the argument of the function.
!
!    Output, real ( kind = 8 ) FX, the value of the function.
!
  implicit none

  integer, parameter :: n_max = 21

  real ( kind = 8 ) fx
  real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ &
    0.0000000000D+00, 0.1124629160D+00, 0.2227025892D+00, 0.3286267595D+00, &
    0.4283923550D+00, 0.5204998778D+00, 0.6038560908D+00, 0.6778011938D+00, &
    0.7421009647D+00, 0.7969082124D+00, 0.8427007929D+00, 0.8802050696D+00, &
    0.9103139782D+00, 0.9340079449D+00, 0.9522851198D+00, 0.9661051465D+00, &
    0.9763483833D+00, 0.9837904586D+00, 0.9890905016D+00, 0.9927904292D+00, &
    0.9953222650D+00 /)
  integer n_data
  real ( kind = 8 ) x
  real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ &
    0.0D+00, 0.1D+00, 0.2D+00, 0.3D+00, &
    0.4D+00, 0.5D+00, 0.6D+00, 0.7D+00, &
    0.8D+00, 0.9D+00, 1.0D+00, 1.1D+00, &
    1.2D+00, 1.3D+00, 1.4D+00, 1.5D+00, &
    1.6D+00, 1.7D+00, 1.8D+00, 1.9D+00, &
    2.0D+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
    x = 0.0D+00
    fx = 0.0D+00
  else
    x = x_vec(n_data)
    fx = fx_vec(n_data)
  end if

  return
end
function error_f ( x )

!*****************************************************************************80
!
!! ERROR_F evaluates the error function.
!
!  Discussion:
!
!    Since some compilers already supply a routine named ERF which evaluates
!    the error function, this routine has been given a distinct, if
!    somewhat unnatural, name.
!    
!    The function is defined by:
!
!      ERF(X) = ( 2 / sqrt ( PI ) ) * Integral ( 0 <= T <= X ) EXP ( - T**2 ) dT.
!
!    Properties of the function include:
!
!      Limit ( X -> -Infinity ) ERF(X) =          -1.0;
!                               ERF(0) =           0.0;
!                               ERF(0.476936...) = 0.5;
!      Limit ( X -> +Infinity ) ERF(X) =          +1.0.
!
!      0.5D+00 * ( ERF(X/sqrt(2)) + 1 ) = Normal_01_CDF(X)
!
!  Modified:
!
!    17 November 2006
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  Reference:
!
!    Armido DiDinato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, 1993, pages 360-373.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, the argument.
!
!    Output, real ( kind = 8 ) ERF, the value of the error function at X.
!
  implicit none

  real ( kind = 8 ), parameter, dimension ( 5 ) :: a = (/ &
    0.771058495001320D-04, &
    -0.133733772997339D-02, &
    0.323076579225834D-01, &
    0.479137145607681D-01, &
    0.128379167095513D+00 /)
  real ( kind = 8 ) ax
  real ( kind = 8 ), parameter, dimension ( 3 ) :: b = (/ &
    0.301048631703895D-02, &
    0.538971687740286D-01, &
    0.375795757275549D+00 /)
  real ( kind = 8 ) bot
  real ( kind = 8 ), parameter :: c = 0.564189583547756D+00
  real ( kind = 8 ) error_f
  real ( kind = 8 ), dimension ( 8 ) :: p = (/   &
   -1.36864857382717D-07, 5.64195517478974D-01, &
    7.21175825088309D+00, 4.31622272220567D+01, &
    1.52989285046940D+02, 3.39320816734344D+02, &
    4.51918953711873D+02, 3.00459261020162D+02 /)
  real ( kind = 8 ), dimension ( 8 ) :: q = (/ &
    1.00000000000000D+00, 1.27827273196294D+01, &
    7.70001529352295D+01, 2.77585444743988D+02, &
    6.38980264465631D+02, 9.31354094850610D+02, & 
    7.90950925327898D+02, 3.00459260956983D+02 /)
  real ( kind = 8 ), dimension ( 5 ) :: r = (/ &
    2.10144126479064D+00, 2.62370141675169D+01, &
    2.13688200555087D+01, 4.65807828718470D+00, &
    2.82094791773523D-01 /)
  real ( kind = 8 ), parameter, dimension ( 4 ) :: s = (/ &
    9.41537750555460D+01, 1.87114811799590D+02, &
    9.90191814623914D+01, 1.80124575948747D+02 /)
  real ( kind = 8 ) t
  real ( kind = 8 ) top
  real ( kind = 8 ) x
  real ( kind = 8 ) x2

  ax = abs ( x )

  if ( ax <= 0.5D+00 ) then

    t = x * x

    top = (((( a(1)   * t &
             + a(2) ) * t &
             + a(3) ) * t &
             + a(4) ) * t &
             + a(5) ) + 1.0D+00

    bot = (( b(1) * t + b(2) ) * t + b(3) ) * t + 1.0D+00
    error_f = ax * ( top / bot )

  else if ( ax <= 4.0D+00 ) 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)

    error_f = 0.5D+00 &
      + ( 0.5D+00 - exp ( - x * x ) * top / bot )

  else if ( ax < 5.8D+00 ) then

    x2 = x * x
    t = 1.0D+00 / 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.0D+00

    error_f = ( c - top / ( x2 * bot )) / ax
    error_f = 0.5D+00 &
      + ( 0.5D+00 - exp ( - x2 ) * error_f )

  else

    error_f = 1.0D+00

  end if

  if ( x < 0.0D+00 ) then
    error_f = -error_f
  end if

  return
end
function error_fc ( ind, x )

!*****************************************************************************80
!
!! ERROR_FC evaluates the complementary error function.
!
!  Modified:
!
!    09 December 1999
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  Reference:
!
!    Armido DiDinato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, 1993, pages 360-373.
!
!  Parameters:
!
!    Input, integer IND, chooses the scaling.
!    If IND is nonzero, then the value returned has been multiplied by
!    EXP(X*X).
!
!    Input, real ( kind = 8 ) X, the argument of the function.
!
!    Output, real ( kind = 8 ) ERROR_FC, the value of the complementary 
!    error function.
!
  implicit none

  real ( kind = 8 ), dimension ( 5 ) :: a = (/ &
     0.771058495001320D-04,  -0.133733772997339D-02, &
     0.323076579225834D-01,   0.479137145607681D-01, &
     0.128379167095513D+00 /)
  real ( kind = 8 ) ax
  real ( kind = 8 ), dimension(3) :: b = (/ &
    0.301048631703895D-02, &
    0.538971687740286D-01, &
    0.375795757275549D+00 /)
  real ( kind = 8 ) bot
  real ( kind = 8 ), parameter :: c = 0.564189583547756D+00
  real ( kind = 8 ) e
  real ( kind = 8 ) error_fc
  real ( kind = 8 ) exparg
  integer ind
  real ( kind = 8 ), dimension ( 8 ) :: p = (/ &
    -1.36864857382717D-07, 5.64195517478974D-01, &
     7.21175825088309D+00, 4.31622272220567D+01, &
     1.52989285046940D+02, 3.39320816734344D+02, &
     4.51918953711873D+02, 3.00459261020162D+02 /)
  real ( kind = 8 ), dimension ( 8 ) :: q = (/  &
    1.00000000000000D+00, 1.27827273196294D+01, &
    7.70001529352295D+01, 2.77585444743988D+02, &
    6.38980264465631D+02, 9.31354094850610D+02, &
    7.90950925327898D+02, 3.00459260956983D+02 /)
  real ( kind = 8 ), dimension ( 5 ) :: r = (/ &
    2.10144126479064D+00, 2.62370141675169D+01, &
    2.13688200555087D+01, 4.65807828718470D+00, &
    2.82094791773523D-01 /)
  real ( kind = 8 ), dimension ( 4 ) :: s = (/ &
    9.41537750555460D+01, 1.87114811799590D+02, &
    9.90191814623914D+01, 1.80124575948747D+02 /)
  real ( kind = 8 ) t
  real ( kind = 8 ) top
  real ( kind = 8 ) w
  real ( kind = 8 ) x
!
!  ABS ( X ) <= 0.5
!
  ax = abs ( x )

  if ( ax <= 0.5D+00 ) then

    t = x * x

    top = (((( a(1) * t + a(2) ) * t + a(3) ) * t + a(4) ) * t + a(5) ) &
      + 1.0D+00

    bot = (( b(1) * t + b(2) ) * t + b(3) ) * t + 1.0D+00

    error_fc = 0.5D+00 + ( 0.5D+00 &
      - x * ( top / bot ) )

    if ( ind /= 0 ) then
      error_fc = exp ( t ) * error_fc
    end if

    return

  end if
!
!  0.5 < abs ( X ) <= 4
!
  if ( ax <= 4.0D+00 ) 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)

    error_fc = top / bot
!
!  4 < ABS ( X )
!
  else

    if ( x <= -5.6D+00 ) then

      if ( ind == 0 ) then
        error_fc =  2.0D+00 
      else
        error_fc =  2.0D+00  * exp ( x * x )
      end if

      return

    end if

    if ( ind == 0 ) then

      if ( 100.0D+00 < x ) then
        error_fc = 0.0D+00
        return
      end if

      if ( -exparg ( 1 ) < x * x ) then
        error_fc = 0.0D+00
        return
      end if

    end if

    t = ( 1.0D+00 / 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.0D+00

    error_fc = ( c - t * top / bot ) / ax

  end if
!
!  Final assembly.
!
  if ( ind /= 0 ) then

    if ( x < 0.0D+00 ) then
      error_fc =  2.0D+00  * exp ( x * x ) - error_fc
    end if

  else

    w = x * x
    t = w
    e = w - t
    error_fc = (( 0.5D+00 &
      + ( 0.5D+00 - e ) ) * exp ( - t ) ) * error_fc

    if ( x < 0.0D+00 ) then
      error_fc =  2.0D+00  - error_fc
    end if

  end if

  return
end
function esum ( mu, x )

subroutine f_cdf_values ( n_data, a, b, x, fx )

!*****************************************************************************80
!
!! F_CDF_VALUES returns some values of the F CDF test function.
!
!  Discussion:
!
!    The value of F_CDF ( DFN, DFD, X ) can be evaluated in Mathematica by
!    commands like:
!
!      Needs["Statistics`ContinuousDistributions`"]
!      CDF[FRatioDistribution[ DFN, DFD ], X ]
!
!  Modified:
!
!    11 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964.
!
!    Stephen Wolfram,
!    The Mathematica Book,
!    Fourth Edition,
!    Wolfram Media / Cambridge University Press, 1999.
!
!  Parameters:
!
!    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
!    first call.  On each call, the routine increments N_DATA by 1, and
!    returns the corresponding data; when there is no more data, the
!    output value of N_DATA will be 0 again.
!
!    Output, integer A, integer B, real ( kind = 8 ) X, the arguments
!    of the function.
!
!    Output, real ( kind = 8 ) FX, the value of the function.
!
  implicit none

  integer, parameter :: n_max = 20

  integer a
  integer, save, dimension ( n_max ) :: a_vec = (/ &
    1, 1, 5, 1, &
    2, 4, 1, 6, &
    8, 1, 3, 6, &
    1, 1, 1, 1, &
    2, 3, 4, 5 /)
  integer b
  integer, save, dimension ( n_max ) :: b_vec = (/ &
     1,  5,  1,  5, &
    10, 20,  5,  6, &
    16,  5, 10, 12, &
     5,  5,  5,  5, &
     5,  5,  5,  5 /)
  real ( kind = 8 ) fx
  real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ &
    0.500000D+00, 0.499971D+00, 0.499603D+00, 0.749699D+00, &
    0.750466D+00, 0.751416D+00, 0.899987D+00, 0.899713D+00, &
    0.900285D+00, 0.950025D+00, 0.950057D+00, 0.950193D+00, &
    0.975013D+00, 0.990002D+00, 0.994998D+00, 0.999000D+00, &
    0.568799D+00, 0.535145D+00, 0.514343D+00, 0.500000D+00 /)
  integer n_data
  real ( kind = 8 ) x
  real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ &
    1.00D+00,  0.528D+00, 1.89D+00,  1.69D+00, &
    1.60D+00,  1.47D+00,  4.06D+00,  3.05D+00, &
    2.09D+00,  6.61D+00,  3.71D+00,  3.00D+00, &
   10.01D+00, 16.26D+00, 22.78D+00, 47.18D+00, &
    1.00D+00,  1.00D+00,  1.00D+00,  1.00D+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
    b = 0
    x = 0.0D+00
    fx = 0.0D+00
  else
    a = a_vec(n_data)
    b = b_vec(n_data)
    x = x_vec(n_data)
    fx = fx_vec(n_data)
  end if

  return
end
subroutine f_noncentral_cdf_values ( n_data, a, b, lambda, x, fx )

!*****************************************************************************80
!
!! F_NONCENTRAL_CDF_VALUES returns some values of the F CDF test function.
!
!  Discussion:
!
!    The value of NONCENTRAL_F_CDF ( DFN, DFD, LAMDA, X ) can be evaluated 
!    in Mathematica by commands like:
!
!      Needs["Statistics`ContinuousDistributions`"]
!      CDF[NoncentralFRatioDistribution[ DFN, DFD, LAMBDA ], X ]
!
!  Modified:
!
!    12 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964.
!
!    Stephen Wolfram,
!    The Mathematica Book,
!    Fourth Edition,
!    Wolfram Media / Cambridge University Press, 1999.
!
!  Parameters:
!
!    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
!    first call.  On each call, the routine increments N_DATA by 1, and
!    returns the corresponding data; when there is no more data, the
!    output value of N_DATA will be 0 again.
!
!    Output, integer A, integer B, real ( kind = 8 ) LAMBDA, the 
!    parameters of the function.
!
!    Output, real ( kind = 8 ) X, the argument of the function.
!
!    Output, real ( kind = 8 ) FX, the value of the function.
!
  implicit none

  integer, parameter :: n_max = 22

  integer a
  integer, save, dimension ( n_max ) :: a_vec = (/ &
     1,  1,  1,  1, &
     1,  1,  1,  1, &
     1,  1,  2,  2, &
     3,  3,  4,  4, &
     5,  5,  6,  6, &
     8, 16 /)
  integer b
  integer, save, dimension ( n_max ) :: b_vec = (/ &
     1,  5,  5,  5, &
     5,  5,  5,  5, &
     5,  5,  5, 10, &
     5,  5,  5,  5, &
     1,  5,  6, 12, &
    16,  8 /)
  real ( kind = 8 ) fx
  real, save, dimension ( n_max ) :: fx_vec = (/ &
    0.500000D+00, 0.636783D+00, 0.584092D+00, 0.323443D+00, &
    0.450119D+00, 0.607888D+00, 0.705928D+00, 0.772178D+00, &
    0.819105D+00, 0.317035D+00, 0.432722D+00, 0.450270D+00, &
    0.426188D+00, 0.337744D+00, 0.422911D+00, 0.692767D+00, &
    0.363217D+00, 0.421005D+00, 0.426667D+00, 0.446402D+00, &
    0.844589D+00, 0.816368D+00 /)
  real ( kind = 8 ) lambda
  real, save, dimension ( n_max ) :: lambda_vec = (/ &
    0.00D+00,  0.000D+00, 0.25D+00,  1.00D+00, &
    1.00D+00,  1.00D+00,  1.00D+00,  1.00D+00, &
    1.00D+00,  2.00D+00,  1.00D+00,  1.00D+00, &
    1.00D+00,  2.00D+00,  1.00D+00,  1.00D+00, &
    0.00D+00,  1.00D+00,  1.00D+00,  1.00D+00, &
    1.00D+00,  1.00D+00 /)
  integer n_data
  real ( kind = 8 ) x
  real, save, dimension ( n_max ) :: x_vec = (/ &
    1.00D+00,  1.00D+00, 1.00D+00,  0.50D+00, &
    1.00D+00,  2.00D+00, 3.00D+00,  4.00D+00, &
    5.00D+00,  1.00D+00, 1.00D+00,  1.00D+00, &
    1.00D+00,  1.00D+00, 1.00D+00,  2.00D+00, &
    1.00D+00,  1.00D+00, 1.00D+00,  1.00D+00, &
    2.00D+00,  2.00D+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
    b = 0
    lambda = 0.0D+00
    x = 0.0D+00
    fx = 0.0D+00
  else
    a = a_vec(n_data)
    b = b_vec(n_data)
    lambda = lambda_vec(n_data)
    x = x_vec(n_data)
    fx = fx_vec(n_data)
  end if

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

!*****************************************************************************80
!
!! GAM1 computes 1 / GAMMA(A+1) - 1 for -0.5 <= A <= 1.5
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  Reference:
!
!    Armido DiDinato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, 1993, pages 360-373.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) A, forms the argument of the Gamma function.
!
!    Output, real ( kind = 8 ) GAM1, the value of 1 / GAMMA ( A + 1 ) - 1.
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) bot
  real ( kind = 8 ) d
  real ( kind = 8 ) gam1
  real ( kind = 8 ), parameter, dimension ( 7 ) :: p = (/ &
     0.577215664901533D+00, -0.409078193005776D+00, &
    -0.230975380857675D+00,  0.597275330452234D-01, &
     0.766968181649490D-02, -0.514889771323592D-02, &
     0.589597428611429D-03 /)
  real ( kind = 8 ), dimension ( 5 ) :: q = (/ &
    0.100000000000000D+01, 0.427569613095214D+00, &
    0.158451672430138D+00, 0.261132021441447D-01, &
    0.423244297896961D-02 /)
  real ( kind = 8 ), dimension ( 9 ) :: r = (/ &
    -0.422784335098468D+00, -0.771330383816272D+00, &
    -0.244757765222226D+00,  0.118378989872749D+00, &
     0.930357293360349D-03, -0.118290993445146D-01, &
     0.223047661158249D-02,  0.266505979058923D-03, &
    -0.132674909766242D-03 /)
  real ( kind = 8 ), parameter :: s1 = 0.273076135303957D+00
  real ( kind = 8 ), parameter :: s2 = 0.559398236957378D-01
  real ( kind = 8 ) t
  real ( kind = 8 ) top
  real ( kind = 8 ) w

  d = a - 0.5D+00

  if ( 0.0D+00 < d ) then
    t = d - 0.5D+00
  else
    t = a
  end if

  if ( t == 0.0D+00 ) then

    gam1 = 0.0D+00

  else if ( 0.0D+00 < t ) then

    top = (((((    &
            p(7)   &
      * t + p(6) ) &
      * t + p(5) ) &
      * t + p(4) ) &
      * t + p(3) ) &
      * t + p(2) ) &
      * t + p(1)

    bot = ((( q(5) * t + q(4) ) * t + q(3) ) * t + q(2) ) * t &
      + 1.0D+00

    w = top / bot

    if ( d <= 0.0D+00 ) then
      gam1 = a * w
    else
      gam1 = ( t / a ) * ( ( w - 0.5D+00 ) &
        - 0.5D+00 )
    end if

  else if ( t < 0.0D+00 ) then

    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.0D+00
    w = top / bot

    if ( d <= 0.0D+00 ) then
      gam1 = a * ( ( w + 0.5D+00 ) + 0.5D+00 )
    else
      gam1 = t * w / a
    end if

  end if

  return
end
function gamma ( a )

!*****************************************************************************80
!
!! GAMMA evaluates the gamma function.
!
!  Author:
!
!    Alfred Morris,
!    Naval Surface Weapons Center,
!    Dahlgren, Virginia.
!
!  Reference:
!
!    Armido DiDinato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, 1993, pages 360-373.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) A, the argument of the Gamma function.
!
!    Output, real ( kind = 8 ) GAMMA, the value of the Gamma function.
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) bot
  real ( kind = 8 ), parameter :: d = 0.41893853320467274178D+00
  real ( kind = 8 ) exparg
  real ( kind = 8 ) g
  real ( kind = 8 ) gamma
  integer i
  integer j
  real ( kind = 8 ) lnx
  integer m
  integer n
  real ( kind = 8 ), dimension ( 7 ) :: p = (/ &
    0.539637273585445D-03, 0.261939260042690D-02, &
    0.204493667594920D-01, 0.730981088720487D-01, &
    0.279648642639792D+00, 0.553413866010467D+00, &
    1.0D+00 /)
  real ( kind = 8 ), parameter :: pi = 3.1415926535898D+00
  real ( kind = 8 ), dimension ( 7 ) :: q = (/ &
    -0.832979206704073D-03,  0.470059485860584D-02, &
     0.225211131035340D-01, -0.170458969313360D+00, &
    -0.567902761974940D-01,  0.113062953091122D+01, &
     1.0D+00 /)
  real ( kind = 8 ), parameter :: r1 =  0.820756370353826D-03
  real ( kind = 8 ), parameter :: r2 = -0.595156336428591D-03
  real ( kind = 8 ), parameter :: r3 =  0.793650663183693D-03
  real ( kind = 8 ), parameter :: r4 = -0.277777777770481D-02
  real ( kind = 8 ), parameter :: r5 =  0.833333333333333D-01
  real ( kind = 8 ) s
  real ( kind = 8 ) t
  real ( kind = 8 ) top
  real ( kind = 8 ) w
  real ( kind = 8 ) x
  real ( kind = 8 ) z

  gamma = 0.0D+00
  x = a

  if ( abs ( a ) < 15.0D+00 ) then
!
!  Evaluation of GAMMA(A) for |A| < 15
!
    t = 1.0D+00
    m = int ( a ) - 1
!
!  Let T be the product of A-J when 2 <= A.
!
    if ( 0 <= m ) then

      do j = 1, m
        x = x - 1.0D+00
        t = x * t
      end do

      x = x - 1.0D+00
!
!  Let T be the product of A+J WHEN A < 1
!
    else

      t = a

      if ( a <= 0.0D+00 ) then

        m = - m - 1

        do j = 1, m
          x = x + 1.0D+00
          t = x * t
        end do

        x = ( x + 0.5D+00 ) + 0.5D+00
        t = x * t
        if ( t == 0.0D+00 ) then
          return
        end if

      end if
!
!  Check if 1/T can overflow. 
!
      if ( abs ( t ) < 1.0D-30 ) then
        if ( 1.0001D+00 < abs ( t ) * huge ( t ) ) then
          gamma = 1.0D+00 / t
        end if
        return
      end if

    end if
!
!  Compute Gamma(1 + X) for 0 <= X < 1.
!
    top = p(1)
    bot = q(1)
    do i = 2, 7
      top = top * x + p(i)
      bot = bot * x + q(i)
    end do

    gamma = top / bot
!
!  Termination.
!
    if ( 1.0D+00 <= a ) then
      gamma = gamma * t
    else
      gamma = gamma / t
    end if
!
!  Evaluation of Gamma(A) FOR 15 <= ABS ( A ).
!
  else

    if ( 1000.0D+00 <= abs ( a ) ) then
      return
    end if

    if ( a <= 0.0D+00 ) then

      x = -a
      n = x
      t = x - n

      if ( 0.9D+00 < t ) then
        t = 1.0D+00 - t
      end if

      s = sin ( pi * t ) / pi

      if ( mod ( n, 2 ) == 0 ) then
        s = -s
      end if

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

    end if
!
!  Compute the modified asymptotic sum.
!
    t = 1.0D+00 / ( x * x )

    g = (((( r1 * t + r2 ) * t + r3 ) * t + r4 ) * t + r5 ) / x

    lnx = log ( x )
!
!  Final assembly.
!
    z = x
    g = ( d + g ) + ( z - 0.5D+00 ) &
      * ( lnx - 1.0D+00 )
    w = g
    t = g - real ( w, kind = 8 )
  
    if ( 0.99999D+00 * exparg ( 0 ) < w ) then
      return
    end if

    gamma = exp ( w )* ( 1.0D+00 + t )

    if ( a < 0.0D+00 ) then
      gamma = ( 1.0D+00 / ( gamma * s ) ) / x
    end if

  end if

  return
end
subroutine gamma_inc ( a, x, ans, qans, ind )

!*****************************************************************************80
!
!! GAMMA_INC evaluates the incomplete gamma ratio functions P(A,X) and Q(A,X).
!
!  Modified:
!
!    16 April 2005
!
!  Author:
!
!    Alfred Morris,
!    Naval Surface Weapons Center,
!    Dahlgren, Virginia.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) A, X, the arguments of the incomplete
!    gamma ratio.  A and X must be nonnegative.  A and X cannot
!    both be zero.
!
!    Output, real ( kind = 8 ) ANS, QANS.  On normal output,
!    ANS = P(A,X) and QANS = Q(A,X).  However, ANS is set to 2 if
!    A or X is negative, or both are 0, or when the answer is
!    computationally indeterminate because A is extremely large
!    and X is very close to A.
!
!    Input, integer IND, indicates the accuracy request:
!    0, as much accuracy as possible.
!    1, to within 1 unit of the 6-th significant digit, 
!    otherwise, to within 1 unit of the 3rd significant digit.
!
!  Local Parameters:
!
!     ALOG10 = LN(10)
!     RT2PIN = 1/SQRT(2*PI)
!     RTPI   = SQRT(PI)
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) a2n
  real ( kind = 8 ) a2nm1
  real ( kind = 8 ) acc
  real ( kind = 8 ), dimension ( 3 ) :: acc0 = (/ &
    5.0D-15, 5.0D-07, 5.0D-04 /)
  real ( kind = 8 ), parameter :: alog10 = 2.30258509299405D+00
  real ( kind = 8 ) am0
  real ( kind = 8 ) amn
  real ( kind = 8 ) an
  real ( kind = 8 ) an0
  real ( kind = 8 ) ans
  real ( kind = 8 ) apn
  real ( kind = 8 ) b2n
  real ( kind = 8 ) b2nm1
  real ( kind = 8 ) big(3)
  real ( kind = 8 ) c
  real ( kind = 8 ) c0
  real ( kind = 8 ) c1
  real ( kind = 8 ) c2
  real ( kind = 8 ) c3
  real ( kind = 8 ) c4
  real ( kind = 8 ) c5
  real ( kind = 8 ) c6
  real ( kind = 8 ) cma
  real ( kind = 8 ) d0(13)
  real ( kind = 8 ) d1(12)
  real ( kind = 8 ) d2(10)
  real ( kind = 8 ) d3(8)
  real ( kind = 8 ) d4(6)
  real ( kind = 8 ) d5(4)
  real ( kind = 8 ) d6(2)
  real ( kind = 8 ) d10
  real ( kind = 8 ) d20
  real ( kind = 8 ) d30
  real ( kind = 8 ) d40
  real ( kind = 8 ) d50
  real ( kind = 8 ) d60
  real ( kind = 8 ) d70
  real ( kind = 8 ) e
  real ( kind = 8 ) e0
  real ( kind = 8 ) e00(3)
  real ( kind = 8 ) error_f
  real ( kind = 8 ) error_fc
  real ( kind = 8 ) g
  real ( kind = 8 ) gam1
  real ( kind = 8 ) gamma
  real ( kind = 8 ) h
  integer i
  integer ind
  integer iop
  real ( kind = 8 ) j
  real ( kind = 8 ) l
  integer m
  integer n
  integer n_max
  real ( kind = 8 ) qans
  real ( kind = 8 ) r
  real ( kind = 8 ) rexp
  real ( kind = 8 ) rlog
  real ( kind = 8 ), parameter :: rt2pin = 0.398942280401433D+00
  real ( kind = 8 ) rta
  real ( kind = 8 ), parameter :: rtpi = 1.77245385090552D+00
  real ( kind = 8 ) rtx
  real ( kind = 8 ) s
  real ( kind = 8 ) sum1
  real ( kind = 8 ) t
  real ( kind = 8 ) t1
  real ( kind = 8 ) tol
  real ( kind = 8 ) twoa
  real ( kind = 8 ) u
  real ( kind = 8 ) w
  real ( kind = 8 ) wk(20)
  real ( kind = 8 ) x
  real ( kind = 8 ) x0
  real ( kind = 8 ) x00(3)
  real ( kind = 8 ) y
  real ( kind = 8 ) z

  data big(1)/20.0D+00/,big(2)/14.0D+00/,big(3)/10.0D+00/
  data e00(1)/0.25D-03/,e00(2)/0.25D-01/,e00(3)/0.14D+00/
  data x00(1)/31.0D+00/,x00(2)/17.0D+00/,x00(3)/9.7D+00/
  data d0(1)/0.833333333333333D-01/
  data d0(2)/-0.148148148148148D-01/
  data d0(3)/0.115740740740741D-02/,d0(4)/0.352733686067019D-03/
  data d0(5)/-0.178755144032922D-03/,d0(6)/0.391926317852244D-04/
  data d0(7)/-0.218544851067999D-05/,d0(8)/-0.185406221071516D-05/
  data d0(9)/0.829671134095309D-06/,d0(10)/-0.176659527368261D-06/
  data d0(11)/0.670785354340150D-08/,d0(12)/0.102618097842403D-07/
  data d0(13)/-0.438203601845335D-08/
  data d10/-0.185185185185185D-02/,d1(1)/-0.347222222222222D-02/
  data d1(2)/0.264550264550265D-02/,d1(3)/-0.990226337448560D-03/
  data d1(4)/0.205761316872428D-03/,d1(5)/-0.401877572016461D-06/
  data d1(6)/-0.180985503344900D-04/,d1(7)/0.764916091608111D-05/
  data d1(8)/-0.161209008945634D-05/,d1(9)/0.464712780280743D-08/
  data d1(10)/0.137863344691572D-06/,d1(11)/-0.575254560351770D-07/
  data d1(12)/0.119516285997781D-07/
  data d20/0.413359788359788D-02/,d2(1)/-0.268132716049383D-02/
  data d2(2)/0.771604938271605D-03/,d2(3)/0.200938786008230D-05/
  data d2(4)/-0.107366532263652D-03/,d2(5)/0.529234488291201D-04/
  data d2(6)/-0.127606351886187D-04/,d2(7)/0.342357873409614D-07/
  data d2(8)/0.137219573090629D-05/,d2(9)/-0.629899213838006D-06/
  data d2(10)/0.142806142060642D-06/
  data d30/0.649434156378601D-03/,d3(1)/0.229472093621399D-03/
  data d3(2)/-0.469189494395256D-03/,d3(3)/0.267720632062839D-03/
  data d3(4)/-0.756180167188398D-04/,d3(5)/-0.239650511386730D-06/
  data d3(6)/0.110826541153473D-04/,d3(7)/-0.567495282699160D-05/
  data d3(8)/0.142309007324359D-05/
  data d40/-0.861888290916712D-03/,d4(1)/0.784039221720067D-03/
  data d4(2)/-0.299072480303190D-03/,d4(3)/-0.146384525788434D-05/
  data d4(4)/0.664149821546512D-04/,d4(5)/-0.396836504717943D-04/
  data d4(6)/0.113757269706784D-04/
  data d50/-0.336798553366358D-03/,d5(1)/-0.697281375836586D-04/
  data d5(2)/0.277275324495939D-03/,d5(3)/-0.199325705161888D-03/
  data d5(4)/0.679778047793721D-04/
  data d60/0.531307936463992D-03/,d6(1)/-0.592166437353694D-03/
  data d6(2)/0.270878209671804D-03/
  data d70 / 0.344367606892378D-03/

  e = epsilon ( 1.0D+00 )

  if ( a < 0.0D+00 .or. x < 0.0D+00 ) then
    ans = 2.0D+00
    return
  end if

  if ( a == 0.0D+00 .and. x == 0.0D+00 ) then
    ans = 2.0D+00
    return
  end if

  if ( a * x == 0.0D+00 ) then
    if ( x <= a ) then
      ans = 0.0D+00
      qans = 1.0D+00
    else
      ans = 1.0D+00
      qans = 0.0D+00
    end if
    return
  end if

  iop = ind + 1
  if ( iop /= 1 .and. iop /= 2 ) iop = 3
  acc = max ( acc0(iop), e )
  e0 = e00(iop)
  x0 = x00(iop)
!
!  Select the appropriate algorithm.
!
  if ( 1.0D+00 <= a ) then
    go to 10
  end if

  if ( a == 0.5D+00 ) then
    go to 390
  end if

  if ( x < 1.1D+00 ) then
    go to 160
  end if

  t1 = a * log ( x ) - x
  u = a * exp ( t1 )

  if ( u == 0.0D+00 ) then
    ans = 1.0D+00
    qans = 0.0D+00
    return
  end if

  r = u * ( 1.0D+00 + gam1 ( a ) )
  go to 250

   10 continue

  if ( big(iop) <= a ) then
    go to 30
  end if

  if ( x < a .or. x0 <= x ) then
    go to 20
  end if

  twoa = a + a
  m = int ( twoa )

  if ( twoa == real ( m, kind = 8 ) ) then
    i = m / 2
    if ( a == real ( i, kind = 8 ) ) then
      go to 210
    end if
    go to 220
  end if

   20 continue

  t1 = a * log ( x ) - x
  r = exp ( t1 ) / gamma ( a )
  go to 40

   30 continue

  l = x / a

  if ( l == 0.0D+00 ) then
    ans = 0.0D+00
    qans = 1.0D+00
    return
  end if

  s = 0.5D+00 + ( 0.5D+00 - l )
  z = rlog ( l )
  if ( 700.0D+00 / a <= z ) then
    go to 410
  end if

  y = a * z
  rta = sqrt ( a )

  if ( abs ( s ) <= e0 / rta ) then
    go to 330
  end if

  if ( abs ( s ) <= 0.4D+00 ) then
    go to 270
  end if

  t = ( 1.0D+00 / a )**2
  t1 = ((( 0.75D+00 * t - 1.0D+00 ) * t + 3.5D+00 ) &
    * t - 105.0D+00 ) / ( a * 1260.0D+00 )
  t1 = t1 - y
  r = rt2pin * rta * exp ( t1 )

40    continue

  if ( r == 0.0D+00 ) then
    if ( x <= a ) then
      ans = 0.0D+00
      qans = 1.0D+00
    else
      ans = 1.0D+00
      qans = 0.0D+00
    end if
    return
  end if

  if ( x <= max ( a, alog10 ) ) then
    go to 50
  end if

  if ( x < x0 ) then
    go to 250
  end if

  go to 100
!
!  Taylor series for P/R.
!
50    continue

  apn = a + 1.0D+00
  t = x / apn
  wk(1) = t

  n = 20

  do i = 2, 20
    apn = apn + 1.0D+00
    t = t * ( x / apn )
    if ( t <= 1.0D-03 ) then
      n = i
      exit
    end if
    wk(i) = t
  end do

  sum1 = t

  tol = 0.5D+00 * acc

  do

    apn = apn + 1.0D+00
    t = t * ( x / apn )
    sum1 = sum1 + t

    if ( t <= tol ) then
      exit
    end if

  end do

  n_max = n - 1
  do m = 1, n_max
    n = n - 1
    sum1 = sum1 + wk(n)
  end do

  ans = ( r / a ) * ( 1.0D+00 + sum1 )
  qans = 0.5D+00 + ( 0.5D+00 - ans )
  return
!
!  Asymptotic expansion.
!
  100 continue

  amn = a - 1.0D+00
  t = amn / x
  wk(1) = t

  n = 20

  do i = 2, 20
    amn = amn - 1.0D+00
    t = t * ( amn / x )
    if ( abs ( t ) <= 1.0D-03 ) then
      n = i
      exit
    end if
    wk(i) = t
  end do

  sum1 = t

  do

    if ( abs ( t ) <= acc ) then
      exit
    end if

    amn = amn - 1.0D+00
    t = t * ( amn / x )
    sum1 = sum1 + t
 
  end do

  n_max = n - 1
  do m = 1, n_max
    n = n - 1
    sum1 = sum1 + wk(n)
  end do
  qans = ( r / x ) * ( 1.0D+00 + sum1 )
  ans = 0.5D+00 + ( 0.5D+00 - qans )
  return
!
!  Taylor series for P(A,X)/X**A
!
  160 continue

  an = 3.0D+00
  c = x
  sum1 = x / ( a + 3.0D+00 )
  tol = 3.0D+00 * acc / ( a + 1.0D+00 )

  do

    an = an + 1.0D+00
    c = -c * ( x / an )
    t = c / ( a + an )
    sum1 = sum1 + t

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

  end do

  j = a * x * ( ( sum1 / 6.0D+00 - 0.5D+00 / &
    ( a +  2.0D+00  ) ) * x + 1.0D+00 &
    / ( a + 1.0D+00 ) )

  z = a * log ( x )
  h = gam1 ( a )
  g = 1.0D+00 + h

  if ( x < 0.25D+00 ) then
    go to 180
  end if

  if ( a < x / 2.59D+00 ) then
    go to 200
  end if

  go to 190

  180 continue

  if ( -0.13394D+00 < z ) then
    go to 200
  end if

  190 continue

  w = exp ( z )
  ans = w * g * ( 0.5D+00 + ( 0.5D+00 - j ))
  qans = 0.5D+00 + ( 0.5D+00 - ans )
  return

200   continue

  l = rexp ( z )
  w = 0.5D+00 + ( 0.5D+00 + l )
  qans = ( w * j - l ) * g - h

  if ( qans < 0.0D+00 ) then
    ans = 1.0D+00
    qans = 0.0D+00
    return
  end if

  ans = 0.5D+00 + ( 0.5D+00 - qans )
  return
!
!  Finite sums for Q when 1 <= A and 2*A is an integer.
!
210   continue

  sum1 = exp ( - x )
  t = sum1
  n = 1
  c = 0.0D+00
  go to 230

220   continue

  rtx = sqrt ( x )
  sum1 = error_fc ( 0, rtx )
  t = exp ( -x ) / ( rtpi * rtx )
  n = 0
  c = -0.5D+00

  230 continue

  do while ( n /= i )
    n = n + 1
    c = c + 1.0D+00
    t = ( x * t ) / c
    sum1 = sum1 + t
  end do

  240 continue

  qans = sum1
  ans = 0.5D+00 + ( 0.5D+00 - qans )
  return
!
!  Continued fraction expansion.
!
250   continue

  tol = max ( 5.0D+00 * e, acc )
  a2nm1 = 1.0D+00
  a2n = 1.0D+00
  b2nm1 = x
  b2n = x + ( 1.0D+00 - a )
  c = 1.0D+00

  do

    a2nm1 = x * a2n + c * a2nm1
    b2nm1 = x * b2n + c * b2nm1
    am0 = a2nm1 / b2nm1
    c = c + 1.0D+00
    cma = c - a
    a2n = a2nm1 + cma * a2n
    b2n = b2nm1 + cma * b2n
    an0 = a2n / b2n

    if ( abs ( an0 - am0 ) < tol * an0 ) then
      exit
    end if

  end do

  qans = r * an0
  ans = 0.5D+00 + ( 0.5D+00 - qans )
  return
!
!  General Temme expansion.
!
270   continue

  if ( abs ( s ) <= 2.0D+00 * e .and. 3.28D-03 < a * e * e ) then
    ans =  2.0D+00 
    return
  end if

  c = exp ( - y )
  w = 0.5D+00 * error_fc ( 1, sqrt ( y ) )
  u = 1.0D+00 / a
  z = sqrt ( z + z )

  if ( l < 1.0D+00 ) then
    z = -z
  end if

  if ( iop < 2 ) then

    if ( abs ( s ) <= 1.0D-03 ) then

      c0 = ((((((     &
              d0(7)   &
        * z + d0(6) ) &
        * z + d0(5) ) &
        * z + d0(4) ) &
        * z + d0(3) ) &
        * z + d0(2) ) &
        * z + d0(1) ) &
        * z - 1.0D+00 / 3.0D+00

      c1 = (((((      &
              d1(6)   &
        * z + d1(5) ) &
        * z + d1(4) ) &
        * z + d1(3) ) &
        * z + d1(2) ) &
        * z + d1(1) ) &
        * z + d10

      c2 = ((((d2(5)*z+d2(4))*z+d2(3))*z+d2(2))*z+d2(1))*z + d20

      c3 = (((d3(4)*z+d3(3))*z+d3(2))*z+d3(1))*z + d30

      c4 = ( d4(2) * z + d4(1) ) * z + d40
      c5 = ( d5(2) * z + d5(1) ) * z + d50
      c6 = d6(1) * z + d60

      t = (((((( d70 &
        * u + c6 ) &
        * u + c5 ) &
        * u + c4 ) &
        * u + c3 ) &
        * u + c2 ) &
        * u + c1 ) &
        * u + c0

    else

      c0 = (((((((((((( &
              d0(13)   &
        * z + d0(12) ) &
        * z + d0(11) ) &
        * z + d0(10) ) &
        * z + d0(9)  ) &
        * z + d0(8)  ) &
        * z + d0(7)  ) &
        * z + d0(6)  ) &
        * z + d0(5)  ) &
        * z + d0(4)  ) &
        * z + d0(3)  ) &
        * z + d0(2)  ) &
        * z + d0(1)  ) &
        * z - 1.0D+00 / 3.0D+00

      c1 = ((((((((((( &
                d1(12) &
          * z + d1(11) &
        ) * z + d1(10) &
        ) * z + d1(9)  &
        ) * z + d1(8)  &
        ) * z + d1(7)  &
        ) * z + d1(6)  &
        ) * z + d1(5)  &
        ) * z + d1(4)  & 
        ) * z + d1(3)  &
        ) * z + d1(2)  &
        ) * z + d1(1)  &
        ) * z + d10

      c2 = ((((((((( &
                d2(10) &
          * z + d2(9) &
        ) * z + d2(8) &
        ) * z + d2(7) &
        ) * z + d2(6) &
        ) * z + d2(5) &
        ) * z + d2(4) &
        ) * z + d2(3) &
        ) * z + d2(2) &
        ) * z + d2(1) &
        ) * z + d20 

      c3 = ((((((( &
                d3(8) &
          * z + d3(7) &
        ) * z + d3(6) &
        ) * z + d3(5) &
        ) * z + d3(4) &
        ) * z + d3(3) &
        ) * z + d3(2) &
        ) * z + d3(1) &
        ) * z + d30

      c4 = ((((( d4(6)*z+d4(5))*z+d4(4))*z+d4(3))*z+d4(2))*z+d4(1))*z + d40

      c5 = (((d5(4)*z+d5(3))*z+d5(2))*z+d5(1))*z + d50

      c6 = ( d6(2) * z + d6(1) ) * z + d60

      t = ((((((   &
            d70    &
        * u + c6 ) &
        * u + c5 ) &
        * u + c4 ) &
        * u + c3 ) &
        * u + c2 ) &
        * u + c1 ) &
        * u + c0

    end if

  else if ( iop == 2 ) then

    c0 = (((((      &
            d0(6)   &
      * z + d0(5) ) &
      * z + d0(4) ) &
      * z + d0(3) ) &
      * z + d0(2) ) &
      * z + d0(1) ) &
      * z - 1.0D+00 / 3.0D+00

    c1 = ((( d1(4) * z + d1(3) ) * z + d1(2) ) * z + d1(1) ) * z + d10
    c2 = d2(1) * z + d20
    t = ( c2 * u + c1 ) * u + c0

  else if ( 2 < iop ) then

    t = (( d0(3) * z + d0(2) ) * z + d0(1) ) * z - 1.0D+00 / 3.0D+00

  end if

310   continue

  if ( 1.0D+00 <= l ) then
    qans = c * ( w + rt2pin * t / rta )
    ans = 0.5D+00 + ( 0.5D+00 - qans )
  else
    ans = c * ( w - rt2pin * t / rta )
    qans = 0.5D+00 + ( 0.5D+00 - ans )
  end if

  return
!
!  Temme expansion for L = 1
!
  330 continue

  if ( 3.28D-03 < a * e * e ) then
    ans =  2.0D+00 
    return
  end if

  c = 0.5D+00 + ( 0.5D+00 - y )
  w = ( 0.5D+00 - sqrt ( y ) &
    * ( 0.5D+00 &
    + ( 0.5D+00 - y / 3.0D+00 ) ) / rtpi ) / c
  u = 1.0D+00 / a
  z = sqrt ( z + z )

  if ( l < 1.0D+00 ) then
    z = -z
  end if

  if ( iop < 2 ) then

    c0 = ((((((     &
            d0(7)   &
      * z + d0(6) ) &
      * z + d0(5) ) &
      * z + d0(4) ) &
      * z + d0(3) ) &
      * z + d0(2) ) &
      * z + d0(1) ) &
      * z - 1.0D+00 / 3.0D+00

    c1 = (((((      &
            d1(6)   &
      * z + d1(5) ) &
      * z + d1(4) ) &
      * z + d1(3) ) &
      * z + d1(2) ) &
      * z + d1(1) ) &
      * z + d10

    c2 = ((((d2(5)*z+d2(4))*z+d2(3))*z+d2(2))*z+d2(1))*z + d20

    c3 = (((d3(4)*z+d3(3))*z+d3(2))*z+d3(1))*z + d30

    c4 = ( d4(2) * z + d4(1) ) * z + d40
    c5 = ( d5(2) * z + d5(1) ) * z + d50
    c6 = d6(1) * z + d60

    t = (((((( d70 &
      * u + c6 ) &
      * u + c5 ) &
      * u + c4 ) &
      * u + c3 ) &
      * u + c2 ) &
      * u + c1 ) &
      * u + c0

  else if ( iop == 2 ) then

    c0 = ( d0(2) * z + d0(1) ) * z - 1.0D+00 / 3.0D+00
    c1 = d1(1) * z + d10
    t = ( d20 * u + c1 ) * u + c0
  
  else if ( 2 < iop ) then

    t = d0(1) * z - 1.0D+00 / 3.0D+00

  end if

  go to 310
!
!  Special cases
!
  390 continue

  if ( x < 0.25D+00 ) then
    ans = error_f ( sqrt ( x ) )
    qans = 0.5D+00 + ( 0.5D+00 - ans )
  else
    qans = error_fc ( 0, sqrt ( x ) )
    ans = 0.5D+00 + ( 0.5D+00 - qans )
  end if

  return

  410 continue

  if ( abs ( s ) <= 2.0D+00 * e ) then
    ans =  2.0D+00 
    return
  end if

  if ( x <= a ) then
    ans = 0.0D+00
    qans = 1.0D+00
  else
    ans = 1.0D+00
    qans = 0.0D+00
  end if

  return
end
subroutine gamma_inc_inv ( a, x, x0, p, q, ierr )

!*****************************************************************************80
!
!! GAMMA_INC_INV computes the inverse incomplete gamma ratio function.
!
!  Discussion:
!
!    The routine is given positive A, and nonnegative P and Q where P + Q = 1.
!    The value X is computed with the property that P(A,X) = P and Q(A,X) = Q.  
!    Schroder iteration is employed.  The routine attempts to compute X
!    to 10 significant digits if this is possible for the particular computer 
!    arithmetic being used.
!
!  Author:
!
!    Alfred Morris,
!    Naval Surface Weapons Center,
!    Dahlgren, Virginia.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) A, the parameter in the incomplete gamma
!    ratio.  A must be positive.
!
!    Output, real ( kind = 8 ) X, the computed point for which the
!    incomplete gamma functions have the values P and Q.
!
!    Input, real ( kind = 8 ) X0, an optional initial approximation
!    for the solution X.  If the user does not want to supply an
!    initial approximation, then X0 should be set to 0, or a negative
!    value.
!
!    Input, real ( kind = 8 ) P, Q, the values of the incomplete gamma
!    functions, for which the corresponding argument is desired.
!
!    Output, integer IERR, error flag.
!    0, the solution was obtained. Iteration was not used.
!    0 < K, The solution was obtained. IERR iterations were performed.
!    -2, A <= 0
!    -3, No solution was obtained. The ratio Q/A is too large.
!    -4, P + Q /= 1
!    -6, 20 iterations were performed. The most recent value obtained 
!        for X is given.  This cannot occur if X0 <= 0.
!    -7, Iteration failed. No value is given for X.
!        This may occur when X is approximately 0.
!    -8, A value for X has been obtained, but the routine is not certain
!        of its accuracy.  Iteration cannot be performed in this
!        case. If X0 <= 0, this can occur only when P or Q is 
!        approximately 0. If X0 is positive then this can occur when A is
!        exceedingly close to X and A is extremely large (say A .GE. 1.E20).
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ), parameter :: a0 = 3.31125922108741D+00
  real ( kind = 8 ), parameter :: a1 = 11.6616720288968D+00
  real ( kind = 8 ), parameter :: a2 = 4.28342155967104D+00
  real ( kind = 8 ), parameter :: a3 = 0.213623493715853D+00
  real ( kind = 8 ) alnrel
  real ( kind = 8 ) am1
  real ( kind = 8 ) amax
  real ( kind = 8 ), dimension(2) :: amin = (/ &
    500.0D+00, 100.0D+00 /)
  real ( kind = 8 ) ap1
  real ( kind = 8 ) ap2
  real ( kind = 8 ) ap3
  real ( kind = 8 ) apn
  real ( kind = 8 ) b
  real ( kind = 8 ), parameter :: b1 = 6.61053765625462D+00
  real ( kind = 8 ), parameter :: b2 = 6.40691597760039D+00
  real ( kind = 8 ), parameter :: b3 = 1.27364489782223D+00
  real ( kind = 8 ), parameter :: b4 = .036117081018842D+00
  real ( kind = 8 ), dimension ( 2 ) :: bmin = (/ &
    1.0D-28, 1.0D-13 /)
  real ( kind = 8 ), parameter :: c = 0.577215664901533D+00
  real ( kind = 8 ) c1
  real ( kind = 8 ) c2
  real ( kind = 8 ) c3
  real ( kind = 8 ) c4
  real ( kind = 8 ) c5
  real ( kind = 8 ) d
  real ( kind = 8 ), dimension ( 2 ) :: dmin = (/ &
    1.0D-06, 1.0D-04 /)
  real ( kind = 8 ) e
  real ( kind = 8 ) e2
  real ( kind = 8 ), dimension ( 2 ) :: emin = (/ &
    2.0D-03, 6.0D-03 /)
  real ( kind = 8 ) eps
  real ( kind = 8 ), dimension ( 2 ) :: eps0 = (/ &
    1.0D-10, 1.0D-08 /)
  real ( kind = 8 ) g
  real ( kind = 8 ) gamma_log
  real ( kind = 8 ) gamma_ln1
  real ( kind = 8 ) gamma
  real ( kind = 8 ) h
  real ( kind = 8 ), parameter :: half = 0.5D+00
  integer ierr
  integer iop
  real ( kind = 8 ), parameter :: ln10 = 2.302585D+00
  real ( kind = 8 ) p
  real ( kind = 8 ) pn
  real ( kind = 8 ) q
  real ( kind = 8 ) qg
  real ( kind = 8 ) qn
  real ( kind = 8 ) r
  real ( kind = 8 ) rcomp
  real ( kind = 8 ) rta
  real ( kind = 8 ) s
  real ( kind = 8 ) s2
  real ( kind = 8 ) sum1
  real ( kind = 8 ) t
  real ( kind = 8 ), parameter :: tol = 1.0D-05
  real ( kind = 8 ), parameter :: two =  2.0D+00 
  real ( kind = 8 ) u
  real ( kind = 8 ) w
  real ( kind = 8 ) x
  real ( kind = 8 ) x0
  real ( kind = 8 ) xn
  real ( kind = 8 ) y
  real ( kind = 8 ) z

  e = epsilon ( e )

  x = 0.0D+00

  if ( a <= 0.0D+00 ) then
    ierr = -2
    return
  end if

  t = p + q - 1.0D+00

  if ( e < abs ( t ) ) then
    ierr = -4
    return
  end if

  ierr = 0

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

  if ( q == 0.0D+00 ) then
    x = huge ( x )
    return
  end if

  if ( a == 1.0D+00 ) then
    if ( 0.9D+00 <= q ) then
      x = -alnrel ( - p )
    else
      x = -log ( q )
    end if
    return
  end if

  e2 = two * e
  amax = 0.4D-10 / ( e * e )

  if ( 1.0D-10 < e ) then
    iop = 2
  else
    iop = 1
  end if

  eps = eps0(iop)
  xn = x0

  if ( 0.0D+00 < x0 ) then
    go to 160
  end if
!
!  Selection of the initial approximation XN of X when A < 1.
!
  if ( 1.0D+00 < a ) then
    go to 80
  end if

  g = gamma ( a + 1.0D+00 )
  qg = q * g

  if ( qg == 0.0D+00 ) then
    x = huge ( x )
    ierr = -8
    return
  end if

  b = qg / a

  if ( 0.6D+00 * a < qg ) then
    go to 40
  end if

  if ( a < 0.30D+00 .and. 0.35D+00 <= b ) then
    t = exp ( - ( b + c ) )
    u = t * exp ( t )
    xn = t * exp ( u )
    go to 160
  end if

  if ( 0.45D+00 <= b ) then
    go to 40
  end if

  if ( b == 0.0D+00 ) then
    x = huge ( x )
    ierr = -8
    return
  end if

  y = -log ( b ) 
  s = half + ( half - a )
  z = log ( y ) 
  t = y - s * z

  if ( 0.15D+00 <= b ) then
    xn = y - s * log ( t ) - log ( 1.0D+00 + s / ( t + 1.0D+00 ) )
    go to 220
  end if

  if ( 0.01D+00 < b ) then
    u = ( ( t + two * ( 3.0D+00 - a ) ) * t &
      + ( two - a ) * ( 3.0D+00 - a )) / &
      ( ( t + ( 5.0D+00 - a ) ) * t + two )
    xn = y - s * log ( t ) - log ( u )
    go to 220
  end if

30 continue

  c1 = -s * z
  c2 = -s * ( 1.0D+00 + c1 )

  c3 = s * (( half * c1 &
    + ( two - a ) ) * c1 + ( 2.5D+00 - 1.5D+00 * a ) )

  c4 = -s * ((( c1 / 3.0D+00 + ( 2.5D+00 - 1.5D+00 * a ) ) * c1 &
    + ( ( a - 6.0D+00 ) * a + 7.0D+00 ) ) &
    * c1 + ( ( 11.0D+00 * a - 46.0D+00 ) * a + 47.0D+00 ) / 6.0D+00 )

  c5 = -s * (((( - c1 / 4.0D+00 + ( 11.0D+00 * a - 17.0D+00 ) / 6.0D+00 ) * c1 &
     + ( ( -3.0D+00 * a + 13.0D+00 ) * a - 13.0D+00 ) ) * c1 &
     + half &
     * ( ( ( two * a - 25.0D+00 ) * a + 72.0D+00 ) &
     * a - 61.0D+00 ) ) * c1 &
     + ( ( ( 25.0D+00 * a - 195.0D+00 ) * a &
     + 477.0D+00 ) * a - 379.0D+00 ) / 12.0D+00 )

  xn = (((( c5 / y + c4 ) / y + c3 ) / y + c2 ) / y + c1 ) + y

  if ( 1.0D+00 < a ) then
    go to 220
  end if

  if ( bmin(iop) < b ) then
    go to 220
  end if

  x = xn
  return

   40 continue

  if ( b * q <= 1.0D-08 ) then
    xn = exp ( - ( q / a + c ))
  else if ( 0.9D+00 < p ) then
    xn = exp ( ( alnrel ( - q ) + gamma_ln1 ( a )) / a )
  else
    xn = exp ( log ( p * g ) / a )
  end if

  if ( xn == 0.0D+00 ) then
    ierr = -3
    return
  end if

  t = half + ( half - xn / ( a + 1.0D+00 ))
  xn = xn / t
  go to 160
!
!  Selection of the initial approximation XN of X when 1 < A.
!
   80 continue

  if ( 0.5D+00 < q ) then
    w = log ( p )
  else
    w = log ( q )
  end if

  t = sqrt ( - two * w )

  s = t - ((( a3 * t + a2 ) * t + a1 ) * t + a0 ) / (((( &
    b4 * t + b3 ) * t + b2 ) * t + b1 ) * t + 1.0D+00 )

  if ( 0.5D+00 < q ) then
    s = -s
  end if

  rta = sqrt ( a )
  s2 = s * s

  xn = a + s * rta + ( s2 - 1.0D+00 ) / 3.0D+00 + s * ( s2 - 7.0D+00 ) &
    / ( 36.0D+00 * rta ) - ( ( 3.0D+00 * s2 + 7.0D+00 ) * s2 - 16.0D+00 ) &
    / ( 810.0D+00 * a ) + s * (( 9.0D+00 * s2 + 256.0D+00 ) * s2 - 433.0D+00 ) &
    / ( 38880.0D+00 * a * rta )

  xn = max ( xn, 0.0D+00 )

  if ( amin(iop) <= a ) then

    x = xn
    d = half + ( half - x / a )

    if ( abs ( d ) <= dmin(iop) ) then
      return
    end if

  end if

  110 continue

  if ( p <= 0.5D+00 ) then
    go to 130
  end if

  if ( xn < 3.0D+00 * a ) then
    go to 220
  end if

  y = - ( w + gamma_log ( a ) )
  d = max ( two, a * ( a - 1.0D+00 ) )

  if ( ln10 * d <= y ) then
    s = 1.0D+00 - a
    z = log ( y )
    go to 30
  end if

  120 continue

  t = a - 1.0D+00
  xn = y + t * log ( xn ) - alnrel ( -t / ( xn + 1.0D+00 ) )
  xn = y + t * log ( xn ) - alnrel ( -t / ( xn + 1.0D+00 ) )
  go to 220

  130 continue

  ap1 = a + 1.0D+00

  if ( 0.70D+00 * ap1 < xn ) then
    go to 170
  end if

  w = w + gamma_log ( ap1 )

  if ( xn <= 0.15 * ap1 ) then
    ap2 = a + two
    ap3 = a + 3.0D+00
    x = exp ( ( w + x ) / a )
    x = exp ( ( w + x - log ( 1.0D+00 + ( x / ap1 ) &
      * ( 1.0D+00 + x / ap2 ) ) ) / a )
    x = exp ( ( w + x - log ( 1.0D+00 + ( x / ap1 ) &
      * ( 1.0D+00 + x / ap2 ) ) ) / a )
    x = exp ( ( w + x - log ( 1.0D+00 + ( x / ap1 ) &
      * ( 1.0D+00 + ( x / ap2 ) &
      * ( 1.0D+00 + x / ap3 ) ) ) ) / a )
    xn = x

    if ( xn <= 1.0D-02 * ap1 ) then
      if ( xn <= emin(iop) * ap1 ) then
        return
      end if
      go to 170
    end if

  end if

  apn = ap1
  t = xn / apn
  sum1 = 1.0D+00 + t

  do

    apn = apn + 1.0D+00
    t = t * ( xn / apn )
    sum1 = sum1 + t

    if ( t <= 1.0D-04 ) then
      exit
    end if

  end do

  t = w - log ( sum1 )
  xn = exp ( ( xn + t ) / a )
  xn = xn * ( 1.0D+00 - ( a * log ( xn ) - xn - t ) / ( a - xn ) )
  go to 170
!
!  Schroder iteration using P.
!
  160 continue

  if ( 0.5D+00 < p ) then
    go to 220
  end if

  170 continue

  if ( p <= 1.0D+10 * tiny ( p ) ) then
    x = xn
    ierr = -8
    return
  end if

  am1 = ( a - half ) - half

  180 continue

  if ( amax < a ) then
    d = half + ( half - xn / a )
    if ( abs ( d ) <= e2 ) then
      x = xn
      ierr = -8
      return
    end if
  end if

  190 continue

  if ( 20 <= ierr ) then
    ierr = -6
    return
  end if

  ierr = ierr + 1
  call gamma_inc ( a, xn, pn, qn, 0 )

  if ( pn == 0.0D+00 .or. qn == 0.0D+00 ) then
    x = xn
    ierr = -8
    return
  end if

  r = rcomp ( a, xn )

  if ( r == 0.0D+00 ) then
    x = xn
    ierr = -8
    return
  end if

  t = ( pn - p ) / r
  w = half * ( am1 - xn )

  if ( abs ( t ) <= 0.1D+00 .and. abs ( w * t ) <= 0.1D+00 ) then
    go to 200
  end if

  x = xn * ( 1.0D+00 - t )

  if ( x <= 0.0D+00 ) then
    ierr = -7
    return
  end if

  d = abs ( t )
  go to 210

  200 continue

  h = t * ( 1.0D+00 + w * t )
  x = xn * ( 1.0D+00 - h )

  if ( x <= 0.0D+00 ) then
    ierr = -7
    return
  end if

  if ( 1.0D+00 <= abs ( w ) .and. abs ( w ) * t * t <= eps ) then
    return
  end if

  d = abs ( h )

  210 continue

  xn = x

  if ( d <= tol ) then

    if ( d <= eps ) then
      return
    end if

    if ( abs ( p - pn ) <= tol * p ) then
      return
    end if

  end if

  go to 180
!
!  Schroder iteration using Q.
!
  220 continue

  if ( q <= 1.0D+10 * tiny ( q ) ) then
    x = xn
    ierr = -8
    return
  end if

  am1 = ( a - half ) - half

  230 continue

  if ( amax < a ) then
    d = half + ( half - xn / a )
    if ( abs ( d ) <= e2 ) then
      x = xn
      ierr = -8
      return
    end if
  end if

  if ( 20 <= ierr ) then
    ierr = -6
    return
  end if

  ierr = ierr + 1
  call gamma_inc ( a, xn, pn, qn, 0 )

  if ( pn == 0.0D+00 .or. qn == 0.0D+00 ) then
    x = xn
    ierr = -8
    return
  end if

  r = rcomp ( a, xn )

  if ( r == 0.0D+00 ) then
    x = xn
    ierr = -8
    return
  end if

  t = ( q - qn ) / r
  w = half * ( am1 - xn )

  if ( abs ( t ) <= 0.1 .and. abs ( w * t ) <= 0.1 ) then
    go to 250
  end if

  x = xn * ( 1.0D+00 - t )

  if ( x <= 0.0D+00 ) then
    ierr = -7
    return
  end if

  d = abs ( t )
  go to 260

  250 continue

  h = t * ( 1.0D+00 + w * t )
  x = xn * ( 1.0D+00 - h )

  if ( x <= 0.0D+00 ) then
    ierr = -7
    return
  end if

  if ( 1.0D+00 <= abs ( w ) .and. abs ( w ) * t * t <= eps ) then
    return
  end if

  d = abs ( h )

  260 continue

  xn = x

  if ( tol < d ) then
    go to 230
  end if

  if ( d <= eps ) then
    return
  end if

  if ( abs ( q - qn ) <= tol * q ) then
    return
  end if

  go to 230
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
!
!    Mathematica can compute this value as
!
!      1 - GammaRegularized[A,X]
!
!  Modified:
!
!    08 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964.
!
!  Parameters:
!
!    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
!    first call.  On each call, the routine increments N_DATA by 1, and
!    returns the corresponding data; when there is no more data, the
!    output value of N_DATA will be 0 again.
!
!    Output, real ( kind = 8 ) A, X, the arguments of the function.
!
!    Output, real ( kind = 8 ) FX, the value of the function.
!
  implicit none

  integer, parameter :: n_max = 20

  real ( kind = 8 ) a
  real ( kind = 8 ), save, dimension ( n_max ) :: a_vec = (/ &
    0.1D+00,  0.1D+00,  0.1D+00,  0.5D+00, &
    0.5D+00,  0.5D+00,  1.0D+00,  1.0D+00, &
    1.0D+00,  1.1D+00,  1.1D+00,  1.1D+00, &
    2.0D+00,  2.0D+00,  2.0D+00,  6.0D+00, &
    6.0D+00, 11.0D+00, 26.0D+00, 41.0D+00 /)
  real ( kind = 8 ) fx
  real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ &
    0.7420263D+00, 0.9119753D+00, 0.9898955D+00, 0.2931279D+00, &
    0.7656418D+00, 0.9921661D+00, 0.0951626D+00, 0.6321206D+00, &
    0.9932621D+00, 0.0757471D+00, 0.6076457D+00, 0.9933425D+00, &
    0.0091054D+00, 0.4130643D+00, 0.9931450D+00, 0.0387318D+00, &
    0.9825937D+00, 0.9404267D+00, 0.4863866D+00, 0.7359709D+00 /)
  integer n_data
  real ( kind = 8 ) x
  real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ &
    3.1622777D-02, 3.1622777D-01, 1.5811388D+00, 7.0710678D-02, &
    7.0710678D-01, 3.5355339D+00, 0.1000000D+00, 1.0000000D+00, &
    5.0000000D+00, 1.0488088D-01, 1.0488088D+00, 5.2440442D+00, &
    1.4142136D-01, 1.4142136D+00, 7.0710678D+00, 2.4494897D+00, &
    1.2247449D+01, 1.6583124D+01, 2.5495098D+01, 4.4821870D+01 /)

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

  n_data = n_data + 1

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

  return
end
function gamma_ln1 ( a )

!*****************************************************************************80
!
!! GAMMA_LN1 evaluates ln ( Gamma ( 1 + A ) ), for -0.2 <= A <= 1.25.
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  Reference:
!
!    Armido DiDinato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, 1993, pages 360-373.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) A, defines the argument of the function.
!
!    Output, real ( kind = 8 ) GAMMA_LN1, the value of ln ( Gamma ( 1 + A ) ).
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) bot
  real ( kind = 8 ) gamma_ln1
  real ( kind = 8 ), parameter :: p0 =  0.577215664901533D+00
  real ( kind = 8 ), parameter :: p1 =  0.844203922187225D+00
  real ( kind = 8 ), parameter :: p2 = -0.168860593646662D+00
  real ( kind = 8 ), parameter :: p3 = -0.780427615533591D+00
  real ( kind = 8 ), parameter :: p4 = -0.402055799310489D+00
  real ( kind = 8 ), parameter :: p5 = -0.673562214325671D-01
  real ( kind = 8 ), parameter :: p6 = -0.271935708322958D-02
  real ( kind = 8 ), parameter :: q1 =  0.288743195473681D+01
  real ( kind = 8 ), parameter :: q2 =  0.312755088914843D+01
  real ( kind = 8 ), parameter :: q3 =  0.156875193295039D+01
  real ( kind = 8 ), parameter :: q4 =  0.361951990101499D+00
  real ( kind = 8 ), parameter :: q5 =  0.325038868253937D-01
  real ( kind = 8 ), parameter :: q6 =  0.667465618796164D-03
  real ( kind = 8 ), parameter :: r0 = 0.422784335098467D+00
  real ( kind = 8 ), parameter :: r1 = 0.848044614534529D+00
  real ( kind = 8 ), parameter :: r2 = 0.565221050691933D+00
  real ( kind = 8 ), parameter :: r3 = 0.156513060486551D+00
  real ( kind = 8 ), parameter :: r4 = 0.170502484022650D-01
  real ( kind = 8 ), parameter :: r5 = 0.497958207639485D-03
  real ( kind = 8 ), parameter :: s1 = 0.124313399877507D+01
  real ( kind = 8 ), parameter :: s2 = 0.548042109832463D+00
  real ( kind = 8 ), parameter :: s3 = 0.101552187439830D+00
  real ( kind = 8 ), parameter :: s4 = 0.713309612391000D-02
  real ( kind = 8 ), parameter :: s5 = 0.116165475989616D-03
  real ( kind = 8 ) top
  real ( kind = 8 ) x

  if ( a < 0.6D+00 ) then

    top = (((((  &
            p6   &
      * a + p5 ) &
      * a + p4 ) &
      * a + p3 ) &
      * a + p2 ) &
      * a + p1 ) &
      * a + p0  

    bot = (((((  &
            q6   &
      * a + q5 ) &
      * a + q4 ) &
      * a + q3 ) &
      * a + q2 ) &
      * a + q1 ) &
      * a + 1.0D+00

    gamma_ln1 = -a * ( top / bot )

  else

    x = ( a - 0.5D+00 ) - 0.5D+00

    top = ((((( r5 * x + r4 ) * x + r3 ) * x + r2 ) * x + r1 ) * x + r0 ) 

    bot = ((((( s5 * x + s4 ) * x + s3 ) * x + s2 ) * x + s1 ) * x + 1.0D+00 )

    gamma_ln1 = x * ( top / bot )

  end if

  return
end
function gamma_log ( a )

!*****************************************************************************80
!
!! GAMMA_LOG evaluates ln ( Gamma ( A ) ) for positive A.
!
!  Author:
!
!    Alfred Morris,
!    Naval Surface Weapons Center,
!    Dahlgren, Virginia.
!
!  Reference:
!
!    Armido DiDinato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, 1993, pages 360-373.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) A, the argument of the function.
!    A should be positive.
!
!    Output, real ( kind = 8 ), GAMMA_LOG, the value of ln ( Gamma ( A ) ).
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ), parameter :: c0 =  0.833333333333333D-01
  real ( kind = 8 ), parameter :: c1 = -0.277777777760991D-02
  real ( kind = 8 ), parameter :: c2 =  0.793650666825390D-03
  real ( kind = 8 ), parameter :: c3 = -0.595202931351870D-03
  real ( kind = 8 ), parameter :: c4 =  0.837308034031215D-03
  real ( kind = 8 ), parameter :: c5 = -0.165322962780713D-02
  real ( kind = 8 ), parameter :: d  =  0.418938533204673D+00
  real ( kind = 8 ) gamma_log
  real ( kind = 8 ) gamma_ln1
  integer i
  integer n
  real ( kind = 8 ) t
  real ( kind = 8 ) w

  if ( a <= 0.8D+00 ) then

    gamma_log = gamma_ln1 ( a ) - log ( a )

  else if ( a <= 2.25D+00 ) then

    t = ( a - 0.5D+00 ) - 0.5D+00
    gamma_log = gamma_ln1 ( t )

  else if ( a < 10.0D+00 ) then

    n = a - 1.25D+00
    t = a
    w = 1.0D+00
    do i = 1, n
      t = t - 1.0D+00
      w = t * w
    end do

    gamma_log = gamma_ln1 ( t - 1.0D+00 ) + log ( w )

  else

    t = ( 1.0D+00 / a )**2

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

    gamma_log = ( d + w ) + ( a - 0.5D+00 ) &
      * ( log ( a ) - 1.0D+00 )

  end if

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


!*****************************************************************************80
!
!! GAMMA_VALUES returns some values of the Gamma function.
!
!  Definition:
!
!    Gamma(Z) = Integral ( 0 <= T < Infinity) T**(Z-1) exp(-T) dT
!
!  Recursion:
!
!    Gamma(X+1) = X * Gamma(X)
!
!  Restrictions:
!
!    0 < X ( a software restriction).
!
!  Special values:
!
!    GAMMA(0.5) = sqrt(PI)
!
!    For N a positive integer, GAMMA(N+1) = N!, the standard factorial.
!
!  Modified:
!
!    17 April 2001
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964.
!
!  Parameters:
!
!    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
!    first call.  On each call, the routine increments N_DATA by 1, and
!    returns the corresponding data; when there is no more data, the
!    output value of N_DATA will be 0 again.
!
!    Output, real ( kind = 8 ) X, the argument of the function.
!
!    Output, real ( kind = 8 ) FX, the value of the function.
!
  implicit none

  integer, parameter :: n_max = 18

  real ( kind = 8 ) fx
  real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ &
    4.590845D+00,     2.218160D+00,     1.489192D+00,     1.164230D+00, &
    1.0000000000D+00, 0.9513507699D+00, 0.9181687424D+00, 0.8974706963D+00, &
    0.8872638175D+00, 0.8862269255D+00, 0.8935153493D+00, 0.9086387329D+00, &
    0.9313837710D+00, 0.9617658319D+00, 1.0000000000D+00, 3.6288000D+05, &
    1.2164510D+17,    8.8417620D+30 /)
  integer n_data
  real ( kind = 8 ) x
  real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ &
    0.2D+00,  0.4D+00,  0.6D+00,  0.8D+00, &
    1.0D+00,  1.1D+00,  1.2D+00,  1.3D+00, &
    1.4D+00,  1.5D+00,  1.6D+00,  1.7D+00, &
    1.8D+00,  1.9D+00,  2.0D+00, 10.0D+00, &
   20.0D+00, 30.0D+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
    x = 0.0D+00
    fx = 0.0D+00
  else
    x = x_vec(n_data)
    fx = fx_vec(n_data)
  end if

  return
end
function gsumln ( a, b )


!*****************************************************************************80
!
!! IPMPAR returns integer machine constants. 
!
!  Discussion:
!
!    Input arguments 1 through 3 are queries about integer arithmetic.
!    We assume integers are represented in the N-digit, base A form
!
!      sign * ( X(N-1)*A**(N-1) + ... + X(1)*A + X(0) )
!
!    where 0 <= X(0:N-1) < A.
!
!    Then:
!
!      IPMPAR(1) = A, the base of integer arithmetic;
!      IPMPAR(2) = N, the number of base A digits;
!      IPMPAR(3) = A**N - 1, the largest magnitude.
!
!    It is assumed that the single and real ( kind = 8 ) 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(1:M) is one of { 0, 1,..., B-1 }, and 1 <= X(1) and
!    EMIN <= E <= EMAX.
!
!    Input argument 4 is a query about the base of real arithmetic:
!
!      IPMPAR(4) = B, the base of single and real ( kind = 8 ) arithmetic.
!
!    Input arguments 5 through 7 are queries about single precision
!    floating point arithmetic:
!
!     IPMPAR(5) = M, the number of base B digits for single precision.
!     IPMPAR(6) = EMIN, the smallest exponent E for single precision.
!     IPMPAR(7) = EMAX, the largest exponent E for single precision.
!
!    Input arguments 8 through 10 are queries about real ( kind = 8 )
!    floating point arithmetic:
!
!     IPMPAR(8) = M, the number of base B digits for real ( kind = 8 ).
!     IPMPAR(9) = EMIN, the smallest exponent E for real ( kind = 8 ).
!     IPMPAR(10) = EMAX, the largest exponent E for real ( kind = 8 ).
!
!  Reference:
!
!    Phyllis Fox, Andrew Hall, Norman Schryer,
!    Algorithm 528:
!    Framework for a Portable FORTRAN Subroutine Library,
!    ACM Transactions on Mathematical Software,
!    Volume 4, 1978, pages 176-188.
!
!  Parameters:
!
!    Input, integer I, the index of the desired constant.
!
!    Output, integer IPMPAR, the value of the desired constant.
!
  implicit none

  integer i
  integer imach(10)
  integer ipmpar
!
!     MACHINE CONSTANTS FOR AMDAHL MACHINES.
!
!     data imach( 1) /   2 /
!     data imach( 2) /  31 /
!     data imach( 3) / 2147483647 /
!     data imach( 4) /  16 /
!     data imach( 5) /   6 /
!     data imach( 6) / -64 /
!     data imach( 7) /  63 /
!     data imach( 8) /  14 /
!     data imach( 9) / -64 /
!     data imach(10) /  63 /
!
!     Machine constants for the AT&T 3B SERIES, AT&T
!     PC 7300, AND AT&T 6300.
!
!     data imach( 1) /     2 /
!     data imach( 2) /    31 /
!     data imach( 3) / 2147483647 /
!     data imach( 4) /     2 /
!     data imach( 5) /    24 /
!     data imach( 6) /  -125 /
!     data imach( 7) /   128 /
!     data imach( 8) /    53 /
!     data imach( 9) / -1021 /
!     data imach(10) /  1024 /
!
!     Machine constants for the BURROUGHS 1700 SYSTEM.
!
!     data imach( 1) /    2 /
!     data imach( 2) /   33 /
!     data imach( 3) / 8589934591 /
!     data imach( 4) /    2 /
!     data imach( 5) /   24 /
!     data imach( 6) / -256 /
!     data imach( 7) /  255 /
!     data imach( 8) /   60 /
!     data imach( 9) / -256 /
!     data imach(10) /  255 /
!
!     Machine constants for the BURROUGHS 5700 SYSTEM.
!
!     data imach( 1) /    2 /
!     data imach( 2) /   39 /
!     data imach( 3) / 549755813887 /
!     data imach( 4) /    8 /
!     data imach( 5) /   13 /
!     data imach( 6) /  -50 /
!     data imach( 7) /   76 /
!     data imach( 8) /   26 /
!     data imach( 9) /  -50 /
!     data imach(10) /   76 /
!
!     Machine constants for the BURROUGHS 6700/7700 SYSTEMS.
!
!     data imach( 1) /      2 /
!     data imach( 2) /     39 /
!     data imach( 3) / 549755813887 /
!     data imach( 4) /      8 /
!     data imach( 5) /     13 /
!     data imach( 6) /    -50 /
!     data imach( 7) /     76 /
!     data imach( 8) /     26 /
!     data imach( 9) / -32754 /
!     data imach(10) /  32780 /
!
!     Machine constants for the CDC 6000/7000 SERIES
!     60 BIT ARITHMETIC, AND THE CDC CYBER 995 64 BIT
!     ARITHMETIC (NOS OPERATING SYSTEM).
!
!     data imach( 1) /    2 /
!     data imach( 2) /   48 /
!     data imach( 3) / 281474976710655 /
!     data imach( 4) /    2 /
!     data imach( 5) /   48 /
!     data imach( 6) / -974 /
!     data imach( 7) / 1070 /
!     data imach( 8) /   95 /
!     data imach( 9) / -926 /
!     data imach(10) / 1070 /
!
!     Machine constants for the CDC CYBER 995 64 BIT
!     ARITHMETIC (NOS/VE OPERATING SYSTEM).
!
!     data imach( 1) /     2 /
!     data imach( 2) /    63 /
!     data imach( 3) / 9223372036854775807 /
!     data imach( 4) /     2 /
!     data imach( 5) /    48 /
!     data imach( 6) / -4096 /
!     data imach( 7) /  4095 /
!     data imach( 8) /    96 /
!     data imach( 9) / -4096 /
!     data imach(10) /  4095 /
!
!     Machine constants for the CRAY 1, XMP, 2, AND 3.
!
!     data imach( 1) /     2 /
!     data imach( 2) /    63 /
!     data imach( 3) / 9223372036854775807 /
!     data imach( 4) /     2 /
!     data imach( 5) /    47 /
!     data imach( 6) / -8189 /
!     data imach( 7) /  8190 /
!     data imach( 8) /    94 /
!     data imach( 9) / -8099 /
!     data imach(10) /  8190 /
!
!     Machine constants for the data GENERAL ECLIPSE S/200.
!
!     data imach( 1) /    2 /
!     data imach( 2) /   15 /
!     data imach( 3) / 32767 /
!     data imach( 4) /   16 /
!     data imach( 5) /    6 /
!     data imach( 6) /  -64 /
!     data imach( 7) /   63 /
!     data imach( 8) /   14 /
!     data imach( 9) /  -64 /
!     data imach(10) /   63 /
!
!     Machine constants for the HARRIS 220.
!
!     data imach( 1) /    2 /
!     data imach( 2) /   23 /
!     data imach( 3) / 8388607 /
!     data imach( 4) /    2 /
!     data imach( 5) /   23 /
!     data imach( 6) / -127 /
!     data imach( 7) /  127 /
!     data imach( 8) /   38 /
!     data imach( 9) / -127 /
!     data imach(10) /  127 /
!
!     Machine constants for the HONEYWELL 600/6000
!     AND DPS 8/70 SERIES.
!
!     data imach( 1) /    2 /
!     data imach( 2) /   35 /
!     data imach( 3) / 34359738367 /
!     data imach( 4) /    2 /
!     data imach( 5) /   27 /
!     data imach( 6) / -127 /
!     data imach( 7) /  127 /
!     data imach( 8) /   63 /
!     data imach( 9) / -127 /
!     data imach(10) /  127 /
!
!     Machine constants for the HP 2100
!     3 WORD real ( kind = 8 ) 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 real ( kind = 8 ) 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-11 FORTRAN SUPPORTING
!     32-BIT integer ARITHMETIC.
!
!     data imach( 1) /    2 /
!     data imach( 2) /   31 /
!     data imach( 3) / 2147483647 /
!     data imach( 4) /    2 /
!     data imach( 5) /   24 /
!     data imach( 6) / -127 /
!     data imach( 7) /  127 /
!     data imach( 8) /   56 /
!     data imach( 9) / -127 /
!     data imach(10) /  127 /
!
!     Machine constants for the SEQUENT BALANCE 8000.
!
!     data imach( 1) /     2 /
!     data imach( 2) /    31 /
!     data imach( 3) / 2147483647 /
!     data imach( 4) /     2 /
!     data imach( 5) /    24 /
!     data imach( 6) /  -125 /
!     data imach( 7) /   128 /
!     data imach( 8) /    53 /
!     data imach( 9) / -1021 /
!     data imach(10) /  1024 /
!
!     Machine constants for the SILICON GRAPHICS IRIS-4D
!     SERIES (MIPS R3000 PROCESSOR).
!
!     data imach( 1) /     2 /
!     data imach( 2) /    31 /
!     data imach( 3) / 2147483647 /
!     data imach( 4) /     2 /
!     data imach( 5) /    24 /
!     data imach( 6) /  -125 /
!     data imach( 7) /   128 /
!     data imach( 8) /    53 /
!     data imach( 9) / -1021 /
!     data imach(10) /  1024 /
!
!     MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T
!     3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T
!     PC 7300), AND 8087 BASED MICROS (E.G. IBM PC 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 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
subroutine negative_binomial_cdf_values ( n_data, f, s, p, cdf )

!*****************************************************************************80
!
!! NEGATIVE_BINOMIAL_CDF_VALUES returns values of the negative binomial CDF.
!
!  Discussion:
!
!    Assume that a coin has a probability P of coming up heads on
!    any one trial.  Suppose that we plan to flip the coin until we
!    achieve a total of S heads.  If we let F represent the number of
!    tails that occur in this process, then the value of F satisfies
!    a negative binomial PDF:
!
!      PDF(F,S,P) = Choose ( F from F+S-1 ) * P**S * (1-P)**F
!
!    The negative binomial CDF is the probability that there are F or
!    fewer failures upon the attainment of the S-th success.  Thus,
!
!      CDF(F,S,P) = sum ( 0 <= G <= F ) PDF(G,S,P)
!
!  Modified:
!
!    07 June 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    FC Powell,
!    Statistical Tables for Sociology, Biology and Physical Sciences,
!    Cambridge University Press, 1982.
!
!  Parameters:
!
!    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
!    first call.  On each call, the routine increments N_DATA by 1, and
!    returns the corresponding data; when there is no more data, the
!    output value of N_DATA will be 0 again.
!
!    Output, integer F, the maximum number of failures.
!
!    Output, integer S, the number of successes.
!
!    Output, real ( kind = 8 ) P, the probability of a success on one trial.
!
!    Output, real ( kind = 8 ) CDF, the probability of at most F failures 
!    before the S-th success.
!
  implicit none

  integer, parameter :: n_max = 27

  real ( kind = 8 ) cdf
  real ( kind = 8 ), save, dimension ( n_max ) :: cdf_vec = (/ &
    0.6367D+00, 0.3633D+00, 0.1445D+00, &
    0.5000D+00, 0.2266D+00, 0.0625D+00, &
    0.3438D+00, 0.1094D+00, 0.0156D+00, &
    0.1792D+00, 0.0410D+00, 0.0041D+00, &
    0.0705D+00, 0.0109D+00, 0.0007D+00, &
    0.9862D+00, 0.9150D+00, 0.7472D+00, &
    0.8499D+00, 0.5497D+00, 0.2662D+00, &
    0.6513D+00, 0.2639D+00, 0.0702D+00, &
    1.0000D+00, 0.0199D+00, 0.0001D+00 /)
  integer f
  integer, save, dimension ( n_max ) :: f_vec = (/ &
     4,  3,  2, &
     3,  2,  1, &
     2,  1,  0, &
     2,  1,  0, &
     2,  1,  0, &
    11, 10,  9, &
    17, 16, 15, &
     9,  8,  7, &
     2,  1,  0 /)
  integer n_data
  real ( kind = 8 ) p
  real ( kind = 8 ), save, dimension ( n_max ) :: p_vec = (/ &
    0.50D+00, 0.50D+00, 0.50D+00, &
    0.50D+00, 0.50D+00, 0.50D+00, &
    0.50D+00, 0.50D+00, 0.50D+00, &
    0.40D+00, 0.40D+00, 0.40D+00, &
    0.30D+00, 0.30D+00, 0.30D+00, &
    0.30D+00, 0.30D+00, 0.30D+00, &
    0.10D+00, 0.10D+00, 0.10D+00, &
    0.10D+00, 0.10D+00, 0.10D+00, &
    0.01D+00, 0.01D+00, 0.01D+00 /)
  integer s
  integer, save, dimension ( n_max ) :: s_vec = (/ &
    4, 5, 6, &
    4, 5, 6, &
    4, 5, 6, &
    4, 5, 6, &
    4, 5, 6, &
    1, 2, 3, &
    1, 2, 3, &
    1, 2, 3, &
    0, 1, 2 /)

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

  n_data = n_data + 1

  if ( n_max < n_data ) then
    n_data = 0
    f = 0
    s = 0
    p = 0.0D+00
    cdf = 0.0D+00
  else
    f = f_vec(n_data)
    s = s_vec(n_data)
    p = p_vec(n_data)
    cdf = cdf_vec(n_data)
  end if

  return
end
subroutine normal_01_cdf_values ( n_data, x, fx )

!*****************************************************************************80
!
!! NORMAL_01_CDF_VALUES returns some values of the Normal 01 CDF.
!
!  Discussion:
!
!    In Mathematica, the function can be evaluated by:
!
!      Needs["Statistics`ContinuousDistributions`"]
!      dist = NormalDistribution [ 0, 1 ]
!      CDF [ dist, x ]
!
!  Modified:
!
!    28 August 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964.
!
!    Stephen Wolfram,
!    The Mathematica Book,
!    Fourth Edition,
!    Wolfram Media / Cambridge University Press, 1999.
!
!  Parameters:
!
!    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
!    first call.  On each call, the routine increments N_DATA by 1, and
!    returns the corresponding data; when there is no more data, the
!    output value of N_DATA will be 0 again.
!
!    Output, real ( kind = 8 ) X, the argument of the function.
!
!    Output, real ( kind = 8 ) FX, the value of the function.
!
  implicit none

  integer, parameter :: n_max = 17

  real ( kind = 8 ) fx
  real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ &
    0.5000000000000000D+00, &
    0.5398278372770290D+00, &
    0.5792597094391030D+00, &
    0.6179114221889526D+00, &
    0.6554217416103242D+00, &
    0.6914624612740131D+00, &
    0.7257468822499270D+00, &
    0.7580363477769270D+00, &
    0.7881446014166033D+00, &
    0.8159398746532405D+00, &
    0.8413447460685429D+00, &
    0.9331927987311419D+00, &
    0.9772498680518208D+00, &
    0.9937903346742239D+00, &
    0.9986501019683699D+00, &
    0.9997673709209645D+00, &
    0.9999683287581669D+00 /)
  integer n_data
  real ( kind = 8 ) x
  real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ &
    0.0000000000000000D+00, &
    0.1000000000000000D+00, &
    0.2000000000000000D+00, &
    0.3000000000000000D+00, &
    0.4000000000000000D+00, &
    0.5000000000000000D+00, &
    0.6000000000000000D+00, &
    0.7000000000000000D+00, &
    0.8000000000000000D+00, &
    0.9000000000000000D+00, &
    0.1000000000000000D+01, &
    0.1500000000000000D+01, &
    0.2000000000000000D+01, &
    0.2500000000000000D+01, &
    0.3000000000000000D+01, &
    0.3500000000000000D+01, &
    0.4000000000000000D+01 /)

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

  n_data = n_data + 1

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

  return
end
subroutine normal_cdf_values ( n_data, mu, sigma, x, fx )

!*****************************************************************************80
!
!! NORMAL_CDF_VALUES returns some values of the Normal CDF.
!
!  Discussion:
!
!    In Mathematica, the function can be evaluated by:
!
!      Needs["Statistics`ContinuousDistributions`"]
!      dist = NormalDistribution [ mu, sigma ]
!      CDF [ dist, x ]
!
!  Modified:
!
!    05 August 2004
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964.
!
!    Stephen Wolfram,
!    The Mathematica Book,
!    Fourth Edition,
!    Wolfram Media / Cambridge University Press, 1999.
!
!  Parameters:
!
!    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
!    first call.  On each call, the routine increments N_DATA by 1, and
!    returns the corresponding data; when there is no more data, the
!    output value of N_DATA will be 0 again.
!
!    Output, real ( kind = 8 ) MU, the mean of the distribution.
!
!    Output, real ( kind = 8 ) SIGMA, the variance of the distribution.
!
!    Output, real ( kind = 8 ) X, the argument of the function.
!
!    Output, real ( kind = 8 ) FX, the value of the function.
!
  implicit none

  integer, parameter :: n_max = 12

  real ( kind = 8 ) fx
  real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ &
    0.5000000000000000D+00, &
    0.9772498680518208D+00, &
    0.9999683287581669D+00, &
    0.9999999990134124D+00, &
    0.6914624612740131D+00, &
    0.6305586598182364D+00, &
    0.5987063256829237D+00, &
    0.5792597094391030D+00, &
    0.6914624612740131D+00, &
    0.5000000000000000D+00, &
    0.3085375387259869D+00, &
    0.1586552539314571D+00 /)
  real ( kind = 8 ) mu
  real ( kind = 8 ), save, dimension ( n_max ) :: mu_vec = (/ &
    0.1000000000000000D+01, &
    0.1000000000000000D+01, &
    0.1000000000000000D+01, &
    0.1000000000000000D+01, &
    0.1000000000000000D+01, &
    0.1000000000000000D+01, &
    0.1000000000000000D+01, &
    0.1000000000000000D+01, &
    0.2000000000000000D+01, &
    0.3000000000000000D+01, &
    0.4000000000000000D+01, &
    0.5000000000000000D+01 /)
  integer n_data
  real ( kind = 8 ) sigma
  real ( kind = 8 ), save, dimension ( n_max ) :: sigma_vec = (/ &
    0.5000000000000000D+00, &
    0.5000000000000000D+00, &
    0.5000000000000000D+00, &
    0.5000000000000000D+00, &
    0.2000000000000000D+01, &
    0.3000000000000000D+01, &
    0.4000000000000000D+01, &
    0.5000000000000000D+01, &
    0.2000000000000000D+01, &
    0.2000000000000000D+01, &
    0.2000000000000000D+01, &
    0.2000000000000000D+01 /)
  real ( kind = 8 ) x
  real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ &
    0.1000000000000000D+01, &
    0.2000000000000000D+01, &
    0.3000000000000000D+01, &
    0.4000000000000000D+01, &
    0.2000000000000000D+01, &
    0.2000000000000000D+01, &
    0.2000000000000000D+01, &
    0.2000000000000000D+01, &
    0.3000000000000000D+01, &
    0.3000000000000000D+01, &
    0.3000000000000000D+01, &
    0.3000000000000000D+01 /)

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

  n_data = n_data + 1

  if ( n_max < n_data ) then
    n_data = 0
    mu = 0.0D+00
    sigma = 0.0D+00
    x = 0.0D+00
    fx = 0.0D+00
  else
    mu = mu_vec(n_data)
    sigma = sigma_vec(n_data)
    x = x_vec(n_data)
    fx = fx_vec(n_data)
  end if

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

!*****************************************************************************80
!
!! POISSON_CDF_VALUES returns some values of the Poisson CDF.
!
!  Discussion:
!
!    CDF(X)(A) is the probability of at most X successes in unit time,
!    given that the expected mean number of successes is A.
!
!  Modified:
!
!    28 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964.
!
!    Daniel Zwillinger,
!    CRC Standard Mathematical Tables and Formulae,
!    30th Edition, CRC Press, 1996, pages 653-658.
!
!  Parameters:
!
!    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
!    first call.  On each call, the routine increments N_DATA by 1, and
!    returns the corresponding data; when there is no more data, the
!    output value of N_DATA will be 0 again.
!
!    Output, real ( kind = 8 ) A, integer X, the arguments of the function.
!
!    Output, real ( kind = 8 ) FX, the value of the function.
!
  implicit none

  integer, parameter :: n_max = 21

  real ( kind = 8 ) a
  real ( kind = 8 ), save, dimension ( n_max ) :: a_vec = (/ &
    0.02D+00, 0.10D+00, 0.10D+00, 0.50D+00, &
    0.50D+00, 0.50D+00, 1.00D+00, 1.00D+00, &
    1.00D+00, 1.00D+00, 2.00D+00, 2.00D+00, &
    2.00D+00, 2.00D+00, 5.00D+00, 5.00D+00, &
    5.00D+00, 5.00D+00, 5.00D+00, 5.00D+00, &
    5.00D+00 /)
  real ( kind = 8 ) fx
  real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ &
    0.980D+00, 0.905D+00, 0.995D+00, 0.607D+00, &
    0.910D+00, 0.986D+00, 0.368D+00, 0.736D+00, &
    0.920D+00, 0.981D+00, 0.135D+00, 0.406D+00, &
    0.677D+00, 0.857D+00, 0.007D+00, 0.040D+00, &
    0.125D+00, 0.265D+00, 0.441D+00, 0.616D+00, &
    0.762D+00 /)
  integer n_data
  integer x
  integer, save, dimension ( n_max ) :: x_vec = (/ &
     0, 0, 1, 0, &
     1, 2, 0, 1, &
     2, 3, 0, 1, &
     2, 3, 0, 1, &
     2, 3, 4, 5, &
     6 /)

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

  n_data = n_data + 1

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

  return
end
function psi ( xx )

!*****************************************************************************80
!
!! PSI evaluates the psi or digamma function, d/dx ln(gamma(x)).
!
!  Discussion:
!
!    The main computation involves evaluation of rational Chebyshev
!    approximations.  PSI was written at Argonne National Laboratory 
!    for FUNPACK, and subsequently modified by A. H. Morris of NSWC.
!
!  Reference:
!
!    William Cody, Anthony Strecok, Henry Thacher,
!    Chebyshev Approximations for the Psi Function,
!    Mathematics of Computation,
!    Volume 27, 1973, pages 123-127.
!
!    Armido DiDinato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, 1993, pages 360-373.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) XX, the argument of the psi function.
!
!    Output, real ( kind = 8 ) PSI, the value of the psi function.  PSI 
!    is assigned the value 0 when the psi function is undefined.
!
  implicit none

  real ( kind = 8 ) aug
  real ( kind = 8 ) den
  real ( kind = 8 ), parameter :: dx0 = &
    1.461632144968362341262659542325721325D+00
  integer i
  integer ipmpar
  integer m
  integer n
  integer nq
  real ( kind = 8 ), parameter, dimension ( 7 ) :: p1 = (/ &
   0.895385022981970D-02, &
   0.477762828042627D+01, &
   0.142441585084029D+03, &
   0.118645200713425D+04, &
   0.363351846806499D+04, &
   0.413810161269013D+04, &
   0.130560269827897D+04/)
  real ( kind = 8 ), dimension ( 4 ) :: p2 = (/ &
    -0.212940445131011D+01, &
    -0.701677227766759D+01, &
    -0.448616543918019D+01, &
    -0.648157123766197D+00 /)
  real ( kind = 8 ), parameter :: piov4 = 0.785398163397448D+00
  real ( kind = 8 ) psi
!
!  Coefficients for rational approximation of
!  PSI(X) / (X - X0),  0.5D+00 <= X <= 3.0D+00
!
  real ( kind = 8 ), dimension ( 6 ) :: q1 = (/ &
    0.448452573429826D+02, &
    0.520752771467162D+03, &
    0.221000799247830D+04, &
    0.364127349079381D+04, &
    0.190831076596300D+04, &
    0.691091682714533D-05 /)
  real ( kind = 8 ), dimension ( 4 ) :: q2 = (/ &
    0.322703493791143D+02, &
    0.892920700481861D+02, &
    0.546117738103215D+02, &
    0.777788548522962D+01 /)
  real ( kind = 8 ) sgn
  real ( kind = 8 ) upper
  real ( kind = 8 ) w
  real ( kind = 8 ) x
  real ( kind = 8 ) xmax1
  real ( kind = 8 ) xmx0
  real ( kind = 8 ) xsmall
  real ( kind = 8 ) xx
  real ( kind = 8 ) z
!
!  XMAX1 is the largest positive floating point constant with entirely 
!  integer representation.  It is 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).
!
  xmax1 = real ( ipmpar(3), kind = 8 )
  xmax1 = min ( xmax1, 1.0D+00 / epsilon ( xmax1 ) )
!
!  XSMALL is the absolute argument below which PI*COTAN(PI*X)
!  may be represented by 1/X.
!
  xsmall = 1.0D-09

  x = xx
  aug = 0.0D+00

  if ( x == 0.0D+00 ) then
    psi = 0.0D+00
    return
  end if
!
!  X < 0.5,  Use reflection formula PSI(1-X) = PSI(X) + PI * COTAN(PI*X)
!
  if ( x < 0.5D+00 ) then
!
!  0 < ABS ( X ) <= XSMALL.  Use 1/X as a substitute for PI*COTAN(PI*X)
!
    if ( abs ( x ) <= xsmall ) then
      aug = -1.0D+00 / x
      go to 40
    end if
!
!  Reduction of argument for cotangent.
!
    w = -x
    sgn = piov4

    if ( w <= 0.0D+00 ) then
      w = -w
      sgn = -sgn
    end if
!
!  Make an error exit if X <= -XMAX1
!
    if ( xmax1 <= w ) then
      psi = 0.0D+00
      return
    end if

    nq = int ( w )
    w = w - real ( nq, kind = 8 )
    nq = int ( w * 4.0D+00 )
    w = 4.0D+00 * ( w - real ( nq, kind = 8 ) * 0.25D+00 )
!
!  W is now related to the fractional part of 4.0D+00 * X.
!  Adjust argument to correspond to values in first
!  quadrant and determine sign.
!
    n = nq / 2
    if ( n + n /= nq ) then
      w = 1.0D+00 - 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 ) then

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

      aug = 4.0D+00 * sgn * ( cos(z) / sin(z) )

    else

      aug = 4.0D+00 * sgn * ( sin(z) / cos(z) )

    end if

   40   continue

    x = 1.0D+00 - x

  end if
!
!  0.5 <= X <= 3 
!
  if ( x <= 3.0D+00 ) then

    den = x
    upper = p1(1) * x

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

    den = ( upper + p1(7) ) / ( den + q1(6) )
    xmx0 = real ( x, kind = 8 ) - dx0
    psi = den * xmx0 + aug
!
!  3 < X < XMAX1
!
  else if ( x < xmax1 ) then

    w = 1.0D+00 / x**2
    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.5D+00 / x + aug
    psi = aug + log ( x )
!
!  XMAX1 <= X
!
  else

    psi = aug + log ( x )

  end if

  return
end
subroutine psi_values ( n_data, x, fx )

!*****************************************************************************80
!
!! PSI_VALUES returns some values of the Psi or Digamma function.
!
!  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
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964.
!
!  Parameters:
!
!    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
!    first call.  On each call, the routine increments N_DATA by 1, and
!    returns the corresponding data; when there is no more data, the
!    output value of N_DATA will be 0 again.
!
!    Output, real ( kind = 8 ) X, the argument of the function.
!
!    Output, real ( kind = 8 ) FX, the value of the function.
!
  implicit none

  integer, parameter :: n_max = 11

  real ( kind = 8 ) fx
  real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ &
    -0.5772156649D+00, -0.4237549404D+00, -0.2890398966D+00, &
    -0.1691908889D+00, -0.0613845446D+00, -0.0364899740D+00, &
     0.1260474528D+00,  0.2085478749D+00,  0.2849914333D+00, &
     0.3561841612D+00,  0.4227843351D+00 /)
  integer n_data
  real ( kind = 8 ) x
  real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ &
    1.0D+00,  1.1D+00,  1.2D+00,  &
    1.3D+00,  1.4D+00,  1.5D+00,  &
    1.6D+00,  1.7D+00,  1.8D+00,  &
    1.9D+00,  2.0D+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
    x = 0.0D+00
    fx = 0.0D+00
  else
    x = x_vec(n_data)
    fx = fx_vec(n_data)
  end if

  return
end
subroutine r8_swap ( x, y )

function rcomp ( a, x )

!*****************************************************************************80
!
!! REXP evaluates the function EXP(X) - 1.
!
!  Modified:
!
!    09 December 1999
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  Reference:
!
!    Armido DiDinato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, 1993, pages 360-373.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, the argument of the function.
!
!    Output, real ( kind = 8 ) REXP, the value of EXP(X)-1.
!
  implicit none

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

  if ( abs ( x ) <= 0.15D+00 ) then

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

  else

    w = exp ( x )

    if ( x <= 0.0D+00 ) then
      rexp = ( w - 0.5D+00 ) - 0.5D+00
    else
      rexp = w * ( 0.5D+00 + ( 0.5D+00 - 1.0D+00 / w ) )
    end if

  end if

  return
end
function rlog ( x )

!*****************************************************************************80
!
!! RLOG computes X - 1 - LN(X).
!
!  Modified:
!
!    06 August 2004
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  Reference:
!
!    Armido DiDinato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, 1993, pages 360-373.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, the argument of the function.
!
!    Output, real ( kind = 8 ) RLOG, the value of the function.
!
  implicit none

  real ( kind = 8 ), parameter :: a  =  0.566749439387324D-01
  real ( kind = 8 ), parameter :: b  =  0.456512608815524D-01
  real ( kind = 8 ), parameter :: half = 0.5D+00
  real ( kind = 8 ), parameter :: p0 =  0.333333333333333D+00
  real ( kind = 8 ), parameter :: p1 = -0.224696413112536D+00
  real ( kind = 8 ), parameter :: p2 =  0.620886815375787D-02
  real ( kind = 8 ), parameter :: q1 = -0.127408923933623D+01
  real ( kind = 8 ), parameter :: q2 =  0.354508718369557D+00
  real ( kind = 8 ) r
  real ( kind = 8 ) rlog
  real ( kind = 8 ) t
  real ( kind = 8 ), parameter :: two =  2.0D+00 
  real ( kind = 8 ) u
  real ( kind = 8 ) w
  real ( kind = 8 ) w1
  real ( kind = 8 ) x

  if ( x < 0.61D+00 ) then

    r = ( x - 0.5D+00 ) - 0.5D+00
    rlog = r - log ( x )

  else if ( x < 1.57D+00 ) then

    if ( x < 0.82D+00 ) then

      u = x - 0.7D+00
      u = u / 0.7D+00
      w1 = a - u * 0.3D+00

    else if ( x < 1.18D+00 ) then

      u = ( x - half ) - half
      w1 = 0.0D+00

    else if ( x < 1.57D+00 ) then

      u = 0.75D+00 * x - 1.0D+00
      w1 = b + u / 3.0D+00

    end if

    r = u / ( u + two )
    t = r * r
    w = ( ( p2 * t + p1 ) * t + p0 ) / ( ( q2 * t + q1 ) * t + 1.0D+00 )
    rlog = two * t * ( 1.0D+00 / ( 1.0D+00 - r ) - r * w ) + w1

  else if ( 1.57D+00 <= x ) then

    r = ( x - half ) - half
    rlog = r - log ( x )

  end if

  return
end
function rlog1 ( x )

!*****************************************************************************80
!
!! RLOG1 evaluates the function X - ln ( 1 + X ).
!
!  Author:
!
!    Armido DiDinato, Alfred Morris
!
!  Reference:
!
!    Armido DiDinato, Alfred Morris,
!    Algorithm 708: 
!    Significant Digit Computation of the Incomplete Beta Function Ratios,
!    ACM Transactions on Mathematical Software,
!    Volume 18, 1993, pages 360-373.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, the argument.
!
!    Output, real ( kind = 8 ) RLOG1, the value of X - ln ( 1 + X ).
!
  implicit none

  real ( kind = 8 ), parameter :: a = 0.566749439387324D-01
  real ( kind = 8 ), parameter :: b = 0.456512608815524D-01
  real ( kind = 8 ) h
  real ( kind = 8 ), parameter :: half = 0.5D+00
  real ( kind = 8 ), parameter :: p0 = 0.333333333333333D+00
  real ( kind = 8 ), parameter :: p1 = -0.224696413112536D+00
  real ( kind = 8 ), parameter :: p2 = 0.620886815375787D-02
  real ( kind = 8 ), parameter :: q1 = -0.127408923933623D+01
  real ( kind = 8 ), parameter :: q2 = 0.354508718369557D+00
  real ( kind = 8 ) r
  real ( kind = 8 ) rlog1
  real ( kind = 8 ) t
  real ( kind = 8 ), parameter :: two =  2.0D+00 
  real ( kind = 8 ) w
  real ( kind = 8 ) w1
  real ( kind = 8 ) x

  if ( x < -0.39D+00 ) then

    w = ( x + half ) + half
    rlog1 = x - log ( w )

  else if ( x < -0.18D+00 ) then

    h = x + 0.3D+00
    h = h / 0.7D+00
    w1 = a - h * 0.3D+00

    r = h / ( h + 2.0D+00 )
    t = r * r
    w = ( ( p2 * t + p1 ) * t + p0 ) / ( ( q2 * t + q1 ) * t + 1.0D+00 )
    rlog1 = two * t * ( 1.0D+00 / ( 1.0D+00 - r ) - r * w ) + w1

  else if ( x <= 0.18D+00 ) then

    h = x
    w1 = 0.0D+00

    r = h / ( h + two )
    t = r * r
    w = ( ( p2 * t + p1 ) * t + p0 ) / ( ( q2 * t + q1 ) * t + 1.0D+00 )
    rlog1 = two * t * ( 1.0D+00 / ( 1.0D+00 - r ) - r * w ) + w1

  else if ( x <= 0.57D+00 ) then

    h = 0.75D+00 * x - 0.25D+00
    w1 = b + h / 3.0D+00

    r = h / ( h + 2.0D+00 )
    t = r * r
    w = ( ( p2 * t + p1 ) * t + p0 ) / ( ( q2 * t + q1 ) * t + 1.0D+00 )
    rlog1 = two * t * ( 1.0D+00 / ( 1.0D+00 - r ) - r * w ) + w1

  else 

    w = ( x + half ) + half
    rlog1 = x - log ( w )

  end if

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

!*****************************************************************************80
!
!! STUDENT_CDF_VALUES returns some values of the Student CDF.
!
!  Modified:
!
!    02 June 2001
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Milton Abramowitz, Irene Stegun,
!    Handbook of Mathematical Functions,
!    US Department of Commerce, 1964.
!
!  Parameters:
!
!    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
!    first call.  On each call, the routine increments N_DATA by 1, and
!    returns the corresponding data; when there is no more data, the
!    output value of N_DATA will be 0 again.
!
!    Output, integer A, real ( kind = 8 ) X, the arguments of the function.
!
!    Output, real ( kind = 8 ) FX, the value of the function.
!
  implicit none

  integer, parameter :: n_max = 13

  integer a
  integer, save, dimension ( n_max ) :: a_vec = (/ &
    1, 2, 3, 4, &
    5, 2, 5, 2, &
    5, 2, 3, 4, &
    5 /)
  real ( kind = 8 ) fx
  real ( kind = 8 ), save, dimension ( n_max ) :: fx_vec = (/ &
    0.60D+00, 0.60D+00, 0.60D+00, 0.60D+00, &
    0.60D+00, 0.75D+00, 0.75D+00, 0.95D+00, &
    0.95D+00, 0.99D+00, 0.99D+00, 0.99D+00, &
    0.99D+00 /)
  integer n_data
  real ( kind = 8 ) x
  real ( kind = 8 ), save, dimension ( n_max ) :: x_vec = (/ &
    0.325D+00, 0.289D+00, 0.277D+00, 0.271D+00, &
    0.267D+00, 0.816D+00, 0.727D+00, 2.920D+00, &
    2.015D+00, 6.965D+00, 4.541D+00, 3.747D+00, &
    3.365D+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
    x = 0.0D+00
    fx = 0.0D+00
  else
    a = a_vec(n_data)
    x = x_vec(n_data)
    fx = fx_vec(n_data)
  end if

  return
end
function stvaln ( p )

!*****************************************************************************80
!
!! STVALN provides starting values for the inverse of the normal distribution.
!
!  Discussion:
!
!    The routine returns an X for which it is approximately true that 
!      P = CUMNOR(X),  
!    that is, 
!      P = Integral ( -infinity < U <= X ) exp(-U*U/2)/sqrt(2*PI) dU.
!
!  Reference:
!
!    William Kennedy, James Gentle,
!    Statistical Computing, 
!    Marcel Dekker, NY, 1980, page 95,
!    QA276.4 K46
!
!  Parameters:
!
!    Input, real ( kind = 8 ) P, the probability whose normal deviate 
!    is sought.
!
!    Output, real ( kind = 8 ) STVALN, the normal deviate whose probability
!    is approximately P.
!
  implicit none

  real ( kind = 8 ) eval_pol
  real ( kind = 8 ) p
  real ( kind = 8 ) sgn
  real ( kind = 8 ) stvaln
  real ( kind = 8 ), parameter, dimension(0:4) :: xden = (/ &
    0.993484626060D-01, &
    0.588581570495D+00, &
    0.531103462366D+00, &
    0.103537752850D+00, &
    0.38560700634D-02 /)
  real ( kind = 8 ), parameter, dimension(0:4) :: xnum = (/ &
    -0.322232431088D+00, &
    -1.000000000000D+00, &
    -0.342242088547D+00, &
    -0.204231210245D-01, &
    -0.453642210148D-04 /)
  real ( kind = 8 ) y
  real ( kind = 8 ) z

  if ( p <= 0.5D+00 ) then

    sgn = -1.0D+00
    z = p

  else

    sgn = 1.0D+00
    z = 1.0D+00 - p

  end if

  y = sqrt ( -2.0D+00 * log ( z ) )
  stvaln = y + eval_pol ( xnum, 4, y ) / eval_pol ( xden, 4, y )
  stvaln = sgn * stvaln

  return
end
subroutine timestamp ( )

C     ALGORITHM 654, COLLECTED ALGORITHMS FROM ACM.
C     THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C     VOL. 13, NO. 3, P. 318.
C     PROGRAM GTST(OUTPUT,TAPE6=OUTPUT)
C ----------------------------------------------------------------------
C     SAMPLE PROGRAM EMPLOYING GRATIO AND GAMINV. GIVEN A AND X,
C     P AND Q ARE COMPUTED BY GRATIO. THEN FOR A, THE INVERSE OF
C     P AND Q, DENOTED BY XN, IS OBTAINED BY GAMINV. D IS THE
C     RELATIVE DIFFERENCE BETWEEN X AND XN.
C
C     NO DATA IS READ. THE OUTPUT FOR THE PROGRAM IS WRITTEN ON
C     UNIT 6. THE FIRST STATMENT OF THIS TEXT MAY BE USED TO
C     BEGIN THE PROGRAM FOR THE CDC 6000-7000 SERIES COMPUTERS.
C
C     THE LAST FUNCTION IN THIS PACKAGE, SPMPAR, MUST BE DEFINED
C     FOR THE PARTICULAR COMPUTER BEING USED. FOR DETAILS SEE THE
C     IN-LINE DOCUMENTATION OF SPMPAR.
C ----------------------------------------------------------------------
      WRITE (6,1)
    1 FORMAT(11H1   A     X,10X,1HP,15X,1HQ,15X,2HXN,13X,
     *       10HD     IERR)
    2 FORMAT(2F6.2,3E16.6,E12.2,I5)
    3 FORMAT(1H0)
C
      A = 0.1
      X0 = 0.0
      DO 20 L = 1,10
      WRITE (6,3)
      X = 0.1
         DO 10 I = 1,10
         CALL GRATIO (A,X,P,Q,0)
         CALL GAMINV (A,XN,X0,P,Q,IERR)
         D = ABS((X - XN)/X)
         WRITE (6,2) A,X,P,Q,XN,D,IERR
   10    X = X + 0.1
   20 A = A + 0.1
      STOP
      END
      SUBROUTINE GRATIO (A, X, ANS, QANS, IND)
C ----------------------------------------------------------------------
C        EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS
C                      P(A,X) AND Q(A,X)
C
C                        ----------
C
C     IT IS ASSUMED THAT A AND X ARE NONNEGATIVE, WHERE A AND X
C     ARE NOT BOTH 0.
C
C     ANS AND QANS ARE VARIABLES. GRATIO ASSIGNS ANS THE VALUE
C     P(A,X) AND QANS THE VALUE Q(A,X). IND MAY BE ANY INTEGER.
C     IF IND = 0 THEN THE USER IS REQUESTING AS MUCH ACCURACY AS
C     POSSIBLE (UP TO 14 SIGNIFICANT DIGITS). OTHERWISE, IF
C     IND = 1 THEN ACCURACY IS REQUESTED TO WITHIN 1 UNIT OF THE
C     6-TH SIGNIFICANT DIGIT, AND IF IND .NE. 0,1 THEN ACCURACY
C     IS REQUESTED TO WITHIN 1 UNIT OF THE 3RD SIGNIFICANT DIGIT.
C
C     ERROR RETURN ...
C        ANS IS ASSIGNED THE VALUE 2 WHEN A OR X IS NEGATIVE,
C     WHEN A*X = 0, OR WHEN P(A,X) AND Q(A,X) ARE INDETERMINANT.
C     P(A,X) AND Q(A,X) ARE COMPUTATIONALLY INDETERMINANT WHEN
C     X IS EXCEEDINGLY CLOSE TO A AND A IS EXTREMELY LARGE.
C ----------------------------------------------------------------------
C     WRITTEN BY ALFRED H. MORRIS, JR.
C        NAVAL SURFACE WEAPONS CENTER
C        DAHLGREN, VIRGINIA
C     --------------------
      REAL J, L, ACC0(3), BIG(3), E00(3), X00(3), WK(20)
      REAL D0(13), D1(12), D2(10), D3(8), D4(6), D5(4), D6(2)
C     --------------------
      DATA ACC0(1)/5.E-15/, ACC0(2)/5.E-7/, ACC0(3)/5.E-4/
      DATA BIG(1)/20.0/, BIG(2)/14.0/, BIG(3)/10.0/
      DATA E00(1)/.25E-3/, E00(2)/.25E-1/, E00(3)/.14/
      DATA X00(1)/31.0/, X00(2)/17.0/, X00(3)/9.7/
C     --------------------
C     ALOG10 = LN(10)
C     RT2PIN = 1/SQRT(2*PI)
C     RTPI   = SQRT(PI)
C     --------------------
      DATA ALOG10/2.30258509299405/
      DATA RT2PIN/.398942280401433/
      DATA RTPI  /1.77245385090552/
      DATA THIRD /.333333333333333/
C     --------------------
      DATA D0(1) / .833333333333333E-01/, D0(2) /-.148148148148148E-01/,
     *     D0(3) / .115740740740741E-02/, D0(4) / .352733686067019E-03/,
     *     D0(5) /-.178755144032922E-03/, D0(6) / .391926317852244E-04/,
     *     D0(7) /-.218544851067999E-05/, D0(8) /-.185406221071516E-05/,
     *     D0(9) / .829671134095309E-06/, D0(10)/-.176659527368261E-06/,
     *     D0(11)/ .670785354340150E-08/, D0(12)/ .102618097842403E-07/,
     *     D0(13)/-.438203601845335E-08/
C     --------------------
      DATA D10   /-.185185185185185E-02/, D1(1) /-.347222222222222E-02/,
     *     D1(2) / .264550264550265E-02/, D1(3) /-.990226337448560E-03/,
     *     D1(4) / .205761316872428E-03/, D1(5) /-.401877572016461E-06/,
     *     D1(6) /-.180985503344900E-04/, D1(7) / .764916091608111E-05/,
     *     D1(8) /-.161209008945634E-05/, D1(9) / .464712780280743E-08/,
     *     D1(10)/ .137863344691572E-06/, D1(11)/-.575254560351770E-07/,
     *     D1(12)/ .119516285997781E-07/
C     --------------------
      DATA D20   / .413359788359788E-02/, D2(1) /-.268132716049383E-02/,
     *     D2(2) / .771604938271605E-03/, D2(3) / .200938786008230E-05/,
     *     D2(4) /-.107366532263652E-03/, D2(5) / .529234488291201E-04/,
     *     D2(6) /-.127606351886187E-04/, D2(7) / .342357873409614E-07/,
     *     D2(8) / .137219573090629E-05/, D2(9) /-.629899213838006E-06/,
     *     D2(10)/ .142806142060642E-06/
C     --------------------
      DATA D30   / .649434156378601E-03/, D3(1) / .229472093621399E-03/,
     *     D3(2) /-.469189494395256E-03/, D3(3) / .267720632062839E-03/,
     *     D3(4) /-.756180167188398E-04/, D3(5) /-.239650511386730E-06/,
     *     D3(6) / .110826541153473E-04/, D3(7) /-.567495282699160E-05/,
     *     D3(8) / .142309007324359E-05/
C     --------------------
      DATA D40   /-.861888290916712E-03/, D4(1) / .784039221720067E-03/,
     *     D4(2) /-.299072480303190E-03/, D4(3) /-.146384525788434E-05/,
     *     D4(4) / .664149821546512E-04/, D4(5) /-.396836504717943E-04/,
     *     D4(6) / .113757269706784E-04/
C     --------------------
      DATA D50   /-.336798553366358E-03/, D5(1) /-.697281375836586E-04/,
     *     D5(2) / .277275324495939E-03/, D5(3) /-.199325705161888E-03/,
     *     D5(4) / .679778047793721E-04/
C     --------------------
      DATA D60   / .531307936463992E-03/, D6(1) /-.592166437353694E-03/,
     *     D6(2) / .270878209671804E-03/
C     --------------------
      DATA D70   / .344367606892378E-03/
C     --------------------
C     ****** E IS A MACHINE DEPENDENT CONSTANT. E IS THE SMALLEST
C            FLOATING POINT NUMBER FOR WHICH 1.0 + E .GT. 1.0 .
C
                    E = SPMPAR(1)
C
C     --------------------
      IF (A .LT. 0.0 .OR. X .LT. 0.0) GO TO 400
      IF (A .EQ. 0.0 .AND. X .EQ. 0.0) GO TO 400
      IF (A*X .EQ. 0.0) GO TO 331
C
      IOP = IND + 1
      IF (IOP .NE. 1 .AND. IOP .NE. 2) IOP = 3
      ACC = AMAX1(ACC0(IOP),E)
      E0 = E00(IOP)
      X0 = X00(IOP)
C
C            SELECT THE APPROPRIATE ALGORITHM
C
      IF (A .GE. 1.0) GO TO 10
      IF (A .EQ. 0.5) GO TO 320
      IF (X .LT. 1.1) GO TO 110
      T1 = A*ALOG(X) - X
      U = A*EXP(T1)
      IF (U .EQ. 0.0) GO TO 310
      R = U*(1.0 + GAM1(A))
      GO TO 170
C
   10 IF (A .GE. BIG(IOP)) GO TO 20
      IF (A .GT. X .OR. X .GE. X0) GO TO 11
         TWOA = A + A
         M = INT(TWOA)
         IF (TWOA .NE. FLOAT(M)) GO TO 11
         I = M/2
         IF (A .EQ. FLOAT(I)) GO TO 140
         GO TO 150
   11 T1 = A*ALOG(X) - X
      R = EXP(T1)/GAMMA(A)
      GO TO 30
C
   20 L = X/A
      IF (L .EQ. 0.0) GO TO 300
      S = 0.5 + (0.5 - L)
      Z = RLOG(L)
      IF (Z .GE. 700.0/A) GO TO 330
      Y = A*Z
      RTA = SQRT(A)
      IF (ABS(S) .LE. E0/RTA) GO TO 250
      IF (ABS(S) .LE. 0.4) GO TO 200
C
      T = (1.0/A)**2
      T1 = (((0.75*T - 1.0)*T + 3.5)*T - 105.0)/(A*1260.0)
      T1 = T1 - Y
      R = RT2PIN*RTA*EXP(T1)
C
   30 IF (R .EQ. 0.0) GO TO 331
      IF (X .LE. AMAX1(A,ALOG10)) GO TO 50
      IF (X .LT. X0) GO TO 170
      GO TO 80
C
C                 TAYLOR SERIES FOR P/R
C
   50 APN = A + 1.0
      T = X/APN
      WK(1) = T
      DO 51 N = 2,20
         APN = APN + 1.0
         T = T*(X/APN)
         IF (T .LE. 1.E-3) GO TO 60
   51 WK(N) = T
      N = 20
C
   60 SUM = T
      TOL = 0.5*ACC
   61 APN = APN + 1.0
      T = T*(X/APN)
      SUM = SUM + T
      IF (T .GT. TOL) GO TO 61
C
      MAX = N - 1
      DO 70 M = 1,MAX
         N = N - 1
   70 SUM = SUM + WK(N)
      ANS = (R/A)*(1.0 + SUM)
      QANS = 0.5 + (0.5 - ANS)
      RETURN
C
C                 ASYMPTOTIC EXPANSION
C
   80 AMN = A - 1.0
      T = AMN/X
      WK(1) = T
      DO 81 N = 2,20
         AMN = AMN - 1.0
         T = T*(AMN/X)
         IF (ABS(T) .LE. 1.E-3) GO TO 90
   81 WK(N) = T
      N = 20
C
   90 SUM = T
   91 IF (ABS(T) .LE. ACC) GO TO 100
      AMN = AMN - 1.0
      T = T*(AMN/X)
      SUM = SUM + T
      GO TO 91
C
  100 MAX = N - 1
      DO 101 M = 1,MAX
         N = N - 1
  101 SUM = SUM + WK(N)
      QANS = (R/X)*(1.0 + SUM)
      ANS = 0.5 + (0.5 - QANS)
      RETURN
C
C             TAYLOR SERIES FOR P(A,X)/X**A
C
  110 AN = 3.0
      C = X
      SUM = X/(A + 3.0)
      TOL = 3.0*ACC/(A + 1.0)
  111    AN = AN + 1.0
         C = -C*(X/AN)
         T = C/(A + AN)
         SUM = SUM + T
         IF (ABS(T) .GT. TOL) GO TO 111
      J = A*X*((SUM/6.0 - 0.5/(A + 2.0))*X + 1.0/(A + 1.0))
C
      Z = A*ALOG(X)
      H = GAM1(A)
      G = 1.0 + H
      IF (X .LT. 0.25) GO TO 120
         IF (A .LT. X/2.59) GO TO 135
         GO TO 130
  120 IF (Z .GT. -.13394) GO TO 135
C
  130 W = EXP(Z)
      ANS = W*G*(0.5 + (0.5 - J))
      QANS = 0.5 + (0.5 - ANS)
      RETURN
C
  135 L = REXP(Z)
      W = 0.5 + (0.5 + L)
      QANS = (W*J - L)*G - H
      IF (QANS .LT. 0.0) GO TO 310
      ANS = 0.5 + (0.5 - QANS)
      RETURN
C
C             FINITE SUMS FOR Q WHEN A .GE. 1
C                 AND 2*A IS AN INTEGER
C
  140 SUM = EXP(-X)
      T = SUM
      N = 1
      C = 0.0
      GO TO 160
C
  150 RTX = SQRT(X)
      SUM = ERFC1(0,RTX)
      T = EXP(-X)/(RTPI*RTX)
      N = 0
      C = -0.5
C
  160 IF (N .EQ. I) GO TO 161
         N = N + 1
         C = C + 1.0
         T = (X*T)/C
         SUM = SUM + T
         GO TO 160
  161 QANS = SUM
      ANS = 0.5 + (0.5 - QANS)
      RETURN
C
C              CONTINUED FRACTION EXPANSION
C
  170 TOL = AMAX1(5.0*E,ACC)
      A2NM1 = 1.0
      A2N = 1.0
      B2NM1 = X
      B2N = X + (1.0 - A)
      C = 1.0
  171 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. TOL*AN0) GO TO 171
C
      QANS = R*AN0
      ANS = 0.5 + (0.5 - QANS)
      RETURN
C
C                GENERAL TEMME EXPANSION
C
  200 IF (ABS(S) .LE. 2.0*E .AND. A*E*E .GT. 3.28E-3) GO TO 400
      C = EXP(-Y)
      W = 0.5*ERFC1(1,SQRT(Y))
      U = 1.0/A
      Z = SQRT(Z + Z)
      IF (L .LT. 1.0) Z = -Z
      IF (IOP - 2) 210,220,230
C
  210 IF (ABS(S) .LE. 1.E-3) GO TO 260
      C0 = ((((((((((((D0(13) * Z + D0(12)) * Z + D0(11)) * Z
     *     + D0(10)) * Z + D0(9)) * Z + D0(8)) * Z + D0(7)) * Z
     *     + D0(6)) * Z + D0(5)) * Z + D0(4)) * Z + D0(3)) * Z
     *     + D0(2)) * Z + D0(1)) * Z - THIRD
      C1 = (((((((((((D1(12) * Z + D1(11)) * Z + D1(10)) * Z
     *     + D1(9)) * Z + D1(8)) * Z + D1(7)) * Z + D1(6)) * Z
     *     + D1(5)) * Z + D1(4)) * Z + D1(3)) * Z + D1(2)) * Z
     *     + D1(1)) * Z + D10
      C2 = (((((((((D2(10) * Z + D2(9)) * Z + D2(8)) * Z
     *     + D2(7)) * Z + D2(6)) * Z + D2(5)) * Z + D2(4)) * Z
     *     + D2(3)) * Z + D2(2)) * Z + D2(1)) * Z + D20
      C3 = (((((((D3(8) * Z + D3(7)) * Z + D3(6)) * Z
     *     + D3(5)) * Z + D3(4)) * Z + D3(3)) * Z + D3(2)) * Z
     *     + D3(1)) * Z + D30
      C4 = (((((D4(6) * Z + D4(5)) * Z + D4(4)) * Z + D4(3)) * Z
     *     + D4(2)) * Z + D4(1)) * Z + D40
      C5 = (((D5(4) * Z + D5(3)) * Z + D5(2)) * Z + D5(1)) * Z
     *     + D50
      C6 = (D6(2) * Z + D6(1)) * Z + D60
      T  = ((((((D70*U + C6)*U + C5)*U + C4)*U + C3)*U + C2)*U
     *                 + C1)*U + C0
      GO TO 240
C
  220 C0 = (((((D0(6) * Z + D0(5)) * Z + D0(4)) * Z + D0(3)) * Z
     *     + D0(2)) * Z + D0(1)) * Z - THIRD
      C1 = (((D1(4) * Z + D1(3)) * Z + D1(2)) * Z + D1(1)) * Z
     *     + D10
      C2 = D2(1) * Z + D20
      T  = (C2*U + C1)*U + C0
      GO TO 240
C
  230 T  = ((D0(3) * Z + D0(2)) * Z + D0(1)) * Z - THIRD
C
  240 IF (L .LT. 1.0) GO TO 241
      QANS = C*(W + RT2PIN*T/RTA)
      ANS = 0.5 + (0.5 - QANS)
      RETURN
  241 ANS = C*(W - RT2PIN*T/RTA)
      QANS = 0.5 + (0.5 - ANS)
      RETURN
C
C               TEMME EXPANSION FOR L = 1
C
  250 IF (A*E*E .GT. 3.28E-3) GO TO 400
      C = 0.5 + (0.5 - Y)
      W = (0.5 - SQRT(Y)*(0.5 + (0.5 - Y/3.0))/RTPI)/C
      U = 1.0/A
      Z = SQRT(Z + Z)
      IF (L .LT. 1.0) Z = -Z
      IF (IOP - 2) 260,270,280
C
  260 C0 = ((((((D0(7) * Z + D0(6)) * Z + D0(5)) * Z + D0(4)) * Z
     *     + D0(3)) * Z + D0(2)) * Z + D0(1)) * Z - THIRD
      C1 = (((((D1(6) * Z + D1(5)) * Z + D1(4)) * Z + D1(3)) * Z
     *     + D1(2)) * Z + D1(1)) * Z + D10
      C2 = ((((D2(5) * Z + D2(4)) * Z + D2(3)) * Z + D2(2)) * Z
     *     + D2(1)) * Z + D20
      C3 = (((D3(4) * Z + D3(3)) * Z + D3(2)) * Z + D3(1)) * Z
     *     + D30
      C4 = (D4(2) * Z + D4(1)) * Z + D40
      C5 = (D5(2) * Z + D5(1)) * Z + D50
      C6 = D6(1) * Z + D60
      T  = ((((((D70*U + C6)*U + C5)*U + C4)*U + C3)*U + C2)*U
     *                 + C1)*U + C0
      GO TO 240
C
  270 C0 = (D0(2) * Z + D0(1)) * Z - THIRD
      C1 = D1(1) * Z + D10
      T  = (D20*U + C1)*U + C0
      GO TO 240
C
  280 T  = D0(1) * Z - THIRD
      GO TO 240
C
C                     SPECIAL CASES
C
  300 ANS = 0.0
      QANS = 1.0
      RETURN
C
  310 ANS = 1.0
      QANS = 0.0
      RETURN
C
  320 IF (X .GE. 0.25) GO TO 321
      ANS = ERF(SQRT(X))
      QANS = 0.5 + (0.5 - ANS)
      RETURN
  321 QANS = ERFC1(0,SQRT(X))
      ANS = 0.5 + (0.5 - QANS)
      RETURN
C
  330 IF (ABS(S) .LE. 2.0*E) GO TO 400
  331 IF (X .LE. A) GO TO 300
      GO TO 310
C
C                     ERROR RETURN
C
  400 ANS = 2.0
      RETURN
      END
      SUBROUTINE GAMINV (A, X, X0, P, Q, IERR)
C ----------------------------------------------------------------------
C            INVERSE INCOMPLETE GAMMA RATIO FUNCTION
C
C     GIVEN POSITIVE A, AND NONEGATIVE P AND Q WHERE P + Q = 1.
C     THEN X IS COMPUTED WHERE P(A,X) = P AND Q(A,X) = Q. SCHRODER
C     ITERATION IS EMPLOYED. THE ROUTINE ATTEMPTS TO COMPUTE X
C     TO 10 SIGNIFICANT DIGITS IF THIS IS POSSIBLE FOR THE
C     PARTICULAR COMPUTER ARITHMETIC BEING USED.
C
C                      ------------
C
C     X IS A VARIABLE. IF P = 0 THEN X IS ASSIGNED THE VALUE 0,
C     AND IF Q = 0 THEN X IS SET TO THE LARGEST FLOATING POINT
C     NUMBER AVAILABLE. OTHERWISE, GAMINV ATTEMPTS TO OBTAIN
C     A SOLUTION FOR P(A,X) = P AND Q(A,X) = Q. IF THE ROUTINE
C     IS SUCCESSFUL THEN THE SOLUTION IS STORED IN X.
C
C     X0 IS AN OPTIONAL INITIAL APPROXIMATION FOR X. IF THE USER
C     DOES NOT WISH TO SUPPLY AN INITIAL APPROXIMATION, THEN SET
C     X0 .LE. 0.
C
C     IERR IS A VARIABLE THAT REPORTS THE STATUS OF THE RESULTS.
C     WHEN THE ROUTINE TERMINATES, IERR HAS ONE OF THE FOLLOWING
C     VALUES ...
C
C       IERR =  0    THE SOLUTION WAS OBTAINED. ITERATION WAS
C                    NOT USED.
C       IERR.GT.0    THE SOLUTION WAS OBTAINED. IERR ITERATIONS
C                    WERE PERFORMED.
C       IERR = -2    (INPUT ERROR) A .LE. 0
C       IERR = -3    NO SOLUTION WAS OBTAINED. THE RATIO Q/A
C                    IS TOO LARGE.
C       IERR = -4    (INPUT ERROR) P + Q .NE. 1
C       IERR = -6    20 ITERATIONS WERE PERFORMED. THE MOST
C                    RECENT VALUE OBTAINED FOR X IS GIVEN.
C                    THIS CANNOT OCCUR IF X0 .LE. 0.
C       IERR = -7    ITERATION FAILED. NO VALUE IS GIVEN FOR X.
C                    THIS MAY OCCUR WHEN X IS APPROXIMATELY 0.
C       IERR = -8    A VALUE FOR X HAS BEEN OBTAINED, BUT THE
C                    ROUTINE IS NOT CERTAIN OF ITS ACCURACY.
C                    ITERATION CANNOT BE PERFORMED IN THIS
C                    CASE. IF X0 .LE. 0, THIS CAN OCCUR ONLY
C                    WHEN P OR Q IS APPROXIMATELY 0. IF X0 IS
C                    POSITIVE THEN THIS CAN OCCUR WHEN A IS
C                    EXCEEDINGLY CLOSE TO X AND A IS EXTREMELY
C                    LARGE (SAY A .GE. 1.E20).
C ----------------------------------------------------------------------
C     WRITTEN BY ALFRED H. MORRIS, JR.
C        NAVAL SURFACE WEAPONS CENTER
C        DAHLGREN, VIRGINIA
C     -------------------
      REAL LN10, EPS0(2), AMIN(2), BMIN(2), DMIN(2), EMIN(2)
C     -------------------
C     LN10 = LN(10)
C     C = EULER CONSTANT
C     -------------------
      DATA LN10 /2.302585/
      DATA C  /.577215664901533/
C     -------------------
      DATA A0 /3.31125922108741/, A1 /11.6616720288968/,
     *     A2 /4.28342155967104/, A3 /.213623493715853/
      DATA B1 /6.61053765625462/, B2 /6.40691597760039/,
     *     B3 /1.27364489782223/, B4 /.036117081018842/
C     -------------------
      DATA EPS0(1) /1.E-10/, EPS0(2) /1.E-08/
      DATA AMIN(1) / 500.0/, AMIN(2) / 100.0/
      DATA BMIN(1) /1.E-28/, BMIN(2) /1.E-13/
      DATA DMIN(1) /1.E-06/, DMIN(2) /1.E-04/
      DATA EMIN(1) /2.E-03/, EMIN(2) /6.E-03/
C     -------------------
      DATA TOL /1.E-5/
C     -------------------
C     ****** E, XMIN, AND XMAX ARE MACHINE DEPENDENT CONSTANTS.
C            E IS THE SMALLEST NUMBER FOR WHICH 1.0 + E .GT. 1.0.
C            XMIN IS THE SMALLEST POSITIVE NUMBER AND XMAX IS THE
C            LARGEST POSITIVE NUMBER.
C
                   E = SPMPAR(1)
                   XMIN = SPMPAR(2)
                   XMAX = SPMPAR(3)
C     -------------------
      X = 0.0
      IF (A .LE. 0.0) GO TO 500
      T = DBLE(P) + DBLE(Q) - 1.D0
      IF (ABS(T) .GT. E) GO TO 520
C
      IERR = 0
      IF (P .EQ. 0.0) RETURN
      IF (Q .EQ. 0.0) GO TO 400
      IF (A .EQ. 1.0) GO TO 410
C
      E2 = 2.0*E
      AMAX = 0.4E-10/(E*E)
      IOP = 1
      IF (E .GT. 1.E-10) IOP = 2
      EPS = EPS0(IOP)
      XN = X0
      IF (X0 .GT. 0.0) GO TO 100
C
C        SELECTION OF THE INITIAL APPROXIMATION XN OF X
C                       WHEN A .LT. 1
C
      IF (A .GT. 1.0) GO TO 30
      G = GAMMA(A + 1.0)
      QG = Q*G
      IF (QG .EQ. 0.0) GO TO 560
      B = QG/A
      IF (QG .GT. 0.6*A) GO TO 20
      IF (A .GE. 0.30 .OR. B .LT. 0.35) GO TO 10
         T = EXP(-(B + C))
         U = T*EXP(T)
         XN = T*EXP(U)
         GO TO 100
C
   10 IF (B .GE. 0.45) GO TO 20
      IF (B .EQ. 0.0) GO TO 560
      Y = -ALOG(B)
      S = 0.5 + (0.5 - A)
      Z = ALOG(Y)
      T = Y - S*Z
      IF (B .LT. 0.15) GO TO 11
         XN = Y - S*ALOG(T) - ALOG(1.0 + S/(T + 1.0))
         GO TO 200
   11 IF (B .LE. 0.01) GO TO 12
         U = ((T + 2.0*(3.0 - A))*T + (2.0 - A)*(3.0 - A))/
     *           ((T + (5.0 - A))*T + 2.0)
         XN = Y - S*ALOG(T) - ALOG(U)
         GO TO 200
   12 C1 = -S*Z
      C2 = -S*(1.0 + C1)
      C3 =  S*((0.5*C1 + (2.0 - A))*C1 + (2.5 - 1.5*A))
      C4 = -S*(((C1/3.0 + (2.5 - 1.5*A))*C1 + ((A - 6.0)*A + 7.0))*C1
     *           + ((11.0*A - 46)*A + 47.0)/6.0)
      C5 = -S*((((-C1/4.0 + (11.0*A - 17.0)/6.0)*C1
     *           + ((-3.0*A + 13.0)*A - 13.0))*C1
     *           + 0.5*(((2.0*A - 25.0)*A + 72.0)*A - 61.0))*C1
     *           + (((25.0*A - 195.0)*A + 477.0)*A -379.0)/12.0)
      XN = ((((C5/Y + C4)/Y + C3)/Y + C2)/Y + C1) + Y
      IF (A .GT. 1.0) GO TO 200
      IF (B .GT. BMIN(IOP)) GO TO 200
      X = XN
      RETURN
C
   20 IF (B*Q .GT. 1.E-8) GO TO 21
         XN = EXP(-(Q/A + C))
         GO TO 23
   21 IF (P .LE. 0.9) GO TO 22
         XN = EXP((ALNREL(-Q) + GAMLN1(A))/A)
         GO TO 23
   22 XN = EXP(ALOG(P*G)/A)
   23 IF (XN .EQ. 0.0) GO TO 510
      T = 0.5 + (0.5 - XN/(A + 1.0))
      XN = XN/T
      GO TO 100
C
C        SELECTION OF THE INITIAL APPROXIMATION XN OF X
C                       WHEN A .GT. 1
C
   30 IF (Q .LE. 0.5) GO TO 31
         W = ALOG(P)
         GO TO 32
   31 W = ALOG(Q)
   32 T = SQRT(-2.0*W)
      S = T - (((A3*T + A2)*T + A1)*T + A0)/((((B4*T + B3)*T
     *                + B2)*T + B1)*T + 1.0)
      IF (Q .GT. 0.5) S = -S
C
      RTA = SQRT(A)
      S2 = S*S
      XN = A + S*RTA + (S2 - 1.0)/3.0 + S*(S2 - 7.0)/(36.0*RTA)
     1       - ((3.0*S2 + 7.0)*S2 - 16.0)/(810.0*A)
     2       + S*((9.0*S2 + 256.0)*S2 - 433.0)/(38880.0*A*RTA)
      XN = AMAX1(XN, 0.0)
      IF (A .LT. AMIN(IOP)) GO TO 40
      X = XN
      D = 0.5 + (0.5 - X/A)
      IF (ABS(D) .LE. DMIN(IOP)) RETURN
C
   40 IF (P .LE. 0.5) GO TO 50
      IF (XN .LT. 3.0*A) GO TO 200
      Y = -(W + GAMLN(A))
      D = AMAX1(2.0, A*(A - 1.0))
      IF (Y .LT. LN10*D) GO TO 41
         S = 1.0 - A
         Z = ALOG(Y)
         GO TO 12
   41 T = A - 1.0
      XN = Y + T*ALOG(XN) - ALNREL(-T/(XN + 1.0))
      XN = Y + T*ALOG(XN) - ALNREL(-T/(XN + 1.0))
      GO TO 200
C
   50 AP1 = A + 1.0
      IF (XN .GT. 0.70*AP1) GO TO 101
      W = W + GAMLN(AP1)
      IF (XN .GT. 0.15*AP1) GO TO 60
         AP2 = A + 2.0
         AP3 = A + 3.0
         X = EXP((W + X)/A)
         X = EXP((W + X - ALOG(1.0 + (X/AP1)*(1.0 + X/AP2)))/A)
         X = EXP((W + X - ALOG(1.0 + (X/AP1)*(1.0 + X/AP2)))/A)
         X = EXP((W + X - ALOG(1.0 + (X/AP1)*(1.0 + (X/AP2)*(1.0
     *                      + X/AP3))))/A)
         XN = X
         IF (XN .GT. 1.E-2*AP1) GO TO 60
         IF (XN .LE. EMIN(IOP)*AP1) RETURN
         GO TO 101
C
   60 APN = AP1
      T = XN/APN
      SUM = 1.0 + T
   61    APN = APN + 1.0
         T = T*(XN/APN)
         SUM = SUM + T
         IF (T .GT. 1.E-4) GO TO 61
      T = W - ALOG(SUM)
      XN = EXP((XN + T)/A)
      XN = XN*(1.0 - (A*ALOG(XN) - XN -T)/(A - XN))
      GO TO 101
C
C                 SCHRODER ITERATION USING P
C
  100 IF (P .GT. 0.5) GO TO 200
  101 IF (P .LE. 1.E10*XMIN) GO TO 550
      AM1 = (A - 0.5) - 0.5
  102 IF (A .LE. AMAX) GO TO 110
      D = 0.5 + (0.5 - XN/A)
      IF (ABS(D) .LE. E2) GO TO 550
C
  110 IF (IERR .GE. 20) GO TO 530
      IERR = IERR + 1
      CALL GRATIO(A,XN,PN,QN,0)
      IF (PN .EQ. 0.0 .OR. QN .EQ. 0.0) GO TO 550
      R = RCOMP(A,XN)
      IF (R .EQ. 0.0) GO TO 550
      T = (PN - P)/R
      W = 0.5*(AM1 - XN)
      IF (ABS(T) .LE. 0.1 .AND. ABS(W*T) .LE. 0.1) GO TO 120
         X = XN*(1.0 - T)
         IF (X .LE. 0.0) GO TO 540
         D = ABS(T)
         GO TO 121
C
  120 H = T*(1.0 + W*T)
      X = XN*(1.0 - H)
      IF (X .LE. 0.0) GO TO 540
      IF (ABS(W) .GE. 1.0 .AND. ABS(W)*T*T .LE. EPS) RETURN
      D = ABS(H)
  121 XN = X
      IF (D .GT. TOL) GO TO 102
      IF (D .LE. EPS) RETURN
      IF (ABS(P - PN) .LE. TOL*P) RETURN
      GO TO 102
C
C                 SCHRODER ITERATION USING Q
C
  200 IF (Q .LE. 1.E10*XMIN) GO TO 550
      AM1 = (A - 0.5) - 0.5
  201 IF (A .LE. AMAX) GO TO 210
      D = 0.5 + (0.5 - XN/A)
      IF (ABS(D) .LE. E2) GO TO 550
C
  210 IF (IERR .GE. 20) GO TO 530
      IERR = IERR + 1
      CALL GRATIO(A,XN,PN,QN,0)
      IF (PN .EQ. 0.0 .OR. QN .EQ. 0.0) GO TO 550
      R = RCOMP(A,XN)
      IF (R .EQ. 0.0) GO TO 550
      T = (Q - QN)/R
      W = 0.5*(AM1 - XN)
      IF (ABS(T) .LE. 0.1 .AND. ABS(W*T) .LE. 0.1) GO TO 220
         X = XN*(1.0 - T)
         IF (X .LE. 0.0) GO TO 540
         D = ABS(T)
         GO TO 221
C
  220 H = T*(1.0 + W*T)
      X = XN*(1.0 - H)
      IF (X .LE. 0.0) GO TO 540
      IF (ABS(W) .GE. 1.0 .AND. ABS(W)*T*T .LE. EPS) RETURN
      D = ABS(H)
  221 XN = X
      IF (D .GT. TOL) GO TO 201
      IF (D .LE. EPS) RETURN
      IF (ABS(Q - QN) .LE. TOL*Q) RETURN
      GO TO 201
C
C                       SPECIAL CASES
C
  400 X = XMAX
      RETURN
C
  410 IF (Q .LT. 0.9) GO TO 411
         X = -ALNREL(-P)
         RETURN
  411 X = -ALOG(Q)
      RETURN
C
C                       ERROR RETURN
C
  500 IERR = -2
      RETURN
C
  510 IERR = -3
      RETURN
C
  520 IERR = -4
      RETURN
C
  530 IERR = -6
      RETURN
C
  540 IERR = -7
      RETURN
C
  550 X = XN
      IERR = -8
      RETURN
C
  560 X = XMAX
      IERR = -8
      RETURN
      END
      FUNCTION ERF(X)
C     ******************************************************************
C     EVALUATION OF THE REAL ERROR FUNCTION
C     ******************************************************************
      DIMENSION A(4),B(4),P(8),Q(8),R(5),S(5)
      DATA A(1)/-1.65581836870402E-4/, A(2)/3.25324098357738E-2/,
     1     A(3)/1.02201136918406E-1/,  A(4)/1.12837916709552E00/
      DATA B(1)/4.64988945913179E-3/,  B(2)/7.01333417158511E-2/,
     1     B(3)/4.23906732683201E-1/,  B(4)/1.00000000000000E00/
      DATA P(1)/-1.36864857382717E-7/, P(2)/5.64195517478974E-1/,
     1     P(3)/7.21175825088309E00/,  P(4)/4.31622272220567E01/,
     2     P(5)/1.52989285046940E02/,  P(6)/3.39320816734344E02/,
     3     P(7)/4.51918953711873E02/,  P(8)/3.00459261020162E02/
      DATA Q(1)/1.00000000000000E00/,  Q(2)/1.27827273196294E01/,
     1     Q(3)/7.70001529352295E01/,  Q(4)/2.77585444743988E02/,
     2     Q(5)/6.38980264465631E02/,  Q(6)/9.31354094850610E02/,
     3     Q(7)/7.90950925327898E02/,  Q(8)/3.00459260956983E02/
      DATA R(1)/2.10144126479064E00/,  R(2)/2.62370141675169E01/,
     1     R(3)/2.13688200555087E01/,  R(4)/4.65807828718470E00/,
     2     R(5)/2.82094791773523E-1/
      DATA S(1)/9.41537750555460E01/,  S(2)/1.87114811799590E02/,
     1     S(3)/9.90191814623914E01/,  S(4)/1.80124575948747E01/,
     2     S(5)/1.00000000000000E00/
      DATA C/5.64189583547756E-1/
C     -------------------
      AX=ABS(X)
      X2=AX*AX
      IF (AX.GE.0.5) GO TO 10
      TOP=((A(1)*X2+A(2))*X2+A(3))*X2+A(4)
      BOT=((B(1)*X2+B(2))*X2+B(3))*X2+B(4)
      ERF=X*TOP/BOT
      RETURN
C
   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)
      ERF=1.0-EXP(-X2)*TOP/BOT
      IF (X.LT.0.0) ERF=-ERF
      RETURN
C
   20 ERF=1.0
      IF (AX.GE.5.54) GO TO 21
      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+S(5)
      ERF=C-TOP/(X2*BOT)
      ERF=1.0-EXP(-X2)*ERF/AX
   21 IF (X.LT.0.0) ERF=-ERF
      RETURN
      END
      REAL FUNCTION ERFC1(IND,X)
C ----------------------------------------------------------------------
C     EVALUATION OF THE REAL COMPLEMENTARY ERROR FUNCTION
C
C        ERFC1(IND,X) = ERFC(X)            IF IND = 0
C        ERFC1(IND,X) = EXP(X*X)*ERFC(X)   OTHERWISE
C ----------------------------------------------------------------------
      DIMENSION A(4),B(4),P(8),Q(8),R(5),S(5)
      DATA A(1)/-1.65581836870402E-4/, A(2)/3.25324098357738E-2/,
     1     A(3)/1.02201136918406E-1/,  A(4)/1.12837916709552E00/
      DATA B(1)/4.64988945913179E-3/,  B(2)/7.01333417158511E-2/,
     1     B(3)/4.23906732683201E-1/,  B(4)/1.00000000000000E00/
      DATA P(1)/-1.36864857382717E-7/, P(2)/5.64195517478974E-1/,
     1     P(3)/7.21175825088309E00/,  P(4)/4.31622272220567E01/,
     2     P(5)/1.52989285046940E02/,  P(6)/3.39320816734344E02/,
     3     P(7)/4.51918953711873E02/,  P(8)/3.00459261020162E02/
      DATA Q(1)/1.00000000000000E00/,  Q(2)/1.27827273196294E01/,
     1     Q(3)/7.70001529352295E01/,  Q(4)/2.77585444743988E02/,
     2     Q(5)/6.38980264465631E02/,  Q(6)/9.31354094850610E02/,
     3     Q(7)/7.90950925327898E02/,  Q(8)/3.00459260956983E02/
      DATA R(1)/2.10144126479064E00/,  R(2)/2.62370141675169E01/,
     1     R(3)/2.13688200555087E01/,  R(4)/4.65807828718470E00/,
     2     R(5)/2.82094791773523E-1/
      DATA S(1)/9.41537750555460E01/,  S(2)/1.87114811799590E02/,
     1     S(3)/9.90191814623914E01/,  S(4)/1.80124575948747E01/,
     2     S(5)/1.00000000000000E00/
      DATA C/5.64189583547756E-1/
C     -------------------
      AX=ABS(X)
      X2=AX*AX
      IF (AX.GE.0.47) GO TO 10
      TOP=((A(1)*X2+A(2))*X2+A(3))*X2+A(4)
      BOT=((B(1)*X2+B(2))*X2+B(3))*X2+B(4)
      ERFC1=1.0-X*TOP/BOT
      IF (IND.NE.0) ERFC1=EXP(X2)*ERFC1
      RETURN
C
   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
      IF (IND.EQ.0) GO TO 11
      IF (X.LT.0.0) ERFC1=2.0*EXP(X2)-ERFC1
      RETURN
   11 ERFC1=EXP(-X2)*ERFC1
      IF (X.LT.0.0) ERFC1=2.0-ERFC1
      RETURN
C
   20 IF (X.LE.-5.33) GO TO 30
      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+S(5)
      ERFC1=(C-TOP/(X2*BOT))/AX
      IF (IND.EQ.0) GO TO 11
      IF (X.LT.0.0) ERFC1=2.0*EXP(X2)-ERFC1
      RETURN
C
   30 ERFC1=2.0
      IF (IND.NE.0) ERFC1=EXP(X2)*ERFC1
      RETURN
      END
      REAL FUNCTION REXP(X)
C     ------------------------------------------------------------------
C     COMPUTATION OF EXP(X) - 1
C     ------------------------------------------------------------------
      DATA P1/ .914041914819518E-09/, P2/ .238082361044469E-01/,
     *     Q1/-.499999999085958E+00/, Q2/ .107141568980644E+00/,
     *     Q3/-.119041179760821E-01/, Q4/ .595130811860248E-03/
C     ------------------
      IF (ABS(X) .GT. 0.15) GO TO 10
      REXP = X*(((P2*X + P1)*X + 1.0)/((((Q4*X + Q3)*X + Q2)*X
     *                 + Q1)*X + 1.0))
      RETURN
C
   10 W = EXP(X)
      IF (X .GT. 0.0) GO TO 20
         REXP = (W - 0.5) - 0.5
         RETURN
   20 REXP = W*(0.5 + (0.5 - 1.0/W))
      RETURN
      END
      REAL FUNCTION ALNREL(A)
C     ------------------------------------------------------------------
C     EVALUATION OF THE FUNCTION LN(1 + A)
C     ------------------------------------------------------------------
      DATA P1/-.129418923021993E+01/, P2/.405303492862024E+00/,
     *     P3/-.178874546012214E-01/
      DATA Q1/-.162752256355323E+01/, Q2/.747811014037616E+00/,
     *     Q3/-.845104217945565E-01/
C     ---------------------
      IF (ABS(A) .GT. 0.375) GO TO 10
      T = A/(A + 2.0)
      T2 = T*T
      W = (((P3*T2 + P2)*T2 + P1)*T2 + 1.0)/
     *    (((Q3*T2 + Q2)*T2 + Q1)*T2 + 1.0)
      ALNREL = 2.0*T*W
      RETURN
C
   10 X = 1.D0 + DBLE(A)
      ALNREL = ALOG(X)
      RETURN
      END
      REAL FUNCTION RLOG(X)
C     -------------------
C     COMPUTATION OF  X - 1 - LN(X)
C     -------------------
      DATA A/.566749439387324E-01/
      DATA B/.456512608815524E-01/
C     -------------------
      DATA P0/ .333333333333333E+00/, P1/-.224696413112536E+00/,
     *     P2/ .620886815375787E-02/
      DATA Q1/-.127408923933623E+01/, Q2/ .354508718369557E+00/
C     -------------------
      IF (X .LT. 0.61 .OR. X .GT. 1.57) GO TO 100
      IF (X .LT. 0.82) GO TO 10
      IF (X .GT. 1.18) GO TO 20
C
C              ARGUMENT REDUCTION
C
      U = (X - 0.5) - 0.5
      W1 = 0.0
      GO TO 30
C
   10 U = DBLE(X) - 0.7D0
      U = U/0.7
      W1 = A - U*0.3
      GO TO 30
C
   20 U = 0.75D0*DBLE(X) - 1.D0
      W1 = B + U/3.0
C
C               SERIES EXPANSION
C
   30 R = U/(U + 2.0)
      T = R*R
      W = ((P2*T + P1)*T + P0)/((Q2*T + Q1)*T + 1.0)
      RLOG = 2.0*T*(1.0/(1.0 - R) - R*W) + W1
      RETURN
C
C
  100 R = (X - 0.5) - 0.5
      RLOG = R - ALOG(X)
      RETURN
      END
      REAL FUNCTION GAMMA(A)
C-----------------------------------------------------------------------
C
C         EVALUATION OF THE GAMMA FUNCTION FOR REAL ARGUMENTS
C
C                           -----------
C
C     GAMMA(A) IS ASSIGNED THE VALUE 0 WHEN THE GAMMA FUNCTION CANNOT
C     BE COMPUTED.
C
C-----------------------------------------------------------------------
C     WRITTEN BY ALFRED H. MORRIS, JR.
C          NAVAL SURFACE WEAPONS CENTER
C          DAHLGREN, VIRGINIA
C-----------------------------------------------------------------------
      REAL P(7), Q(7)
      DOUBLE PRECISION D, G, Z, LNX, GLOG
C--------------------------
C     D = 0.5*(LN(2*PI) - 1)
C--------------------------
      DATA PI /3.1415926535898/
      DATA D /.41893853320467274178D0/
C--------------------------
      DATA P(1)/ .539637273585445E-03/,  P(2)/ .261939260042690E-02/,
     1     P(3)/ .204493667594920E-01/,  P(4)/ .730981088720487E-01/,
     2     P(5)/ .279648642639792E+00/,  P(6)/ .553413866010467E+00/,
     3     P(7)/ 1.0/
      DATA Q(1)/-.832979206704073E-03/,  Q(2)/ .470059485860584E-02/,
     1     Q(3)/ .225211131035340E-01/,  Q(4)/-.170458969313360E+00/,
     2     Q(5)/-.567902761974940E-01/,  Q(6)/ .113062953091122E+01/,
     3     Q(7)/ 1.0/
C--------------------------
      DATA R1/.820756370353826E-03/, R2/-.595156336428591E-03/,
     1     R3/.793650663183693E-03/, R4/-.277777777770481E-02/,
     2     R5/.833333333333333E-01/
C--------------------------
      GAMMA = 0.0
      X = A
      IF (ABS(A) .GE. 15.0) GO TO 60
C-----------------------------------------------------------------------
C            EVALUATION OF GAMMA(A) FOR ABS(A) .LT. 15
C-----------------------------------------------------------------------
      T = 1.0
      M = INT(A) - 1
C
C     LET T BE THE PRODUCT OF A-J WHEN A .GE. 2
C
      IF (M) 20,12,10
   10 DO 11 J = 1,M
        X = X - 1.0
   11   T = X*T
   12 X = X - 1.0
      GO TO 40
C
C     LET T BE THE PRODUCT OF A+J WHEN A .LT. 1
C
   20 T = A
      IF (A .GT. 0.0) GO TO 30
      M = - M - 1
      IF (M .EQ. 0) GO TO 22
         DO 21 J = 1,M
         X = X + 1.0
   21    T = X*T
   22 X = (X + 0.5) + 0.5
      T = X*T
      IF (T .EQ. 0.0) RETURN
C
   30 CONTINUE
C
C     THE FOLLOWING CODE CHECKS IF 1/T CAN OVERFLOW. THIS
C     CODE MAY BE OMITTED IF DESIRED.
C
      IF (ABS(T) .GE. 1.E-30) GO TO 40
      IF (ABS(T)*SPMPAR(3) .LE. 1.0001) RETURN
      GAMMA = 1.0/T
      RETURN
C
C     COMPUTE GAMMA(1 + X) FOR  0 .LE. X .LT. 1
C
   40 TOP = P(1)
      BOT = Q(1)
      DO 41 I = 2,7
         TOP = P(I) + X*TOP
   41    BOT = Q(I) + X*BOT
      GAMMA = TOP/BOT
C
C     TERMINATION
C
      IF (A .LT. 1.0) GO TO 50
      GAMMA = GAMMA*T
      RETURN
   50 GAMMA = GAMMA/T
      RETURN
C-----------------------------------------------------------------------
C            EVALUATION OF GAMMA(A) FOR ABS(A) .GE. 15
C-----------------------------------------------------------------------
   60 IF (ABS(A) .GE. 1.E3) RETURN
      IF (A .GT. 0.0) GO TO 70
      X = -A
      N = X
      T = X - N
      IF (T .GT. 0.9) T = 1.0 - T
      S = SIN(PI*T)/PI
      IF (MOD(N,2) .EQ. 0) S = -S
      IF (S .EQ. 0.0) RETURN
C
C     COMPUTE THE MODIFIED ASYMPTOTIC SUM
C
   70 T = 1.0/(X*X)
      G = ((((R1*T + R2)*T + R3)*T + R4)*T + R5)/X
C
C     ONE MAY REPLACE THE NEXT STATEMENT WITH  LNX = ALOG(X)
C     BUT LESS ACCURACY WILL NORMALLY BE OBTAINED.
C
      LNX = GLOG(X)
C
C     FINAL ASSEMBLY
C
      Z = X
      G = (D + G) + (Z - 0.5D0)*(LNX - 1.D0)
      W = G
      T = G - DBLE(W)
      IF (W .GT. 0.99999*EXPARG(0)) RETURN
      GAMMA = EXP(W)*(1.0 + T)
      IF (A .LT. 0.0) GAMMA = (1.0/(GAMMA*S))/X
      RETURN
      END
      DOUBLE PRECISION FUNCTION GLOG(X)
C     -------------------
C     EVALUATION OF LN(X) FOR X .GE. 15
C     -------------------
      REAL X
      DOUBLE PRECISION Z, W(163)
C     -------------------
      DATA C1/.286228750476730/, C2/.399999628131494/,
     1     C3/.666666666752663/
C     -------------------
C     W(J) = LN(J + 14) FOR EACH J
C     -------------------
      DATA W(1) /.270805020110221007D+01/,
     1 W(2) /.277258872223978124D+01/, W(3) /.283321334405621608D+01/,
     2 W(4) /.289037175789616469D+01/, W(5) /.294443897916644046D+01/,
     3 W(6) /.299573227355399099D+01/, W(7) /.304452243772342300D+01/,
     4 W(8) /.309104245335831585D+01/, W(9) /.313549421592914969D+01/,
     5 W(10)/.317805383034794562D+01/, W(11)/.321887582486820075D+01/,
     6 W(12)/.325809653802148205D+01/, W(13)/.329583686600432907D+01/,
     7 W(14)/.333220451017520392D+01/, W(15)/.336729582998647403D+01/,
     8 W(16)/.340119738166215538D+01/, W(17)/.343398720448514625D+01/,
     9 W(18)/.346573590279972655D+01/, W(19)/.349650756146648024D+01/,
     1 W(20)/.352636052461616139D+01/, W(21)/.355534806148941368D+01/,
     2 W(22)/.358351893845611000D+01/, W(23)/.361091791264422444D+01/,
     3 W(24)/.363758615972638577D+01/, W(25)/.366356164612964643D+01/,
     4 W(26)/.368887945411393630D+01/, W(27)/.371357206670430780D+01/,
     5 W(28)/.373766961828336831D+01/, W(29)/.376120011569356242D+01/,
     6 W(30)/.378418963391826116D+01/
      DATA W(31)/.380666248977031976D+01/,
     1 W(32)/.382864139648909500D+01/, W(33)/.385014760171005859D+01/,
     2 W(34)/.387120101090789093D+01/, W(35)/.389182029811062661D+01/,
     3 W(36)/.391202300542814606D+01/, W(37)/.393182563272432577D+01/,
     4 W(38)/.395124371858142735D+01/, W(39)/.397029191355212183D+01/,
     5 W(40)/.398898404656427438D+01/, W(41)/.400733318523247092D+01/,
     6 W(42)/.402535169073514923D+01/, W(43)/.404305126783455015D+01/,
     7 W(44)/.406044301054641934D+01/, W(45)/.407753744390571945D+01/,
     8 W(46)/.409434456222210068D+01/, W(47)/.411087386417331125D+01/,
     9 W(48)/.412713438504509156D+01/, W(49)/.414313472639153269D+01/,
     1 W(50)/.415888308335967186D+01/, W(51)/.417438726989563711D+01/,
     2 W(52)/.418965474202642554D+01/, W(53)/.420469261939096606D+01/,
     3 W(54)/.421950770517610670D+01/, W(55)/.423410650459725938D+01/,
     4 W(56)/.424849524204935899D+01/, W(57)/.426267987704131542D+01/,
     5 W(58)/.427666611901605531D+01/, W(59)/.429045944114839113D+01/,
     6 W(60)/.430406509320416975D+01/
      DATA W(61)/.431748811353631044D+01/,
     1 W(62)/.433073334028633108D+01/, W(63)/.434380542185368385D+01/,
     2 W(64)/.435670882668959174D+01/, W(65)/.436944785246702149D+01/,
     3 W(66)/.438202663467388161D+01/, W(67)/.439444915467243877D+01/,
     4 W(68)/.440671924726425311D+01/, W(69)/.441884060779659792D+01/,
     5 W(70)/.443081679884331362D+01/, W(71)/.444265125649031645D+01/,
     6 W(72)/.445434729625350773D+01/, W(73)/.446590811865458372D+01/,
     7 W(74)/.447733681447820647D+01/, W(75)/.448863636973213984D+01/,
     8 W(76)/.449980967033026507D+01/, W(77)/.451085950651685004D+01/,
     9 W(78)/.452178857704904031D+01/, W(79)/.453259949315325594D+01/,
     1 W(80)/.454329478227000390D+01/, W(81)/.455387689160054083D+01/,
     2 W(82)/.456434819146783624D+01/, W(83)/.457471097850338282D+01/,
     3 W(84)/.458496747867057192D+01/, W(85)/.459511985013458993D+01/,
     4 W(86)/.460517018598809137D+01/, W(87)/.461512051684125945D+01/,
     5 W(88)/.462497281328427108D+01/, W(89)/.463472898822963577D+01/,
     6 W(90)/.464439089914137266D+01/
      DATA W(91) /.465396035015752337D+01/,
     1 W(92) /.466343909411206714D+01/, W(93) /.467282883446190617D+01/,
     2 W(94) /.468213122712421969D+01/, W(95) /.469134788222914370D+01/,
     3 W(96) /.470048036579241623D+01/, W(97) /.470953020131233414D+01/,
     4 W(98) /.471849887129509454D+01/, W(99) /.472738781871234057D+01/,
     5 W(100)/.473619844839449546D+01/, W(101)/.474493212836325007D+01/,
     6 W(102)/.475359019110636465D+01/, W(103)/.476217393479775612D+01/,
     7 W(104)/.477068462446566476D+01/, W(105)/.477912349311152939D+01/,
     8 W(106)/.478749174278204599D+01/, W(107)/.479579054559674109D+01/,
     9 W(108)/.480402104473325656D+01/, W(109)/.481218435537241750D+01/,
     1 W(110)/.482028156560503686D+01/, W(111)/.482831373730230112D+01/,
     2 W(112)/.483628190695147800D+01/, W(113)/.484418708645859127D+01/,
     3 W(114)/.485203026391961717D+01/, W(115)/.485981240436167211D+01/,
     4 W(116)/.486753445045558242D+01/, W(117)/.487519732320115154D+01/,
     5 W(118)/.488280192258637085D+01/, W(119)/.489034912822175377D+01/,
     6 W(120)/.489783979995091137D+01/
      DATA W(121)/.490527477843842945D+01/,
     1 W(122)/.491265488573605201D+01/, W(123)/.491998092582812492D+01/,
     2 W(124)/.492725368515720469D+01/, W(125)/.493447393313069176D+01/,
     3 W(126)/.494164242260930430D+01/, W(127)/.494875989037816828D+01/,
     4 W(128)/.495582705760126073D+01/, W(129)/.496284463025990728D+01/,
     5 W(130)/.496981329957600062D+01/, W(131)/.497673374242057440D+01/,
     6 W(132)/.498360662170833644D+01/, W(133)/.499043258677873630D+01/,
     7 W(134)/.499721227376411506D+01/, W(135)/.500394630594545914D+01/,
     8 W(136)/.501063529409625575D+01/, W(137)/.501727983681492433D+01/,
     9 W(138)/.502388052084627639D+01/, W(139)/.503043792139243546D+01/,
     1 W(140)/.503695260241362916D+01/, W(141)/.504342511691924662D+01/,
     2 W(142)/.504985600724953705D+01/, W(143)/.505624580534830806D+01/,
     3 W(144)/.506259503302696680D+01/, W(145)/.506890420222023153D+01/,
     4 W(146)/.507517381523382692D+01/, W(147)/.508140436498446300D+01/,
     5 W(148)/.508759633523238407D+01/, W(149)/.509375020080676233D+01/,
     6 W(150)/.509986642782419842D+01/
      DATA W(151)/.510594547390058061D+01/,
     1 W(152)/.511198778835654323D+01/, W(153)/.511799381241675511D+01/,
     2 W(154)/.512396397940325892D+01/, W(155)/.512989871492307347D+01/,
     3 W(156)/.513579843705026176D+01/, W(157)/.514166355650265984D+01/,
     4 W(158)/.514749447681345304D+01/, W(159)/.515329159449777895D+01/,
     5 W(160)/.515905529921452903D+01/, W(161)/.516478597392351405D+01/,
     6 W(162)/.517048399503815178D+01/, W(163)/.517614973257382914D+01/
C
      IF (X .GE. 178.0) GO TO 10
      N = X
      T = (X - N)/(X + N)
      T2 = T*T
      Z = (((C1*T2 + C2)*T2 + C3)*T2 + 2.0)*T
      GLOG = W(N - 14) + Z
      RETURN
C
   10 GLOG = ALOG(X)
      RETURN
      END
      REAL FUNCTION EXPARG (IDUMMY)
C--------------------------------------------------------------------
C     COMPUTATION OF THE LARGEST ARGUMENT W FOR WHICH EXP(W)
C     MAY BE COMPUTED. (ONLY AN APPROXIMATE VALUE IS NEEDED.)
C--------------------------------------------------------------------
      EXPARG = 0.99999*ALOG(SPMPAR(3))
      RETURN
      END
      REAL FUNCTION GAM1(A)
C     ------------------------------------------------------------------
C     COMPUTATION OF 1/GAMMA(A+1) - 1  FOR -0.5 .LE. A .LE. 1.5
C     ------------------------------------------------------------------
      REAL P(7), Q(5), R(9)
C     -------------------
      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/
C     -------------------
      DATA Q(1)/ .100000000000000E+01/, Q(2)/ .427569613095214E+00/,
     *     Q(3)/ .158451672430138E+00/, Q(4)/ .261132021441447E-01/,
     *     Q(5)/ .423244297896961E-02/
C     -------------------
      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/
C     -------------------
      DATA S1  / .273076135303957E+00/, S2  / .559398236957378E-01/
C     -------------------
      T = A
      D = A - 0.5
      IF (D .GT. 0.0) T = D - 0.5
      IF (T) 30,10,20
C
   10 GAM1 = 0.0
      RETURN
C
   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
C
   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 (D .GT. 0.0) GO TO 31
         GAM1 = A*((W + 0.5) + 0.5)
         RETURN
   31 GAM1 = T*W/A
      RETURN
      END
      REAL FUNCTION GAMLN(A)
C     ******************************************************************
C     EVALUATION OF LN(GAMMA(A)) FOR POSITIVE A
C     ******************************************************************
C     WRITTEN BY ALFRED H. MORRIS
C          NAVAL SURFACE WEAPONS CENTER
C          DAHLGREN, VIRGINIA
C     ---------------------
C     D = 0.5*(LN(2*PI) - 1)
C     ---------------------
      DATA D/.418938533204673/
C     ---------------------
      DATA C0/.833333333333333E-01/, C1/-.277777777770481E-02/,
     1     C2/.793650663183693E-03/, C3/-.595156336428591E-03/,
     2     C4/.820756370353826E-03/
C     ------------------------------------------------------------------
      IF (A .GT. 0.8) GO TO 10
         GAMLN = GAMLN1(A) - ALOG(A)
         RETURN
   10 IF (A .GT. 2.25) GO TO 20
         X = DBLE(A) - 1.D0
         GAMLN = GAMLN1(X)
         RETURN
C
   20 IF (A .GE. 15.0) GO TO 30
      N = A - 1.25
      X = A
      W = 1.0
      DO 21 I = 1,N
         X = X - 1.0
   21    W = X*W
      GAMLN = GAMLN1(X - 1.0) + ALOG(W)
      RETURN
C
   30 X = (1.0/A)**2
      W = ((((C4*X + C3)*X + C2)*X + C1)*X + C0)/A
      GAMLN = (D + W) + (A - 0.5)*(ALOG(A) - 1.0)
      END
      REAL FUNCTION GAMLN1(A)
C     ------------------------------------------------------------------
C     EVALUATION OF LN(GAMMA(1 + A)) FOR -0.2 .LE. A .LE. 1.25
C     ------------------------------------------------------------------
      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/
C     -----------------
      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/
C     -----------------
      IF (A .GE. 0.6) GO TO 10
      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
      RETURN
C
   10 X = DBLE(A) - 1.D0
      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
      RETURN
      END
      REAL FUNCTION RCOMP(A,X)
C     -------------------
C     EVALUATION OF EXP(-X)*X**A/GAMMA(A)
C     -------------------
C     RT2PIN = 1/SQRT(2*PI)
C     -------------------
      DATA RT2PIN/.398942280401433/
C     -------------------
      RCOMP = 0.0
      IF (A .GE. 20.0) GO TO 20
      T = A*ALOG(X) - X
      IF (A .GE. 1.0) GO TO 10
         RCOMP = (A*EXP(T))*(1.0 + GAM1(A))
         RETURN
   10 RCOMP = EXP(T)/GAMMA(A)
      RETURN
C
   20 U = X/A
      IF (U .EQ. 0.0) RETURN
      T = (1.0/A)**2
      T1 = (((0.75*T - 1.0)*T + 3.5)*T - 105.0)/(A*1260.0)
      T1 = T1 - A*RLOG(U)
      RCOMP = RT2PIN*SQRT(A)*EXP(T1)
      RETURN
      END
