# include <cstdlib>
# include <iostream>
# include <iomanip>
# include <cmath>
# include <ctime>

using namespace std;

# include "test_nint.H"

//****************************************************************************80

char ch_cap ( char ch )

//****************************************************************************80
//
//  Purpose:
//
//    CH_CAP capitalizes a single character.
//
//  Discussion:
//
//    This routine should be equivalent to the library "toupper" function.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    19 July 1998
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char CH, the character to capitalize.
//
//    Output, char CH_CAP, the capitalized character.
//
{
  if ( 97 <= ch && ch <= 122 ) 
  {
    ch = ch - 32;
  }   

  return ch;
}
//****************************************************************************80

double csevl ( double x, double cs[], int n )

//****************************************************************************80
//
//  Purpose:
//
//    CSEVL evaluates an N term Chebyshev series.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    FORTRAN77 original version by Roger Broucke,
//    FORTRAN90 version by John Burkardt.
//
//  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, double X, the argument at which the series is 
//    to be evaluated.
//    X must satisfy -1.0 <= X <= 1.0.
//
//    Input, double CS(N), the array of N terms of a Chebyshev series.
//    In evaluating CS, only half the first coefficient is summed.
//
//    Input, int N, the number of terms in array CS.
//    N must be at least 1, and no more than 1000.
//
//    Output, double CSEVL, the value of the Chebyshev series.
//
{
  double b0;
  double b1;
  double b2;
  int i;
  double value;

  if ( n < 1 )
  {
    cout << "\n";
    cout << "CSEVL - Fatal error!\n";
    cout << "  Number of terms N is less than 1.\n";
    exit ( 1 );
  }

  if ( 1000 < n )
  {
    cout << "\n";
    cout << "CSEVL - Fatal error!\n";
    cout << "  The number of terms is more than 1000.\n";
    exit ( 1 );
  }

  if ( x < -1.0 || 1.0 < x )
  {
    cout << "\n";
    cout << "CSEVL - Fatal error!\n";
    cout << "  The input argument X is outside the interval [-1,1].\n";
    exit ( 1 );
  }

  b1 = 0.0;
  b0 = 0.0;

  for ( i = n; 1 <= i; i-- )
  {
    b2 = b1;
    b1 = b0;
    b0 = 2.0 * x * b1 - b2 + cs[i-1];
  }

  value = 0.5 * ( b0 - b2 );

  return value;
}
//****************************************************************************80

double error_f ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    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
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    FORTRAN77 original version by Wayne Fullerton,
//    FORTRAN90 version by John Burkardt.
//
//  Reference:
//
//    David Kahaner, Cleve Moler, Steven Nash,
//    Numerical Methods and Software,
//    Prentice Hall, 1989,
//    ISBN: 0-13-627258-4,
//    LC: TA345.K34.
//
//  Parameters:
//
//    Input, double X, the argument of the error function.
//
//    Output, double ERROR_F, the value of the error function at X.
//
{
  double erfcs[13] = {
    -0.049046121234691808, -0.14226120510371364, 
     0.010035582187599796, -0.000576876469976748, 
     0.000027419931252196, -0.000001104317550734, 
     0.000000038488755420, -0.000000001180858253, 
     0.000000000032334215, -0.000000000000799101, 
     0.000000000000017990, -0.000000000000000371, 
     0.000000000000000007 };
  static int nterf = 0;
  static double sqeps = 0.0;
  double sqrtpi = 1.7724538509055160;
  double value;
  static double xbig = 0.0;
  double y;
//
//  Initialize the Chebyshev series.
//
  if ( nterf == 0 )
  {
    nterf = inits ( erfcs, 13, 0.1 * r8_epsilon ( ) );
    xbig = sqrt ( - log ( sqrtpi * r8_epsilon ( ) ) );
    sqeps = sqrt ( 2.0 * r8_epsilon ( ) );
  }

  y = r8_abs ( x );

  if ( y <= sqeps )
  {
    value = 2.0 * x / sqrtpi;
  }
  else if ( y <= 1.0 )
  {
    value = x * ( 1.0 + csevl ( 2.0 * x * x - 1.0, erfcs, nterf ) );
  }
  else if ( y <= xbig )
  {
    value = ( 1.0 - error_fc ( y ) ) * r8_sign ( x );
  }
  else
  {
    value = r8_sign ( x );
  }

  return value;
}
//****************************************************************************80

double error_fc ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    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.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    FORTRAN77 original version by Wayne Fullerton,
//    FORTRAN90 version by John Burkardt.
//
//  Reference:
//
//    David Kahaner, Cleve Moler, Steven Nash,
//    Numerical Methods and Software,
//    Prentice Hall, 1989,
//    ISBN: 0-13-627258-4,
//    LC: TA345.K34.
//
//  Parameters:
//
//    Input, double X, the argument of the error function.
//
//    Output, double ERROR_FC, the value of the complementary 
//    error function at X.
//
{
  double erfcs[13] = {
   -0.049046121234691808, 
   -0.14226120510371364,  
    0.010035582187599796, 
   -0.000576876469976748, 
    0.000027419931252196, 
   -0.000001104317550734, 
    0.000000038488755420, 
   -0.000000001180858253, 
    0.000000000032334215, 
   -0.000000000000799101, 
    0.000000000000017990, 
   -0.000000000000000371, 
    0.000000000000000007 };
  double erfccs[24] = {
    0.0715179310202925, 
   -0.026532434337606719, 
    0.001711153977920853, 
   -0.000163751663458512, 
    0.000019871293500549, 
   -0.000002843712412769, 
    0.000000460616130901, 
   -0.000000082277530261, 
    0.000000015921418724, 
   -0.000000003295071356, 
    0.000000000722343973, 
   -0.000000000166485584, 
    0.000000000040103931, 
   -0.000000000010048164, 
    0.000000000002608272, 
   -0.000000000000699105, 
    0.000000000000192946, 
   -0.000000000000054704, 
    0.000000000000015901, 
   -0.000000000000004729, 
    0.000000000000001432, 
   -0.000000000000000439,
    0.000000000000000138, 
   -0.000000000000000048 };
  double erc2cs[23] = {
   -0.069601346602309501, 
   -0.041101339362620893, 
    0.003914495866689626, 
   -0.000490639565054897, 
    0.000071574790013770, 
   -0.000011530716341312, 
    0.000001994670590201, 
   -0.000000364266647159, 
    0.000000069443726100, 
   -0.000000013712209021, 
    0.000000002788389661, 
   -0.000000000581416472, 
    0.000000000123892049, 
   -0.000000000026906391, 
    0.000000000005942614, 
   -0.000000000001332386, 
    0.000000000000302804, 
   -0.000000000000069666, 
    0.000000000000016208, 
   -0.000000000000003809, 
    0.000000000000000904, 
   -0.000000000000000216, 
    0.000000000000000052 };
  double eta;
  static int nterc2 = 0;
  static int nterf = 0;
  static int nterfc = 0;
  static double sqeps = 0.0;
  double sqrtpi = 1.7724538509055160;
  double value;
  static double xmax = 0.0;
  static double xsml = 0.0;
  double y;

  if ( nterf == 0 )
  {
    eta = 0.1 * r8_epsilon ( );
    nterf = inits ( erfcs, 13, eta );
    nterfc = inits ( erfccs, 24, eta );
    nterc2 = inits ( erc2cs, 23, eta );

    xsml = -sqrt ( - log ( sqrtpi * r8_epsilon ( ) ) );
    xmax = sqrt ( - log ( sqrtpi * r8_tiny ( ) ) );
    xmax = xmax - 0.5 * log ( xmax ) / xmax - 0.01;
    sqeps = sqrt ( 2.0 * r8_epsilon ( ) );
  }

  if ( x <= xsml )
  {
    value = 2.0;
    return value;
  }
//
//  X so big that ERFC will underflow.
//
  if ( xmax < x )
  {
    value = 0.0;
    return value;
  }

  y = r8_abs ( x );
//
//  erfc(x) = 1.0 - erf(x) for -1 <= x <= 1.
//
  if ( y <= 1.0 )
  {
    if ( y < sqeps )
    {
      value = 1.0 - 2.0 * x / sqrtpi;
    }
    else if ( sqeps <= y )
    {
      value = 1.0 - x * ( 1.0 + csevl ( 2.0 * x * x - 1.0, erfcs, nterf ) );
    }
    return value;
  }
//
//  For 1 < |x| <= xmax, erfc(x) = 1.0 - erf(x)
//
  y = y * y;

  if ( y <= 4.0 )
  {
    value = exp ( -y ) / r8_abs ( x ) * ( 0.5 
      + csevl ( ( 8.0 / y - 5.0 ) / 3.0, erc2cs, nterc2 ) );
  }
  else
  {
    value = exp ( -y ) / r8_abs ( x ) * ( 0.5 
      + csevl ( 8.0 / y - 1.0, erfccs, nterfc ) );
  }

  if ( x < 0.0 )
  {
    value = 2.0 - value;
  }

  return value;
}
//****************************************************************************80

