function csevl ( x, cs, n )

!*****************************************************************************80
!
!! CSEVL evaluates an N term Chebyshev series.
!
!  Modified:
!
!    15 April 2003
!
!  Reference:
!
!    Roger Broucke,
!    Algorithm 446:
!    Ten Subroutines for the Manipulation of Chebyshev Series,
!    Communications of the ACM,
!    Volume 16, 1973, pages 254-256.
!
!    Leslie Fox, Ian Parker,
!    Chebyshev Polynomials in Numerical Analysis,
!    Oxford Press, 1968,
!    LC: QA297.F65.
!
!    David Kahaner, Cleve Moler, Steven Nash,
!    Numerical Methods and Software,
!    Prentice Hall, 1989,
!    ISBN: 0-13-627258-4,
!    LC: TA345.K34.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, the argument at which the series is 
!    to be evaluated.
!    X must satisfy -1.0 <= X <= 1.0.
!
!    Input, real ( kind = 8 ) CS(N), the array of N terms of a Chebyshev series.
!    In evaluating CS, only half the first coefficient is summed.
!
!    Input, integer ( kind = 4 ) N, the number of terms in array CS.
!    N must be at least 1, and no more than 1000.
!
!    Output, real ( kind = 8 ) CSEVL, the value of the Chebyshev series.
!
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) b0
  real ( kind = 8 ) b1
  real ( kind = 8 ) b2
  real ( kind = 8 ) cs(n)
  real ( kind = 8 ) csevl
  integer ( kind = 4 ) i
  real ( kind = 8 ) x

  if ( n < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CSEVL - Fatal error!'
    write ( *, '(a)' ) '  Number of terms N is less than 1.'
    stop
  end if

  if ( 1000 < n ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CSEVL - Fatal error!'
    write ( *, '(a)' ) '  The number of terms is more than 1000.'
    stop
  end if

  if ( x < -1.0D+00 .or. 1.0D+00 < x ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CSEVL - Fatal error!'
    write ( *, '(a)' ) '  The input argument X is outside the interval [-1,1].'
    stop
  end if

  b1 = 0.0D+00
  b0 = 0.0D+00

  do i = n, 1, -1
    b2 = b1
    b1 = b0
    b0 = 2.0D+00 * x * b1 - b2 + cs(i)
  end do

  csevl = 0.5D+00 * ( b0 - b2 )

  return
end
function error_f ( x )

!*****************************************************************************80
!
!! ERROR_F computes the error function.
!
!  Discussion:
!
!    This function was renamed "ERROR_F" from "ERF", to avoid a conflict
!    with the name of a corresponding routine often, but not always,
!    supplied as part of the math support library.
!
!    The definition of the error function is:
!
!      ERF(X) = ( 2 / SQRT ( PI ) ) * Integral ( 0 <= T <= X ) EXP ( -T**2 ) dT
!
!  Modified:
!
!    15 April 2003
!
!  Author:
!
!    Wayne Fullerton
!
!  Reference:
!
!    David Kahaner, Cleve Moler, Steven Nash,
!    Numerical Methods and Software,
!    Prentice Hall, 1989,
!    ISBN: 0-13-627258-4,
!    LC: TA345.K34.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, the argument of the error function.
!
!    Output, real ( kind = 8 ) ERROR_F, the value of the error function at X.
!
  implicit none

  real ( kind = 8 ) csevl
  real ( kind = 8 ) error_f
  real ( kind = 8 ) error_fc
  real ( kind = 8 ), parameter, dimension ( 13 ) :: erfcs = (/ &
    -0.049046121234691808D+00, -0.14226120510371364D+00, &
     0.010035582187599796D+00, -0.000576876469976748D+00, &
     0.000027419931252196D+00, -0.000001104317550734D+00, &
     0.000000038488755420D+00, -0.000000001180858253D+00, &
     0.000000000032334215D+00, -0.000000000000799101D+00, &
     0.000000000000017990D+00, -0.000000000000000371D+00, &
     0.000000000000000007D+00 /)
  integer ( kind = 4 ) inits
  integer ( kind = 4 ), save :: nterf = 0
  real ( kind = 8 ), save :: sqeps = 0.0D+00
  real ( kind = 8 ), parameter :: sqrtpi = 1.7724538509055160D+00
  real ( kind = 8 ) x
  real ( kind = 8 ), save :: xbig = 0.0D+00
  real ( kind = 8 ) y
!
!  Initialize the Chebyshev series.
!
  if ( nterf == 0 ) then
    nterf = inits ( erfcs, 13, 0.1D+00 * epsilon ( erfcs ) )
    xbig = sqrt ( - log ( sqrtpi * epsilon ( xbig ) ) )
    sqeps = sqrt ( 2.0D+00 * epsilon ( sqeps ) )
  end if

  y = abs ( x )

  if ( y <= sqeps ) then
    error_f = 2.0D+00 * x / sqrtpi
  else if ( y <= 1.0D+00 ) then
    error_f = x * ( 1.0D+00 + csevl ( 2.0D+00 * x**2 - 1.0D+00, erfcs, nterf ) )
  else if ( y <= xbig ) then
    error_f = sign ( 1.0D+00 - error_fc ( y ), x )
  else
    error_f = sign ( 1.0D+00, x )
  end if

  return
end
function error_fc ( x )

!*****************************************************************************80
!
!! ERROR_FC computes the complementary error function.
!
!  Discussion:
!
!    This function was renamed "ERROR_FC" from "ERFC", to avoid a conflict
!    with the name of a corresponding routine often, but not always,
!    supplied as part of the math support library.
!
!    The definition of the complementary error function is:
!
!      ERFC(X) = 1 - ERF(X)
!
!    where ERF(X) is the error function.
!
!  Modified:
!
!    26 August 2001
!
!  Author:
!
!    Wayne Fullerton
!
!  Reference:
!
!    David Kahaner, Cleve Moler, Steven Nash,
!    Numerical Methods and Software,
!    Prentice Hall, 1989,
!    ISBN: 0-13-627258-4,
!    LC: TA345.K34.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, the argument of the error function.
!
!    Output, real ( kind = 8 ) ERROR_FC, the value of the complementary 
!    error function at X.
!
  implicit none

  real ( kind = 8 ) csevl
  real ( kind = 8 ) error_fc
  real ( kind = 8 ), parameter, dimension ( 13 ) :: erfcs = (/ &
   -0.049046121234691808D+00, &
   -0.14226120510371364D+00,  &
    0.010035582187599796D+00, &
   -0.000576876469976748D+00, &
    0.000027419931252196D+00, &
   -0.000001104317550734D+00, &
    0.000000038488755420D+00, &
   -0.000000001180858253D+00, &
    0.000000000032334215D+00, &
   -0.000000000000799101D+00, &
    0.000000000000017990D+00, &
   -0.000000000000000371D+00, &
    0.000000000000000007D+00 /)
  real ( kind = 8 ), parameter, dimension ( 24 ) :: erfccs = (/ &
    0.0715179310202925D+00, &
   -0.026532434337606719D+00, &
    0.001711153977920853D+00, &
   -0.000163751663458512D+00, &
    0.000019871293500549D+00, &
   -0.000002843712412769D+00, &
    0.000000460616130901D+00, &
   -0.000000082277530261D+00, &
    0.000000015921418724D+00, &
   -0.000000003295071356D+00, &
    0.000000000722343973D+00, &
   -0.000000000166485584D+00, &
    0.000000000040103931D+00, &
   -0.000000000010048164D+00, &
    0.000000000002608272D+00, &
   -0.000000000000699105D+00, &
    0.000000000000192946D+00, &
   -0.000000000000054704D+00, &
    0.000000000000015901D+00, &
   -0.000000000000004729D+00, &
    0.000000000000001432D+00, &
   -0.000000000000000439D+00, &
    0.000000000000000138D+00, &
   -0.000000000000000048D+00 /)
  real ( kind = 8 ), parameter, dimension ( 23 ) :: erc2cs = (/ &
   -0.069601346602309501D+00, &
   -0.041101339362620893D+00, &
    0.003914495866689626D+00, &
   -0.000490639565054897D+00, &
    0.000071574790013770D+00, &
   -0.000011530716341312D+00, &
    0.000001994670590201D+00, &
   -0.000000364266647159D+00, &
    0.000000069443726100D+00, &
   -0.000000013712209021D+00, &
    0.000000002788389661D+00, &
   -0.000000000581416472D+00, &
    0.000000000123892049D+00, &
   -0.000000000026906391D+00, &
    0.000000000005942614D+00, &
   -0.000000000001332386D+00, &
    0.000000000000302804D+00, &
   -0.000000000000069666D+00, &
    0.000000000000016208D+00, &
   -0.000000000000003809D+00, &
    0.000000000000000904D+00, &
   -0.000000000000000216D+00, &
    0.000000000000000052D+00 /)
  real ( kind = 8 ) eta
  integer ( kind = 4 ) inits
  integer ( kind = 4 ), save :: nterc2 = 0
  integer ( kind = 4 ), save :: nterf = 0
  integer ( kind = 4 ), save :: nterfc = 0
  real ( kind = 8 ), save :: sqeps = 0.0D+00
  real ( kind = 8 ), parameter :: sqrtpi = 1.7724538509055160D+00
  real ( kind = 8 ) x
  real ( kind = 8 ), save :: xmax = 0.0D+00
  real ( kind = 8 ), save :: xsml = 0.0D+00
  real ( kind = 8 ) y

  if ( nterf == 0 ) then

    eta = 0.1D+00 * epsilon ( eta )
    nterf = inits ( erfcs, 13, eta )
    nterfc = inits ( erfccs, 24, eta )
    nterc2 = inits ( erc2cs, 23, eta )

    xsml = -sqrt ( - log ( sqrtpi * epsilon ( xsml ) ) )
    xmax = sqrt ( - log ( sqrtpi * tiny ( xmax ) ) )
    xmax = xmax - 0.5D+00 * log ( xmax ) / xmax - 0.01D+00
    sqeps = sqrt ( 2.0D+00 * epsilon ( sqeps ) )

  end if

  if ( x <= xsml ) then
    error_fc = 2.0D+00
    return
  end if
!
!  X so big that ERFC will underflow.
!
  if ( xmax < x ) then
    error_fc = 0.0D+00
    return
  end if

  y = abs ( x )
!
!  erfc(x) = 1.0D+00 - erf(x) for -1 <= x <= 1.
!
  if ( y <= 1.0D+00 ) then

    if ( y < sqeps ) then
      error_fc = 1.0D+00 - 2.0D+00 * x / sqrtpi
    else if ( sqeps <= y ) then
      error_fc = 1.0D+00 - x * ( 1.0D+00 + &
        csevl ( 2.0D+00 * x * x - 1.0D+00, erfcs, nterf ) )
    end if

    return

  end if
!
!  For 1 < |x| <= xmax, erfc(x) = 1.0D+00 - erf(x)
!
  y = y * y

  if ( y <= 4.0D+00 ) then
    error_fc = exp ( -y ) / abs ( x ) * ( 0.5D+00 &
      + csevl ( ( 8.0D+00 / y - 5.0D+00 ) / 3.0D+00, erc2cs, nterc2 ) )
  else
    error_fc = exp ( -y ) / abs ( x ) * ( 0.5D+00 &
      + csevl ( 8.0D+00 / y - 1.0D+00, erfccs, nterfc ) )
  end if

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

  return
end
function gamma_log ( x )

!*****************************************************************************80
!
!! GAMMA_LOG calculates the natural logarithm of GAMMA ( X ) for positive X.
!
!  Discussion:
!
!    Computation is based on an algorithm outlined in references 1 and 2.
!    The program uses rational functions that theoretically approximate
!    log ( GAMMA(X) ) to at least 18 significant decimal digits.  The
!    approximation for 12 < X is from Hart, while approximations
!    for X < 12.0 are similar to those in reference Cody and Hillstrom,
!    but are unpublished.
!
!    The accuracy achieved depends on the arithmetic system, the compiler,
!    intrinsic functions, and proper selection of the machine-dependent
!    constants.
!
!  Modified:
!
!    16 June 1999
!
!  Authors:
!
!    William Cody and Laura Stoltz
!    Argonne National Laboratory
!
!  Reference:
!
!    William Cody, Kenneth Hillstrom,
!    Chebyshev Approximations for the Natural Logarithm of the Gamma Function,
!    Mathematics of Computation,
!    Volume 21, Number 98, April 1967, pages 198-203.
!
!    Kenneth Hillstrom,
!    ANL/AMD Program ANLC366S, DGAMMA/DLGAMA,
!    May 1969.
!
!    John Hart, Ward Cheney, Charles Lawson, Hans Maely, Charles Mesztenyi,
!    John Rice, Henry Thatcher, Christop Witzgall,
!    Computer Approximations,
!    Wiley, 1968,
!    LC: QA297.C64.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, the argument of the Gamma function.  
!    X must be positive.
!
!    Output, real ( kind = 8 ) GAMMA_LOG, the logarithm of the Gamma 
!    function of X.
!    If X <= 0.0, or if overflow would occur, the program returns the
!    value HUGE().
!
!  Machine-dependent constants:
!
!    BETA   - radix for the floating-point representation.
!
!    MAXEXP - the smallest positive power of BETA that overflows.
!
!    XBIG   - largest argument for which LN(GAMMA(X)) is representable
!             in the machine, i.e., the solution to the equation
!             LN(GAMMA(XBIG)) = BETA**MAXEXP.
!
!    XINF   - largest machine representable floating-point number;
!             approximately BETA**MAXEXP.
!
!    FRTBIG - Rough estimate of the fourth root of XBIG
!
!
!    Approximate values for some important machines are:
!
!                              BETA      MAXEXP         XBIG
!
!    CRAY-1        (S.P.)        2        8191       9.62D+2461
!    Cyber 180/855
!      under NOS   (S.P.)        2        1070       1.72D+319
!    IEEE (IBM/XT,
!      SUN, etc.)  (S.P.)        2         128       4.08D+36
!    IEEE (IBM/XT,
!    SUN, etc.)  (D.P.)        2        1024       2.55D+305
!    IBM 3033      (D.P.)       16          63       4.29D+73
!    VAX D-Format  (D.P.)        2         127       2.05D+36
!    VAX G-Format  (D.P.)        2        1023       1.28D+305
!
!
!                            FRTBIG
!
!    CRAY-1        (S.P.)   3.13D+615
!    Cyber 180/855
!      under NOS   (S.P.)   6.44D+79
!    IEEE (IBM/XT,
!      SUN, etc.)  (S.P.)   1.42D+9
!    IEEE (IBM/XT,
!      SUN, etc.)  (D.P.)   2.25D+76
!    IBM 3033      (D.P.)   2.56D+18
!    VAX D-Format  (D.P.)   1.20D+9
!    VAX G-Format  (D.P.)   1.89D+76
!
  implicit none

  real ( kind = 8 ), parameter, dimension ( 7 ) :: c = (/ &
    -1.910444077728D-03, &
     8.4171387781295D-04, &
    -5.952379913043012D-04, &
     7.93650793500350248D-04, &
    -2.777777777777681622553D-03, &
     8.333333333333333331554247D-02, &
     5.7083835261D-03 /)
  real ( kind = 8 ) corr
  real ( kind = 8 ), parameter :: d1 = - 5.772156649015328605195174D-01
  real ( kind = 8 ), parameter :: d2 =   4.227843350984671393993777D-01
  real ( kind = 8 ), parameter :: d4 =   1.791759469228055000094023D+00
  real ( kind = 8 ), parameter :: frtbig = 1.42D+09
  integer ( kind = 4 ) i
  real ( kind = 8 ) gamma_log
  real ( kind = 8 ), parameter, dimension ( 8 ) :: p1 = (/ &
    4.945235359296727046734888D+00, &
    2.018112620856775083915565D+02, &
    2.290838373831346393026739D+03, &
    1.131967205903380828685045D+04, &
    2.855724635671635335736389D+04, &
    3.848496228443793359990269D+04, &
    2.637748787624195437963534D+04, &
    7.225813979700288197698961D+03 /)
  real ( kind = 8 ), parameter, dimension ( 8 ) :: p2 = (/ &
    4.974607845568932035012064D+00, &
    5.424138599891070494101986D+02, &
    1.550693864978364947665077D+04, &
    1.847932904445632425417223D+05, &
    1.088204769468828767498470D+06, &
    3.338152967987029735917223D+06, &
    5.106661678927352456275255D+06, &
    3.074109054850539556250927D+06 /)
  real ( kind = 8 ), parameter, dimension ( 8 ) :: p4 = (/ &
    1.474502166059939948905062D+04, &
    2.426813369486704502836312D+06, &
    1.214755574045093227939592D+08, &
    2.663432449630976949898078D+09, &
    2.940378956634553899906876D+10, &
    1.702665737765398868392998D+11, &
    4.926125793377430887588120D+11, &
    5.606251856223951465078242D+11 /)
  real ( kind = 8 ), parameter :: pnt68 = 0.6796875D+00
  real ( kind = 8 ), parameter, dimension ( 8 ) :: q1 = (/ &
    6.748212550303777196073036D+01, &
    1.113332393857199323513008D+03, &
    7.738757056935398733233834D+03, &
    2.763987074403340708898585D+04, &
    5.499310206226157329794414D+04, &
    6.161122180066002127833352D+04, &
    3.635127591501940507276287D+04, &
    8.785536302431013170870835D+03 /)
  real ( kind = 8 ), parameter, dimension ( 8 ) :: q2 = (/ &
    1.830328399370592604055942D+02, &
    7.765049321445005871323047D+03, &
    1.331903827966074194402448D+05, &
    1.136705821321969608938755D+06, &
    5.267964117437946917577538D+06, &
    1.346701454311101692290052D+07, &
    1.782736530353274213975932D+07, &
    9.533095591844353613395747D+06 /)
  real ( kind = 8 ), parameter, dimension ( 8 ) :: q4 = (/ &
    2.690530175870899333379843D+03, &
    6.393885654300092398984238D+05, &
    4.135599930241388052042842D+07, &
    1.120872109616147941376570D+09, &
    1.488613728678813811542398D+10, &
    1.016803586272438228077304D+11, &
    3.417476345507377132798597D+11, &
    4.463158187419713286462081D+11 /)
  real ( kind = 8 ) res
  real ( kind = 8 ) roundoff
  real ( kind = 8 ), parameter :: sqrtpi = 0.9189385332046727417803297D+00
  real ( kind = 8 ) x
  real ( kind = 8 ), parameter :: xbig = 4.08D+36
  real ( kind = 8 ) xden
  real ( kind = 8 ) xm1
  real ( kind = 8 ) xm2
  real ( kind = 8 ) xm4
  real ( kind = 8 ) xnum
  real ( kind = 8 ) xsq
!
!  Return immediately if the argument is out of range.
!
  if ( x <= 0.0D+00 .or. xbig < x ) then
    gamma_log = huge ( gamma_log )
    return
  end if

  roundoff = epsilon ( roundoff )

  if ( x <= roundoff ) then

    res = - log ( x )

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

    if ( x < pnt68 ) then
      corr = - log ( x )
      xm1 = x
    else
      corr = 0.0D+00
      xm1 = ( x - 0.5D+00 ) - 0.5D+00
    end if

    if ( x <= 0.5D+00 .or. pnt68 <= x ) then

      xden = 1.0D+00
      xnum = 0.0D+00

      do i = 1, 8
        xnum = xnum * xm1 + p1(i)
        xden = xden * xm1 + q1(i)
      end do

      res = corr + ( xm1 * ( d1 + xm1 * ( xnum / xden ) ) )

    else

      xm2 = ( x - 0.5D+00 ) - 0.5D+00
      xden = 1.0D+00
      xnum = 0.0D+00
      do i = 1, 8
        xnum = xnum * xm2 + p2(i)
        xden = xden * xm2 + q2(i)
      end do

      res = corr + xm2 * ( d2 + xm2 * ( xnum / xden ) )

    end if

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

    xm2 = x - 2.0D+00
    xden = 1.0D+00
    xnum = 0.0D+00
    do i = 1, 8
      xnum = xnum * xm2 + p2(i)
      xden = xden * xm2 + q2(i)
    end do

    res = xm2 * ( d2 + xm2 * ( xnum / xden ) )

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

    xm4 = x - 4.0D+00
    xden = - 1.0D+00
    xnum = 0.0D+00
    do i = 1, 8
      xnum = xnum * xm4 + p4(i)
      xden = xden * xm4 + q4(i)
    end do

    res = d4 + xm4 * ( xnum / xden )

  else

    res = 0.0D+00

    if ( x <= frtbig ) then

      res = c(7)
      xsq = x * x

      do i = 1, 6
        res = res / xsq + c(i)
      end do

    end if

    res = res / x
    corr = log ( x )
    res = res + sqrtpi - 0.5D+00 * corr
    res = res + x * ( corr - 1.0D+00 )

  end if

  gamma_log = res

  return
end
subroutine get_problem_num ( problem_num )

!*****************************************************************************80
!
!! GET_PROBLEM_NUM returns the number of test integration problems.
!
!  Modified:
!
!    06 June 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer ( kind = 4 ) PROBLEM_NUM, the number of test problems.
!
  implicit none

  integer ( kind = 4 ) problem_num

  problem_num = 32

  return
end
subroutine get_seed ( seed )

!*****************************************************************************80
!
!! GET_SEED returns a seed for the random number generator.
!
!  Discussion:
!
!    The seed depends on the current time, and ought to be (slightly)
!    different every millisecond.  Once the seed is obtained, a random
!    number generator should be called a few times to further process
!    the seed.
!
!  Modified:
!
!    02 August 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer ( kind = 4 ) SEED, a pseudorandom seed value.
!
  implicit none

  integer ( kind = 4 ) i4_huge
  integer ( kind = 4 ) seed
  real    ( kind = 8 ) temp
  character ( len = 10 ) time
  character ( len = 8 ) today
  integer ( kind = 4 ) values(8)
  character ( len = 5 ) zone

  call date_and_time ( today, time, zone, values )

  temp = 0.0D+00

  temp = temp + real ( values(2) - 1, kind = 8 ) /  11.0D+00
  temp = temp + real ( values(3) - 1, kind = 8 ) /  30.0D+00
  temp = temp + real ( values(5),     kind = 8 ) /  23.0D+00
  temp = temp + real ( values(6),     kind = 8 ) /  59.0D+00
  temp = temp + real ( values(7),     kind = 8 ) /  59.0D+00
  temp = temp + real ( values(8),     kind = 8 ) / 999.0D+00
  temp = temp                                    /   6.0D+00

  do while ( temp <= 0.0D+00 )
    temp = temp + 1.0D+00
  end do

  do while ( 1.0D+00 < temp )
    temp = temp - 1.0D+00
  end do

  seed = int ( real ( i4_huge ( ), kind = 8 ) * temp )
!
!  Never use a seed of 0 or maximum integer.
!
  if ( seed == 0 ) then
    seed = 1
  end if

  if ( seed == i4_huge ( ) ) then
    seed = seed - 1
  end if

  return
end
function i4_huge ( )

!*****************************************************************************80
!
!! I4_HUGE returns a "huge" I4.
!
!  Discussion:
!
!    On an IEEE 32 bit machine, I4_HUGE should be 2**31 - 1, and its
!    bit pattern should be
!
!     01111111111111111111111111111111
!
!    In this case, its numerical value is 2147483647.
!
!    Using the Dec/Compaq/HP Alpha FORTRAN compiler FORT, I could
!    use I4_HUGE() and HUGE interchangeably.
!
!    However, when using the G95, the values returned by HUGE were
!    not equal to 2147483647, apparently, and were causing severe
!    and obscure errors in my random number generator, which needs to
!    add I4_HUGE to the seed whenever the seed is negative.  So I
!    am backing away from invoking HUGE, whereas I4_HUGE is under
!    my control.
!
!    Explanation: because under G95 the default integer type is 64 bits!
!    So HUGE ( 1 ) = a very very huge integer indeed, whereas
!    I4_HUGE ( ) = the same old 32 bit big value.
!
!    An I4 is an integer ( kind = 4 ) value.
!
!  Modified:
!
!    26 January 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer ( kind = 4 ) I4_HUGE, a "huge" I4.
!
  implicit none

  integer ( kind = 4 ) i4
  integer ( kind = 4 ) i4_huge

  i4_huge = 2147483647

  return
end
function inits ( os, nos, eta )

!*****************************************************************************80
!
!! INITS estimates the order of an orthogonal series for a given accuracy.
!
!  Discussion:
!
!    Because this is a function, it is not possible to print out
!    warning error messages.  Therefore, if an error condition is
!    detected, a bogus value of INITS is returned.
!
!  Modified:
!
!    15 April 2003
!
!  Reference:
!
!    David Kahaner, Cleve Moler, Steven Nash,
!    Numerical Methods and Software,
!    Prentice Hall, 1989,
!    ISBN: 0-13-627258-4,
!    LC: TA345.K34.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) OS(NOS), the coefficients in the series.
!
!    Input, integer ( kind = 4 ) NOS, the number of coefficients.  NOS must be
!    at least 1, or an error condition arises.
!
!    Input, real ( kind = 8 ) ETA, the requested accuracy of the series.
!    Ordinarily, ETA will be chosen to be one-tenth machine precision.
!
!    Output, integer ( kind = 4 ) INITS, the order of the series guaranteeing 
!    the given accuracy.  However, on error, INITS will be returned
!    as a negative number.
!
  implicit none

  integer ( kind = 4 ) nos

  real ( kind = 8 ) err
  real ( kind = 8 ) eta
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ii
  integer ( kind = 4 ) inits
  real ( kind = 8 ) os(nos)
!
!  Fatal error.  Number of coefficients less than 1.
!  But an error message cannot be printed out because this is a function.
!
  if ( nos < 1 ) then
    inits = - huge ( i )
    return
  end if

  err = 0.0D+00

  i = 0

  do ii = 1, nos

    i = nos + 1 - ii
    err = err + abs ( os(i) )

    if ( eta < err ) then
      i = nos + 1 - ii
      exit
    end if

  end do
!
!  Eta may be too small.  Accuracy cannot be guaranteed.
!  But an error message cannot be printed out because this is a function.
!
  if ( i == 0 ) then
    i = - nos
  end if

  inits = i

  return
end
subroutine normal_01_cdf_inv ( cdf, x )

!*****************************************************************************80
!
!! NORMAL_01_CDF_INV inverts the Normal 01 CDF.
!
!  Discussion:
!
!    Modified to handle case where R = 0 would otherwise require
!    evaluation of LOG(0).
!
!  Modified:
!
!    05 June 2007
!
!  Reference:
!
!    JD Beasley, SG Springer,
!    Algorithm AS 111:
!    The Percentage Points of the Normal Distribution,
!    Applied Statistics,
!    Volume 26, 1977, pages 118-121.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) CDF, the value of the CDF.
!    0.0 <= CDF <= 1.0.
!
!    Input, real ( kind = 8 ) X, the corresponding argument.
!
  implicit none

  real ( kind = 8 ), parameter :: a0 =   2.50662823884D+00
  real ( kind = 8 ), parameter :: a1 = -18.61500062529D+00
  real ( kind = 8 ), parameter :: a2 =  41.39119773534D+00
  real ( kind = 8 ), parameter :: a3 = -25.44106049637D+00
  real ( kind = 8 ), parameter :: b1 =  -8.47351093090D+00
  real ( kind = 8 ), parameter :: b2 =  23.08336743743D+00
  real ( kind = 8 ), parameter :: b3 = -21.06224101826D+00
  real ( kind = 8 ), parameter :: b4 =   3.13082909833D+00
  real ( kind = 8 ), parameter :: c0 =  -2.78718931138D+00
  real ( kind = 8 ), parameter :: c1 =  -2.29796479134D+00
  real ( kind = 8 ), parameter :: c2 =   4.85014127135D+00
  real ( kind = 8 ), parameter :: c3 =   2.32121276858D+00
  real ( kind = 8 ) cdf
  real ( kind = 8 ), parameter :: d1 =   3.54388924762D+00
  real ( kind = 8 ), parameter :: d2 =   1.63706781897D+00
  real ( kind = 8 ) q
  real ( kind = 8 ) r
  real ( kind = 8 ) r8_huge
  real ( kind = 8 ) x

  if ( cdf < 0.0D+00 .or. 1.0D+00 < cdf ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'NORMAL_01_CDF_INV - Fatal error!'
    write ( *, '(a)' ) '  CDF < 0 or 1 < CDF.'
    stop
  end if

  q = cdf - 0.5D+00

  q = min ( q, 0.5D+00 )
  q = max ( q, -0.5D+00 )
!
!  0.08 < CDF < 0.92
!
  if ( abs ( q ) <= 0.42D+00 ) then

    r = q * q

    x = q * ( ( ( &
           a3   * r &
         + a2 ) * r &
         + a1 ) * r &
         + a0 ) / ( ( ( ( &
           b4   * r &
         + b3 ) * r &
         + b2 ) * r &
         + b1 ) * r + 1.0D+00 )
!
!  CDF < 0.08 or 0.92 < CDF.
!
  else

    r = min ( cdf, 1.0D+00 - cdf )
    r = max ( r, 0.0D+00 )
    r = min ( r, 1.0D+00 )

    if ( r <= 0.0D+00 ) then

      x = r8_huge ( )

    else

      r = sqrt ( - log ( r ) )

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

    end if

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

  end if

  return
end
subroutine p00_box_gl05 ( problem, dim_num, sub_num, result )

!*****************************************************************************80
!
!! P00_BOX_GL05 uses Gauss-Legendre quadrature in an N-dimensional box. 
!
!  Discussion:
!
!    The rule is the product of 5-point Gauss-Legendre rules 
!    in each spatial dimension.
!
!  Modified:
!
!    03 June 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) PROBLEM, the problem number.
!
!    Input, integer ( kind = 4 ) DIM_NUM, the spatial dimension.
!
!    Input, integer ( kind = 4 ) SUB_NUM, the number of subdivisions 
!    in each dimension.
!
!    Output, real ( kind = 8 ) RESULT, the approximate integral.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ), parameter :: norder = 5
  integer ( kind = 4 ), parameter :: point_num = 1

  real ( kind = 8 ) a(dim_num)
  real ( kind = 8 ) a_sub(dim_num)
  real ( kind = 8 ) b(dim_num)
  real ( kind = 8 ) b_sub(dim_num)
  integer ( kind = 4 ) dim
  integer ( kind = 4 ) indx(dim_num)
  integer ( kind = 4 ) problem
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  real ( kind = 8 ) result
  integer ( kind = 4 ) sub_indx(dim_num)
  integer ( kind = 4 ) sub_num
  real ( kind = 8 ) value(point_num)
  real ( kind = 8 ) volume
  real ( kind = 8 ) w
  real ( kind = 8 ), save, dimension ( norder ) :: weight = (/ &
    0.236926885056189087514264040720D+00, &
    0.478628670499366468041291514836D+00, &
    0.568888888888888888888888888889D+00, &
    0.478628670499366468041291514836D+00, &
    0.236926885056189087514264040720D+00 /)
  real ( kind = 8 ) x(dim_num,point_num)
  real ( kind = 8 ), save, dimension ( norder ) :: xtab = (/ &
    - 0.906179845938663992797626878299D+00, &
    - 0.538469310105683091036314420700D+00, &
      0.0D+00, &
      0.538469310105683091036314420700D+00, &
      0.906179845938663992797626878299D+00 /)