double gamma_log ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    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 reference 3, while approximations
//    for X < 12.0 are similar to those in reference 1, but are unpublished.
//
//    The accuracy achieved depends on the arithmetic system, the compiler,
//    intrinsic functions, and proper selection of the machine-dependent
//    constants.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    28 May 2003
//
//  Author:
//
//    FORTRAN77 original version by William Cody, Laura Stoltz,
//    FORTRAN90 version by John Burkardt.
//
//  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, double X, the argument of the Gamma function.  X must be positive.
//
//    Output, double GAMMA_LOG, the logarithm of the Gamma function of X.
//    If X <= 0.0, or if overflow would occur, the program returns the
//    value XINF, the largest representable double precision number.
//
//
//  Explanation of machine-dependent constants
//
//  BETA   - radix for the real number 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.
//
//  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.62E+2461
//  Cyber 180/855
//    under NOS   (S.P.)        2        1070       1.72E+319
//  IEEE (IBM/XT,
//    SUN, etc.)  (S.P.)        2         128       4.08E+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.13E+615
//  Cyber 180/855
//    under NOS   (S.P.)   6.44E+79
//  IEEE (IBM/XT,
//    SUN, etc.)  (S.P.)   1.42E+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
//
{
  double c[7] = {
    -1.910444077728E-03, 
     8.4171387781295E-04, 
    -5.952379913043012E-04, 
     7.93650793500350248E-04, 
    -2.777777777777681622553E-03, 
     8.333333333333333331554247E-02, 
     5.7083835261E-03 };
  double corr;
  double d1 = - 5.772156649015328605195174E-01;
  double d2 =   4.227843350984671393993777E-01;
  double d4 =   1.791759469228055000094023E+00;
  double frtbig = 1.42E+09;
  int i;
  double p1[8] = {
    4.945235359296727046734888E+00, 
    2.018112620856775083915565E+02, 
    2.290838373831346393026739E+03, 
    1.131967205903380828685045E+04, 
    2.855724635671635335736389E+04, 
    3.848496228443793359990269E+04, 
    2.637748787624195437963534E+04, 
    7.225813979700288197698961E+03 };
  double p2[8] = {
    4.974607845568932035012064E+00, 
    5.424138599891070494101986E+02, 
    1.550693864978364947665077E+04, 
    1.847932904445632425417223E+05, 
    1.088204769468828767498470E+06, 
    3.338152967987029735917223E+06, 
    5.106661678927352456275255E+06, 
    3.074109054850539556250927E+06 };
  double p4[8] = {
    1.474502166059939948905062E+04, 
    2.426813369486704502836312E+06, 
    1.214755574045093227939592E+08, 
    2.663432449630976949898078E+09, 
    2.940378956634553899906876E+010,
    1.702665737765398868392998E+011,
    4.926125793377430887588120E+011, 
    5.606251856223951465078242E+011 };
  double pnt68 = 0.6796875E+00;
  double q1[8] = {
    6.748212550303777196073036E+01, 
    1.113332393857199323513008E+03, 
    7.738757056935398733233834E+03, 
    2.763987074403340708898585E+04, 
    5.499310206226157329794414E+04, 
    6.161122180066002127833352E+04, 
    3.635127591501940507276287E+04, 
    8.785536302431013170870835E+03 };
  double q2[8] = {
    1.830328399370592604055942E+02, 
    7.765049321445005871323047E+03, 
    1.331903827966074194402448E+05, 
    1.136705821321969608938755E+06, 
    5.267964117437946917577538E+06, 
    1.346701454311101692290052E+07, 
    1.782736530353274213975932E+07, 
    9.533095591844353613395747E+06 };
  double q4[8] = {
    2.690530175870899333379843E+03, 
    6.393885654300092398984238E+05, 
    4.135599930241388052042842E+07, 
    1.120872109616147941376570E+09, 
    1.488613728678813811542398E+010, 
    1.016803586272438228077304E+011, 
    3.417476345507377132798597E+011, 
    4.463158187419713286462081E+011 };
  double res;
  double sqrtpi = 0.9189385332046727417803297;
  double xbig = 4.08E+36;
  double xden;
  double xm1;
  double xm2;
  double xm4;
  double xnum;
  double xsq;
//
//  Return immediately if the argument is out of range.
//
  if ( x <= 0.0 || xbig < x )
  {
    return r8_huge ( );
  }

  if ( x <= r8_epsilon ( ) )
  {
    res = -log ( x );
  }
  else if ( x <= 1.5 )
  {
    if ( x < pnt68 )
    {
      corr = - log ( x );
      xm1 = x;
    }
    else
    {
      corr = 0.0;
      xm1 = ( x - 0.5 ) - 0.5;
    }

    if ( x <= 0.5 || pnt68 <= x )
    {
      xden = 1.0;
      xnum = 0.0;

      for ( i = 0; i < 8; i++ )
      {
        xnum = xnum * xm1 + p1[i];
        xden = xden * xm1 + q1[i];
      }

      res = corr + ( xm1 * ( d1 + xm1 * ( xnum / xden ) ) );
    }
    else
    {
      xm2 = ( x - 0.5 ) - 0.5;
      xden = 1.0;
      xnum = 0.0;
      for ( i = 0; i < 8; i++ )
      {
        xnum = xnum * xm2 + p2[i];
        xden = xden * xm2 + q2[i];
      }
      res = corr + xm2 * ( d2 + xm2 * ( xnum / xden ) );
    }
  }
  else if ( x <= 4.0 )
  {
    xm2 = x - 2.0;
    xden = 1.0;
    xnum = 0.0;
    for ( i = 0; i < 8; i++ )
    {
      xnum = xnum * xm2 + p2[i];
      xden = xden * xm2 + q2[i];
    }

    res = xm2 * ( d2 + xm2 * ( xnum / xden ) );
  }
  else if ( x <= 12.0 )
  {
    xm4 = x - 4.0;
    xden = - 1.0;
    xnum = 0.0;
    for ( i = 0; i < 8; i++ )
    {
      xnum = xnum * xm4 + p4[i];
      xden = xden * xm4 + q4[i];
    }
    res = d4 + xm4 * ( xnum / xden );
  }
  else
  {
    res = 0.0;

    if ( x <= frtbig )
    {
      res = c[6];
      xsq = x * x;

      for ( i = 0; i < 6; i++ )
      {
        res = res / xsq + c[i];
      }
    }
    res = res / x;
    corr = log ( x );
    res = res + sqrtpi - 0.5 * corr;
    res = res + x * ( corr - 1.0 );
  }

  return res;
}
//****************************************************************************80

int get_problem_num ( void )

//****************************************************************************80
//
//  Purpose:
//
//    GET_PROBLEM_NUM returns the number of test integration problems.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, int PROBLEM_NUM, the number of test integration problems.
//
{
  int problem_num;

  problem_num = 32;

  return problem_num;
}
//****************************************************************************80

int i4_huge ( void )

//****************************************************************************80
//
//  Purpose:
//
//    I4_HUGE returns a "huge" I4.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    16 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, int I4_HUGE, a "huge" I4.
//
{
  return 2147483647;
}
//****************************************************************************80

int i4_power ( int i, int j )

//****************************************************************************80
//
//  Purpose:
//
//    I4_POWER returns the value of I^J.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    01 April 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int I, J, the base and the power.  J should be nonnegative.
//
//    Output, int I4_POWER, the value of I^J.
//
{
  int k;
  int value;

  if ( j < 0 )
  {
    if ( i == 1 )
    {
      value = 1;
    }
    else if ( i == 0 )
    {
      cout << "\n";
      cout << "I4_POWER - Fatal error!\n";
      cout << "  I^J requested, with I = 0 and J negative.\n";
      exit ( 1 );
    }
    else
    {
      value = 0;
    }
  }
  else if ( j == 0 )
  {
    if ( i == 0 )
    {
      cout << "\n";
      cout << "I4_POWER - Fatal error!\n";
      cout << "  I^J requested, with I = 0 and J = 0.\n";
      exit ( 1 );
    }
    else
    {
      value = 1;
    }
  }
  else if ( j == 1 )
  {
    value = i;
  }
  else
  {
    value = 1;
    for ( k = 1; k <= j; k++ )
    {
      value = value * i;
    }
  }
  return value;
}
//****************************************************************************80

void i4vec_copy ( int n, int a1[], int a2[] )

//****************************************************************************80
//
//  Purpose:
//
//    I4VEC_COPY copies an I4VEC.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    29 June 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of entries in the vectors.
//
//    Input, int A1[N], the vector to be copied.
//
//    Input, int A2[N], the copy of A1.
//
{
  int i;

  for ( i = 0; i < n; i++ )
  {
    a2[i] = a2[i];
  }
  return;
}
//****************************************************************************80

int i4vec_sum ( int n, int a[] )

//****************************************************************************80
//
//  Purpose:
//
//    I4VEC_SUM sums the entries of an I4VEC.
//
//  Example:
//
//    Input:
//
//      A = ( 1, 2, 3, 4 )
//
//    Output:
//
//      I4VEC_SUM = 10
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    26 May 1999
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of entries in the vector.
//
//    Input, int A[N], the vector to be summed.
//
//    Output, int I4VEC_SUM, the sum of the entries of A.
//
{
  int i;
  int sum;

  sum = 0;
  for ( i = 0; i < n; i++ )
  {
    sum = sum + a[i];
  }

  return sum;
}
//****************************************************************************80

int inits ( double os[], int nos, double eta )

//****************************************************************************80
//
//  Purpose:
//
//    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.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    FORTRAN77 original version by Roger Broucke,
//    FORTRAN90 version by John Burkardt.
//
//  Reference:
//
//    Roger Broucke,
//    Algorithm 446:
//    Ten Subroutines for the Manipulation of Chebyshev Series,
//    Communications of the ACM,
//    Volume 16, 1973, pages 254-256.
//
//    David Kahaner, Cleve Moler, Steven Nash,
//    Numerical Methods and Software,
//    Prentice Hall, 1989,
//    ISBN: 0-13-627258-4,
//    LC: TA345.K34.
//
//  Parameters:
//
//    Input, double OS[NOS], the coefficients in the series.
//
//    Input, int NOS, the number of coefficients.  NOS must be
//    at least 1, or an error condition arises.
//
//    Input, double ETA, the requested accuracy of the series.
//    Ordinarily, ETA will be chosen to be one-tenth machine precision.
//
//    Output, int INITS, the order of the series guaranteeing the
//    given accuracy.  However, on error, INITS will be returned
//    as a negative number.
//
{
  double err;
  int i;
  int ii;
  int value;
//
//  Fatal error.  Number of coefficients less than 1.
//  But an error message cannot be printed out because this is a function.
//
  if ( nos < 1 )
  {
    value = - i4_huge ( );
    return value;
  }

  err = 0.0;

  i = 0;

  for ( ii = 1; ii <= nos; ii++ )
  {
    i = nos + 1 - ii;
    err = err + r8_abs ( os[i-1] );

    if ( eta < err )
    {
      i = nos + 1 - ii;
      break;
    }
  }
//
//  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 )
  {
    i = - nos;
  }

  value = i;

  return value;
}
//****************************************************************************80

double normal_01_cdf_inv ( double cdf )

//****************************************************************************80
//
//  Purpose:
//
//    NORMAL_01_CDF_INV inverts the Normal 01 CDF.
//
//  Discussion:
//
//    Modified to handle cases where R = 0 would result in LOG(0).
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    05 June 2007
//
//  Author:
//
//    FORTRAN77 original version by JD Beasley, SG Springer,
//    FORTRAN90 version by John Burkardt.
//
//  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, double CDF, the value of the CDF.
//    0.0 <= CDF <= 1.0.
//
//    Input, double NORMAL_01_CDF_INV, the corresponding argument.
//
{
  double a0 =   2.50662823884;
  double a1 = -18.61500062529;
  double a2 =  41.39119773534;
  double a3 = -25.44106049637;
  double b1 =  -8.47351093090;
  double b2 =  23.08336743743;
  double b3 = -21.06224101826;
  double b4 =   3.13082909833;
  double c0 =  -2.78718931138;
  double c1 =  -2.29796479134;
  double c2 =   4.85014127135;
  double c3 =   2.32121276858;
  double d1 =   3.54388924762;
  double d2 =   1.63706781897;
  double q;
  double r;
  double x;

  if ( cdf < 0.0 || 1.0 < cdf )
  {
    cout << "\n";
    cout << "NORMAL_01_CDF_INV - Fatal error!\n";
    cout << "  CDF < 0 or 1 < CDF.\n";
    exit ( 1 );
  }

  q = cdf - 0.5;

  q = r8_min ( q,  0.5 );
  q = r8_max ( q, -0.5 );
//
//  0.08 < CDF < 0.92
//
  if ( r8_abs ( q ) <= 0.42 )
  {
    r = q * q;

    x = q * ( ( ( 
           a3   * r 
         + a2 ) * r 
         + a1 ) * r 
         + a0 ) / ( ( ( ( 
           b4   * r 
         + b3 ) * r 
         + b2 ) * r 
         + b1 ) * r + 1.0 );
  }
//
//  CDF < 0.08 or 0.92 < CDF.
//
  else
  {
    r = r8_min ( cdf, 1.0 - cdf );
    r = r8_max ( r, 0.0 );
    r = r8_min ( r, 1.0 );

    if ( r <= 0.0 )
    {
      x = r8_huge ( );
    }
    else
    {
      r = sqrt ( - log ( r ) );

      x = ( ( ( 
             c3   * r 
           + c2 ) * r 
           + c1 ) * r 
           + c0 ) / ( ( 
             d2   * r 
           + d1 ) * r + 1.0 );
    }

    if ( q < 0.0 )
    {
      x = - x;
    }
  }

  return x;
}
//****************************************************************************80

double p00_box_gl05 ( int problem, int dim_num, int sub_num )

//****************************************************************************80
//
//  Purpose:
//
//    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.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    03 June 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROBLEM, the problem number.
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Input, int SUB_NUM, the number of subdivisions in each dimension.
//
//    Output, double P00_BOX_GL05, the approximate integral.
//
{
# define NORDER 5

  double *a;
  double *a_sub;
  double *b;
  double *b_sub;
  int dim;
  int *indx;
  int j;
  int k;
  int point;
  int point_num = 1;
  double result;
  int *sub_indx;
  double *value;
  double volume;
  double w;
  double weight[NORDER] = {
    0.236926885056189087514264040720, 
    0.478628670499366468041291514836, 
    0.568888888888888888888888888889, 
    0.478628670499366468041291514836, 
    0.236926885056189087514264040720 };
  double *x;
  double xtab[NORDER] = {
    - 0.906179845938663992797626878299, 
    - 0.538469310105683091036314420700, 
      0.0, 
      0.538469310105683091036314420700, 
      0.906179845938663992797626878299 };
//
//  Get the integration limits, and the weight adjustment factor.
//
  a = new double[dim_num];
  a_sub = new double[dim_num];
  b = new double[dim_num];
  b_sub = new double[dim_num];
  indx = new int[dim_num];
  sub_indx = new int[dim_num];
  x = new double[dim_num*point_num];

  p00_lim ( problem, dim_num, a, b );
//
//  Carry out the product rule.
//
  result = 0.0;
//
//  In the outer loop, we pick a sub-box.
//
  j = 0;

  for ( ; ; )
  {
    tuple_next ( 1, sub_num, dim_num, &j, sub_indx );

    if ( j == 0 )
    {
      break;
    }

    for ( dim = 0; dim < dim_num; dim++ )
    {
      a_sub[dim] = ( ( double ) ( sub_num + 1 - sub_indx[dim] ) * a[dim] 
                   + ( double ) (         - 1 + sub_indx[dim] ) * b[dim] ) 
                   / ( double ) ( sub_num );

      b_sub[dim] = ( ( double ) ( sub_num     - sub_indx[dim] ) * a[dim] 
                   + ( double ) (               sub_indx[dim] ) * b[dim] ) 
                   / ( double ) ( sub_num );
    }
//
//  In the inner loop, we go through all the points in the
//  cross product of the quadrature rule.
//
    k = 0;

    for ( ; ; )
    {
      tuple_next ( 1, NORDER, dim_num, &k, indx );

      if ( k == 0  )
      {
        break;
      }

      w = 1.0;
      for ( dim = 0; dim < dim_num; dim++ )
      {
      w = w * weight[indx[dim]-1];
      }

      point = 0;
      for ( dim = 0; dim < dim_num; dim++ )
      {
        x[dim+point*dim_num] = 0.5 * ( 
            a_sub[dim] * ( 1.0 - xtab[indx[dim]-1] ) 
          + b_sub[dim] * ( 1.0 + xtab[indx[dim]-1] ) );
      }
      value = p00_f ( problem, dim_num, point_num, x );

      for ( point = 0; point < point_num; point++ )
      {
        result = result + w * value[point];
      }
      delete [] value;
    }
  }
//
//  Get the volume.
//
  volume = p00_volume ( problem, dim_num );

  result = result * volume / ( double ) ( i4_power ( 2 * sub_num, dim_num ) );

  delete [] a;
  delete [] a_sub;
  delete [] b;
  delete [] b_sub;
  delete [] indx;
  delete [] sub_indx;
  delete [] x;

  return result;
# undef NORDER
}
//****************************************************************************80

double p00_box_mc ( int problem, int dim_num, int point_num )

//****************************************************************************80
//
//  Purpose:
//
//    P00_BOX_MC integrates over an multi-dimensional box using Monte Carlo.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    03 June 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROBLEM, the problem number.
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Input, int POINT_NUM, the number of points to use.
//
//    Output, double P00_BOX_MC, the approximate integral.
//
{
  double *a;
  double *b;
  int dim;
  int point;
  double quad;
  double result;
  int seed = 123456789;
  double *value;
  double volume;
  double *x;
//
//  Get the volume.
//
  volume = p00_volume ( problem, dim_num );
//
//  Get the integration limits.
//
  a = new double[dim_num];
  b = new double[dim_num];

  p00_lim ( problem, dim_num, a, b );
//
//  Get the random values.
//
  x = r8vec_uniform_01 ( dim_num*point_num, &seed );
//
//  Map them to the integration retion.
//
  for ( point = 0; point < point_num; point++ )
  {
    for ( dim = 0; dim < dim_num; dim++ )
    {
      x[dim+point*dim_num] = ( a[dim] * ( 1.0 - x[dim+point*dim_num] ) 
                             + b[dim] * (     + x[dim+point*dim_num] ) );
    }
  }
//
//  Evaluate the function.
//
  value = p00_f ( problem, dim_num, point_num, x );
//
//  Estimate the integral.
//
  result = r8vec_sum ( point_num, value ) * volume / ( double ) ( point_num );

  delete [] a;
  delete [] b;
  delete [] value;
  delete [] x;

  return result;
}
//****************************************************************************80