!
!  Get the integration limits, and the weight adjustment factor.
!
  call p00_lim ( problem, dim_num, a, b )
!
!  Carry out the product rule.
!
  result = 0.0D+00
!
!  In the outer loop, we pick a sub-box.
!
  j = 0

  do

    call tuple_next ( 1, sub_num, dim_num, j, sub_indx )

    if ( j == 0 ) then
      exit
    end if

    do dim = 1, dim_num

      a_sub(dim) = ( real ( sub_num + 1 - sub_indx(dim), kind = 8 ) * a(dim)   &
                   + real (         - 1 + sub_indx(dim), kind = 8 ) * b(dim) ) &
                   / real ( sub_num,                     kind = 8 )

      b_sub(dim) = ( real ( sub_num     - sub_indx(dim), kind = 8 ) * a(dim)   &
                   + real (               sub_indx(dim), kind = 8 ) * b(dim) ) &
                   / real ( sub_num,                     kind = 8 )

    end do
!
!  In the inner loop, we go through all the points in the
!  cross product of the quadrature rule.
!
    k = 0

    do

      call tuple_next ( 1, norder, dim_num, k, indx )

      if ( k == 0  ) then
        exit
      end if

      w = product ( weight(indx(1:dim_num)) )

      x(1:dim_num,1) = 0.5D+00 * ( &
          a_sub(1:dim_num) * ( 1.0D+00 - xtab(indx(1:dim_num)) ) &
        + b_sub(1:dim_num) * ( 1.0D+00 + xtab(indx(1:dim_num)) ) )

      call p00_f ( problem, dim_num, point_num, x, value )

      result = result + w * value(1)

    end do

  end do
!
!  Get the volume.
!
  call p00_volume ( problem, dim_num, volume )

  result = result * volume / real ( ( 2 * sub_num )**dim_num, kind = 8 )

  return
end
subroutine p00_box_mc ( problem, dim_num, point_num, result )

!*****************************************************************************80
!
!! P00_BOX_MC integrates over an multi-dimensional box using Monte Carlo.
!
!  Modified:
!
!    03 June 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) PROBLEM, the problem number.
!
!    Input, integer ( kind = 4 ) DIM_NUM, the spatial dimension.
!
!    Input, integer ( kind = 4 ) POINT_NUM, the number of points to use.
!
!    Output, real ( kind = 8 ) RESULT, the approximate integral.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) point_num

  real ( kind = 8 ) a(dim_num)
  real ( kind = 8 ) b(dim_num)
  integer ( kind = 4 ) dim
  integer ( kind = 4 ) problem
  real ( kind = 8 ) quad
  real ( kind = 8 ) result
  real ( kind = 8 ), allocatable, dimension ( : ) :: value
  real ( kind = 8 ) volume
  real ( kind = 8 ), allocatable, dimension ( :, : ) :: x

  allocate ( value(1:point_num) )
  allocate ( x(1:dim_num,1:point_num) )
!
!  Get the volume.
!
  call p00_volume ( problem, dim_num, volume )
!
!  Get the integration limits.
!
  call p00_lim ( problem, dim_num, a, b )
!
!  Get the random values.
!
  call random_number ( harvest = x(1:dim_num,1:point_num) )
!
!  Map them to the domain.
!
  do dim = 1, dim_num
    x(dim,1:point_num) = ( ( 1.0D+00 - x(dim,1:point_num) ) * a(dim) &
                         + (         + x(dim,1:point_num) ) * b(dim) )
  end do
!
!  Evaluate the function.
!
  call p00_f ( problem, dim_num, point_num, x, value )

  result = sum ( value(1:point_num) ) * volume / real ( point_num, kind = 8 )

  deallocate ( value )
  deallocate ( x )

  return
end
subroutine p00_default ( problem, dim_num )

!*****************************************************************************80
!
!! P00_DEFAULT sets a problem to a default state.
!
!  Discussion:
!
!    If a problem uses vector parameters, then the spatial dimension
!    DIM_NUM is needed as input, so that the vector parameters can be
!    properly dimensioned.
!
!    Every problem keeps a count of the number of function calls.  Calling
!    this routine causes that count to be reset to 0.
!
!  Modified:
!
!    10 March 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) PROBLEM, the index of the desired problem.
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the problem.
!
  implicit none

  integer ( kind = 4 ) problem
  integer ( kind = 4 ) dim_num

  if ( problem == 1 ) then
    call p01_default ( dim_num )
  else if ( problem == 2 ) then
    call p02_default ( dim_num )
  else if ( problem == 3 ) then
    call p03_default ( dim_num )
  else if ( problem == 4 ) then
    call p04_default ( dim_num )
  else if ( problem == 5 ) then
    call p05_default ( dim_num )
  else if ( problem == 6 ) then
    call p06_default ( dim_num )
  else if ( problem == 7 ) then
    call p07_default ( dim_num )
  else if ( problem == 8 ) then
    call p08_default ( dim_num )
  else if ( problem == 9 ) then
    call p09_default ( dim_num )
  else if ( problem == 10 ) then
    call p10_default ( dim_num )
  else if ( problem == 11 ) then
    call p11_default ( dim_num )
  else if ( problem == 12 ) then
    call p12_default ( dim_num )
  else if ( problem == 13 ) then
    call p13_default ( dim_num )
  else if ( problem == 14 ) then
    call p14_default ( dim_num )
  else if ( problem == 15 ) then
    call p15_default ( dim_num )
  else if ( problem == 16 ) then
    call p16_default ( dim_num )
  else if ( problem == 17 ) then
    call p17_default ( dim_num )
  else if ( problem == 18 ) then
    call p18_default ( dim_num )
  else if ( problem == 19 ) then
    call p19_default ( dim_num )
  else if ( problem == 20 ) then
    call p20_default ( dim_num )
  else if ( problem == 21 ) then
    call p21_default ( dim_num )
  else if ( problem == 22 ) then
    call p22_default ( dim_num )
  else if ( problem == 23 ) then
    call p23_default ( dim_num )
  else if ( problem == 24 ) then
    call p24_default ( dim_num )
  else if ( problem == 25 ) then
    call p25_default ( dim_num )
  else if ( problem == 26 ) then
    call p26_default ( dim_num )
  else if ( problem == 27 ) then
    call p27_default ( dim_num )
  else if ( problem == 28 ) then
    call p28_default ( dim_num )
  else if ( problem == 29 ) then
    call p29_default ( dim_num )
  else if ( problem == 30 ) then
    call p30_default ( dim_num )
  else if ( problem == 31 ) then
    call p31_default ( dim_num )
  else if ( problem == 32 ) then
    call p32_default ( dim_num )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_DEFAULT - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal problem number = ', problem
    stop
  end if

  return
end
subroutine p00_exact ( problem, dim_num, exact )

!*****************************************************************************80
!
!! P00_EXACT returns the exact integral for any problem.
!
!  Modified:
!
!    26 February 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) PROBLEM, the index of the desired problem.
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the problem.
!
!    Output, real ( kind = 8 ) EXACT, the exact value of the integral.
!
  implicit none

  real ( kind = 8 ) exact
  integer ( kind = 4 ) problem
  integer ( kind = 4 ) dim_num

  if ( problem == 1 ) then
    call p01_exact ( dim_num, exact )
  else if ( problem == 2 ) then
    call p02_exact ( dim_num, exact )
  else if ( problem == 3 ) then
    call p03_exact ( dim_num, exact )
  else if ( problem == 4 ) then
    call p04_exact ( dim_num, exact )
  else if ( problem == 5 ) then
    call p05_exact ( dim_num, exact )
  else if ( problem == 6 ) then
    call p06_exact ( dim_num, exact )
  else if ( problem == 7 ) then
    call p07_exact ( dim_num, exact )
  else if ( problem == 8 ) then
    call p08_exact ( dim_num, exact )
  else if ( problem == 9 ) then
    call p09_exact ( dim_num, exact )
  else if ( problem == 10 ) then
    call p10_exact ( dim_num, exact )
  else if ( problem == 11 ) then
    call p11_exact ( dim_num, exact )
  else if ( problem == 12 ) then
    call p12_exact ( dim_num, exact )
  else if ( problem == 13 ) then
    call p13_exact ( dim_num, exact )
  else if ( problem == 14 ) then
    call p14_exact ( dim_num, exact )
  else if ( problem == 15 ) then
    call p15_exact ( dim_num, exact )
  else if ( problem == 16 ) then
    call p16_exact ( dim_num, exact )
  else if ( problem == 17 ) then
    call p17_exact ( dim_num, exact )
  else if ( problem == 18 ) then
    call p18_exact ( dim_num, exact )
  else if ( problem == 19 ) then
    call p19_exact ( dim_num, exact )
  else if ( problem == 20 ) then
    call p20_exact ( dim_num, exact )
  else if ( problem == 21 ) then
    call p21_exact ( dim_num, exact )
  else if ( problem == 22 ) then
    call p22_exact ( dim_num, exact )
  else if ( problem == 23 ) then
    call p23_exact ( dim_num, exact )
  else if ( problem == 24 ) then
    call p24_exact ( dim_num, exact )
  else if ( problem == 25 ) then
    call p25_exact ( dim_num, exact )
  else if ( problem == 26 ) then
    call p26_exact ( dim_num, exact )
  else if ( problem == 27 ) then
    call p27_exact ( dim_num, exact )
  else if ( problem == 28 ) then
    call p28_exact ( dim_num, exact )
  else if ( problem == 29 ) then
    call p29_exact ( dim_num, exact )
  else if ( problem == 30 ) then
    call p30_exact ( dim_num, exact )
  else if ( problem == 31 ) then
    call p31_exact ( dim_num, exact )
  else if ( problem == 32 ) then
    call p32_exact ( dim_num, exact )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_EXACT - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal problem number = ', problem
    stop
  end if

  return
end
subroutine p00_f ( problem, dim_num, point_num, x, value )

!*****************************************************************************80
!
!! P00_F evaluates the integrand for any problem.
!
!  Modified:
!
!    03 June 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) PROBLEM, the index of the desired problem.
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Input, integer ( kind = 4 ) POINT_NUM, the number of points.
!
!    Input, real ( kind = 8 ) X(DIM_NUM,POINT_NUM), the evaluation points.
!
!    Output, real ( kind = 8 ) VALUE(POINT_NUM), the function values.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) point_num

  integer ( kind = 4 ) problem
  real ( kind = 8 ) value(point_num)
  real ( kind = 8 ) x(dim_num,point_num)

  if ( problem == 1 ) then
    call p01_f ( dim_num, point_num, x, value )
  else if ( problem == 2 ) then
    call p02_f ( dim_num, point_num, x, value )
  else if ( problem == 3 ) then
    call p03_f ( dim_num, point_num, x, value )
  else if ( problem == 4 ) then
    call p04_f ( dim_num, point_num, x, value )
  else if ( problem == 5 ) then
    call p05_f ( dim_num, point_num, x, value )
  else if ( problem == 6 ) then
    call p06_f ( dim_num, point_num, x, value )
  else if ( problem == 7 ) then
    call p07_f ( dim_num, point_num, x, value )
  else if ( problem == 8 ) then
    call p08_f ( dim_num, point_num, x, value )
  else if ( problem == 9 ) then
    call p09_f ( dim_num, point_num, x, value )
  else if ( problem == 10 ) then
    call p10_f ( dim_num, point_num, x, value )
  else if ( problem == 11 ) then
    call p11_f ( dim_num, point_num, x, value )
  else if ( problem == 12 ) then
    call p12_f ( dim_num, point_num, x, value )
  else if ( problem == 13 ) then
    call p13_f ( dim_num, point_num, x, value )
  else if ( problem == 14 ) then
    call p14_f ( dim_num, point_num, x, value )
  else if ( problem == 15 ) then
    call p15_f ( dim_num, point_num, x, value )
  else if ( problem == 16 ) then
    call p16_f ( dim_num, point_num, x, value )
  else if ( problem == 17 ) then
    call p17_f ( dim_num, point_num, x, value )
  else if ( problem == 18 ) then
    call p18_f ( dim_num, point_num, x, value )
  else if ( problem == 19 ) then
    call p19_f ( dim_num, point_num, x, value )
  else if ( problem == 20 ) then
    call p20_f ( dim_num, point_num, x, value )
  else if ( problem == 21 ) then
    call p21_f ( dim_num, point_num, x, value )
  else if ( problem == 22 ) then
    call p22_f ( dim_num, point_num, x, value )
  else if ( problem == 23 ) then
    call p23_f ( dim_num, point_num, x, value )
  else if ( problem == 24 ) then
    call p24_f ( dim_num, point_num, x, value )
  else if ( problem == 25 ) then
    call p25_f ( dim_num, point_num, x, value )
  else if ( problem == 26 ) then
    call p26_f ( dim_num, point_num, x, value )
  else if ( problem == 27 ) then
    call p27_f ( dim_num, point_num, x, value )
  else if ( problem == 28 ) then
    call p28_f ( dim_num, point_num, x, value )
  else if ( problem == 29 ) then
    call p29_f ( dim_num, point_num, x, value )
  else if ( problem == 30 ) then
    call p30_f ( dim_num, point_num, x, value )
  else if ( problem == 31 ) then
    call p31_f ( dim_num, point_num, x, value )
  else if ( problem == 32 ) then
    call p32_f ( dim_num, point_num, x, value )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_F - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal problem number = ', problem
    stop
  end if

  return
end
subroutine p00_i4 ( problem, action, name, value )