void p00_default ( int problem, int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    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.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    10 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROBLEM, the number of the desired test problem.
//
//    Input, int DIM_NUM, the dimension of the problem.
//
{
  if ( problem == 1 )
  {
    p01_default ( dim_num );
  }
  else if ( problem == 2 )
  {
    p02_default ( dim_num );
  }
  else if ( problem == 3 )
  {
    p03_default ( dim_num );
  }
  else if ( problem == 4 )
  {
    p04_default ( dim_num );
  }
  else if ( problem == 5 )
  {
    p05_default ( dim_num );
  }
  else if ( problem == 6 )
  {
    p06_default ( dim_num );
  }
  else if ( problem == 7 )
  {
    p07_default ( dim_num );
  }
  else if ( problem == 8 )
  {
    p08_default ( dim_num );
  }
  else if ( problem == 9 )
  {
    p09_default ( dim_num );
  }
  else if ( problem == 10 )
  {
    p10_default ( dim_num );
  }
  else if ( problem == 11 )
  {
    p11_default ( dim_num );
  }
  else if ( problem == 12 )
  {
    p12_default ( dim_num );
  }
  else if ( problem == 13 )
  {
    p13_default ( dim_num );
  }
  else if ( problem == 14 )
  {
    p14_default ( dim_num );
  }
  else if ( problem == 15 )
  {
    p15_default ( dim_num );
  }
  else if ( problem == 16 )
  {
    p16_default ( dim_num );
  }
  else if ( problem == 17 )
  {
    p17_default ( dim_num );
  }
  else if ( problem == 18 )
  {
    p18_default ( dim_num );
  }
  else if ( problem == 19 )
  {
    p19_default ( dim_num );
  }
  else if ( problem == 20 )
  {
    p20_default ( dim_num );
  }
  else if ( problem == 21 )
  {
    p21_default ( dim_num );
  }
  else if ( problem == 22 )
  {
    p22_default ( dim_num );
  }
  else if ( problem == 23 )
  {
    p23_default ( dim_num );
  }
  else if ( problem == 24 )
  {
    p24_default ( dim_num );
  }
  else if ( problem == 25 )
  {
    p25_default ( dim_num );
  }
  else if ( problem == 26 )
  {
    p26_default ( dim_num );
  }
  else if ( problem == 27 )
  {
    p27_default ( dim_num );
  }
  else if ( problem == 28 )
  {
    p28_default ( dim_num );
  }
  else if ( problem == 29 )
  {
    p29_default ( dim_num );
  }
  else if ( problem == 30 )
  {
    p30_default ( dim_num );
  }
  else if ( problem == 31 )
  {
    p31_default ( dim_num );
  }
  else if ( problem == 32 )
  {
    p32_default ( dim_num );
  }
  else
  {
    cout << "\n";
    cout << "P00_DEFAULT - Fatal error!\n";
    cout << "  Illegal problem number = " << problem << "\n";
    exit ( 1 );
  }
  return;
}
//****************************************************************************80

double p00_exact ( int problem, int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    P00_EXACT returns the exact integral for any problem.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROBLEM, the number of the desired test problem.
//
//    Input, int DIM_NUM, the dimension of the problem.
//
//    Output, double P00_EXACT, the exact value of the integral.
//
{
  double exact;

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

  return exact;
}
//****************************************************************************80

double *p00_f ( int problem, int dim_num, int point_num, double x[] )

//****************************************************************************80
//
//  Purpose:
//
//    P00_F evaluates the integrand for any problem.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    03 June 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROBLEM, the number of the desired test problem.
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Input, int POINT_NUM, the number of points.
//
//    Input, double X[DIM_NUM*POINT_NUM], the evaluation points.
//
//    Output, double P00_F[POINT_NUM], the integrand values.
//
{
  double *value;

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

  return value;
}
//****************************************************************************80

void p00_i4 ( int problem, char action, char name, int *value )

//****************************************************************************80
//
//  Purpose:
//
//    P00_I4 sets or gets I4 parameters for any problem.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    09 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROBLEM, the number of the desired test problem.
//
//    Input, char ACTION,
//    'D' to set a parameter to its default value;
//    'S' to set a parameter to VALUE,
//    'G' to set VALUE to the parameter's value,
//    'I' to increment a parameter.
//
//    Input, char NAME, the name of the variable.
//    '#' is the number of calls to the integrand routine.
//
//    Input/output, int 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.
//
{
  if ( problem == 1 )
  {
    p01_i4 ( action, name, value );
  }
  else if ( problem == 2 )
  {
    p02_i4 ( action, name, value );
  }
  else if ( problem == 3 )
  {
    p03_i4 ( action, name, value );
  }
  else if ( problem == 4 )
  {
    p04_i4 ( action, name, value );
  }
  else if ( problem == 5 )
  {
    p05_i4 ( action, name, value );
  }
  else if ( problem == 6 )
  {
    p06_i4 ( action, name, value );
  }
  else if ( problem == 7 )
  {
    p07_i4 ( action, name, value );
  }
  else if ( problem == 8 )
  {
    p08_i4 ( action, name, value );
  }
  else if ( problem == 9 )
  {
    p09_i4 ( action, name, value );
  }
  else if ( problem == 10 )
  {
    p10_i4 ( action, name, value );
  }
  else if ( problem == 11 )
  {
    p11_i4 ( action, name, value );
  }
  else if ( problem == 12 )
  {
    p12_i4 ( action, name, value );
  }
  else if ( problem == 13 )
  {
    p13_i4 ( action, name, value );
  }
  else if ( problem == 14 )
  {
    p14_i4 ( action, name, value );
  }
  else if ( problem == 15 )
  {
    p15_i4 ( action, name, value );
  }
  else if ( problem == 16 )
  {
    p16_i4 ( action, name, value );
  }
  else if ( problem == 17 )
  {
    p17_i4 ( action, name, value );
  }
  else if ( problem == 18 )
  {
    p18_i4 ( action, name, value );
  }
  else if ( problem == 19 )
  {
    p19_i4 ( action, name, value );
  }
  else if ( problem == 20 )
  {
    p20_i4 ( action, name, value );
  }
  else if ( problem == 21 )
  {
    p21_i4 ( action, name, value );
  }
  else if ( problem == 22 )
  {
    p22_i4 ( action, name, value );
  }
  else if ( problem == 23 )
  {
    p23_i4 ( action, name, value );
  }
  else if ( problem == 24 )
  {
    p24_i4 ( action, name, value );
  }
  else if ( problem == 25 )
  {
    p25_i4 ( action, name, value );
  }
  else if ( problem == 26 )
  {
    p26_i4 ( action, name, value );
  }
  else if ( problem == 27 )
  {
    p27_i4 ( action, name, value );
  }
  else if ( problem == 28 )
  {
    p28_i4 ( action, name, value );
  }
  else if ( problem == 29 )
  {
    p29_i4 ( action, name, value );
  }
  else if ( problem == 30 )
  {
    p30_i4 ( action, name, value );
  }
  else if ( problem == 31 )
  {
    p31_i4 ( action, name, value );
  }
  else if ( problem == 32 )
  {
    p32_i4 ( action, name, value );
  }
  else
  {
    cout << "\n";
    cout << "P00_I4 - Fatal error!\n";
    cout << "  Illegal problem number = " << problem << "\n";
    exit ( 1 );
  }

  return;
}
//****************************************************************************80

void p00_lim ( int problem, int dim_num, double a[], double b[] )

//****************************************************************************80
//
//  Purpose:
//
//    P00_LIM returns the integration limits for any problem.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROBLEM, the number of the desired test problem.
//
//    Input, int DIM_NUM, the dimension of the argument.
//
//    Output, double 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.
//
{
  if ( problem == 1 )
  {
    p01_lim ( dim_num, a, b );
  }
  else if ( problem == 2 )
  {
    p02_lim ( dim_num, a, b );
  }
  else if ( problem == 3 )
  {
    p03_lim ( dim_num, a, b );
  }
  else if ( problem == 4 )
  {
    p04_lim ( dim_num, a, b );
  }
  else if ( problem == 5 )
  {
    p05_lim ( dim_num, a, b );
  }
  else if ( problem == 6 )
  {
    p06_lim ( dim_num, a, b );
  }
  else if ( problem == 7 )
  {
    p07_lim ( dim_num, a, b );
  }
  else if ( problem == 8 )
  {
    p08_lim ( dim_num, a, b );
  }
  else if ( problem == 9 )
  {
    p09_lim ( dim_num, a, b );
  }
  else if ( problem == 10 )
  {
    p10_lim ( dim_num, a, b );
  }
  else if ( problem == 11 )
  {
    p11_lim ( dim_num, a, b );
  }
  else if ( problem == 12 )
  {
    p12_lim ( dim_num, a, b );
  }
  else if ( problem == 13 )
  {
    p13_lim ( dim_num, a, b );
  }
  else if ( problem == 14 )
  {
    p14_lim ( dim_num, a, b );
  }
  else if ( problem == 15 )
  {
    p15_lim ( dim_num, a, b );
  }
  else if ( problem == 16 )
  {
    p16_lim ( dim_num, a, b );
  }
  else if ( problem == 17 )
  {
    p17_lim ( dim_num, a, b );
  }
  else if ( problem == 18 )
  {
    p18_lim ( dim_num, a, b );
  }
  else if ( problem == 19 )
  {
    p19_lim ( dim_num, a, b );
  }
  else if ( problem == 20 )
  {
    p20_lim ( dim_num, a, b );
  }
  else if ( problem == 21 )
  {
    p21_lim ( dim_num, a, b );
  }
  else if ( problem == 22 )
  {
    p22_lim ( dim_num, a, b );
  }
  else if ( problem == 23 )
  {
    p23_lim ( dim_num, a, b );
  }
  else if ( problem == 24 )
  {
    p24_lim ( dim_num, a, b );
  }
  else if ( problem == 25 )
  {
    p25_lim ( dim_num, a, b );
  }
  else if ( problem == 26 )
  {
    p26_lim ( dim_num, a, b );
  }
  else if ( problem == 27 )
  {
    p27_lim ( dim_num, a, b );
  }
  else if ( problem == 28 )
  {
    p28_lim ( dim_num, a, b );
  }
  else if ( problem == 29 )
  {
    p29_lim ( dim_num, a, b );
  }
  else if ( problem == 30 )
  {
    p30_lim ( dim_num, a, b );
  }
  else if ( problem == 31 )
  {
    p31_lim ( dim_num, a, b );
  }
  else if ( problem == 32 )
  {
    p32_lim ( dim_num, a, b );
  }
  else
  {
    cout << "\n";
    cout << "P00_LIM - Fatal error!\n";
    cout << "  Illegal problem number = " << problem << "\n";
    exit ( 1 );
  }
  return;
}
//****************************************************************************80

char *p00_name ( int problem )

//****************************************************************************80
//
//  Purpose:
//
//    P00_NAME returns the name of the problem.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROBLEM, the number of the desired test problem.
//
//    Output, char *P00_NAME, the name of the problem.
//
{
  char *value;

  if ( problem == 1 )
  {
    value = p01_name ( );
  }
  else if ( problem == 2 )
  {
    value = p02_name ( );
  }
  else if ( problem == 3 )
  {
    value = p03_name ( );
  }
  else if ( problem == 4 )
  {
    value = p04_name ( );
  }
  else if ( problem == 5 )
  {
    value = p05_name ( );
  }
  else if ( problem == 6 )
  {
    value = p06_name ( );
  }
  else if ( problem == 7 )
  {
    value = p07_name ( );
  }
  else if ( problem == 8 )
  {
    value = p08_name ( );
  }
  else if ( problem == 9 )
  {
    value = p09_name ( );
  }
  else if ( problem == 10 )
  {
    value = p10_name ( );
  }
  else if ( problem == 11 )
  {
    value = p11_name ( );
  }
  else if ( problem == 12 )
  {
    value = p12_name ( );
  }
  else if ( problem == 13 )
  {
    value = p13_name ( );
  }
  else if ( problem == 14 )
  {
    value = p14_name ( );
  }
  else if ( problem == 15 )
  {
    value = p15_name ( );
  }
  else if ( problem == 16 )
  {
    value = p16_name ( );
  }
  else if ( problem == 17 )
  {
    value = p17_name ( );
  }
  else if ( problem == 18 )
  {
    value = p18_name ( );
  }
  else if ( problem == 19 )
  {
    value = p19_name ( );
  }
  else if ( problem == 20 )
  {
    value = p20_name ( );
  }
  else if ( problem == 21 )
  {
    value = p21_name ( );
  }
  else if ( problem == 22 )
  {
    value = p22_name ( );
  }
  else if ( problem == 23 )
  {
    value = p23_name ( );
  }
  else if ( problem == 24 )
  {
    value = p24_name ( );
  }
  else if ( problem == 25 )
  {
    value = p25_name ( );
  }
  else if ( problem == 26 )
  {
    value = p26_name ( );
  }
  else if ( problem == 27 )
  {
    value = p27_name ( );
  }
  else if ( problem == 28 )
  {
    value = p28_name ( );
  }
  else if ( problem == 29 )
  {
    value = p29_name ( );
  }
  else if ( problem == 30 )
  {
    value = p30_name ( );
  }
  else if ( problem == 31 )
  {
    value = p31_name ( );
  }
  else if ( problem == 32 )
  {
    value = p32_name ( );
  }
  else
  {
    cout << "\n";
    cout << "P00_NAME - Fatal error!\n";
    cout << "  Illegal problem number = " << problem << "\n";
    exit ( 1 );
  }
  return value;
}
//****************************************************************************80

void p00_r8vec ( int problem, char action, char name, int dim_num, 
  double value[] )

//****************************************************************************80
//
//  Purpose:
//
//    P00_R8VEC sets or gets R8VEC parameters for any problem.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    09 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROBLEM, the number of the desired test problem.
//
//    Input, char 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 should be returned.
//    'R' means randomize the object and return it.
//    'S' means the input value of the object should be stored.
//
//    Input, char NAME, the name of the variable.
//
//    Input, int DIM_NUM, the dimension of the variable.
//
//    Input/output, double VALUE[DIM_NUM], the value of the variable.
//
{
  if ( problem == 1 )
  {
  }
  else if ( problem == 2 )
  {
  }
  else if ( problem == 3 )
  {
  }
  else if ( problem == 4 )
  {
  }
  else if ( problem == 5 )
  {
  }
  else if ( problem == 6 )
  {
  }
  else if ( problem == 7 )
  {
  }
  else if ( problem == 8 )
  {
  }
  else if ( problem == 9 )
  {
    p09_r8vec ( action, name, dim_num, value );
  }
  else if ( problem == 10 )
  {
  }
  else if ( problem == 11 )
  {
  }
  else if ( problem == 12 )
  {
  }
  else if ( problem == 13 )
  {
  }
  else if ( problem == 14 )
  {
  }
  else if ( problem == 15 )
  {
  }
  else if ( problem == 16 )
  {
    p16_r8vec ( action, name, dim_num, value );
  }
  else if ( problem == 17 )
  {
    p17_r8vec ( action, name, dim_num, value );
  }
  else if ( problem == 18 )
  {
    p18_r8vec ( action, name, dim_num, value );
  }
  else if ( problem == 19 )
  {
    p19_r8vec ( action, name, dim_num, value );
  }
  else if ( problem == 20 )
  {
  }
  else if ( problem == 21 )
  {
  }
  else if ( problem == 22 )
  {
  }
  else if ( problem == 23 )
  {
  }
  else if ( problem == 24 )
  {
    p24_r8vec ( action, name, dim_num, value );
  }
  else if ( problem == 25 )
  {
  }
  else if ( problem == 26 )
  {
    p26_r8vec ( action, name, dim_num, value );
  }
  else if ( problem == 27 )
  {
    p27_r8vec ( action, name, dim_num, value );
  }
  else if ( problem == 28 )
  {
    p28_r8vec ( action, name, dim_num, value );
  }
  else if ( problem == 29 )
  {
    p29_r8vec ( action, name, dim_num, value );
  }
  else if ( problem == 30 )
  {
    p30_r8vec ( action, name, dim_num, value );
  }
  else if ( problem == 31 )
  {
    p31_r8vec ( action, name, dim_num, value );
  }
  else if ( problem == 32 )
  {
    p32_r8vec ( action, name, dim_num, value );
  }
  else
  {
    cout << "\n";
    cout << "P00_R8VEC - Fatal error!\n";
    cout << "  Illegal problem number = " << problem << "\n";
    exit ( 1 );
  }

  return;
}
//****************************************************************************80

char *p00_region ( int problem )

//****************************************************************************80
//
//  Purpose:
//
//    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.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    09 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROBLEM, the number of the desired test problem.
//
//    Output, char *REGION, the name of the integration region.
//
{
  char *value;

  if ( problem == 1 )
  {
    value = p01_region ( );
  }
  else if ( problem == 2 )
  {
    value = p02_region ( );
  }
  else if ( problem == 3 )
  {
    value = p03_region ( );
  }
  else if ( problem == 4 )
  {
    value = p04_region ( );
  }
  else if ( problem == 5 )
  {
    value = p05_region ( );
  }
  else if ( problem == 6 )
  {
    value = p06_region ( );
  }
  else if ( problem == 7 )
  {
    value = p07_region ( );
  }
  else if ( problem == 8 )
  {
    value = p08_region ( );
  }
  else if ( problem == 9 )
  {
    value = p09_region ( );
  }
  else if ( problem == 10 )
  {
    value = p10_region ( );
  }
  else if ( problem == 11 )
  {
    value = p11_region ( );
  }
  else if ( problem == 12 )
  {
    value = p12_region ( );
  }
  else if ( problem == 13 )
  {
    value = p13_region ( );
  }
  else if ( problem == 14 )
  {
    value = p14_region ( );
  }
  else if ( problem == 15 )
  {
    value = p15_region ( );
  }
  else if ( problem == 16 )
  {
    value = p16_region ( );
  }
  else if ( problem == 17 )
  {
    value = p17_region ( );
  }
  else if ( problem == 18 )
  {
    value = p18_region ( );
  }
  else if ( problem == 19 )
  {
    value = p19_region ( );
  }
  else if ( problem == 20 )
  {
    value = p20_region ( );
  }
  else if ( problem == 21 )
  {
    value = p21_region ( );
  }
  else if ( problem == 22 )
  {
    value = p22_region ( );
  }
  else if ( problem == 23 )
  {
    value = p23_region ( );
  }
  else if ( problem == 24 )
  {
    value = p24_region ( );
  }
  else if ( problem == 25 )
  {
    value = p25_region ( );
  }
  else if ( problem == 26 )
  {
    value = p26_region ( );
  }
  else if ( problem == 27 )
  {
    value = p27_region ( );
  }
  else if ( problem == 28 )
  {
    value = p28_region ( );
  }
  else if ( problem == 29 )
  {
    value = p29_region ( );
  }
  else if ( problem == 30 )
  {
    value = p30_region ( );
  }
  else if ( problem == 31 )
  {
    value = p31_region ( );
  }
  else if ( problem == 32 )
  {
    value = p32_region ( );
  }
  else
  {
    cout << "\n";
    cout << "P00_REGION - Fatal error!\n";
    cout << "  Illegal problem number = " << problem << "\n";
    exit ( 1 );
  }

  return value;
}
//****************************************************************************80

void p00_remap01 ( int problem, int dim_num, int point_num, double x01[], 
  double xab[] )

//****************************************************************************80
//
//  Purpose:
//
//    P00_REMAP01 remaps points in [0,1]^DIM_NUM into [A(1:DIM_NUM),B(1:DIM_NUM)].
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    08 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROBLEM, the problem number.
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Input, int POINT_NUM, the number of points.
//
//    Input, double X01[DIM_NUM*N], the points, in [0,1], 
//    to be transformed.
//
//    Output, double XAB[DIM_NUM*N], the transformed points.
//
{
  double *a;
  double *b;
  int dim;
  int point;

  a = new double[dim_num];
  b = new double[dim_num];

  p00_lim ( problem, dim_num, a, b );

  for ( dim = 0; dim < dim_num; dim++ )
  {
    for ( point = 0; point < point_num; point++ )
    {
      xab[dim+point*dim_num] = 
        a[dim] + ( b[dim] - a[dim] ) * x01[dim+point*dim_num];
    }
  }

  delete [] a;
  delete [] b;

  return;
}
//****************************************************************************80

void p00_title ( int problem )

//****************************************************************************80
//
//  Purpose:
//
//    P00_TITLE prints a title for any problem.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROBLEM, the number of the desired test problem.
//
{
  if ( problem == 1 )
  {
    p01_title ( );
  }
  else if ( problem == 2 )
  {
    p02_title ( );
  }
  else if ( problem == 3 )
  {
    p03_title ( );
  }
  else if ( problem == 4 )
  {
    p04_title ( );
  }
  else if ( problem == 5 )
  {
    p05_title ( );
  }
  else if ( problem == 6 )
  {
    p06_title ( );
  }
  else if ( problem == 7 )
  {
    p07_title ( );
  }
  else if ( problem == 8 )
  {
    p08_title ( );
  }
  else if ( problem == 9 )
  {
    p09_title ( );
  }
  else if ( problem == 10 )
  {
    p10_title ( );
  }
  else if ( problem == 11 )
  {
    p11_title ( );
  }
  else if ( problem == 12 )
  {
    p12_title ( );
  }
  else if ( problem == 13 )
  {
    p13_title ( );
  }
  else if ( problem == 14 )
  {
    p14_title ( );
  }
  else if ( problem == 15 )
  {
    p15_title ( );
  }
  else if ( problem == 16 )
  {
    p16_title ( );
  }
  else if ( problem == 17 )
  {
    p17_title ( );
  }
  else if ( problem == 18 )
  {
    p18_title ( );
  }
  else if ( problem == 19 )
  {
    p19_title ( );
  }
  else if ( problem == 20 )
  {
    p20_title ( );
  }
  else if ( problem == 21 )
  {
    p21_title ( );
  }
  else if ( problem == 22 )
  {
    p22_title ( );
  }
  else if ( problem == 23 )
  {
    p23_title ( );
  }
  else if ( problem == 24 )
  {
    p24_title ( );
  }
  else if ( problem == 25 )
  {
    p25_title ( );
  }
  else if ( problem == 26 )
  {
    p26_title ( );
  }
  else if ( problem == 27 )
  {
    p27_title ( );
  }
  else if ( problem == 28 )
  {
    p28_title ( );
  }
  else if ( problem == 29 )
  {
    p29_title ( );
  }
  else if ( problem == 30 )
  {
    p30_title ( );
  }
  else if ( problem == 31 )
  {
    p31_title ( );
  }
  else if ( problem == 32 )
  {
    p32_title ( );
  }
  else
  {
    cout << "\n";
    cout << "P00_TITLE - Fatal error!\n";
    cout << "  Illegal problem number = " << problem << "\n";
    exit ( 1 );
  }
  return;
}
//****************************************************************************80

double p00_volume ( int problem, int dim_num  )

//****************************************************************************80
//
//  Purpose:
//
//    P00_VOLUME returns the volume of the integration region.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROBLEM, the number of the desired test problem.
//
//    Input, int DIM_NUM, the dimension of the problem.
//
//    Output, double P00_VOLUME, the volume of the integration region.
//
{
  double *a;
  double *b;
  int dim;
  double r;
  double volume;

  if ( 1 <= problem && problem <= 20 )
  {
    a = new double[dim_num];
    b = new double[dim_num];

    p00_lim ( problem, dim_num, a, b );

    volume = 1.0;
    for ( dim = 0; dim < dim_num; dim++ )
    {
      volume = volume * ( b[dim] - a[dim] );
    }

    delete [] a;
    delete [] b;
  }
  else if ( problem == 21 )
  {
    volume = sphere_unit_area_nd ( dim_num );
  }
  else if ( problem == 22 )
  {
    p22_r8 ( 'G', 'R', &r );

    volume = sphere_volume_nd ( dim_num, r );
  }
  else if ( problem == 23 )
  {
    volume = simplex_unit_volume_nd ( dim_num );
  }
  else if ( 24 <= problem && problem <= 32 )
  {
    a = new double[dim_num];
    b = new double[dim_num];

    p00_lim ( problem, dim_num, a, b );

    volume = 1.0;
    for ( dim = 0; dim < dim_num; dim++ )
    {
      volume = volume * ( b[dim] - a[dim] );
    }

    delete [] a;
    delete [] b;
  }
  else
  {
    cout << "\n";
    cout << "P00_VOLUME - Fatal error!\n";
    cout << "  Illegal value of PROBLEM.\n";
    exit ( 1 );
  }
  return volume;
}
//****************************************************************************80

void p01_default ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    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.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the problem.
//
{
  int i4;

  p01_i4 ( 'D', '*', &i4 );

  return;
}
//****************************************************************************80

double p01_exact ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    P01_EXACT returns the exact integral for problem 01.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Output, double P01_EXACT, the exact value of the integral.
//
{
  double exact;

  exact = ( double ) ( ( dim_num ) * ( 3 * dim_num + 1 ) ) / 12.0;

  return exact;
}
//****************************************************************************80

double *p01_f ( int dim_num, int point_num, double x[] )

//****************************************************************************80
//
//  Purpose:
//
//    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
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    03 June 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Input, int POINT_NUM, the number of points.
//
//    Input, double X[DIM_NUM*POINT_NUM], the evaluation points.
//
//    Output, double P01_F[POINT_NUM], the integrand values.
//
{
  int dim;
  int inc;
  int point;
  double *value;

  value = new double[point_num];

  for ( point = 0; point < point_num; point++ )
  {
    value[point] = 0.0;
    for ( dim = 0; dim < dim_num; dim++ )
    {
      value[point] = value[point] + x[dim+point*dim_num];
    }
    value[point] = pow ( value[point], 2 );
  }

  inc = point_num;
  p01_i4 ( 'I', '#', &inc );

  return value;
}
//****************************************************************************80

void p01_i4 ( char action, char name, int *value )

//****************************************************************************80
//
//  Purpose:
//
//    P01_I4 sets or gets I4 parameters for problem 01.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    09 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char ACTION,
//    'D' to set a parameter to its default value;
//    'S' to set a parameter to VALUE,
//    'G' to set VALUE to the parameter's value,
//    'I' to increment a parameter.
//
//    Input, char NAME, the name of the variable.
//    '#' is the number of calls to the integrand routine.
//
//    Input/output, int *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.
//
{
  static int calls = 0;

  if ( action == 'D' )
  {
    if ( name == '#' || name == '*' )
    {
      calls = 0;
    }
  }
  else if ( action == 'G' )
  {
    if ( name == '#' )
    {
      *value = calls;
    }
    else
    {
      cout << "\n";
      cout << "P01_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else if ( action == 'I' )
  {
    if ( name == '#' )
    {
      calls = calls + *value;
    }
    else
    {
      cout << "\n";
      cout << "P01_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else if ( action == 'S' )
  {
    if ( name == '#' )
    {
      calls = *value;
    }
    else
    {
      cout << "\n";
      cout << "P01_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else
  {
    cout << "\n";
    cout << "P01_I4 - Fatal error!\n";
    cout << "  Unrecognized action = \"" << action << "\".\n";
    exit ( 1 );
  }
  return;
}
//****************************************************************************80

void p01_lim ( int dim_num, double a[], double b[] )

//****************************************************************************80
//
//  Purpose:
//
//    P01_LIM returns the integration limits for problem 01.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the argument.
//
//    Output, double 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.
//
{
  int dim;

  for ( dim = 0; dim < dim_num; dim++ )
  {
    a[dim] = 0.0;
  }
  for ( dim = 0; dim < dim_num; dim++ )
  {
    b[dim] = 1.0;
  }

  return;
}
//****************************************************************************80

char *p01_name ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P01_NAME returns the name of problem 01.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *P01_NAME, the name of the problem.
//
{
  char *value;

  value = new char[10];

  strcpy ( value, "SquareSum" );

  return value;
}
//****************************************************************************80

char *p01_region ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P01_REGION returns the name of the integration region for problem 01.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *P01_REGION, the name of the integration region.
//
{
  char *value;

  value = new char[4];

  strcpy ( value, "BOX" );

  return value;
}
//****************************************************************************80

void p01_title ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P01_TITLE prints a title for problem 01.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
  cout << "\n";
  cout << "Problem 01\n";
  cout << "  Name:       SquareSum\n";
  cout << "  Region:     0 <= X(i) <= 1\n";
  cout << "  Integrand:  F(X) = ( sum ( X(i) ) )^2\n";

  return;
}
//****************************************************************************80

void p02_default ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    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.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the problem.
//
{
  int i4;

  p02_i4 ( 'D', '*', &i4 );

  return;
}
//****************************************************************************80

double p02_exact ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    P02_EXACT returns the exact integral for problem 02.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Output, double P02_EXACT, the exact value of the integral.
//
{
  double exact;

  exact = ( double ) ( dim_num * ( 5 * dim_num - 2 ) ) / 15.0;

  return exact;
}
//****************************************************************************80

double *p02_f ( int dim_num, int point_num, double x[] )

//****************************************************************************80
//
//  Purpose:
//
//    P02_F evaluates the integrand for problem 02.
//
//  Discussion:
//
//    The spatial dimension DIM_NUM is 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
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  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, int DIM_NUM, the spatial dimension.
//
//    Input, int POINT_NUM, the number of points.
//
//    Input, double X[DIM_NUM*POINT_NUM], the evaluation points.
//
//    Output, double P02_F[POINT_NUM], the integrand values.
//
{
  int dim;
  int inc;
  int point;
  double *value;

  value = new double[point_num];

  for ( point = 0; point < point_num; point++ )
  {
    value[point] = 0.0;
    for ( dim = 0; dim < dim_num; dim++ )
    {
      value[point] = value[point] + ( 2.0 * x[dim+point*dim_num] - 1.0 );
    }
    value[point] = pow ( value[point], 4 );
  }

  inc = point_num;
  p02_i4 ( 'I', '#', &inc );

  return value;
}
//****************************************************************************80

void p02_i4 ( char action, char name, int *value )

//****************************************************************************80
//
//  Purpose:
//
//    P02_I4 sets or gets I4 parameters for problem 02.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    09 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char ACTION,
//    'D' to set a parmater to its default value;
//    'S' to set a parameter to VALUE,
//    'G' to set VALUE to the parameter's value,
//    'I' to increment a parameter.
//
//    Input, char NAME, the name of the variable.
//    '#' is the number of calls to the integrand routine.
//
//    Input/output, int *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.
//
{
  static int calls = 0;

  if ( action == 'D' )
  {
    if ( name == '#' || name == '*' )
    {
      calls = 0;
    }
  }
  else if ( action == 'G' )
  {
    if ( name == '#' )
    {
      *value = calls;
    }
    else
    {
      cout << "\n";
      cout << "P02_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else if ( action == 'I' )
  {
    if ( name == '#' )
    {
      calls = calls + *value;
    }
    else
    {
      cout << "\n";
      cout << "P02_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else if ( action == 'S' )
  {
    if ( name == '#' )
    {
      calls = *value;
    }
    else
    {
      cout << "\n";
      cout << "P02_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else
  {
    cout << "\n";
    cout << "P02_I4 - Fatal error!\n";
    cout << "  Unrecognized action = \"" << action << "\".\n";
    exit ( 1 );
  }
  return;
}
//****************************************************************************80

void p02_lim ( int dim_num, double a[], double b[] )

//****************************************************************************80
//
//  Purpose:
//
//    P02_LIM returns the integration limits for problem 02.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the argument.
//
//    Output, double 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.
//
{
  int dim;

  for ( dim = 0; dim < dim_num; dim++ )
  {
    a[dim] = 0.0;
  }
  for ( dim = 0; dim < dim_num; dim++ )
  {
    b[dim] = 1.0;
  }

  return;
}
//****************************************************************************80

char *p02_name ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P02_NAME returns the name of problem 02.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *P02_NAME, the name of the problem.
//
{
  char *value;

  value = new char[8];

  strcpy ( value, "QuadSum" );

  return value;
}
//****************************************************************************80

char *p02_region ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P02_REGION returns the name of the integration region for problem 02.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *P02_REGION, the name of the integration region.
//
{
  char *value;

  value = new char[4];

  strcpy ( value, "BOX" );

  return value;
}
//****************************************************************************80

void p02_title ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P02_TITLE prints a title for problem 02.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
  cout << "\n";
  cout << "Problem 02\n";
  cout << "  Name:       QuadSum\n";
  cout << "              Davis, Rabinowitz, page 370, #1.\n";
  cout << "  Region:     0 <= X(i) <= 1\n";
  cout << "  Integrand:  F(X) = ( sum ( 2 * X(i) - 1 ) )^4\n";

  return;
}
//****************************************************************************80

void p03_default ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    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.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the problem.
//
{
  int i4;

  p03_i4 ( 'D', '*', &i4 );

  return;
}
//****************************************************************************80

double p03_exact ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    P03_EXACT returns the exact integral for problem 03.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Output, double P03_EXACT, the exact value of the integral.
//
{
  double exact;

  exact = 0.0;

  return exact;
}
//****************************************************************************80

double *p03_f ( int dim_num, int point_num, double x[] )

//****************************************************************************80
//
//  Purpose:
//
//    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
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  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, int DIM_NUM, the spatial dimension.
//
//    Input, int POINT_NUM, the number of points.
//
//    Input, double X[DIM_NUM*POINT_NUM], the evaluation points.
//
//    Output, double P03_F[POINT_NUM], the integrand values.
//
{
  int dim;
  int inc;
  int point;
  double *value;

  value = new double[point_num];

  for ( point = 0; point < point_num; point++ )
  {
    value[point] = 0.0;
    for ( dim = 0; dim < dim_num; dim++ )
    {
      value[point] = value[point] + ( 2.0 * x[dim+point*dim_num] - 1.0 );
    }
    value[point] = pow ( value[point], 5 );
  }

  inc = point_num;
  p03_i4 ( 'I', '#', &inc );

  return value;
}
//****************************************************************************80

void p03_i4 ( char action, char name, int *value )

//****************************************************************************80
//
//  Purpose:
//
//    P03_I4 sets or gets I4 parameters for problem 03.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    09 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char ACTION,
//    'D' to set a parmater to its default value;
//    'S' to set a parameter to VALUE,
//    'G' to set VALUE to the parameter's value,
//    'I' to increment a parameter.
//
//    Input, char NAME, the name of the variable.
//    '#' is the number of calls to the integrand routine.
//
//    Input/output, int *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.
//
{
  static int calls = 0;

  if ( action == 'D' )
  {
    if ( name == '#' || name == '*' )
    {
      calls = 0;
    }
  }
  else if ( action == 'G' )
  {
    if ( name == '#' )
    {
      *value = calls;
    }
    else
    {
      cout << "\n";
      cout << "P03_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else if ( action == 'I' )
  {
    if ( name == '#' )
    {
      calls = calls + *value;
    }
    else
    {
      cout << "\n";
      cout << "P03_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else if ( action == 'S' )
  {
    if ( name == '#' )
    {
      calls = *value;
    }
    else
    {
      cout << "\n";
      cout << "P03_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else
  {
    cout << "\n";
    cout << "P03_I4 - Fatal error!\n";
    cout << "  Unrecognized action = \"" << action << "\".\n";
    exit ( 1 );
  }
  return;
}
//****************************************************************************80

void p03_lim ( int dim_num, double a[], double b[] )

//****************************************************************************80
//
//  Purpose:
//
//    P03_LIM returns the integration limits for problem 03.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the argument.
//
//    Output, double 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.
//
{
  int dim;

  for ( dim = 0; dim < dim_num; dim++ )
  {
    a[dim] = 0.0;
  }
  for ( dim = 0; dim < dim_num; dim++ )
  {
    b[dim] = 1.0;
  }

  return;
}
//****************************************************************************80

char *p03_name ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P03_NAME returns the name of problem 03.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *P03_NAME, the name of the problem.
//
{
  char *value;

  value = new char[9];

  strcpy ( value, "QuintSum" );

  return value;
}
//****************************************************************************80

char *p03_region ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P03_REGION returns the name of the integration region for problem 03.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *P03_REGION, the name of the integration region.
//
{
  char *value;

  value = new char[4];

  strcpy ( value, "BOX" );

  return value;
}
//****************************************************************************80

void p03_title ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P03_TITLE prints a title for problem 03.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
  cout << "\n";
  cout << "Problem 03\n";
  cout << "  Name:       QuintSum\n";
  cout << "              Davis, Rabinowitz, page 370, #3.\n";
  cout << "  Region:     0 <= X(i) <= 1\n";
  cout << "  Integrand:  F(X) = ( sum ( X(i) ) )^5\n";

  return;
}
//****************************************************************************80

void p04_default ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    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.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the problem.
//
{
  int i4;

  p04_i4 ( 'D', '*', &i4 );

  return;
}
//****************************************************************************80

double p04_exact ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    P04_EXACT returns the exact integral for problem 04.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Output, double P04_EXACT, the exact value of the integral.
//
{
  double exact;

  exact = ( double ) ( dim_num * ( 7 * ( dim_num - 1 ) 
    * ( 5 * dim_num - 1 ) + 9 ) ) / 63.0;

  return exact;
}
//****************************************************************************80

double *p04_f ( int dim_num, int point_num, double x[] )

//****************************************************************************80
//
//  Purpose:
//
//    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
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  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, int DIM_NUM, the spatial dimension.
//
//    Input, int POINT_NUM, the number of points.
//
//    Input, double X[DIM_NUM*POINT_NUM], the evaluation points.
//
//    Output, double P04_F[POINT_NUM], the integrand values.
//
{
  int dim;
  int inc;
  int point;
  double *value;

  value = new double[point_num];

  for ( point = 0; point < point_num; point++ )
  {
    value[point] = 0.0;
    for ( dim = 0; dim < dim_num; dim++ )
    {
      value[point] = value[point] + ( 2.0 * x[dim+point*dim_num] - 1.0 );
    }
    value[point] = pow ( value[point], 6 );
  }

  inc = point_num;
  p04_i4 ( 'I', '#', &inc );

  return value;
}
//****************************************************************************80

void p04_i4 ( char action, char name, int *value )

//****************************************************************************80
//
//  Purpose:
//
//    P04_I4 sets or gets I4 parameters for problem 04.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    09 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char ACTION,
//    'D' to set a parmater to its default value;
//    'S' to set a parameter to VALUE,
//    'G' to set VALUE to the parameter's value,
//    'I' to increment a parameter.
//
//    Input, char NAME, the name of the variable.
//    '#' is the number of calls to the integrand routine.
//
//    Input/output, int *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.
//
{
  static int calls = 0;

  if ( action == 'D' )
  {
    if ( name == '#' || name == '*' )
    {
      calls = 0;
    }
  }
  else if ( action == 'G' )
  {
    if ( name == '#' )
    {
      *value = calls;
    }
    else
    {
      cout << "\n";
      cout << "P04_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else if ( action == 'I' )
  {
    if ( name == '#' )
    {
      calls = calls + *value;
    }
    else
    {
      cout << "\n";
      cout << "P04_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else if ( action == 'S' )
  {
    if ( name == '#' )
    {
      calls = *value;
    }
    else
    {
      cout << "\n";
      cout << "P04_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else
  {
    cout << "\n";
    cout << "P04_I4 - Fatal error!\n";
    cout << "  Unrecognized action = \"" << action << "\".\n";
    exit ( 1 );
  }
  return;
}
//****************************************************************************80

void p04_lim ( int dim_num, double a[], double b[] )

//****************************************************************************80
//
//  Purpose:
//
//    P04_LIM returns the integration limits for problem 04.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the dimension of the argument.
//
//    Output, double 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.
//
{
  int dim;

  for ( dim = 0; dim < dim_num; dim++ )
  {
    a[dim] = 0.0;
  }
  for ( dim = 0; dim < dim_num; dim++ )
  {
    b[dim] = 1.0;
  }

  return;
}
//****************************************************************************80

char *p04_name ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P04_NAME returns the name of problem 04.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *P04_NAME, the name of the problem.
//
{
  char *value;

  value = new char[7];

  strcpy ( value, "HexSum" );

  return value;
}
//****************************************************************************80

char *p04_region ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P04_REGION returns the name of the integration region for problem 04.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *P04_REGION, the name of the integration region.
//
{
  char *value;

  value = new char[4];

  strcpy ( value, "BOX" );

  return value;
}
//****************************************************************************80

void p04_title ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P04_TITLE prints a title for problem 04.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
  cout << "\n";
  cout << "Problem 04\n";
  cout << "  Name:       HexSum\n";
  cout << "              Davis, Rabinowitz, page 370, #2.\n";
  cout << "  Region:     0 <= X(i) <= 1\n";
  cout << "  Integrand:  F(X) = ( sum ( 2 * X(i) - 1 ) )^6\n";

  return;
}
//****************************************************************************80

void p05_default ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    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.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the problem.
//
{
  int i4;

  p05_i4 ( 'D', '*', &i4 );

  return;
}
//****************************************************************************80

double p05_exact ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    P05_EXACT returns the exact integral for problem 05.
//
//  Discussion:
//
//    The exact value is given only for DIM_NUM = 1, 2, 3, 4 and 5.
//    For all other dimensions, the value R8_HUGE is returned.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    12 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Output, double P05_EXACT, the exact value of the integral.
//
{
  double exact;

  if ( dim_num == 1 )
  {
    exact = log ( 3.0 );
  }
  else if ( dim_num == 2 )
  {
    exact = 5.0 * log ( 5.0 ) - 6.0 * log ( 3.0 );
  }
  else if ( dim_num == 3 )
  {
    exact = 0.5 * ( 49.0 * log ( 7.0 ) 
      - 75.0 * log ( 5.0 ) + 27.0 * log ( 3.0 ) );
  }
  else if ( dim_num == 4 )
  {
    exact = 225.0 * log ( 3.0 ) + 125.0 * log ( 5.0 ) 
      - 686.0 * log ( 7.0 ) / 3.0;
  }
  else if ( dim_num == 5 )
  {
    exact = ( 
      - 65205.0 * log ( 3.0 ) 
      - 6250.0 * log ( 5.0 ) 
      + 24010.0 * log ( 7.0 ) 
      + 14641.0 * log ( 11.0 ) ) / 24.0;
  }
  else
  {
    exact = r8_huge ( );
  }
  return exact;
}
//****************************************************************************80

double *p05_f ( int dim_num, int point_num, double x[] )

//****************************************************************************80
//
//  Purpose:
//
//    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
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  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, int DIM_NUM, the spatial dimension.
//
//    Input, int POINT_NUM, the number of points.
//
//    Input, double X[DIM_NUM*POINT_NUM], the evaluation points.
//
//    Output, double P05_F[POINT_NUM], the integrand values.
//
{
  double denom;
  int dim;
  int inc;
  int point;
  double *value;

  value = new double[point_num];

  for ( point = 0; point < point_num; point++ )
  {
    denom = 1.0;
    for ( dim = 0; dim < dim_num; dim++ )
    {
      denom = denom + 2.0 * x[dim+point*dim_num];
    }
    value[point] = pow ( 2.0, dim_num ) / denom;
  }

  inc = point_num;
  p05_i4 ( 'I', '#', &inc );

  return value;
}
//****************************************************************************80

void p05_i4 ( char action, char name, int *value )

//****************************************************************************80
//
//  Purpose:
//
//    P05_I4 sets or gets I4 parameters for problem 05.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    09 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char ACTION,
//    'D' to set a parmater to its default value;
//    'S' to set a parameter to VALUE,
//    'G' to set VALUE to the parameter's value,
//    'I' to increment a parameter.
//
//    Input, char NAME, the name of the variable.
//    '#' is the number of calls to the integrand routine.
//
//    Input/output, int *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.
//
{
  static int calls = 0;

  if ( action == 'D' )
  {
    if ( name == '#' || name == '*' )
    {
      calls = 0;
    }
  }
  else if ( action == 'G' )
  {
    if ( name == '#' )
    {
      *value = calls;
    }
    else
    {
      cout << "\n";
      cout << "P05_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else if ( action == 'I' )
  {
    if ( name == '#' )
    {
      calls = calls + *value;
    }
    else
    {
      cout << "\n";
      cout << "P05_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else if ( action == 'S' )
  {
    if ( name == '#' )
    {
      calls = *value;
    }
    else
    {
      cout << "\n";
      cout << "P05_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else
  {
    cout << "\n";
    cout << "P05_I4 - Fatal error!\n";
    cout << "  Unrecognized action = \"" << action << "\".\n";
    exit ( 1 );
  }
  return;
}
//****************************************************************************80

void p05_lim ( int dim_num, double a[], double b[] )

//****************************************************************************80
//
//  Purpose:
//
//    P05_LIM returns the integration limits for problem 05.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the argument.
//
//    Output, double 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.
//
{
  int dim;

  for ( dim = 0; dim < dim_num; dim++ )
  {
    a[dim] = 0.0;
  }
  for ( dim = 0; dim < dim_num; dim++ )
  {
    b[dim] = 1.0;
  }

  return;
}
//****************************************************************************80

char *p05_name ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P05_NAME returns the name of problem 05.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *P05_NAME, the name of the problem.
//
{
  char *value;

  value = new char[5];

  strcpy ( value, "ST04" );

  return value;
}
//****************************************************************************80

char *p05_region ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P05_REGION returns the name of the integration region for problem 05.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *P05_REGION, the name of the integration region.
//
{
  char *value;

  value = new char[4];

  strcpy ( value, "BOX" );

  return value;
}
//****************************************************************************80

void p05_title ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P05_TITLE prints a title for problem 05.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
  cout << "\n";
  cout << "Problem 05\n";
  cout << "  Name:       ST04\n";
  cout << "              Stroud #4, page 26.\n";
  cout << "  Region:     0 <= X(i) <= 1\n";
  cout << "  Integrand:  F(X) = 1 / ( 1 + sum ( 2 * X(i) ) )\n";

  return;
}
//****************************************************************************80

void p06_default ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    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.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the problem.
//
{
  int i4;

  p06_i4 ( 'D', '*', &i4 );

  return;
}
//****************************************************************************80

double p06_exact ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    P06_EXACT returns the exact integral for problem 06.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Output, double P06_EXACT, the exact value of the integral.
//
{
  double exact;

  exact = 1.0;

  return exact;
}
//****************************************************************************80

double *p06_f ( int dim_num, int point_num, double x[] )

//****************************************************************************80
//
//  Purpose:
//
//    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
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  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, int DIM_NUM, the spatial dimension.
//
//    Input, int POINT_NUM, the number of points.
//
//    Input, double X[DIM_NUM*POINT_NUM], the evaluation points.
//
//    Output, double P06_F[POINT_NUM], the integrand values.
//
{
  int dim;
  int inc;
  int point;
  double *value;

  value = new double[point_num];

  for ( point = 0; point < point_num; point++ )
  {
    value[point] = 1.0;
    for ( dim = 0; dim < dim_num; dim++ )
    {
      value[point] = value[point] 
        * 2.0 * r8_abs ( 2.0 * x[dim+point*dim_num] - 1.0 );
    }
  }

  inc = point_num;
  p06_i4 ( 'I', '#', &inc );

  return value;
}
//****************************************************************************80

void p06_i4 ( char action, char name, int *value )

//****************************************************************************80
//
//  Purpose:
//
//    P06_I4 sets or gets I4 parameters for problem 06.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    09 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char ACTION,
//    'D' to set a parmater to its default value;
//    'S' to set a parameter to VALUE,
//    'G' to set VALUE to the parameter's value,
//    'I' to increment a parameter.
//
//    Input, char NAME, the name of the variable.
//    '#' is the number of calls to the integrand routine.
//
//    Input/output, int *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.
//
{
  static int calls = 0;

  if ( action == 'D' )
  {
    if ( name == '#' || name == '*' )
    {
      calls = 0;
    }
  }
  else if ( action == 'G' )
  {
    if ( name == '#' )
    {
      *value = calls;
    }
    else
    {
      cout << "\n";
      cout << "P06_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else if ( action == 'I' )
  {
    if ( name == '#' )
    {
      calls = calls + *value;
    }
    else
    {
      cout << "\n";
      cout << "P06_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else if ( action == 'S' )
  {
    if ( name == '#' )
    {
      calls = *value;
    }
    else
    {
      cout << "\n";
      cout << "P06_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else
  {
    cout << "\n";
    cout << "P06_I4 - Fatal error!\n";
    cout << "  Unrecognized action = \"" << action << "\".\n";
    exit ( 1 );
  }
  return;
}
//****************************************************************************80

void p06_lim ( int dim_num, double a[], double b[] )

//****************************************************************************80
//
//  Purpose:
//
//    P06_LIM returns the integration limits for problem 06.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the argument.
//
//    Output, double 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.
//
{
  int dim;

  for ( dim = 0; dim < dim_num; dim++ )
  {
    a[dim] = 0.0;
  }
  for ( dim = 0; dim < dim_num; dim++ )
  {
    b[dim] = 1.0;
  }

  return;
}
//****************************************************************************80

char *p06_name ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P06_NAME returns the name of problem 06.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *P06_NAME, the name of the problem.
//
{
  char *value;

  value = new char[7];

  strcpy ( value, "DR4061" );

  return value;
}
//****************************************************************************80

char *p06_region ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P06_REGION returns the name of the integration region for problem 06.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *P06_REGION, the name of the integration region.
//
{
  char *value;

  value = new char[4];

  strcpy ( value, "BOX" );

  return value;
}
//****************************************************************************80

void p06_title ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P06_TITLE prints a title for problem 06.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
  cout << "\n";
  cout << "Problem 07\n";
  cout << "  Name:       DR4061\n";
  cout << "              Davis, Rabinowitz, page 406, #1.\n";
  cout << "  Region:     0 <= X(i) <= 1\n";
  cout << "  Integrand:  F(X) = product ( abs ( 4 * X(i) - 2 ) )\n";

  return;
}
//****************************************************************************80

void p07_default ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    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.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the problem.
//
{
  int i4;

  p07_i4 ( 'D', '*', &i4 );

  return;
}
//****************************************************************************80

double p07_exact ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    P07_EXACT returns the exact integral for problem 07.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Output, double P07_EXACT, the exact value of the integral.
//
{
  double exact;

  exact = 1.0;

  return exact;
}
//****************************************************************************80

double *p07_f ( int dim_num, int point_num, double x[] )

//****************************************************************************80
//
//  Purpose:
//
//    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
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  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, int DIM_NUM, the spatial dimension.
//
//    Input, int POINT_NUM, the number of points.
//
//    Input, double X[DIM_NUM*POINT_NUM], the evaluation points.
//
//    Output, double P07_F[POINT_NUM], the integrand values.
//
{
  int dim;
  int inc;
  double pi = 3.141592653589793;
  int point;
  double *value;

  value = new double[point_num];

  for ( point = 0; point < point_num; point++ )
  {
    value[point] = 1.0;
    for ( dim = 0; dim < dim_num; dim++ )
    {
      value[point] = value[point] 
        * ( 0.5 * pi * sin ( pi * x[dim+point*dim_num] ) );
    }
  }

  inc = point_num;
  p07_i4 ( 'I', '#', &inc );

  return value;
}
//****************************************************************************80

void p07_i4 ( char action, char name, int *value )

//****************************************************************************80
//
//  Purpose:
//
//    P07_I4 sets or gets I4 parameters for problem 07.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    09 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char ACTION,
//    'D' to set a parmater to its default value;
//    'S' to set a parameter to VALUE,
//    'G' to set VALUE to the parameter's value,
//    'I' to increment a parameter.
//
//    Input, char NAME, the name of the variable.
//    '#' is the number of calls to the integrand routine.
//
//    Input/output, int *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.
//
{
  static int calls = 0;

  if ( action == 'D' )
  {
    if ( name == '#' || name == '*' )
    {
      calls = 0;
    }
  }
  else if ( action == 'G' )
  {
    if ( name == '#' )
    {
      *value = calls;
    }
    else
    {
      cout << "\n";
      cout << "P07_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else if ( action == 'I' )
  {
    if ( name == '#' )
    {
      calls = calls + *value;
    }
    else
    {
      cout << "\n";
      cout << "P07_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else if ( action == 'S' )
  {
    if ( name == '#' )
    {
      calls = *value;
    }
    else
    {
      cout << "\n";
      cout << "P07_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else
  {
    cout << "\n";
    cout << "P07_I4 - Fatal error!\n";
    cout << "  Unrecognized action = \"" << action << "\".\n";
    exit ( 1 );
  }
  return;
}
//****************************************************************************80

void p07_lim ( int dim_num, double a[], double b[] )

//****************************************************************************80
//
//  Purpose:
//
//    P07_LIM returns the integration limits for problem 07.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the argument.
//
//    Output, double 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.
//
{
  int dim;

  for ( dim = 0; dim < dim_num; dim++ )
  {
    a[dim] = 0.0;
  }
  for ( dim = 0; dim < dim_num; dim++ )
  {
    b[dim] = 1.0;
  }

  return;
}
//****************************************************************************80

char *p07_name ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P07_NAME returns the name of problem 07.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *P07_NAME, the name of the problem.
//
{
  char *value;

  value = new char[7];

  strcpy ( value, "DR4062" );

  return value;
}
//****************************************************************************80

char *p07_region ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P07_REGION returns the name of the integration region for problem 07.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *P07_REGION, the name of the integration region.
//
{
  char *value;

  value = new char[4];

  strcpy ( value, "BOX" );

  return value;
}
//****************************************************************************80

void p07_title ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P07_TITLE prints a title for problem 07.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
  cout << "\n";
  cout << "Problem 07\n";
  cout << "  Name:       DR4062\n";
  cout << "              Davis, Rabinowitz, page 406, #2.\n";
  cout << "  Region:     0 <= X(i) <= 1\n";
  cout << "  Integrand:  F(X) = product ( pi * sin ( pi * X(i) ) / 2 )\n";

  return;
}
//****************************************************************************80

void p08_default ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    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.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the problem.
//
{
  int i4;

  p08_i4 ( 'D', '*', &i4 );

  return;
}
//****************************************************************************80

double p08_exact ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    P08_EXACT returns the exact integral for problem 08.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Output, double P08_EXACT, the exact value of the integral.
//
{
  double exact;
  double pi = 3.141592653589793;

  exact = 0.5 - sqrt ( pow ( 2.0, 3 * dim_num ) ) 
    * cos ( 0.25 * pi * ( double ) ( dim_num ) ) 
    / ( 2.0 * pow ( pi, dim_num ) );

  return exact;
}
//****************************************************************************80

double *p08_f ( int dim_num, int point_num, double x[] )

//****************************************************************************80
//
//  Purpose:
//
//    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 )
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    03 June 2007
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Richard Crandall,
//    Projects in Scientific Computing,
//    Springer, 2000,
//    ISBN: 0387950095,
//    LC: Q183.9.C733.
//
//  Parameters:
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Input, int POINT_NUM, the number of points.
//
//    Input, double X[DIM_NUM*POINT_NUM], the evaluation points.
//
//    Output, double P08_F[POINT_NUM], the integrand values.
//
{
  int inc;
  double pi = 3.141592653589793;
  int point;
  double *value;

  value = new double[point_num];

  for ( point = 0; point < point_num; point++ )
  {
    value[point] = pow ( sin ( pi 
      * r8vec_sum ( dim_num, x+point*dim_num ) / 4.0 ), 2 );
  }

  inc = point_num;
  p08_i4 ( 'I', '#', &inc );

  return value;
}
//****************************************************************************80

void p08_i4 ( char action, char name, int *value )

//****************************************************************************80
//
//  Purpose:
//
//    P08_I4 sets or gets I4 parameters for problem 08.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    09 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char ACTION,
//    'D' to set a parmater to its default value;
//    'S' to set a parameter to VALUE,
//    'G' to set VALUE to the parameter's value,
//    'I' to increment a parameter.
//
//    Input, char NAME, the name of the variable.
//    '#' is the number of calls to the integrand routine.
//
//    Input/output, int *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.
//
{
  static int calls = 0;

  if ( action == 'D' )
  {
    if ( name == '#' || name == '*' )
    {
      calls = 0;
    }
  }
  else if ( action == 'G' )
  {
    if ( name == '#' )
    {
      *value = calls;
    }
    else
    {
      cout << "\n";
      cout << "P08_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else if ( action == 'I' )
  {
    if ( name == '#' )
    {
      calls = calls + *value;
    }
    else
    {
      cout << "\n";
      cout << "P08_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else if ( action == 'S' )
  {
    if ( name == '#' )
    {
      calls = *value;
    }
    else
    {
      cout << "\n";
      cout << "P08_I4 - Fatal error!\n";
      cout << "  Unrecognized name = \"" << name << "\".\n";
      exit ( 1 );
    }
  }
  else
  {
    cout << "\n";
    cout << "P08_I4 - Fatal error!\n";
    cout << "  Unrecognized action = \"" << action << "\".\n";
    exit ( 1 );
  }
  return;
}
//****************************************************************************80

void p08_lim ( int dim_num, double a[], double b[] )

//****************************************************************************80
//
//  Purpose:
//
//    P08_LIM returns the integration limits for problem 08.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the argument.
//
//    Output, double 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.
//
{
  int dim;

  for ( dim = 0; dim < dim_num; dim++ )
  {
    a[dim] = 0.0;
  }
  for ( dim = 0; dim < dim_num; dim++ )
  {
    b[dim] = 1.0;
  }

  return;
}
//****************************************************************************80

char *p08_name ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P08_NAME returns the name of problem 08.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *P08_NAME, the name of the problem.
//
{
  char *value;

  value = new char[5];

  strcpy ( value, "RC01" );

  return value;
}
//****************************************************************************80

char *p08_region ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P08_REGION returns the name of the integration region for problem 08.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *P08_REGION, the name of the integration region.
//
{
  char *value;

  value = new char[4];

  strcpy ( value, "BOX" );

  return value;
}
//****************************************************************************80

void p08_title ( void )

//****************************************************************************80
//
//  Purpose:
//
//    P08_TITLE prints a title for problem 08.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    07 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
  cout << "\n";
  cout << "Problem 08\n";
  cout << "  Name:       RC01\n";
  cout << "              Crandall, page 49, #1\n";
  cout << "  Region:     0 <= X(i) <= 1\n";
  cout << "  Integrand:  F(X) = sin^2 ( pi/4 * sum ( X(i) ) )\n";

  return;
}
//****************************************************************************80

void p09_default ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    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.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    11 March 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the dimension of the problem.
//
{
  int i4;

  p09_i4 ( 'D', '*', &i4 );
  p09_r8vec ( 'D', '*', dim_num, NULL );

  return;
}
//****************************************************************************80

double p09_exact ( int dim_num )

//****************************************************************************80
//
//  Purpose:
//
//    P09_EXACT returns the exact integral for problem 09.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    03 June 2007
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int DIM_NUM, the spatial dimension.
//
//    Output, double P09_EXACT, the exact value of the integral.
//
{
  double *c;
  int dim;
  double exact;

  c = new double[dim_num];