!*****************************************************************************80
!
!! P00_I4 sets or gets I4 parameters for any problem.
!
!  Modified:
!
!    26 February 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) PROBLEM, the index of the desired problem.
!
!    Input, character ( len = * ) ACTION,
!    'D' to set a parameter to its default value.
!    'S' to set a parameter,
!    'G' to get a parameter,
!    'I' to increment a parameter.
!
!    Input, character ( len = * ) NAME, the name of the variable.
!    '#' is the number of calls to the integrand routine.
!
!    Input/output, integer ( kind = 4 ) VALUE.
!    * If ACTION = 'I', then VALUE is an input quantity, and is the
!      new value to be added to NAME.
!    * If ACTION = 'S', then VALUE is an input quantity, and is the
!      new value to be assigned to NAME.
!    * If ACTION = 'G', then VALUE is an output quantity, and is the
!      current value of NAME.
!
  implicit none

  character ( len = * ) action
  integer ( kind = 4 ) problem
  character ( len = * ) name
  integer ( kind = 4 ) value

  if ( problem == 1 ) then
    call p01_i4 ( action, name, value )
  else if ( problem == 2 ) then
    call p02_i4 ( action, name, value )
  else if ( problem == 3 ) then
    call p03_i4 ( action, name, value )
  else if ( problem == 4 ) then
    call p04_i4 ( action, name, value )
  else if ( problem == 5 ) then
    call p05_i4 ( action, name, value )
  else if ( problem == 6 ) then
    call p06_i4 ( action, name, value )
  else if ( problem == 7 ) then
    call p07_i4 ( action, name, value )
  else if ( problem == 8 ) then
    call p08_i4 ( action, name, value )
  else if ( problem == 9 ) then
    call p09_i4 ( action, name, value )
  else if ( problem == 10 ) then
    call p10_i4 ( action, name, value )
  else if ( problem == 11 ) then
    call p11_i4 ( action, name, value )
  else if ( problem == 12 ) then
    call p12_i4 ( action, name, value )
  else if ( problem == 13 ) then
    call p13_i4 ( action, name, value )
  else if ( problem == 14 ) then
    call p14_i4 ( action, name, value )
  else if ( problem == 15 ) then
    call p15_i4 ( action, name, value )
  else if ( problem == 16 ) then
    call p16_i4 ( action, name, value )
  else if ( problem == 17 ) then
    call p17_i4 ( action, name, value )
  else if ( problem == 18 ) then
    call p18_i4 ( action, name, value )
  else if ( problem == 19 ) then
    call p19_i4 ( action, name, value )
  else if ( problem == 20 ) then
    call p20_i4 ( action, name, value )
  else if ( problem == 21 ) then
    call p21_i4 ( action, name, value )
  else if ( problem == 22 ) then
    call p22_i4 ( action, name, value )
  else if ( problem == 23 ) then
    call p23_i4 ( action, name, value )
  else if ( problem == 24 ) then
    call p24_i4 ( action, name, value )
  else if ( problem == 25 ) then
    call p25_i4 ( action, name, value )
  else if ( problem == 26 ) then
    call p26_i4 ( action, name, value )
  else if ( problem == 27 ) then
    call p27_i4 ( action, name, value )
  else if ( problem == 28 ) then
    call p28_i4 ( action, name, value )
  else if ( problem == 29 ) then
    call p29_i4 ( action, name, value )
  else if ( problem == 30 ) then
    call p30_i4 ( action, name, value )
  else if ( problem == 31 ) then
    call p31_i4 ( action, name, value )
  else if ( problem == 32 ) then
    call p32_i4 ( action, name, value )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_I4 - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal problem number = ', problem
    stop
  end if

  return
end
subroutine p00_lim ( problem, dim_num, a, b )

!*****************************************************************************80
!
!! P00_LIM returns the integration limits for any problem.
!
!  Modified:
!
!    26 February 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) PROBLEM, the index of the test problem.
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Output, real ( kind = 8 ) A(DIM_NUM), B(DIM_NUM), the lower and upper 
!    limits of integration.
!    Note that if A = -HUGE(A), the lower limit is
!    actually negative infinity, and if B = HUGE(B), the upper limit
!    is actually infinity.
!
  implicit none

  integer ( kind = 4 ) dim_num

  real ( kind = 8 ) a(dim_num)
  real ( kind = 8 ) b(dim_num)
  integer ( kind = 4 ) problem

  if ( problem == 1 ) then
    call p01_lim ( dim_num, a, b )
  else if ( problem == 2 ) then
    call p02_lim ( dim_num, a, b )
  else if ( problem == 3 ) then
    call p03_lim ( dim_num, a, b )
  else if ( problem == 4 ) then
    call p04_lim ( dim_num, a, b )
  else if ( problem == 5 ) then
    call p05_lim ( dim_num, a, b )
  else if ( problem == 6 ) then
    call p06_lim ( dim_num, a, b )
  else if ( problem == 7 ) then
    call p07_lim ( dim_num, a, b )
  else if ( problem == 8 ) then
    call p08_lim ( dim_num, a, b )
  else if ( problem == 9 ) then
    call p09_lim ( dim_num, a, b )
  else if ( problem == 10 ) then
    call p10_lim ( dim_num, a, b )
  else if ( problem == 11 ) then
    call p11_lim ( dim_num, a, b )
  else if ( problem == 12 ) then
    call p12_lim ( dim_num, a, b )
  else if ( problem == 13 ) then
    call p13_lim ( dim_num, a, b )
  else if ( problem == 14 ) then
    call p14_lim ( dim_num, a, b )
  else if ( problem == 15 ) then
    call p15_lim ( dim_num, a, b )
  else if ( problem == 16 ) then
    call p16_lim ( dim_num, a, b )
  else if ( problem == 17 ) then
    call p17_lim ( dim_num, a, b )
  else if ( problem == 18 ) then
    call p18_lim ( dim_num, a, b )
  else if ( problem == 19 ) then
    call p19_lim ( dim_num, a, b )
  else if ( problem == 20 ) then
    call p20_lim ( dim_num, a, b )
  else if ( problem == 21 ) then
    call p21_lim ( dim_num, a, b )
  else if ( problem == 22 ) then
    call p22_lim ( dim_num, a, b )
  else if ( problem == 23 ) then
    call p23_lim ( dim_num, a, b )
  else if ( problem == 24 ) then
    call p24_lim ( dim_num, a, b )
  else if ( problem == 25 ) then
    call p25_lim ( dim_num, a, b )
  else if ( problem == 26 ) then
    call p26_lim ( dim_num, a, b )
  else if ( problem == 27 ) then
    call p27_lim ( dim_num, a, b )
  else if ( problem == 28 ) then
    call p28_lim ( dim_num, a, b )
  else if ( problem == 29 ) then
    call p29_lim ( dim_num, a, b )
  else if ( problem == 30 ) then
    call p30_lim ( dim_num, a, b )
  else if ( problem == 31 ) then
    call p31_lim ( dim_num, a, b )
  else if ( problem == 32 ) then
    call p32_lim ( dim_num, a, b )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_LIM - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal problem number = ', problem
    stop
  end if

  return
end
subroutine p00_name ( problem, name )

!*****************************************************************************80
!
!! P00_NAME returns the name of the problem.
!
!  Modified:
!
!    26 February 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) PROBLEM, the index of the test problem.
!
!    Output, character ( len = * ) NAME, the name of the problem.
!
  implicit none

  integer ( kind = 4 ) problem
  character ( len = * ) name

  if ( problem == 1 ) then
    call p01_name ( name )
  else if ( problem == 2 ) then
    call p02_name ( name )
  else if ( problem == 3 ) then
    call p03_name ( name )
  else if ( problem == 4 ) then
    call p04_name ( name )
  else if ( problem == 5 ) then
    call p05_name ( name )
  else if ( problem == 6 ) then
    call p06_name ( name )
  else if ( problem == 7 ) then
    call p07_name ( name )
  else if ( problem == 8 ) then
    call p08_name ( name )
  else if ( problem == 9 ) then
    call p09_name ( name )
  else if ( problem == 10 ) then
    call p10_name ( name )
  else if ( problem == 11 ) then
    call p11_name ( name )
  else if ( problem == 12 ) then
    call p12_name ( name )
  else if ( problem == 13 ) then
    call p13_name ( name )
  else if ( problem == 14 ) then
    call p14_name ( name )
  else if ( problem == 15 ) then
    call p15_name ( name )
  else if ( problem == 16 ) then
    call p16_name ( name )
  else if ( problem == 17 ) then
    call p17_name ( name )
  else if ( problem == 18 ) then
    call p18_name ( name )
  else if ( problem == 19 ) then
    call p19_name ( name )
  else if ( problem == 20 ) then
    call p20_name ( name )
  else if ( problem == 21 ) then
    call p21_name ( name )
  else if ( problem == 22 ) then
    call p22_name ( name )
  else if ( problem == 23 ) then
    call p23_name ( name )
  else if ( problem == 24 ) then
    call p24_name ( name )
  else if ( problem == 25 ) then
    call p25_name ( name )
  else if ( problem == 26 ) then
    call p26_name ( name )
  else if ( problem == 27 ) then
    call p27_name ( name )
  else if ( problem == 28 ) then
    call p28_name ( name )
  else if ( problem == 29 ) then
    call p29_name ( name )
  else if ( problem == 30 ) then
    call p30_name ( name )
  else if ( problem == 31 ) then
    call p31_name ( name )
  else if ( problem == 32 ) then
    call p32_name ( name )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_NAME - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal problem number = ', problem
    stop
  end if

  return
end
subroutine p00_r8vec ( problem, action, name, dim_num, value )

!*****************************************************************************80
!
!! P00_R8VEC sets or gets R8VEC parameters for any problem.
!
!  Modified:
!
!    09 March 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) PROBLEM, the number of the test problem.
!
!    Input, character ( len = * ) ACTION, the action.
!    'D' sets the  object to a default value.
!    'G' means the current value of the object should be returned.
!    'R' means randomize the object and return it.
!    'S' means the object should be set to VALUE.
!
!    Input, character ( len = * ) NAME, the name of the variable.
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the variable.
!
!    Input/output, real ( kind = 8 ) VALUE(DIM_NUM), the value of the variable.
!
  implicit none

  integer ( kind = 4 ) dim_num

  character ( len = * ) action
  integer ( kind = 4 ) problem
  character ( len = * ) name
  real ( kind = 8 ) value(dim_num)

  if ( problem == 1 ) then

  else if ( problem == 2 ) then

  else if ( problem == 3 ) then

  else if ( problem == 4 ) then

  else if ( problem == 5 ) then

  else if ( problem == 6 ) then

  else if ( problem == 7 ) then

  else if ( problem == 8 ) then

  else if ( problem == 9 ) then
    call p09_r8vec ( action, name, dim_num, value )
  else if ( problem == 10 ) then

  else if ( problem == 11 ) then

  else if ( problem == 12 ) then

  else if ( problem == 13 ) then

  else if ( problem == 14 ) then

  else if ( problem == 15 ) then

  else if ( problem == 16 ) then
    call p16_r8vec ( action, name, dim_num, value )
  else if ( problem == 17 ) then
    call p17_r8vec ( action, name, dim_num, value )
  else if ( problem == 18 ) then
    call p18_r8vec ( action, name, dim_num, value )
  else if ( problem == 19 ) then
    call p19_r8vec ( action, name, dim_num, value )
  else if ( problem == 20 ) then

  else if ( problem == 21 ) then

  else if ( problem == 22 ) then

  else if ( problem == 23 ) then

  else if ( problem == 24 ) then
    call p24_r8vec ( action, name, dim_num, value )
  else if ( problem == 25 ) then

  else if ( problem == 26 ) then
    call p26_r8vec ( action, name, dim_num, value )
  else if ( problem == 27 ) then
    call p27_r8vec ( action, name, dim_num, value )
  else if ( problem == 28 ) then
    call p28_r8vec ( action, name, dim_num, value )
  else if ( problem == 29 ) then
    call p29_r8vec ( action, name, dim_num, value )
  else if ( problem == 30 ) then
    call p30_r8vec ( action, name, dim_num, value )
  else if ( problem == 31 ) then
    call p31_r8vec ( action, name, dim_num, value )
  else if ( problem == 32 ) then
    call p32_r8vec ( action, name, dim_num, value )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_R8VEC - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal problem number = ', problem
    stop
  end if

  return
end
subroutine p00_region ( problem, region )

!*****************************************************************************80
!
!! P00_REGION returns the name of the integration region for any problem.
!
!  Discussion:
!
!    I thought I was going to use this idea a lot, but most of my test
!    regions are boxes.
!
!    BALL
!      the interior of a 2D circle,
!      the interior of a 3D sphere,
!      the interior of an ND sphere.
!
!    BOX
!      a 1D finite line segment,
!      a 2D finite rectangle,
!      a 3D box,
!      an ND box.
!
!    SIMPLEX 
!      a 2D triangle,
!      a 3D tetrahedron,
!      an ND simplex.
!      The "unit simplex" in ND is the set of nonnegative points X 
!      such that sum ( X(1:N) ) <= 1.
!
!    SPACE
!      a 1D infinite line,
!      a 2D infinite place,
!      a 3D space,
!      an ND space.
!
!    SPHERE
!      the circumference of a 2D circle,
!      the surface of a 3D sphere,
!      the surface of an ND sphere.
!
!  Modified:
!
!    26 February 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) PROBLEM, the number of the test problem.
!
!    Output, character ( len = * ) REGION, the name of the integration region.
!
  implicit none

  integer ( kind = 4 ) problem
  character ( len = * ) region

  if ( problem == 1 ) then
    call p01_region ( region )
  else if ( problem == 2 ) then
    call p02_region ( region )
  else if ( problem == 3 ) then
    call p03_region ( region )
  else if ( problem == 4 ) then
    call p04_region ( region )
  else if ( problem == 5 ) then
    call p05_region ( region )
  else if ( problem == 6 ) then
    call p06_region ( region )
  else if ( problem == 7 ) then
    call p07_region ( region )
  else if ( problem == 8 ) then
    call p08_region ( region )
  else if ( problem == 9 ) then
    call p09_region ( region )
  else if ( problem == 10 ) then
    call p10_region ( region )
  else if ( problem == 11 ) then
    call p11_region ( region )
  else if ( problem == 12 ) then
    call p12_region ( region )
  else if ( problem == 13 ) then
    call p13_region ( region )
  else if ( problem == 14 ) then
    call p14_region ( region )
  else if ( problem == 15 ) then
    call p15_region ( region )
  else if ( problem == 16 ) then
    call p16_region ( region )
  else if ( problem == 17 ) then
    call p17_region ( region )
  else if ( problem == 18 ) then
    call p18_region ( region )
  else if ( problem == 19 ) then
    call p19_region ( region )
  else if ( problem == 20 ) then
    call p20_region ( region )
  else if ( problem == 21 ) then
    call p21_region ( region )
  else if ( problem == 22 ) then
    call p22_region ( region )
  else if ( problem == 23 ) then
    call p23_region ( region )
  else if ( problem == 24 ) then
    call p24_region ( region )
  else if ( problem == 25 ) then
    call p25_region ( region )
  else if ( problem == 26 ) then
    call p26_region ( region )
  else if ( problem == 27 ) then
    call p27_region ( region )
  else if ( problem == 28 ) then
    call p28_region ( region )
  else if ( problem == 29 ) then
    call p29_region ( region )
  else if ( problem == 30 ) then
    call p30_region ( region )
  else if ( problem == 31 ) then
    call p31_region ( region )
  else if ( problem == 32 ) then
    call p32_region ( region )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_REGION - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal problem number = ', problem
    stop
  end if

  return
end
subroutine p00_remap01 ( problem, dim_num, point_num, x01, xab )

!*****************************************************************************80
!
!! P00_REMAP01 remaps points in [0,1]^DIM_NUM into [A(1:DIM_NUM),B(1:DIM_NUM)].
!
!  Modified:
!
!    12 February 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) PROBLEM, the problem number.
!
!    Input, integer ( kind = 4 ) DIM_NUM, the spatial dimension.
!
!    Input, integer ( kind = 4 ) POINT_NUM, the number of points.
!
!    Input, real ( kind = 8 ) X01(DIM_NUM,POINT_NUM), the points, in [0,1], 
!    to be transformed.
!
!    Output, real ( kind = 8 ) XAB(DIM_NUM,POINT_NUM), the transformed points.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) point_num

  real ( kind = 8 ) a(dim_num)
  real ( kind = 8 ) b(dim_num)
  integer ( kind = 4 ) dim
  integer ( kind = 4 ) problem
  real ( kind = 8 ) x01(dim_num,point_num)
  real ( kind = 8 ) xab(dim_num,point_num)

  call p00_lim ( problem, dim_num, a, b )

  do dim = 1, dim_num
    xab(dim,1:point_num) =  a(dim) + ( b(dim) - a(dim) ) * x01(dim,1:point_num)
  end do

  return
end
subroutine p00_title ( problem )

!*****************************************************************************80
!
!! P00_TITLE prints a title for any problem.
!
!  Modified:
!
!    26 February 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) PROBLEM, the number of the test problem.
!
  implicit none

  integer ( kind = 4 ) problem

  if ( problem == 1 ) then
    call p01_title ( )
  else if ( problem == 2 ) then
    call p02_title ( )
  else if ( problem == 3 ) then
    call p03_title ( )
  else if ( problem == 4 ) then
    call p04_title ( )
  else if ( problem == 5 ) then
    call p05_title ( )
  else if ( problem == 6 ) then
    call p06_title ( )
  else if ( problem == 7 ) then
    call p07_title ( )
  else if ( problem == 8 ) then
    call p08_title ( )
  else if ( problem == 9 ) then
    call p09_title ( )
  else if ( problem == 10 ) then
    call p10_title ( )
  else if ( problem == 11 ) then
    call p11_title ( )
  else if ( problem == 12 ) then
    call p12_title ( )
  else if ( problem == 13 ) then
    call p13_title ( )
  else if ( problem == 14 ) then
    call p14_title ( )
  else if ( problem == 15 ) then
    call p15_title ( )
  else if ( problem == 16 ) then
    call p16_title ( )
  else if ( problem == 17 ) then
    call p17_title ( )
  else if ( problem == 18 ) then
    call p18_title ( )
  else if ( problem == 19 ) then
    call p19_title ( )
  else if ( problem == 20 ) then
    call p20_title ( )
  else if ( problem == 21 ) then
    call p21_title ( )
  else if ( problem == 22 ) then
    call p22_title ( )
  else if ( problem == 23 ) then
    call p23_title ( )
  else if ( problem == 24 ) then
    call p24_title ( )
  else if ( problem == 25 ) then
    call p25_title ( )
  else if ( problem == 26 ) then
    call p26_title ( )
  else if ( problem == 27 ) then
    call p27_title ( )
  else if ( problem == 28 ) then
    call p28_title ( )
  else if ( problem == 29 ) then
    call p29_title ( )
  else if ( problem == 30 ) then
    call p30_title ( )
  else if ( problem == 31 ) then
    call p31_title ( )
  else if ( problem == 32 ) then
    call p32_title ( )
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_TITLE - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal problem number = ', problem
    stop
  end if

  return
end
subroutine p00_volume ( problem, dim_num, volume )

!*****************************************************************************80
!
!! P00_VOLUME returns the volume of the integration region.
!
!  Modified:
!
!    26 February 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) PROBLEM, the number of the test problem.
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the problem.
!
!    Output, real ( kind = 8 ) VOLUME, the volume of the integration region.
!
  implicit none

  integer ( kind = 4 ) dim_num

  real ( kind = 8 ) a(dim_num)
  real ( kind = 8 ) b(dim_num)
  integer ( kind = 4 ) problem
  real ( kind = 8 ) r
  real ( kind = 8 ) volume

  if ( 1 <= problem .and. problem <= 20 ) then

    call p00_lim ( problem, dim_num, a, b )

    volume = product ( b(1:dim_num) - a(1:dim_num) )

  else if ( problem == 21 ) then

    call sphere_unit_area_nd ( dim_num, volume )

  else if ( problem == 22 ) then

    call p22_r8 ( 'G', 'R', r )

    call sphere_volume_nd ( dim_num, r, volume )

  else if ( problem == 23 ) then

    call simplex_unit_volume_nd ( dim_num, volume )

  else if ( 24 <= problem .and. problem <= 32 ) then

    call p00_lim ( problem, dim_num, a, b )

    volume = product ( b(1:dim_num) - a(1:dim_num) )

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P00_VOLUME - Fatal error!'
    write ( *, '(a,i8)' ) '  Illegal problem number = ', problem
    stop
  end if

  return
end
subroutine p01_default ( dim_num )

!*****************************************************************************80
!
!! P01_DEFAULT sets default values for problem 01.
!
!  Discussion:
!
!    If a problem uses vector parameters, then the spatial dimension
!    DIM_NUM is needed as input, so that the vector parameters can be
!    properly dimensioned.
!
!    Every problem keeps a count of the number of function calls.  Calling
!    this routine causes that count to be reset to 0.
!
!  Modified:
!
!    10 March 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the problem.
!
  implicit none

  integer ( kind = 4 ) dim_num

  integer ( kind = 4 ) i4

  call p01_i4 ( 'D', '*', i4 )

  return
end
subroutine p01_exact ( dim_num, exact )

!*****************************************************************************80
!
!! P01_EXACT returns the exact integral for problem 01.
!
!  Modified:
!
!    20 February 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the spatial dimension.
!
!    Output, real ( kind = 8 ) EXACT, the exact value of the integral.
!
  implicit none

  real ( kind = 8 ) exact
  integer ( kind = 4 ) dim_num

  exact = real ( ( dim_num ) * ( 3 * dim_num + 1 ), kind = 8 ) / 12.0D+00

  return
end
subroutine p01_f ( dim_num, point_num, x, value )

!*****************************************************************************80
!
!! P01_F evaluates the integrand for problem 01.
!
!  Discussion:
!
!    The spatial dimension DIM_NUM is arbitrary.
!
!  Region:
!
!    0 <= X(1:DIM_NUM) <= 1
!
!  Integrand:
!
!    f(x) = ( sum ( x(1:dim_num) ) )**2
!
!  Exact Integral:
!
!    dim_num * ( 3 * dim_num + 1 ) / 12
!
!  Modified:
!
!    03 June 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Input, integer ( kind = 4 ) POINT_NUM, the number of points.
!
!    Input, real ( kind = 8 ) X(DIM_NUM,POINT_NUM), the evaluation points.
!
!    Output, real ( kind = 8 ) VALUE(POINT_NUM), the function values.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) point_num

  integer ( kind = 4 ) point
  real ( kind = 8 ) value(point_num)
  real ( kind = 8 ) x(dim_num,point_num)

  do point = 1, point_num
    value(point) = ( sum ( x(1:dim_num,point) ) )**2
  end do

  call p01_i4 ( 'I', '#', point_num )

  return
end
subroutine p01_i4 ( action, name, value )

!*****************************************************************************80
!
!! P01_I4 sets or gets I4 parameters for problem 01.
!
!  Modified:
!
!    13 April 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) ACTION,
!    'D' to set a parameter to its default value.
!    'S' to set a parameter,
!    'G' to get a parameter,
!    'I' to increment a parameter.
!
!    Input, character ( len = * ) NAME, the name of the variable.
!    '#' is the number of calls to the integrand routine.
!
!    Input/output, integer ( kind = 4 ) VALUE.
!    * If ACTION = 'I', then VALUE is an input quantity, and is the
!      new value to be added to NAME.
!    * If ACTION = 'S', then VALUE is an input quantity, and is the
!      new value to be assigned to NAME.
!    * If ACTION = 'G', then VALUE is an output quantity, and is the
!      current value of NAME.
!
  implicit none

  character ( len = * ) action
  character ( len = * ) name
  integer ( kind = 4 ), save :: calls = 0
  integer ( kind = 4 ) value

  if ( action(1:1) == 'D' .or. action(1:1) == 'd' ) then

    if ( name == '#' .or. name == '*' ) then
      calls = 0
    end if

  else if ( action(1:1) == 'G' .or. action(1:1) == 'g' ) then

    if ( name == '#' ) then
      value = calls
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P01_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'I' .or. action(1:1) == 'i' ) then

    if ( name == '#' ) then
      calls = calls + value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P01_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'S' .or. action(1:1) == 's' ) then

    if ( name == '#' ) then
      calls = value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P01_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P01_I4 - Fatal error!'
    write ( *, '(a)' ) '  Unrecognized action = "' // trim ( action ) // '".'
    stop

  end if

  return
end
subroutine p01_lim ( dim_num, a, b )

!*****************************************************************************80
!
!! P01_LIM returns the integration limits for problem 01.
!
!  Modified:
!
!    20 February 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Output, real ( kind = 8 ) A(DIM_NUM), B(DIM_NUM), the lower and upper
!    limits of integration.
!    Note that if A = -HUGE(A), the lower limit is
!    actually negative infinity, and if B = HUGE(B), the upper limit
!    is actually infinity.
!
  implicit none

  integer ( kind = 4 ) dim_num

  real ( kind = 8 ) a(dim_num)
  real ( kind = 8 ) b(dim_num)

  a(1:dim_num) = 0.0D+00
  b(1:dim_num) = 1.0D+00

  return
end
subroutine p01_name ( name )

!*****************************************************************************80
!
!! P01_NAME returns the name of problem 01.
!
!  Modified:
!
!    20 February 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) NAME, the name of the problem.
!
  implicit none

  character ( len = * ) name

  name = 'SquareSum'

  return
end
subroutine p01_region ( region )

!*****************************************************************************80
!
!! P01_REGION returns the name of the integration region for problem 01.
!
!  Modified:
!
!    20 February 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) REGION, the name of the integration region.
!
  implicit none

  character ( len = * ) region

  region = 'box'

  return
end
subroutine p01_title ( )

!*****************************************************************************80
!
!! P01_TITLE prints a title for problem 01.
!
!  Modified:
!
!    20 February 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Problem 01'
  write ( *, '(a)' ) '  Name:       SquareSum'
  write ( *, '(a)' ) '  Region:     0 <= X(i) <= 1'
  write ( *, '(a)' ) '  Integrand:  F(X) = ( sum ( X(i) ) )^2'

  return
end
subroutine p02_default ( dim_num )

!*****************************************************************************80
!
!! P02_DEFAULT sets default values for problem 02.
!
!  Discussion:
!
!    If a problem uses vector parameters, then the spatial dimension
!    DIM_NUM is needed as input, so that the vector parameters can be
!    properly dimensioned.
!
!    Every problem keeps a count of the number of function calls.  Calling
!    this routine causes that count to be reset to 0.
!
!  Modified:
!
!    10 March 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the problem.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) i4

  call p02_i4 ( 'D', '*', i4 )

  return
end
subroutine p02_exact ( dim_num, exact )

!*****************************************************************************80
!
!! P02_EXACT returns the exact integral for problem 02.
!
!  Modified:
!
!    25 August 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the spatial dimension.
!
!    Output, real ( kind = 8 ) EXACT, the exact value of the integral.
!
  implicit none

  integer ( kind = 4 ) dim_num
  real ( kind = 8 ) exact

  exact = real ( dim_num * ( 5 * dim_num - 2 ), kind = 8 ) / 15.0D+00

  return
end
subroutine p02_f ( dim_num, point_num, x, value )

!*****************************************************************************80
!
!! P02_F evaluates the integrand for problem 02.
!
!  Discussion:
!
!    The spatial dimension DIM_NUM arbitrary.
!
!  Region:
!
!    0 <= X(1:DIM_NUM) <= 1
!
!  Integrand:
!
!    ( sum ( 2 * x(1:dim_num) - 1 ) )**4
!
!  Exact Integral:
!
!    DIM_NUM * (5*DIM_NUM-2) / 15
!
!  Modified:
!
!    03 June 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Philip Davis, Philip Rabinowitz,
!    Methods of Numerical Integration,
!    Second Edition,
!    Dover, 2007,
!    ISBN: 0486453391,
!    LC: QA299.3.D28.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Input, integer ( kind = 4 ) POINT_NUM, the number of points.
!
!    Input, real ( kind = 8 ) X(DIM_NUM,POINT_NUM), the evaluation points.
!
!    Output, real ( kind = 8 ) VALUE(POINT_NUM), the function values.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) point_num

  integer ( kind = 4 ) point
  real ( kind = 8 ) value(point_num)
  real ( kind = 8 ) x(dim_num,point_num)

  do point = 1, point_num
    value(point) = ( sum ( 2.0D+00 * x(1:dim_num,point) - 1.0D+00 ) )**4
  end do

  call p02_i4 ( 'I', '#', point_num )

  return
end
subroutine p02_i4 ( action, name, value )

!*****************************************************************************80
!
!! P02_I4 sets or gets I4 parameters for problem 02.
!
!  Modified:
!
!    13 April 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) ACTION,
!    'D' to set a parameter to its default value.
!    'S' to set a parameter,
!    'G' to get a parameter,
!    'I' to increment a parameter.
!
!    Input, character ( len = * ) NAME, the name of the variable.
!    '#' is the number of calls to the integrand routine.
!
!    Input/output, integer ( kind = 4 ) VALUE.
!    * If ACTION = 'I', then VALUE is an input quantity, and is the
!      new value to be added to NAME.
!    * If ACTION = 'S', then VALUE is an input quantity, and is the
!      new value to be assigned to NAME.
!    * If ACTION = 'G', then VALUE is an output quantity, and is the
!      current value of NAME.
!
  implicit none

  character ( len = * ) action
  character ( len = * ) name
  integer ( kind = 4 ), save :: calls = 0
  integer ( kind = 4 ) value

  if ( action(1:1) == 'D' .or. action(1:1) == 'd' ) then

    if ( name == '#' .or. name == '*' ) then
      calls = 0
    end if

  else if ( action(1:1) == 'G' .or. action(1:1) == 'g' ) then

    if ( name == '#' ) then
      value = calls
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P02_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'I' .or. action(1:1) == 'i' ) then

    if ( name == '#' ) then
      calls = calls + value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P02_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'S' .or. action(1:1) == 's' ) then

    if ( name == '#' ) then
      calls = value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P02_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P02_I4 - Fatal error!'
    write ( *, '(a)' ) '  Unrecognized action = "' // trim ( action ) // '".'
    stop

  end if

  return
end
subroutine p02_lim ( dim_num, a, b )

!*****************************************************************************80
!
!! P02_LIM returns the integration limits for problem 02.
!
!  Modified:
!
!    25 August 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Output, real ( kind = 8 ) A(DIM_NUM), B(DIM_NUM), the limits
!    of integration.
!    Note that if A = -HUGE(A), the lower limit is
!    actually negative infinity, and if B = HUGE(B), the upper limit
!    is actually infinity.
!
  implicit none

  integer ( kind = 4 ) dim_num

  real ( kind = 8 ) a(dim_num)
  real ( kind = 8 ) b(dim_num)

  a(1:dim_num) =  0.0D+00
  b(1:dim_num) =  1.0D+00

  return
end
subroutine p02_name ( name )

!*****************************************************************************80
!
!! P02_NAME returns the name of problem 02.
!
!  Modified:
!
!    27 October 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) NAME, the name of the problem.
!
  implicit none

  character ( len = * ) name

  name = 'QuadSum'

  return
end
subroutine p02_region ( region )

!*****************************************************************************80
!
!! P02_REGION returns the name of the integration region for problem 02.
!
!  Modified:
!
!    27 October 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) REGION, the name of the integration region.
!
  implicit none

  character ( len = * ) region

  region = 'box'

  return
end
subroutine p02_title ( )

!*****************************************************************************80
!
!! P02_TITLE prints a title for problem 02.
!
!  Modified:
!
!    25 August 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Problem 02'
  write ( *, '(a)' ) '  Name:       QuadSum'
  write ( *, '(a)' ) '              Davis, Rabinowitz, page 370, #1.'
  write ( *, '(a)' ) '  Region:     0 <= X(i) <= 1'
  write ( *, '(a)' ) '  Integrand:  F(X) = ( sum ( 2 * X(i) - 1 ) )^4'

  return
end
subroutine p03_default ( dim_num )

!*****************************************************************************80
!
!! P03_DEFAULT sets default values for problem 03.
!
!  Discussion:
!
!    If a problem uses vector parameters, then the spatial dimension
!    DIM_NUM is needed as input, so that the vector parameters can be
!    properly dimensioned.
!
!    Every problem keeps a count of the number of function calls.  Calling
!    this routine causes that count to be reset to 0.
!
!  Modified:
!
!    10 March 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the problem.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) i4

  call p03_i4 ( 'D', '*', i4 )

  return
end
subroutine p03_exact ( dim_num, exact )

!*****************************************************************************80
!
!! P03_EXACT returns the exact integral for problem 03.
!
!  Modified:
!
!    14 November 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the spatial dimension.
!
!    Output, real ( kind = 8 ) EXACT, the exact value of the integral.
!
  implicit none

  real ( kind = 8 ) exact
  integer ( kind = 4 ) dim_num

  exact = 0.0D+00

  return
end
subroutine p03_f ( dim_num, point_num, x, value )

!*****************************************************************************80
!
!! P03_F evaluates the integrand for problem 03.
!
!  Discussion:
!
!    The spatial dimension DIM_NUM is arbitrary.
!
!  Region:
!
!    0 <= X(1:DIM_NUM) <= 1
!
!  Integrand:
!
!    ( sum ( 2 * x(1:dim_num) - 1 ) )**5
!
!  Exact Integral:
!
!    0
!
!  Modified:
!
!    03 June 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Philip Davis, Philip Rabinowitz,
!    Methods of Numerical Integration,
!    Second Edition,
!    Dover, 2007,
!    ISBN: 0486453391,
!    LC: QA299.3.D28.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Input, integer ( kind = 4 ) POINT_NUM, the number of points.
!
!    Input, real ( kind = 8 ) X(DIM_NUM,POINT_NUM), the evaluation points.
!
!    Output, real ( kind = 8 ) VALUE(POINT_NUM), the function values.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) point_num

  integer ( kind = 4 ) point
  real ( kind = 8 ) value(point_num)
  real ( kind = 8 ) x(dim_num,point_num)

  do point = 1, point_num
    value(point) = ( sum ( 2.0D+00 * x(1:dim_num,point) - 1.0D+00 ) )**5
  end do

  call p03_i4 ( 'I', '#', point_num )

  return
end
subroutine p03_i4 ( action, name, value )

!*****************************************************************************80
!
!! P03_I4 sets or gets I4 parameters for problem 03.
!
!  Modified:
!
!    13 April 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) ACTION,
!    'D' to set a parameter to its default value.
!    'S' to set a parameter,
!    'G' to get a parameter,
!    'I' to increment a parameter.
!
!    Input, character ( len = * ) NAME, the name of the variable.
!    '#' is the number of calls to the integrand routine.
!
!    Input/output, integer ( kind = 4 ) VALUE.
!    * If ACTION = 'I', then VALUE is an input quantity, and is the
!      new value to be added to NAME.
!    * If ACTION = 'S', then VALUE is an input quantity, and is the
!      new value to be assigned to NAME.
!    * If ACTION = 'G', then VALUE is an output quantity, and is the
!      current value of NAME.
!
  implicit none

  character ( len = * ) action
  character ( len = * ) name
  integer ( kind = 4 ), save :: calls = 0
  integer ( kind = 4 ) value

  if ( action(1:1) == 'D' .or. action(1:1) == 'd' ) then

    if ( name == '#' .or. name == '*' ) then
      calls = 0
    end if

  else if ( action(1:1) == 'G' .or. action(1:1) == 'g' ) then

    if ( name == '#' ) then
      value = calls
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P03_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'I' .or. action(1:1) == 'i' ) then

    if ( name == '#' ) then
      calls = calls + value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P03_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'S' .or. action(1:1) == 's' ) then

    if ( name == '#' ) then
      calls = value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P03_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P03_I4 - Fatal error!'
    write ( *, '(a)' ) '  Unrecognized action = "' // trim ( action ) // '".'
    stop

  end if

  return
end
subroutine p03_lim ( dim_num, a, b )

!*****************************************************************************80
!
!! P03_LIM returns the integration limits for problem 03.
!
!  Modified:
!
!    29 October 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Output, real ( kind = 8 ) A(DIM_NUM), B(DIM_NUM), the lower and upper
!    limits of integration.
!    Note that if A = -HUGE(A), the lower limit is
!    actually negative infinity, and if B = HUGE(B), the upper limit
!    is actually infinity.
!
  implicit none

  integer ( kind = 4 ) dim_num

  real ( kind = 8 ) a(dim_num)
  real ( kind = 8 ) b(dim_num)

  a(1:dim_num) =  0.0D+00
  b(1:dim_num) =  1.0D+00

  return
end
subroutine p03_name ( name )

!*****************************************************************************80
!
!! P03_NAME returns the name of problem 03.
!
!  Modified:
!
!    29 October 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) NAME, the name of the problem.
!
  implicit none

  character ( len = * ) name

  name = 'QuintSum'

  return
end
subroutine p03_region ( region )

!*****************************************************************************80
!
!! P03_REGION returns the name of the integration region for problem 03.
!
!  Modified:
!
!    29 October 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) REGION, the name of the integration region.
!
  implicit none

  character ( len = * ) region

  region = 'box'

  return
end
subroutine p03_title ( )

!*****************************************************************************80
!
!! P03_TITLE prints a title for problem 03.
!
!  Modified:
!
!    29 October 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Problem 03'
  write ( *, '(a)' ) '  Name:       QuintSum'
  write ( *, '(a)' ) '              Davis, Rabinowitz, page 370, #3.'
  write ( *, '(a)' ) '  Region:     0 <= X(i) <= 1'
  write ( *, '(a)' ) '  Integrand:  F(X) = ( sum ( X(i) ) )^5'

  return
end
subroutine p04_default ( dim_num )

!*****************************************************************************80
!
!! P04_DEFAULT sets default values for problem 04.
!
!  Discussion:
!
!    If a problem uses vector parameters, then the spatial dimension
!    DIM_NUM is needed as input, so that the vector parameters can be
!    properly dimensioned.
!
!    Every problem keeps a count of the number of function calls.  Calling
!    this routine causes that count to be reset to 0.
!
!  Modified:
!
!    10 March 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the problem.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) i4

  call p04_i4 ( 'D', '*', i4 )

  return
end
subroutine p04_exact ( dim_num, exact )

!*****************************************************************************80
!
!! P04_EXACT returns the exact integral for problem 04.
!
!  Modified:
!
!    25 August 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the spatial dimension.
!
!    Output, real ( kind = 8 ) EXACT, the exact value of the integral.
!
  implicit none

  integer ( kind = 4 ) dim_num
  real ( kind = 8 ) exact

  exact = real ( dim_num * ( 7 * ( dim_num - 1 ) * ( 5 * dim_num - 1 ) + 9 ), &
    kind = 8 ) / 63.0D+00

  return
end
subroutine p04_f ( dim_num, point_num, x, value )

!*****************************************************************************80
!
!! P04_F evaluates the integrand for problem 04.
!
!  Discussion:
!
!    The spatial dimension DIM_NUM is arbitrary.
!
!  Region:
!
!    0 <= X(1:DIM_NUM) <= 1
!
!  Integrand:
!
!    ( sum ( 2 * x(1:dim_num) - 1 ) )**6
!
!  Exact Integral:
!
!    DIM_NUM * ( 7 * (DIM_NUM-1) * (5*DIM_NUM-1) + 9 ) / 63
!
!  Modified:
!
!    03 June 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Philip Davis, Philip Rabinowitz,
!    Methods of Numerical Integration,
!    Second Edition,
!    Dover, 2007,
!    ISBN: 0486453391,
!    LC: QA299.3.D28.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Input, integer ( kind = 4 ) POINT_NUM, the number of points.
!
!    Input, real ( kind = 8 ) X(DIM_NUM,POINT_NUM), the evaluation points.
!
!    Output, real ( kind = 8 ) VALUE(POINT_NUM), the function values.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) point_num

  integer ( kind = 4 ) point
  real ( kind = 8 ) value(point_num)
  real ( kind = 8 ) x(dim_num,point_num)

  do point = 1, point_num
    value(point) = ( sum ( 2.0D+00 * x(1:dim_num,point) - 1.0D+00 ) )**6
  end do

  call p04_i4 ( 'I', '#', point_num )

  return
end
subroutine p04_i4 ( action, name, value )

!*****************************************************************************80
!
!! P04_I4 sets or gets I4 parameters for problem 04.
!
!  Modified:
!
!    13 April 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) ACTION,
!    'D' to set a parameter to its default value.
!    'S' to set a parameter,
!    'G' to get a parameter,
!    'I' to increment a parameter.
!
!    Input, character ( len = * ) NAME, the name of the variable.
!    '#' is the number of calls to the integrand routine.
!
!    Input/output, integer ( kind = 4 ) VALUE.
!    * If ACTION = 'I', then VALUE is an input quantity, and is the
!      new value to be added to NAME.
!    * If ACTION = 'S', then VALUE is an input quantity, and is the
!      new value to be assigned to NAME.
!    * If ACTION = 'G', then VALUE is an output quantity, and is the
!      current value of NAME.
!
  implicit none

  character ( len = * ) action
  character ( len = * ) name
  integer ( kind = 4 ), save :: calls = 0
  integer ( kind = 4 ) value

  if ( action(1:1) == 'D' .or. action(1:1) == 'd' ) then

    if ( name == '#' .or. name == '*' ) then
      calls = 0
    end if

  else if ( action(1:1) == 'G' .or. action(1:1) == 'g' ) then

    if ( name == '#' ) then
      value = calls
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P04_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'I' .or. action(1:1) == 'i' ) then

    if ( name == '#' ) then
      calls = calls + value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P04_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'S' .or. action(1:1) == 's' ) then

    if ( name == '#' ) then
      calls = value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P04_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P04_I4 - Fatal error!'
    write ( *, '(a)' ) '  Unrecognized action = "' // trim ( action ) // '".'
    stop

  end if

  return
end
subroutine p04_lim ( dim_num, a, b )

!*****************************************************************************80
!
!! P04_LIM returns the integration limits for problem 04.
!
!  Modified:
!
!    25 August 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the dimension of the argument.
!
!    Output, real ( kind = 8 ) A(DIM_NUM),B(DIM_NUM), the lower and upper
!    limits of integration.
!    Note that if A = -HUGE(A), the lower limit is
!    actually negative infinity, and if B = HUGE(B), the upper limit
!    is actually infinity.
!
  implicit none

  integer ( kind = 4 ) dim_num

  real ( kind = 8 ) a(dim_num)
  real ( kind = 8 ) b(dim_num)

  a(1:dim_num) = 0.0D+00
  b(1:dim_num) = 1.0D+00

  return
end
subroutine p04_name ( name )

!*****************************************************************************80
!
!! P04_NAME returns the name of problem 04.
!
!  Modified:
!
!    27 October 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) NAME, the name of the problem.
!
  implicit none

  character ( len = * ) name

  name = 'HexSum'

  return
end
subroutine p04_region ( region )

!*****************************************************************************80
!
!! P04_REGION returns the name of the integration region for problem 04.
!
!  Modified:
!
!    27 October 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) REGION, the name of the integration region.
!
  implicit none

  character ( len = * ) region

  region = 'box'

  return
end
subroutine p04_title ( )

!*****************************************************************************80
!
!! P04_TITLE prints a title for problem 04.
!
!  Modified:
!
!    29 October 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Problem 04'
  write ( *, '(a)' ) '  Name:       HexSum'
  write ( *, '(a)' ) '              Davis, Rabinowitz, page 370, #2.'
  write ( *, '(a)' ) '  Region:     0 <= X(i) <= 1'
  write ( *, '(a)' ) '  Integrand:  F(X) = ( sum ( 2 * X(i) - 1 ) )^6'

  return
end
subroutine p05_default ( dim_num )

!*****************************************************************************80
!
!! P05_DEFAULT sets default values for problem 05.
!
!  Discussion:
!
!    If a problem uses vector parameters, then the spatial dimension
!    DIM_NUM is needed as input, so that the vector parameters can be
!    properly dimensioned.
!
!    Every problem keeps a count of the number of function calls.  Calling
!    this routine causes that count to be reset to 0.
!
!  Modified:
!
!    10 March 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the problem.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) i4

  call p05_i4 ( 'D', '*', i4 )

  return
end
subroutine p05_exact ( dim_num, exact )

!*****************************************************************************80
!
!! P05_EXACT returns the exact integral for problem 05.
!
!  Discussion:
!
!    The exact value is given only for DIM_NUM = 1, 2, 3, 4 or 5.
!    For other cases, the value HUGE is returned instead.
!
!  Modified:
!
!    11 March 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the spatial dimension.
!
!    Output, real ( kind = 8 ) EXACT, the exact value of the integral,
!    or HUGE if the exact value is not known.
!
  implicit none

  real ( kind = 8 ) exact
  integer ( kind = 4 ) dim_num

  if ( dim_num == 1 ) then

    exact = log ( 3.0D+00 )

  else if ( dim_num == 2 ) then

    exact = 5.0D+00 * log ( 5.0D+00 ) - 6.0D+00 * log ( 3.0D+00 )

  else if ( dim_num == 3 ) then

    exact = 0.5D+00 * ( 49.0D+00 * log ( 7.0D+00 ) &
      - 75.0D+00 * log ( 5.0D+00 ) + 27.0D+00 * log ( 3.0D+00 ) )

  else if ( dim_num == 4 ) then

    exact = 225.0D+00 * log ( 3.0D+00 ) + 125.0D+00 * log ( 5.0D+00 ) &
      - 686.0D+00 * log ( 7.0D+00 ) / 3.0D+00

  else if ( dim_num == 5 ) then

    exact = ( &
      - 65205.0D+00 * log ( 3.0D+00 ) &
      - 6250.0D+00 * log ( 5.0D+00 ) &
      + 24010.0D+00 * log ( 7.0D+00 ) &
      + 14641.0D+00 * log ( 11.0D+00 ) ) / 24.0D+00

  else

    exact = huge ( exact )

  end if

  return
end
subroutine p05_f ( dim_num, point_num, x, value )

!*****************************************************************************80
!
!! P05_F evaluates the integrand for problem 05.
!
!  Discussion:
!
!    The spatial dimension DIM_NUM is arbitrary.
!
!  Region:
!
!    0 <= X(1:DIM_NUM) <= 1
!
!  Integrand:
!
!    2**DIM_NUM / ( 1 + sum ( 2 * x(1:dim_num) ) )
!
!  Exact Integral:
!
!    For DIM_NUM = 1:
!
!      ln ( 3 )
!
!    For DIM_NUM = 2:
!
!      ln ( 3125 / 729 )
!
!    For DIM_NUM = 3:
!
!      0.5 * ( 49 * ln ( 7 ) - 75 * ln ( 5 ) + 27 * ln ( 3 ) )
! 
!    For DIM_NUM = 4:
!
!      225 * ln ( 3 ) + 125 * ln ( 5 ) - 686 * ln ( 7 ) / 3
!
!    For DIM_NUM = 5:
!
!      ( -65205 * ln ( 3 ) - 6250 * ln ( 5 ) + 24010 * ln ( 7 ) 
!      + 14641 * ln ( 11 ) ) / 24
!
!  Approximate Integral:
!
!    DIM_NUM  VALUE
!
!       1  1.098612289
!       2  1.455514830
!       3  2.152142833
!       4  3.402716587
!       5  5.620255523
!
!  Modified:
!
!    03 June 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Arthur Stroud,
!    Approximate Calculation of Multiple Integrals,
!    Prentice Hall, 1971,
!    ISBN: 0130438936,
!    LC: QA311.S85.
!
!    Arthur Stroud, Don Secrest,
!    Gaussian Quadrature Formulas,
!    Prentice Hall, 1966,
!    LC: QA299.4G3S7.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Input, integer ( kind = 4 ) POINT_NUM, the number of points.
!
!    Input, real ( kind = 8 ) X(DIM_NUM,POINT_NUM), the evaluation points.
!
!    Output, real ( kind = 8 ) VALUE(POINT_NUM), the function values.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) point_num

  integer ( kind = 4 ) point
  real ( kind = 8 ) value(point_num)
  real ( kind = 8 ) x(dim_num,point_num)

  do point = 1, point_num
    value(point) = 2.0D+00**dim_num &
      / ( 1.0D+00 + sum ( 2.0D+00 * x(1:dim_num,point) ) )
  end do

  call p05_i4 ( 'I', '#', point_num )

  return
end
subroutine p05_i4 ( action, name, value )

!*****************************************************************************80
!
!! P05_I4 sets or gets I4 parameters for problem 05.
!
!  Modified:
!
!    13 April 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) ACTION,
!    'D' to set a parameter to its default value.
!    'S' to set a parameter,
!    'G' to get a parameter,
!    'I' to increment a parameter.
!
!    Input, character ( len = * ) NAME, the name of the variable.
!    '#' is the number of calls to the integrand routine.
!
!    Input/output, integer ( kind = 4 ) VALUE.
!    * If ACTION = 'I', then VALUE is an input quantity, and is the
!      new value to be added to NAME.
!    * If ACTION = 'S', then VALUE is an input quantity, and is the
!      new value to be assigned to NAME.
!    * If ACTION = 'G', then VALUE is an output quantity, and is the
!      current value of NAME.
!
  implicit none

  character ( len = * ) action
  character ( len = * ) name
  integer ( kind = 4 ), save :: calls = 0
  integer ( kind = 4 ) value

  if ( action(1:1) == 'D' .or. action(1:1) == 'd' ) then

    if ( name == '#' .or. name == '*' ) then
      calls = 0
    end if

  else if ( action(1:1) == 'G' .or. action(1:1) == 'g' ) then

    if ( name == '#' ) then
      value = calls
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P05_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'I' .or. action(1:1) == 'i' ) then

    if ( name == '#' ) then
      calls = calls + value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P05_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'S' .or. action(1:1) == 's' ) then

    if ( name == '#' ) then
      calls = value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P05_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P05_I4 - Fatal error!'
    write ( *, '(a)' ) '  Unrecognized action = "' // trim ( action ) // '".'
    stop

  end if

  return
end
subroutine p05_lim ( dim_num, a, b )

!*****************************************************************************80
!
!! P05_LIM returns the integration limits for problem 05.
!
!  Modified:
!
!    10 February 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Output, real ( kind = 8 ) A(DIM_NUM), B(DIM_NUM), the lower and upper
!    limits of integration.
!    Note that if A = -HUGE(A), the lower limit is
!    actually negative infinity, and if B = HUGE(B), the upper limit
!    is actually infinity.
!
  implicit none

  integer ( kind = 4 ) dim_num

  real ( kind = 8 ) a(dim_num)
  real ( kind = 8 ) b(dim_num)

  a(1:dim_num) =  0.0D+00
  b(1:dim_num) =  1.0D+00

  return
end
subroutine p05_name ( name )

!*****************************************************************************80
!
!! P05_NAME returns the name of problem 05.
!
!  Modified:
!
!    10 February 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) NAME, the name of the problem.
!
  implicit none

  character ( len = * ) name

  name = 'ST04'

  return
end
subroutine p05_region ( region )

!*****************************************************************************80
!
!! P05_REGION returns the name of the integration region for problem 05.
!
!  Modified:
!
!    10 February 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) REGION, the name of the integration region.
!
  implicit none

  character ( len = * ) region

  region = 'box'

  return
end
subroutine p05_title ( )

!*****************************************************************************80
!
!! P05_TITLE prints a title for problem 05.
!
!  Modified:
!
!    10 February 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Problem 05'
  write ( *, '(a)' ) '  Name:       ST04'
  write ( *, '(a)' ) '              Stroud #4, page 26.'
  write ( *, '(a)' ) '  Region:     0 <= X(i) <= 1'
  write ( *, '(a)' ) '  Integrand:  F(X) = 1 / ( 1 + sum ( 2 * X(i) ) )'

  return
end
subroutine p06_default ( dim_num )

!*****************************************************************************80
!
!! P06_DEFAULT sets default values for problem 06.
!
!  Discussion:
!
!    If a problem uses vector parameters, then the spatial dimension
!    DIM_NUM is needed as input, so that the vector parameters can be
!    properly dimensioned.
!
!    Every problem keeps a count of the number of function calls.  Calling
!    this routine causes that count to be reset to 0.
!
!  Modified:
!
!    10 March 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the problem.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) i4

  call p06_i4 ( 'D', '*', i4 )

  return
end
subroutine p06_exact ( dim_num, exact )

!*****************************************************************************80
!
!! P06_EXACT returns the exact integral for problem 06.
!
!  Modified:
!
!    16 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the spatial dimension.
!
!    Output, real ( kind = 8 ) EXACT, the exact value of the integral.
!
  implicit none

  real ( kind = 8 ) exact
  integer ( kind = 4 ) dim_num

  exact = 1.0D+00

  return
end
subroutine p06_f ( dim_num, point_num, x, value )

!*****************************************************************************80
!
!! P06_F evaluates the integrand for problem 06.
!
!  Discussion:
!
!    The spatial dimension DIM_NUM is arbitrary.
!
!  Region:
!
!    0 <= X(1:DIM_NUM) <= 1
!
!  Integrand:
!
!    product ( 2 * abs ( 2 * x(1:dim_num) - 1 ) )
!
!  Exact Integral:
!
!    1
!
!  Modified:
!
!    03 June 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Paul Bratley, Bennett Fox, Harald Niederreiter,
!    Implementation and Tests of Low-Discrepancy Sequences,
!    ACM Transactions on Modeling and Computer Simulation,
!    Volume 2, Number 3, pages 195-213, 1992.
!
!    Philip Davis, Philip Rabinowitz,
!    Methods of Numerical Integration,
!    Second Edition,
!    Dover, 2007,
!    ISBN: 0486453391,
!    LC: QA299.3.D28.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Input, integer ( kind = 4 ) POINT_NUM, the number of points.
!
!    Input, real ( kind = 8 ) X(DIM_NUM,POINT_NUM), the evaluation points.
!
!    Output, real ( kind = 8 ) VALUE(POINT_NUM), the function values.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) point_num

  integer ( kind = 4 ) point
  real ( kind = 8 ) value(point_num)
  real ( kind = 8 ) x(dim_num,point_num)

  do point = 1, point_num
    value(point) = product ( 2.0D+00 &
      * abs ( 2.0D+00 * x(1:dim_num,point) - 1.0D+00 ) )
  end do

  call p06_i4 ( 'I', '#', point_num )

  return
end
subroutine p06_i4 ( action, name, value )

!*****************************************************************************80
!
!! P06_I4 sets or gets I4 parameters for problem 06.
!
!  Modified:
!
!    13 April 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) ACTION,
!    'D' to set a parameter to its default value.
!    'S' to set a parameter,
!    'G' to get a parameter,
!    'I' to increment a parameter.
!
!    Input, character ( len = * ) NAME, the name of the variable.
!    '#' is the number of calls to the integrand routine.
!
!    Input/output, integer ( kind = 4 ) VALUE.
!    * If ACTION = 'I', then VALUE is an input quantity, and is the
!      new value to be added to NAME.
!    * If ACTION = 'S', then VALUE is an input quantity, and is the
!      new value to be assigned to NAME.
!    * If ACTION = 'G', then VALUE is an output quantity, and is the
!      current value of NAME.
!
  implicit none

  character ( len = * ) action
  character ( len = * ) name
  integer ( kind = 4 ), save :: calls = 0
  integer ( kind = 4 ) value

  if ( action(1:1) == 'D' .or. action(1:1) == 'd' ) then

    if ( name == '#' .or. name == '*' ) then
      calls = 0
    end if

  else if ( action(1:1) == 'G' .or. action(1:1) == 'g' ) then

    if ( name == '#' ) then
      value = calls
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P06_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'I' .or. action(1:1) == 'i' ) then

    if ( name == '#' ) then
      calls = calls + value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P06_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'S' .or. action(1:1) == 's' ) then

    if ( name == '#' ) then
      calls = value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P06_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P06_I4 - Fatal error!'
    write ( *, '(a)' ) '  Unrecognized action = "' // trim ( action ) // '".'
    stop

  end if

  return
end
subroutine p06_lim ( dim_num, a, b )

!*****************************************************************************80
!
!! P06_LIM returns the integration limits for problem 06.
!
!  Modified:
!
!    25 August 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Output, real ( kind = 8 ) A(DIM_NUM), B(DIM_NUM), the lower and upper
!    limits of integration.
!    Note that if A = -HUGE(A), the lower limit is
!    actually negative infinity, and if B = HUGE(B), the upper limit
!    is actually infinity.
!
  implicit none

  integer ( kind = 4 ) dim_num

  real ( kind = 8 ) a(dim_num)
  real ( kind = 8 ) b(dim_num)

  a(1:dim_num) =  0.0D+00
  b(1:dim_num) =  1.0D+00

  return
end
subroutine p06_name ( name )

!*****************************************************************************80
!
!! P06_NAME returns the name of problem 06.
!
!  Modified:
!
!    16 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) NAME, the name of the problem.
!
  implicit none

  character ( len = * ) name

  name = 'DR4061'

  return
end
subroutine p06_region ( region )

!*****************************************************************************80
!
!! P06_REGION returns the name of the integration region for problem 06.
!
!  Modified:
!
!    16 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) REGION, the name of the integration region.
!
  implicit none

  character ( len = * ) region

  region = 'box'

  return
end
subroutine p06_title ( )

!*****************************************************************************80
!
!! P06_TITLE prints a title for problem 06.
!
!  Modified:
!
!    25 August 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Problem 07'
  write ( *, '(a)' ) '  Name:       DR4061'
  write ( *, '(a)' ) '              Davis, Rabinowitz, page 406, #1.'
  write ( *, '(a)' ) '  Region:     0 <= X(i) <= 1'
  write ( *, '(a)' ) '  Integrand:  F(X) = product ( abs ( 4 * X(i) - 2 ) )'

  return
end
subroutine p07_default ( dim_num )

!*****************************************************************************80
!
!! P07_DEFAULT sets default values for problem 07.
!
!  Discussion:
!
!    If a problem uses vector parameters, then the spatial dimension
!    DIM_NUM is needed as input, so that the vector parameters can be
!    properly dimensioned.
!
!    Every problem keeps a count of the number of function calls.  Calling
!    this routine causes that count to be reset to 0.
!
!  Modified:
!
!    10 March 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the problem.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) i4

  call p07_i4 ( 'D', '*', i4 )

  return
end
subroutine p07_exact ( dim_num, exact )

!*****************************************************************************80
!
!! P07_EXACT returns the exact integral for problem 07.
!
!  Modified:
!
!    16 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the spatial dimension.
!
!    Output, real ( kind = 8 ) EXACT, the exact value of the integral.
!
  implicit none

  real ( kind = 8 ) exact
  integer ( kind = 4 ) dim_num

  exact = 1.0D+00

  return
end
subroutine p07_f ( dim_num, point_num, x, value )

!*****************************************************************************80
!
!! P07_F evaluates the integrand for problem 07.
!
!  Discussion:
!
!    The spatial dimension DIM_NUM is arbitrary.
!
!  Region:
!
!    0 <= X(1:DIM_NUM) <= 1
!
!  Integrand:
!
!    product ( pi / 2 ) * sin ( pi * x(1:dim_num) )
!
!  Exact Integral:
!
!    1
!
!  Modified:
!
!    03 June 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Philip Davis, Philip Rabinowitz,
!    Methods of Numerical Integration,
!    Second Edition,
!    Dover, 2007,
!    ISBN: 0486453391,
!    LC: QA299.3.D28.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Input, integer ( kind = 4 ) POINT_NUM, the number of points.
!
!    Input, real ( kind = 8 ) X(DIM_NUM,POINT_NUM), the evaluation points.
!
!    Output, real ( kind = 8 ) VALUE(POINT_NUM), the function values.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) point_num

  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  integer ( kind = 4 ) point
  real ( kind = 8 ) value(point_num)
  real ( kind = 8 ) x(dim_num,point_num)

  do point = 1, point_num
    value(point) = product ( 0.5D+00 * pi * sin ( pi * x(1:dim_num,point) ) )
  end do

  call p07_i4 ( 'I', '#', point_num )

  return
end
subroutine p07_i4 ( action, name, value )

!*****************************************************************************80
!
!! P07_I4 sets or gets I4 parameters for problem 07.
!
!  Modified:
!
!    13 April 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) ACTION,
!    'D' to set a parameter to its default value.
!    'S' to set a parameter,
!    'G' to get a parameter,
!    'I' to increment a parameter.
!
!    Input, character ( len = * ) NAME, the name of the variable.
!    '#' is the number of calls to the integrand routine.
!
!    Input/output, integer ( kind = 4 ) VALUE.
!    * If ACTION = 'I', then VALUE is an input quantity, and is the
!      new value to be added to NAME.
!    * If ACTION = 'S', then VALUE is an input quantity, and is the
!      new value to be assigned to NAME.
!    * If ACTION = 'G', then VALUE is an output quantity, and is the
!      current value of NAME.
!
  implicit none

  character ( len = * ) action
  character ( len = * ) name
  integer ( kind = 4 ), save :: calls = 0
  integer ( kind = 4 ) value

  if ( action(1:1) == 'D' .or. action(1:1) == 'd' ) then

    if ( name == '#' .or. name == '*' ) then
      calls = 0
    end if

  else if ( action(1:1) == 'G' .or. action(1:1) == 'g' ) then

    if ( name == '#' ) then
      value = calls
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P07_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'I' .or. action(1:1) == 'i' ) then

    if ( name == '#' ) then
      calls = calls + value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P07_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'S' .or. action(1:1) == 's' ) then

    if ( name == '#' ) then
      calls = value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P07_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P07_I4 - Fatal error!'
    write ( *, '(a)' ) '  Unrecognized action = "' // trim ( action ) // '".'
    stop

  end if

  return
end
subroutine p07_lim ( dim_num, a, b )

!*****************************************************************************80
!
!! P07_LIM returns the integration limits for problem 07.
!
!  Modified:
!
!    16 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Output, real ( kind = 8 ) A(DIM_NUM), B(DIM_NUM), the lower and upper 
!    limits of integration.
!    Note that if A = -HUGE(A), the lower limit is
!    actually negative infinity, and if B = HUGE(B), the upper limit
!    is actually infinity.
!
  implicit none

  integer ( kind = 4 ) dim_num

  real ( kind = 8 ) a(dim_num)
  real ( kind = 8 ) b(dim_num)

  a(1:dim_num) =  0.0D+00
  b(1:dim_num) =  1.0D+00

  return
end
subroutine p07_name ( name )

!*****************************************************************************80
!
!! P07_NAME returns the name of problem 07.
!
!  Modified:
!
!    16 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) NAME, the name of the problem.
!
  implicit none

  character ( len = * ) name

  name = 'DR4062'

  return
end
subroutine p07_region ( region )

!*****************************************************************************80
!
!! P07_REGION returns the name of the integration region for problem 07.
!
!
!  Modified:
!
!    16 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) REGION, the name of the integration region.
!
  implicit none

  character ( len = * ) region

  region = 'box'

  return
end
subroutine p07_title ( )

!*****************************************************************************80
!
!! P07_TITLE prints a title for problem 07.
!
!  Modified:
!
!    16 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Problem 07'
  write ( *, '(a)' ) '  Name:       DR4062'
  write ( *, '(a)' ) '              Davis, Rabinowitz, page 406, #2.'
  write ( *, '(a)' ) '  Region:     0 <= X(i) <= 1'
  write ( *, '(a)' ) '  Integrand:  F(X) = prod ( pi * sin ( pi * X(i) ) / 2 )'

  return
end
subroutine p08_default ( dim_num )

!*****************************************************************************80
!
!! P08_DEFAULT sets default values for problem 08.
!
!  Discussion:
!
!    If a problem uses vector parameters, then the spatial dimension
!    DIM_NUM is needed as input, so that the vector parameters can be
!    properly dimensioned.
!
!    Every problem keeps a count of the number of function calls.  Calling
!    this routine causes that count to be reset to 0.
!
!  Modified:
!
!    10 March 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the problem.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) i4

  call p08_i4 ( 'D', '*', i4 )

  return
end
subroutine p08_exact ( dim_num, exact )

!*****************************************************************************80
!
!! P08_EXACT returns the exact integral for problem 08.
!
!  Modified:
!
!    02 June 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the spatial dimension.
!
!    Output, real ( kind = 8 ) EXACT, the exact value of the integral.
!
  implicit none

  real ( kind = 8 ) exact
  integer ( kind = 4 ) dim_num
  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00

  exact = 0.5D+00 - sqrt ( 2.0D+00**( 3 * dim_num ) ) &
    * cos ( 0.25D+00 * pi * real ( dim_num, kind = 8 ) ) &
    / ( 2.0D+00 * pi**dim_num )

  return
end
subroutine p08_f ( dim_num, point_num, x, value )

!*****************************************************************************80
!
!! P08_F evaluates the integrand for problem 08.
!
!  Discussion:
!
!    The spatial dimension DIM_NUM is arbitrary.
!
!  Region:
!
!    0 <= X(1:DIM_NUM) <= 1
!
!  Integrand:
!
!    ( sin ( (pi/4) * sum ( x(1:dim_num) ) ) )**2
!
!  Exact Integral:
!
!    1/2 - sqrt ( 2**(3*DIM_NUM) ) * cos ( DIM_NUM * pi / 4 ) ) 
!      / ( 2 * pi**DIM_NUM )
!
!  Modified:
!
!    03 June 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Richard Crandall,
!    Projects in Scientific Computing,
!    Springer, 2000,
!    ISBN: 0387950095,
!    LC: Q183.9.C733.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Input, integer ( kind = 4 ) POINT_NUM, the number of points.
!
!    Input, real ( kind = 8 ) X(DIM_NUM,POINT_NUM), the evaluation points.
!
!    Output, real ( kind = 8 ) VALUE(POINT_NUM), the function values.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) point_num

  real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
  integer ( kind = 4 ) point
  real ( kind = 8 ) value(point_num)
  real ( kind = 8 ) x(dim_num,point_num)

  do point = 1, point_num
    value(point) = ( sin ( pi * sum ( x(1:dim_num,point) ) / 4.0D+00 ) )**2
  end do

  call p08_i4 ( 'I', '#', point_num )

  return
end
subroutine p08_i4 ( action, name, value )

!*****************************************************************************80
!
!! P08_I4 sets or gets I4 parameters for problem 08.
!
!  Modified:
!
!    13 April 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) ACTION,
!    'D' to set a parameter to its default value.
!    'S' to set a parameter,
!    'G' to get a parameter,
!    'I' to increment a parameter.
!
!    Input, character ( len = * ) NAME, the name of the variable.
!    '#' is the number of calls to the integrand routine.
!
!    Input/output, integer ( kind = 4 ) VALUE.
!    * If ACTION = 'I', then VALUE is an input quantity, and is the
!      new value to be added to NAME.
!    * If ACTION = 'S', then VALUE is an input quantity, and is the
!      new value to be assigned to NAME.
!    * If ACTION = 'G', then VALUE is an output quantity, and is the
!      current value of NAME.
!
  implicit none

  character ( len = * ) action
  character ( len = * ) name
  integer ( kind = 4 ), save :: calls = 0
  integer ( kind = 4 ) value

  if ( action(1:1) == 'D' .or. action(1:1) == 'd' ) then

    if ( name == '#' .or. name == '*' ) then
      calls = 0
    end if

  else if ( action(1:1) == 'G' .or. action(1:1) == 'g' ) then

    if ( name == '#' ) then
      value = calls
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P08_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'I' .or. action(1:1) == 'i' ) then

    if ( name == '#' ) then
      calls = calls + value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P08_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'S' .or. action(1:1) == 's' ) then

    if ( name == '#' ) then
      calls = value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P08_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P08_I4 - Fatal error!'
    write ( *, '(a)' ) '  Unrecognized action = "' // trim ( action ) // '".'
    stop

  end if

  return
end
subroutine p08_lim ( dim_num, a, b )

!*****************************************************************************80
!
!! P08_LIM returns the integration limits for problem 08.
!
!  Modified:
!
!    02 June 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Output, real ( kind = 8 ) A(DIM_NUM), B(DIM_NUM), the lower and upper 
!    limits of integration.
!    Note that if A = -HUGE(A), the lower limit is
!    actually negative infinity, and if B = HUGE(B), the upper limit
!    is actually infinity.
!
  implicit none

  integer ( kind = 4 ) dim_num

  real ( kind = 8 ) a(dim_num)
  real ( kind = 8 ) b(dim_num)

  a(1:dim_num) =  0.0D+00
  b(1:dim_num) =  1.0D+00

  return
end
subroutine p08_name ( name )

!*****************************************************************************80
!
!! P08_NAME returns the name of problem 08.
!
!  Modified:
!
!    02 June 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) NAME, the name of the problem.
!
  implicit none

  character ( len = * ) name

  name = 'RC01'

  return
end
subroutine p08_region ( region )

!*****************************************************************************80
!
!! P08_REGION returns the name of the integration region for problem 08.
!
!  Modified:
!
!    02 June 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) REGION, the name of the integration region.
!
  implicit none

  character ( len = * ) region

  region = 'box'

  return
end
subroutine p08_title ( )

!*****************************************************************************80
!
!! P08_TITLE prints a title for problem 08.
!
!  Modified:
!
!    02 June 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Problem 08'
  write ( *, '(a)' ) '  Name:       RC01'
  write ( *, '(a)' ) '              Crandall, page 49, #1'
  write ( *, '(a)' ) '  Region:     0 <= X(i) <= 1'
  write ( *, '(a)' ) '  Integrand:  F(X) = sin^2 ( pi/4 * sum ( X(i) ) )'

  return
end
subroutine p09_default ( dim_num )

!*****************************************************************************80
!
!! P09_DEFAULT sets default values for problem 09.
!
!  Discussion:
!
!    If a problem uses vector parameters, then the spatial dimension
!    DIM_NUM is needed as input, so that the vector parameters can be
!    properly dimensioned.
!
!    Every problem keeps a count of the number of function calls.  Calling
!    this routine causes that count to be reset to 0.
!
!  Modified:
!
!    10 March 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the problem.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) i4
  real ( kind = 8 ) r8vec(1)

  call p09_i4 ( 'D', '*', i4 )
  call p09_r8vec ( 'D', '*', dim_num, r8vec )

  return
end
subroutine p09_exact ( dim_num, exact )

!*****************************************************************************80
!
!! P09_EXACT returns the exact integral for problem 09.
!
!  Modified:
!
!    09 April 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the spatial dimension.
!
!    Output, real ( kind = 8 ) EXACT, the exact value of the integral.
!
  implicit none

  integer ( kind = 4 ) dim_num

  real ( kind = 8 ) c(dim_num)
  real ( kind = 8 ) exact

  call p09_r8vec ( 'G', 'C', dim_num, c )

  exact = product ( ( exp ( c(1:dim_num) ) - 1.0D+00 ) / c(1:dim_num) )

  return
end
subroutine p09_f ( dim_num, point_num, x, value )

!*****************************************************************************80
!
!! P09_F evaluates the integrand for problem 09.
!
!  Discussion:
!
!    The spatial dimension DIM_NUM is arbitrary.
!
!  Region:
!
!    0 <= X(1:DIM_NUM) <= 1
!
!  Integral Parameters:
!
!    The integral depends on a parameter vector C(1:N).
!
!    The reference suggests choosing C at random in [0,1]
!    and then multiplying by the normalizing factor (60/N).
!
!    To get or set C, call P09_R8VEC.
!
!    The default value of C(1:N) is 1/N.
!
!  Integrand:
!
!    exp ( sum ( c(1:dim_num) * x(1:dim_num) ) )
!
!  Exact Integral:
!
!    product ( ( exp ( c(1:n) - 1 ) / c(1:n) )
!
!  Modified:
!
!    03 June 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Thomas Patterson,
!    [Integral #7]
!    On the Construction of a Practical Ermakov-Zolotukhin 
!    Multiple Integrator,
!    in Numerical Integration: Recent Developments, Software
!    and Applications,
!    edited by Patrick Keast and Graeme Fairweather,
!    D. Reidel, 1987, pages 269-290,
!    LC: QA299.3.N38.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Input, integer ( kind = 4 ) POINT_NUM, the number of points.
!
!    Input, real ( kind = 8 ) X(DIM_NUM,POINT_NUM), the evaluation points.
!
!    Output, real ( kind = 8 ) VALUE(POINT_NUM), the function values.
!
  implicit none

  integer ( kind = 4 ) dim_num
  integer ( kind = 4 ) point_num

  real ( kind = 8 ) c(dim_num)
  integer ( kind = 4 ) point
  real ( kind = 8 ) value(point_num)
  real ( kind = 8 ) x(dim_num,point_num)

  call p09_r8vec ( 'G', 'C', dim_num, c )

  do point = 1, point_num
    value(point) = exp ( dot_product ( c(1:dim_num), x(1:dim_num,point) ) )
  end do

  call p09_i4 ( 'I', '#', point_num )

  return
end
subroutine p09_i4 ( action, name, value )

!*****************************************************************************80
!
!! P09_I4 sets or gets I4 parameters for problem 09.
!
!  Modified:
!
!    13 April 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) ACTION,
!    'D' to set a parameter to its default value.
!    'S' to set a parameter,
!    'G' to get a parameter,
!    'I' to increment a parameter.
!
!    Input, character ( len = * ) NAME, the name of the variable.
!    '#' is the number of calls to the integrand routine.
!
!    Input/output, integer ( kind = 4 ) VALUE.
!    * If ACTION = 'I', then VALUE is an input quantity, and is the
!      new value to be added to NAME.
!    * If ACTION = 'S', then VALUE is an input quantity, and is the
!      new value to be assigned to NAME.
!    * If ACTION = 'G', then VALUE is an output quantity, and is the
!      current value of NAME.
!
  implicit none

  character ( len = * ) action
  character ( len = * ) name
  integer ( kind = 4 ), save :: calls = 0
  integer ( kind = 4 ) value

  if ( action(1:1) == 'D' .or. action(1:1) == 'd' ) then

    if ( name == '#' .or. name == '*' ) then
      calls = 0
    end if

  else if ( action(1:1) == 'G' .or. action(1:1) == 'g' ) then

    if ( name == '#' ) then
      value = calls
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P09_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'I' .or. action(1:1) == 'i' ) then

    if ( name == '#' ) then
      calls = calls + value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P09_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'S' .or. action(1:1) == 's' ) then

    if ( name == '#' ) then
      calls = value
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P09_I4 - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P09_I4 - Fatal error!'
    write ( *, '(a)' ) '  Unrecognized action = "' // trim ( action ) // '".'
    stop

  end if

  return
end
subroutine p09_lim ( dim_num, a, b )

!*****************************************************************************80
!
!! P09_LIM returns the integration limits for problem 09.
!
!  Modified:
!
!    18 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the argument.
!
!    Output, real ( kind = 8 ) A(DIM_NUM), B(DIM_NUM), the lower and upper 
!    limits of integration.
!    Note that if A = -HUGE(A), the lower limit is
!    actually negative infinity, and if B = HUGE(B), the upper limit
!    is actually infinity.
!
  implicit none

  integer ( kind = 4 ) dim_num

  real ( kind = 8 ) a(dim_num)
  real ( kind = 8 ) b(dim_num)

  a(1:dim_num) =  0.0D+00
  b(1:dim_num) =  1.0D+00

  return
end
subroutine p09_name ( name )

!*****************************************************************************80
!
!! P09_NAME returns the name of problem 09.
!
!  Modified:
!
!    18 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) NAME, the name of the problem.
!
  implicit none

  character ( len = * ) name

  name = 'Patterson #7'

  return
end
subroutine p09_region ( region )

!*****************************************************************************80
!
!! P09_REGION returns the name of the integration region for problem 09.
!
!  Modified:
!
!    18 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, character ( len = * ) REGION, the name of the integration region.
!
  implicit none

  character ( len = * ) region

  region = 'box'

  return
end
subroutine p09_r8vec ( action, name, dim_num, value )

!*****************************************************************************80
!
!! P09_R8VEC sets or gets R8VEC parameters for problem 09.
!
!  Modified:
!
!    10 March 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) ACTION, the action.
!    'D' sets the internal value of the object to a default value.
!    If NAME = '*', then all variables are defaulted.
!    'G' means the current value of the object should be returned.
!    'R' means randomize the object and return it.
!    'S' means the input values of the object and its dimension should
!    be stored.
!
!    Input, character ( len = * ) NAME, the name of the parameter.
!    'C' is the coefficient vector.
!
!    Input, integer ( kind = 4 ) DIM_NUM, the dimension of the object.
!
!    Input/output, real ( kind = 8 ) VALUE(DIM_NUM), the value of the object.
!
  implicit none

  integer ( kind = 4 ) dim_num

  character ( len = * ) action
  real ( kind = 8 ), allocatable, save, dimension ( : ) :: c
  character ( len = * ) name
  integer ( kind = 4 ), save :: dim_num_save = 0
  real ( kind = 8 ) value(dim_num)

  if ( dim_num_save /= dim_num ) then
    dim_num_save = 0
    if ( allocated ( c ) ) then
      deallocate ( c )
    end if
  end if

  if ( dim_num_save == 0 ) then
    dim_num_save = dim_num
    allocate ( c(1:dim_num) )
  end if

  if ( action(1:1) == 'D' .or. action(1:1) == 'd' ) then

    if ( name == 'C' .or. name == 'c' .or. name == '*' ) then
      c(1:dim_num) = 1.0D+00 / real ( dim_num, kind = 8 )
    end if

  else if ( action(1:1) == 'G' .or. action(1:1) == 'g' ) then

    if ( name == 'C' .or. name == 'c' ) then
      value(1:dim_num) = c(1:dim_num)
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P09_R8VEC - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'R' .or. action(1:1) == 'r' ) then

    if ( name == 'C' .or. name == 'c' ) then
      call random_number ( harvest = c(1:dim_num) )
      c(1:dim_num) = c(1:dim_num) * 60.0D+00 / real ( dim_num, kind = 8 )
      value(1:dim_num) = c(1:dim_num)
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P09_R8VEC - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else if ( action(1:1) == 'S' .or. action(1:1) == 's' ) then

    if ( name == 'C' .or. name == 'c' ) then
      c(1:dim_num) = value(1:dim_num)
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'P09_R8VEC - Fatal error!'
      write ( *, '(a)' ) '  Unrecognized name = "' // trim ( name ) // '".'
      stop
    end if

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P09_R8VEC - Fatal error!'
    write ( *, '(a)' ) '  Unrecognized action = "' // trim ( action ) // '".'
    stop

  end if

  return
end
subroutine p09_title ( )

!*****************************************************************************80
!
!! P09_TITLE prints a title for problem 09.
!
!  Modified:
!
!    18 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Problem 09'
  write ( *, '(a)' ) '  Name:       Patterson #7'
  write ( *, '(a)' ) '  Region: