# include <stdlib.h>
# include <stdio.h>
# include <math.h>
# include <time.h>
# include <string.h>

# include "test_int.h"

/******************************************************************************/

char ch_cap ( char ch )

/******************************************************************************/
/*
  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;
}
/******************************************************************************/

double euler_constant ( )

/******************************************************************************/
/*
  Purpose:

    EULER_CONSTANT returns the value of the Euler-Mascheroni constant.

  Discussion:

    The Euler-Mascheroni constant is often denoted by a lower-case
    Gamma.  Gamma is defined as

      Gamma = limit ( M -> oo ) ( Sum ( 1 <= N <= M ) 1 / N ) - Log ( M )

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    21 December 2011

  Author:

    John Burkardt

  Parameters:

    Output, double EULER_CONSTANT, the value of the Euler-Mascheroni constant.
*/
{
  double value;

  value = 0.577215664901532860606512090082402431042;

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

int get_seed ( void )

/******************************************************************************/
/*
  Purpose:

    GET_SEED returns a random seed for the random number generator.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 November 2004

  Author:

    John Burkardt

  Parameters:

    Output, int GET_SEED, a random seed value.
*/
{
  time_t clock;
  int i4_huge = 2147483647;
  int ihour;
  int imin;
  int isec;
  int seed;
  struct tm *lt;
  time_t tloc;
/*
  If the internal seed is 0, generate a value based on the time.
*/
  clock = time ( &tloc );
  lt = localtime ( &clock );
/*
  Hours is 1, 2, ..., 12.
*/
  ihour = lt->tm_hour;

  if ( 12 < ihour )
  {
    ihour = ihour - 12;
  }
/*
  Move Hours to 0, 1, ..., 11
*/
  ihour = ihour - 1;

  imin = lt->tm_min;

  isec = lt->tm_sec;

  seed = isec + 60 * ( imin + 60 * ihour );
/*
  We want values in [1,43200], not [0,43199].
*/
  seed = seed + 1;
/*
  Remap SEED from [1,43200] to [1,HUGE].
*/
  seed = ( int ) 
    ( ( ( double ) seed )
    * ( ( double ) i4_huge ) / ( 60.0 * 60.0 * 12.0 ) );
/*
  Never use a seed of 0.
*/
  if ( seed == 0 )
  {
    seed = 1;
  }

  return seed;
}
/******************************************************************************/

int i4_huge ( void )

/******************************************************************************/
/*
  Purpose:

    I4_HUGE returns a "huge" I4.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    29 August 2006

  Author:

    John Burkardt

  Parameters:

    Output, int I4_HUGE, a "huge" integer.
*/
{
  static int value = 2147483647;

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

int i4_log_2 ( int i )

/******************************************************************************/
/*
  Purpose:

    I4_LOG_2 returns the integer part of the logarithm base 2 of an I4.

  Example:

        I  I4_LOG_10
    -----  --------
        0    0
        1    0
        2    1
        3    1
        4    2
        5    2
        7    2
        8    3
        9    3
     1000    9
     1024   10

  Discussion:

    I4_LOG_2 ( I ) + 1 is the number of binary digits in I.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    23 October 2007

  Author:

    John Burkardt

  Parameters:

    Input, int I, the number whose logarithm base 2 is desired.

    Output, int I4_LOG_2, the integer part of the logarithm base 2 of
    the absolute value of X.
*/
{
  int i_abs;
  int two_pow;
  int value;

  if ( i == 0 )
  {
    value = 0;
  }
  else
  {
    value = 0;
    two_pow = 2;

    i_abs = abs ( i );

    while ( two_pow <= i_abs )
    {
      value = value + 1;
      two_pow = two_pow * 2;
    }

  }

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

int i4_power ( int i, int j )

/******************************************************************************/
/*
  Purpose:

    I4_POWER returns the value of I^J.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    23 October 2007

  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 )
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "I4_POWER - Fatal error!\n" );
      fprintf ( stderr, "  I^J requested, with I = 0 and J negative.\n" );
      exit ( 1 );
    }
    else
    {
      value = 0;
    }
  }
  else if ( j == 0 )
  {
    if ( i == 0 )
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "I4_POWER - Fatal error!\n" );
      fprintf ( stderr, "  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;
}
/******************************************************************************/

double *i4_to_halton_number_sequence_new ( int seed, int base, int n )

/******************************************************************************/
/*
  Purpose:

    I4_TO_HALTON_NUMBER_SEQUENCE: next N elements of a scalar Halton sequence./

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    25 December 2011

  Author:

    John Burkardt

  Reference:

    John Halton,
    On the efficiency of certain quasi-random sequences of points
    in evaluating multi-dimensional integrals,
    Numerische Mathematik,
    Volume 2, pages 84-90, 1960.
 
  Parameters:

    Input, int SEED, the index of the desired element.
    Only the absolute value of SEED is considered.
    SEED = 0 is allowed, and returns R = 0.

    Input, int BASE, the Halton base, which should 
    be a prime number.  This routine only checks that BASE is greater 
    than 1.

    Input, int N, the number of elements desired.

    Output, double I4_TO_HALTON_NUMBER_SEQUENCE[N], the SEED-th through 
    (SEED+N-1)-th elements of the Halton sequence for base BASE.
*/
{
  double base_inv;
  int digit;
  int i;
  double *r;
  int *seed2;
/*
  Set SEED2 = ( SEED, SEED+1, SEED+2, ..., SEED+N-1 )
*/
  seed2 = i4vec_indicator_new ( n );

  for ( i = 0; i < n; i++ )
  {
    seed2[i] = seed2[i] + abs ( seed ) - 1;
  }

  if ( base <= 1 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "I4_TO_HALTON_NUMBER_SEQUENCE - Fatal error!\n" );
    fprintf ( stderr, "  The input base BASE is <= 1!\n" );
    fprintf ( stderr, "  BASE = %d\n", base );
    exit ( 1 );
  }

  r = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    r[i] = 0.0;
    base_inv = 1.0 / ( double ) ( base );
    while ( seed2[i] != 0 )
    {
      digit = ( seed2[i] % base );
      r[i] = r[i] + ( double ) ( digit ) * base_inv;
      base_inv = base_inv / ( double ) ( base );
      seed2[i] = seed2[i] / base;
    }
  }

  free ( seed2 );

  return r;
}
/******************************************************************************/

int *i4vec_indicator_new ( int n )

/******************************************************************************/
/*
  Purpose:

    I4VEC_INDICATOR_NEW sets an I4VEC to the indicator vector.

  Discussion:

    An I4VEC is a vector of I4's.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    26 August 2008

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of elements of A.

    Output, int I4VEC_INDICATOR_NEW[N], the array.
*/
{
  int *a;
  int i;

  a = ( int * ) malloc ( n * sizeof ( int ) );

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

double p00_even ( int prob, int int_num )

/******************************************************************************/
//
//  Purpose:
//
//    P00_EVEN uses evenly spaced points to integrate a function.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROB, the problem index.
//
//    Input, int INT_NUM, the number of sample points.
//
//    Output, double P00_EVEN, the approximate integral.
//
{
  double a;
  double b;
  double *fx;
  double result;
  double *x;

  p00_lim ( prob, &a, &b );

  x = r8vec_linspace_new ( int_num, a, b );

  fx = p00_fun ( prob, int_num, x );

  result = ( b - a ) * r8vec_sum ( int_num, fx ) / ( double ) ( int_num );

  free ( fx );
  free ( x );

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

double p00_exact ( int prob )

/******************************************************************************/
//
//  Purpose:
//
//    P00_EXACT returns the exact integral for any problem.
//
//  Discussion:
//
//    This routine provides a "generic" interface to the exact integral
//    routines for the various problems, and allows a problem to be called
//    by number (PROB) rather than by name.
//
//    In some cases, the "exact" value of the integral is in fact
//    merely a respectable approximation.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    29 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROB, the number of the desired test problem.
//
//    Output, double P00_EXACT, the exact value of the integral.
//
{
  double exact;

  if ( prob == 1 )
  {
    exact = p01_exact ( );
  }
  else if ( prob == 2 )
  {
    exact = p02_exact ( );
  }
  else if ( prob == 3 )
  {
    exact = p03_exact ( );
  }
  else if ( prob == 4 )
  {
    exact = p04_exact ( );
  }
  else if ( prob == 5 )
  {
    exact = p05_exact ( );
  }
  else if ( prob == 6 )
  {
    exact = p06_exact ( );
  }
  else if ( prob == 7 )
  {
    exact = p07_exact ( );
  }
  else if ( prob == 8 )
  {
    exact = p08_exact ( );
  }
  else if ( prob == 9 )
  {
    exact = p09_exact ( );
  }
  else if ( prob == 10 )
  {
    exact = p10_exact ( );
  }
  else if ( prob == 11 )
  {
    exact = p11_exact ( );
  }
  else if ( prob == 12 )
  {
    exact = p12_exact ( );
  }
  else if ( prob == 13 )
  {
    exact = p13_exact ( );
  }
  else if ( prob == 14 )
  {
    exact = p14_exact ( );
  }
  else if ( prob == 15 )
  {
    exact = p15_exact ( );
  }
  else if ( prob == 16 )
  {
    exact = p16_exact ( );
  }
  else if ( prob == 17 )
  {
    exact = p17_exact ( );
  }
  else if ( prob == 18 )
  {
    exact = p18_exact ( );
  }
  else if ( prob == 19 )
  {
    exact = p19_exact ( );
  }
  else if ( prob == 20 )
  {
    exact = p20_exact ( );
  }
  else if ( prob == 21 )
  {
    exact = p21_exact ( );
  }
  else if ( prob == 22 )
  {
    exact = p22_exact ( );
  }
  else if ( prob == 23 )
  {
    exact = p23_exact ( );
  }
  else if ( prob == 24 )
  {
    exact = p24_exact ( );
  }
  else if ( prob == 25 )
  {
    exact = p25_exact ( );
  }
  else if ( prob == 26 )
  {
    exact = p26_exact ( );
  }
  else if ( prob == 27 )
  {
    exact = p27_exact ( );
  }
  else if ( prob == 28 )
  {
    exact = p28_exact ( );
  }
  else if ( prob == 29 )
  {
    exact = p29_exact ( );
  }
  else if ( prob == 30 )
  {
    exact = p30_exact ( );
  }
  else if ( prob == 31 )
  {
    exact = p31_exact ( );
  }
  else if ( prob == 32 )
  {
    exact = p32_exact ( );
  }
  else if ( prob == 33 )
  {
    exact = p33_exact ( );
  }
  else if ( prob == 34 )
  {
    exact = p34_exact ( );
  }
  else if ( prob == 35 )
  {
    exact = p35_exact ( );
  }
  else if ( prob == 36 )
  {
    exact = p36_exact ( );
  }
  else if ( prob == 37 )
  {
    exact = p37_exact ( );
  }
  else if ( prob == 38 )
  {
    exact = p38_exact ( );
  }
  else if ( prob == 39 )
  {
    exact = p39_exact ( );
  }
  else if ( prob == 40 )
  {
    exact = p40_exact ( );
  }
  else if ( prob == 41 )
  {
    exact = p41_exact ( );
  }
  else if ( prob == 42 )
  {
    exact = p42_exact ( );
  }
  else if ( prob == 43 )
  {
    exact = p43_exact ( );
  }
  else if ( prob == 44 )
  {
    exact = p44_exact ( );
  }
  else if ( prob == 45 )
  {
    exact = p45_exact ( );
  }
  else if ( prob == 46 )
  {
    exact = p46_exact ( );
  }
  else if ( prob == 47 )
  {
    exact = p47_exact ( );
  }
  else if ( prob == 48 )
  {
    exact = p48_exact ( );
  }
  else if ( prob == 49 )
  {
    exact = p49_exact ( );
  }
  else if ( prob == 50 )
  {
    exact = p50_exact ( );
  }
  else if ( prob == 51 )
  {
    exact = p51_exact ( );
  }
  else if ( prob == 52 )
  {
    exact = p52_exact ( );
  }
  else if ( prob == 53 )
  {
    exact = p53_exact ( );
  }
  else if ( prob == 54 )
  {
    exact = p54_exact ( );
  }
  else if ( prob == 55 )
  {
    exact = p55_exact ( );
  }
  else if ( prob == 56 )
  {
    exact = p56_exact ( );
  }
  else if ( prob == 57 )
  {
    exact = p57_exact ( );
  }
  else
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P00_EXACT - Fatal error!\n" );
    fprintf ( stderr, "  Illegal problem number = %d\n", prob );
    exit ( 1 );
  }
  return exact;
}
/******************************************************************************/

double *p00_fun ( int prob, int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P00_FUN evaluates the integrand for any problem.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    29 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROB, the number of the desired test problem.
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double FX[N], the integrand values.
//
{
  double *fx;

  if ( prob == 1 )
  {
    fx = p01_fun ( n, x );
  }
  else if ( prob == 2 )
  {
    fx = p02_fun ( n, x );
  }
  else if ( prob == 3 )
  {
    fx = p03_fun ( n, x );
  }
  else if ( prob == 4 )
  {
    fx = p04_fun ( n, x );
  }
  else if ( prob == 5 )
  {
    fx = p05_fun ( n, x );
  }
  else if ( prob == 6 )
  {
    fx = p06_fun ( n, x );
  }
  else if ( prob == 7 )
  {
    fx = p07_fun ( n, x );
  }
  else if ( prob == 8 )
  {
    fx = p08_fun ( n, x );
  }
  else if ( prob == 9 )
  {
    fx = p09_fun ( n, x );
  }
  else if ( prob == 10 )
  {
    fx = p10_fun ( n, x );
  }
  else if ( prob == 11 )
  {
    fx = p11_fun ( n, x );
  }
  else if ( prob == 12 )
  {
    fx = p12_fun ( n, x );
  }
  else if ( prob == 13 )
  {
    fx = p13_fun ( n, x );
  }
  else if ( prob == 14 )
  {
    fx = p14_fun ( n, x );
  }
  else if ( prob == 15 )
  {
    fx = p15_fun ( n, x );
  }
  else if ( prob == 16 )
  {
    fx = p16_fun ( n, x );
  }
  else if ( prob == 17 )
  {
    fx = p17_fun ( n, x );
  }
  else if ( prob == 18 )
  {
    fx = p18_fun ( n, x );
  }
  else if ( prob == 19 )
  {
    fx = p19_fun ( n, x );
  }
  else if ( prob == 20 )
  {
    fx = p20_fun ( n, x );
  }
  else if ( prob == 21 )
  {
    fx = p21_fun ( n, x );
  }
  else if ( prob == 22 )
  {
    fx = p22_fun ( n, x );
  }
  else if ( prob == 23 )
  {
    fx = p23_fun ( n, x );
  }
  else if ( prob == 24 )
  {
    fx = p24_fun ( n, x );
  }
  else if ( prob == 25 )
  {
    fx = p25_fun ( n, x );
  }
  else if ( prob == 26 )
  {
    fx = p26_fun ( n, x );
  }
  else if ( prob == 27 )
  {
    fx = p27_fun ( n, x );
  }
  else if ( prob == 28 )
  {
    fx = p28_fun ( n, x );
  }
  else if ( prob == 29 )
  {
    fx = p29_fun ( n, x );
  }
  else if ( prob == 30 )
  {
    fx = p30_fun ( n, x );
  }
  else if ( prob == 31 )
  {
    fx = p31_fun ( n, x );
  }
  else if ( prob == 32 )
  {
    fx = p32_fun ( n, x );
  }
  else if ( prob == 33 )
  {
    fx = p33_fun ( n, x );
  }
  else if ( prob == 34 )
  {
    fx = p34_fun ( n, x );
  }
  else if ( prob == 35 )
  {
    fx = p35_fun ( n, x );
  }
  else if ( prob == 36 )
  {
    fx = p36_fun ( n, x );
  }
  else if ( prob == 37 )
  {
    fx = p37_fun ( n, x );
  }
  else if ( prob == 38 )
  {
    fx = p38_fun ( n, x );
  }
  else if ( prob == 39 )
  {
    fx = p39_fun ( n, x );
  }
  else if ( prob == 40 )
  {
    fx = p40_fun ( n, x );
  }
  else if ( prob == 41 )
  {
    fx = p41_fun ( n, x );
  }
  else if ( prob == 42 )
  {
    fx = p42_fun ( n, x );
  }
  else if ( prob == 43 )
  {
    fx = p43_fun ( n, x );
  }
  else if ( prob == 44 )
  {
    fx = p44_fun ( n, x );
  }
  else if ( prob == 45 )
  {
    fx = p45_fun ( n, x );
  }
  else if ( prob == 46 )
  {
    fx = p46_fun ( n, x );
  }
  else if ( prob == 47 )
  {
    fx = p47_fun ( n, x );
  }
  else if ( prob == 48 )
  {
    fx = p48_fun ( n, x );
  }
  else if ( prob == 49 )
  {
    fx = p49_fun ( n, x );
  }
  else if ( prob == 50 )
  {
    fx = p50_fun ( n, x );
  }
  else if ( prob == 51 )
  {
    fx = p51_fun ( n, x );
  }
  else if ( prob == 52 )
  {
    fx = p52_fun ( n, x );
  }
  else if ( prob == 53 )
  {
    fx = p53_fun ( n, x );
  }
  else if ( prob == 54 )
  {
    fx = p54_fun ( n, x );
  }
  else if ( prob == 55 )
  {
    fx = p55_fun ( n, x );
  }
  else if ( prob == 56 )
  {
    fx = p56_fun ( n, x );
  }
  else if ( prob == 57 )
  {
    fx = p57_fun ( n, x );
  }
  else
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P00_FUN - Fatal error!\n" );
    fprintf ( stderr, "  Illegal problem number = %d\n", prob );
    exit ( 1 );
  }
  return fx;
}
/******************************************************************************/

double p00_gauss_legendre ( int prob, int int_num )

/******************************************************************************/
//
//  Purpose:
//
//    P00_GAUSS_LEGENDRE applies a composite Gauss-Legendre rule.
//
//  Discussion:
//
//    A 4 point rule is used.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROB, the problem index.
//
//    Input, int INT_NUM, the number of subintervals.
//
//    Output, double P00_GAUSS_LEGENDRE, the approximate integral.
//
{
# define GAUSS_NUM 4

  double a;
  double a_sub;
  double b;
  double b_sub;
  double *fx;
  double gauss_abs[GAUSS_NUM] = {
    -0.861136311594052575223946488893, 
    -0.339981043584856264802665759103, 
     0.339981043584856264802665759103, 
     0.861136311594052575223946488893 };
  double gauss_weight[GAUSS_NUM] = {
    0.347854845137453857373063949222, 
    0.652145154862546142626936050778, 
    0.652145154862546142626936050778, 
    0.347854845137453857373063949222 };
  double h;
  int i;
  int int_i;
  double result;
  double x[GAUSS_NUM];

  p00_lim ( prob, &a, &b );

  h = ( b - a ) / ( double ) ( int_num );

  result = 0.0;

  for ( int_i = 1; int_i <= int_num; int_i++ )
  {
    a_sub = ( ( double ) ( int_num - int_i + 1 ) * a   
            + ( double ) (           int_i - 1 ) * b ) 
            / ( double ) ( int_num             );

    b_sub = ( ( double ) ( int_num - int_i ) * a   
            + ( double ) (           int_i ) * b ) 
            / ( double ) ( int_num         );

    for ( i = 0; i < GAUSS_NUM; i++ )
    {
      x[i] = 0.5 * ( ( 1.0 - gauss_abs[i] ) * a_sub 
                   + ( 1.0 + gauss_abs[i] ) * b_sub );
    }
    fx = p00_fun ( prob, GAUSS_NUM, x );

    result = result 
      + 0.5 * h * r8vec_dot_product ( GAUSS_NUM, gauss_weight, fx );

    free ( fx );
  }
  return result;
# undef GAUSS_NUM
}
/******************************************************************************/

double p00_halton ( int prob, int int_num )

/******************************************************************************/
//
//  Purpose:
//
//    P00_HALTON applies a Halton sequence rule to integrate a function.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    John Halton,
//    On the efficiency of certain quasi-random sequences of points
//    in evaluating multi-dimensional integrals,
//    Numerische Mathematik,
//    Volume 2, pages 84-90, 1960.
//
//  Parameters:
//
//    Input, int PROB, the problem index.
//
//    Input, int INT_NUM, the number of sample points.
//
//    Output, double P00_HALTON, the approximate integral.
//
{
  double a;
  double b;
  int base;
  double *fx;
  int i;
  double result;
  int seed;
  double *x;

  p00_lim ( prob, &a, &b );

  seed = 1;
  base = 2;
  x = i4_to_halton_number_sequence_new ( seed, base, int_num );

  for ( i = 0; i < int_num; i++ )
  {
    x[i] = a + ( b - a ) * x[i];
  }
  fx = p00_fun ( prob, int_num, x );

  result = ( b - a ) * r8vec_sum ( int_num, fx ) / ( double ) ( int_num );

  free ( fx );
  free ( x );

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

void p00_lim ( int prob, double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P00_LIM returns the integration limits for any problem.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    29 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROB, the number of the desired test problem.
//
//    Output, double *A, *B, the limits of integration.
//
{
  if ( prob == 1 )
  {
    p01_lim ( a, b );
  }
  else if ( prob == 2 )
  {
    p02_lim ( a, b );
  }
  else if ( prob == 3 )
  {
    p03_lim ( a, b );
  }
  else if ( prob == 4 )
  {
    p04_lim ( a, b );
  }
  else if ( prob == 5 )
  {
    p05_lim ( a, b );
  }
  else if ( prob == 6 )
  {
    p06_lim ( a, b );
  }
  else if ( prob == 7 )
  {
    p07_lim ( a, b );
  }
  else if ( prob == 8 )
  {
    p08_lim ( a, b );
  }
  else if ( prob == 9 )
  {
    p09_lim ( a, b );
  }
  else if ( prob == 10 )
  {
    p10_lim ( a, b );
  }
  else if ( prob == 11 )
  {
    p11_lim ( a, b );
  }
  else if ( prob == 12 )
  {
    p12_lim ( a, b );
  }
  else if ( prob == 13 )
  {
    p13_lim ( a, b );
  }
  else if ( prob == 14 )
  {
    p14_lim ( a, b );
  }
  else if ( prob == 15 )
  {
    p15_lim ( a, b );
  }
  else if ( prob == 16 )
  {
    p16_lim ( a, b );
  }
  else if ( prob == 17 )
  {
    p17_lim ( a, b );
  }
  else if ( prob == 18 )
  {
    p18_lim ( a, b );
  }
  else if ( prob == 19 )
  {
    p19_lim ( a, b );
  }
  else if ( prob == 20 )
  {
    p20_lim ( a, b );
  }
  else if ( prob == 21 )
  {
    p21_lim ( a, b );
  }
  else if ( prob == 22 )
  {
    p22_lim ( a, b );
  }
  else if ( prob == 23 )
  {
    p23_lim ( a, b );
  }
  else if ( prob == 24 )
  {
    p24_lim ( a, b );
  }
  else if ( prob == 25 )
  {
    p25_lim ( a, b );
  }
  else if ( prob == 26 )
  {
    p26_lim ( a, b );
  }
  else if ( prob == 27 )
  {
    p27_lim ( a, b );
  }
  else if ( prob == 28 )
  {
    p28_lim ( a, b );
  }
  else if ( prob == 29 )
  {
    p29_lim ( a, b );
  }
  else if ( prob == 30 )
  {
    p30_lim ( a, b );
  }
  else if ( prob == 31 )
  {
    p31_lim ( a, b );
  }
  else if ( prob == 32 )
  {
    p32_lim ( a, b );
  }
  else if ( prob == 33 )
  {
    p33_lim ( a, b );
  }
  else if ( prob == 34 )
  {
    p34_lim ( a, b );
  }
  else if ( prob == 35 )
  {
    p35_lim ( a, b );
  }
  else if ( prob == 36 )
  {
    p36_lim ( a, b );
  }
  else if ( prob == 37 )
  {
    p37_lim ( a, b );
  }
  else if ( prob == 38 )
  {
    p38_lim ( a, b );
  }
  else if ( prob == 39 )
  {
    p39_lim ( a, b );
  }
  else if ( prob == 40 )
  {
    p40_lim ( a, b );
  }
  else if ( prob == 41 )
  {
    p41_lim ( a, b );
  }
  else if ( prob == 42 )
  {
    p42_lim ( a, b );
  }
  else if ( prob == 43 )
  {
    p43_lim ( a, b );
  }
  else if ( prob == 44 )
  {
    p44_lim ( a, b );
  }
  else if ( prob == 45 )
  {
    p45_lim ( a, b );
  }
  else if ( prob == 46 )
  {
    p46_lim ( a, b );
  }
  else if ( prob == 47 )
  {
    p47_lim ( a, b );
  }
  else if ( prob == 48 )
  {
    p48_lim ( a, b );
  }
  else if ( prob == 49 )
  {
    p49_lim ( a, b );
  }
  else if ( prob == 50 )
  {
    p50_lim ( a, b );
  }
  else if ( prob == 51 )
  {
    p51_lim ( a, b );
  }
  else if ( prob == 52 )
  {
    p52_lim ( a, b );
  }
  else if ( prob == 53 )
  {
    p53_lim ( a, b );
  }
  else if ( prob == 54 )
  {
    p54_lim ( a, b );
  }
  else if ( prob == 55 )
  {
    p55_lim ( a, b );
  }
  else if ( prob == 56 )
  {
    p56_lim ( a, b );
  }
  else if ( prob == 57 )
  {
    p57_lim ( a, b );
  }
  else
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P00_LIM - Fatal error!\n" );
    fprintf ( stderr, "  Illegal problem number = %d\n", prob );
    exit ( 1 );
  }
  return;
}
/******************************************************************************/

double p00_midpoint ( int prob, int int_num )

/******************************************************************************/
//
//  Purpose:
//
//    P00_MIDPOINT applies the composite midpoint rule to integrate a function.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROB, the problem index.
//
//    Input, int INT_NUM, the number of subintervals.
//
//    Output, double P00_MIDPOINT, the approximate integral.
//
{
  double a;
  double a_sub;
  double b;
  double b_sub;
  double *fx;
  double h;
  int int_i;
  double result;
  double *x;

  p00_lim ( prob, &a, &b );

  h = b - a;

  x = ( double * ) malloc ( int_num * sizeof ( double ) );

  for ( int_i = 0; int_i < int_num; int_i++ )
  {
    a_sub = ( ( double ) ( int_num - int_i     ) * a   
            + ( double ) (           int_i     ) * b ) 
            / ( double ) ( int_num             );

    b_sub = ( ( double ) ( int_num - int_i - 1 ) * a   
            + ( double ) (           int_i + 1 ) * b ) 
            / ( double ) ( int_num         );

    x[int_i] = 0.5 * ( a_sub + b_sub );
  }

  fx = p00_fun ( prob, int_num, x );

  result = h * r8vec_sum ( int_num, fx ) / ( double ) ( int_num );

  free ( fx );
  free ( x );

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

double p00_montecarlo ( int prob, int int_num )

/******************************************************************************/
//
//  Purpose:
//
//    P00_MONTECARLO applies the Monte Carlo rule to integrate a function.
//
//  Discussion:
//
//    This routine originally used an automatic array for X.  However,
//    under the G95 compiler, this was causing bizarre errors.  Replacing
//    the automatic array by an allocatable array made the problems
//    disappear.  Not an entirely satisfactory conclusion//
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROB, the problem index.
//
//    Input, int INT_NUM, the number of sample points.
//
//    Output, double P00_MONTECARLO, the approximate integral.
//
{
  double a;
  double b;
  double *fx;
  int i;
  double result;
  int seed;
  double *x;

  p00_lim ( prob, &a, &b );

  seed = get_seed ( );

  x = r8vec_uniform_01_new ( int_num, &seed );

  for ( i = 0; i < int_num; i++ )
  {
    x[i] = a + ( b - a ) * x[i];
  }
  fx = p00_fun ( prob, int_num, x );

  result = ( b - a ) * r8vec_sum ( int_num, fx ) / ( double ) ( int_num );

  free ( fx );
  free ( x );

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

int p00_prob_num ( )

/******************************************************************************/
//
//  Purpose:
//
//    P00_PROB_NUM returns the number of test integration problems.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    29 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, int P00_PROB_NUM, the number of test integration
//    problems.
//
{
  int value;

  value = 57;

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

double p00_simpson ( int prob, int int_num )

/******************************************************************************/
//
//  Purpose:
//
//    P00_SIMPSON applies the composite Simpson rule to integrate a function.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROB, the problem index.
//
//    Input, int INT_NUM, the number of subintervals.
//
//    Output, double P00_SIMPSON, the approximate integral.
//
{
  double a;
  double b;
  double *fx;
  double h;
  int i;
  double result;
  double sum2;
  double sum4;
  double *x;

  p00_lim ( prob, &a, &b );

  h = b - a;

  x = r8vec_linspace_new ( 2 * int_num + 1, a, b );
 
  fx = p00_fun ( prob, 2 * int_num + 1, x );

  sum4 = 0.0;
  for ( i = 1; i < 2 * int_num + 1; i = i + 2 )
  {
    sum4 = sum4 + 4.0 * fx[i];
  }
  sum2 = 0.0;
  for ( i = 2; i < 2 * int_num; i = i + 2 )
  {
    sum2 = sum2 + 2.0 * fx[i];
  }
  result = fx[0] + fx[2*int_num] + sum4 + sum2;

  result = h * result / 6.0 / ( double ) int_num;

  free ( fx );
  free ( x );

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

double p00_trapezoid ( int prob, int int_num )

/******************************************************************************/
//
//  Purpose:
//
//    P00_TRAPEZOID applies the composite trapezoid rule to integrate a function.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int PROB, the problem index.
//
//    Input, int INT_NUM, the number of subintervals.
//
//    Output, double P00_TRAPEZOID, the approximate integral.
//
{
  double a;
  double b;
  double *fx;
  double result;
  double *x;

  p00_lim ( prob, &a, &b );

  x = r8vec_linspace_new ( int_num + 1, a, b );

  fx = p00_fun ( prob, int_num + 1, x );

  result = ( b - a ) * ( 
           0.5 * fx[0] 
         + r8vec_sum ( int_num - 1, fx + 1 )
         + 0.5 * fx[int_num] ) / ( double ) ( int_num );

  free ( fx );
  free ( x );

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

double p01_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P01_EXACT returns the exact integral for problem 1.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P01_EXACT, the value of the integral.
//
{
  double exact;

  exact = exp ( 1.0 ) - 1.0;

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

double *p01_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P01_FUN evaluates the integrand for problem 1.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    exp ( x )
//
//  Antiderivative:
//
//    exp ( x )
//
//  Exact Integral:
//
//    exp ( 1 ) - 1
//
//  Approximate Integral (25 digits):
//
//    1.718281828459045235360287...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software, edited by John R Rice,
//    Academic Press, 1971.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double FX[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = exp ( x[i] );
  }
  return fx;
}
/******************************************************************************/

void p01_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P01_LIM returns the integration limits for problem 1.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p02_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P02_EXACT returns the exact integral for problem 2.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P02_EXACT, the value of the integral.
//
{
  double exact;

  exact = 0.7;

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

double *p02_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P02_FUN evaluates the integrand for problem 2.
//
//  Discussion:
//
//    The integrand is discontinuous at X = 0.3.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    if ( x < 0.3 )
//      f(x) = 0
//    else
//      f(x) = 1
//
//  Antiderivative:
//
//    if ( x < 0.3 )
//      g(x) = 0
//    else
//      g(x) = x - 0.3
//
//  Exact Integral:
//
//    0.7
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software, edited by John R Rice,
//    Academic Press, 1971.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P02_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] < 0.3 )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = 1.0;
    }
  }
  return fx;
}
/******************************************************************************/

void p02_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P02_LIM returns the integration limits for problem 2.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p03_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P03_EXACT returns the exact integral for problem 3.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P03_EXACT, the value of the integral.
//
{
  double exact;

  exact = 2.0 / 3.0;

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

double *p03_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P03_FUN evaluates the integrand for problem 3.
//
//  Discussion:
//
//    The integrand is not differentiable at X = 0.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    sqrt ( x )
//
//  Antiderivative:
//
//    ( 2 / 3 ) * x^(3/2)
//
//  Exact Integral:
//
//    2 / 3
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software, edited by John R Rice,
//    Academic Press, 1971.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P03_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = sqrt ( x[i] );
  }
  return fx;
}
/******************************************************************************/

void p03_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P03_LIM returns the integration limits for problem 3.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p04_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P04_EXACT returns the estimated integral for problem 4.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P04_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = 0.47942822668880166736;

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

double *p04_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P04_FUN evaluates the integrand for problem 4.
//
//  Interval:
//
//    -1 <= x <= 1
//
//  Integrand:
//
//    0.92 * cosh ( x ) - cos ( x )
//
//  Antiderivative:
//
//    0.92 * sinh ( x ) - sin ( x )
//
//  Exact Integral:
//
//    1.84 * sinh ( 1 ) - 2 * sin ( 1 )
//
//  Approximate Integral (20 digits):
//
//    0.47942822668880166736...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Charles Clenshaw, Alan Curtis,
//    A Method for Numerical Integration on an Automatic Computer,
//    Numerische Mathematik,
//    Volume 2, Number 1, December 1960, pages 197-205.
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software, edited by John R Rice,
//    Academic Press, 1971.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P04_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = 0.92 * cosh ( x[i] ) - cos ( x[i] );
  }
  return fx;
}
/******************************************************************************/

void p04_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P04_LIM returns the integration limits for problem 4.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = -1.0;
  *b = 1.0;

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

double p05_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P05_EXACT returns the estimated integral for problem 5.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P05_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = 1.5822329637296729331;

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

double *p05_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P05_FUN evaluates the integrand for problem 5.
//
//  Interval:
//
//    -1 <= x <= 1
//
//  Integrand:
//
//    1 / ( x^4 + x^2 + 0.9 )
//
//  Approximate Integral (20 digits):
//
//    1.5822329637296729331...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Charles Clenshaw, Alan Curtis,
//    A Method for Numerical Integration on an Automatic Computer,
//    Numerische Mathematik,
//    Volume 2, Number 1, December 1960, pages 197-205.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P05_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = 1.0 / ( pow ( x[i], 4 ) + pow ( x[i], 2 ) + 0.9 );
  }
  return fx;
}
/******************************************************************************/

void p05_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P05_LIM returns the integration limits for problem 5.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = -1.0;
  *b = 1.0;

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

double p06_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P06_EXACT returns the exact integral for problem 6.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P06_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = 1.460447131787105;

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

double *p06_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P06_FUN evaluates the integrand for problem 6.
//
//  Interval:
//
//    -1 <= x <= 1
//
//  Integrand:
//
//    sqrt ( abs ( x + 0.5 ) )
//
//  Exact Integral:
//
//    ( sqrt ( 2 ) + 3 * sqrt ( 6 ) ) / 6 = 1.460447131787105
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Charles Clenshaw, Alan Curtis,
//    A Method for Numerical Integration on an Automatic Computer,
//    Numerische Mathematik,
//    Volume 2, Number 1, December 1960, pages 197-205.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P06_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = sqrt ( abs ( x[i] + 0.5 ) );
  }
  return fx;
}
/******************************************************************************/

void p06_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P06_LIM returns the integration limits for problem 6.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = -1.0;
  *b = 1.0;

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

double p07_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P07_EXACT returns the exact integral for problem 7.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P07_EXACT, the value of the integral.
//
{
  double exact;

  exact = 2.0;

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

double *p07_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P07_FUN evaluates the integrand for problem 7.
//
//  Discussion:
//
//    The integrand is singular at x = 0.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    1 / sqrt ( x )
//
//  Antiderivative:
//
//    2 * sqrt ( x )
//
//  Exact Integral:
//
//    2
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software, edited by John R Rice,
//    Academic Press, 1971
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P07_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( 0.0 < x[i] )
    {
      fx[i] = 1.0 / sqrt ( x[i] );
    }
    else
    {
      fx[i] = 0.0;
    }
  }
  return fx;
}
/******************************************************************************/

void p07_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P07_LIM returns the integration limits for problem 7.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p08_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P08_EXACT returns the estimated integral for problem 8.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P08_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = 0.86697298733991103757;

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

double *p08_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P08_FUN evaluates the integrand for problem 8.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    1 / ( 1 + x^4 )
//
//  Antiderivative:
//
//    (1/8) * sqrt ( 2 )
//    * ln ( ( x^2 + sqrt ( 2 ) * x + 1 ) / ( x^2 - sqrt ( 2 ) * x + 1 ) )
//    + (1/4) * sqrt ( 2 ) * arctan ( sqrt ( 2 ) * x / ( 1 - x^2 ) )
//
//  Approximate Integral (20 digits):
//
//    0.86697298733991103757...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software, edited by John R Rice,
//    Academic Press, 1971.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P08_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = 1.0 / ( 1.0 + pow ( x[i], 4 ) );
  }
  return fx;
}
/******************************************************************************/

void p08_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P08_LIM returns the integration limits for problem 8.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p09_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P09_EXACT returns the estimated integral for problem 9.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P09_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = 1.1547005383792515290;

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

double *p09_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P09_FUN evaluates the integrand for problem 9.
//
//  Discussion:
//
//    The integrand is oscillatory, going through 5 periods in [0,1].
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    2 / ( 2 + sin ( 10 * pi * x ) )
//
//  Antiderivative:
//
//    1 / ( 5 * pi * sqrt ( 3 ) ) *
//    arctan ( ( 1 + 2 * tan ( 5 * pi * x ) ) / sqrt ( 3 ) )
//
//  Exact Integral:
//
//    2 / sqrt ( 3 )
//
//  Approximate Integral (20 digits):
//
//    1.1547005383792515290...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software, edited by John R Rice,
//    Academic Press, 1971.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P09_FUN[N], the integrand values.
//
{
  double *fx;
  int i;
  const double r8_pi = 3.141592653589793;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = 2.0 / ( 2.0 + sin ( 10.0 * r8_pi * x[i] ) );
  }
  return fx;
}
/******************************************************************************/

void p09_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P09_LIM returns the integration limits for problem 9.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p10_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P10_EXACT returns the estimated integral for problem 10.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P10_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = 0.6931471805599453094172321;

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

double *p10_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P10_FUN evaluates the integrand for problem 10.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    1 / ( 1 + x )
//
//  Antiderivative:
//
//    ln ( 1 + x )
//
//  Exact Integral:
//
//    ln ( 2 )
//
//  Approximate Integral (25 digits):
//
//    0.6931471805599453094172321...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software, edited by John R Rice,
//    Academic Press, 1971.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P10_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = 1.0 / ( 1.0 + x[i] );
  }
  return fx;
}
/******************************************************************************/

void p10_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P10_LIM returns the integration limits for problem 10.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p11_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P11_EXACT returns the estimated integral for problem 11.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P11_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = 0.37988549304172247537;

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

double *p11_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P11_FUN evaluates the integrand for problem 11.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    1 / ( 1 + exp ( x ) )
//
//  Antiderivative:
//
//    ln ( exp ( x ) / ( 1 + exp ( x ) ) )
//
//  Exact Integral:
//
//    ln ( 2 * e / ( 1 + e ) )
//
//  Approximate Integral (20 digits):
//
//    0.37988549304172247537...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software, edited by John R Rice,
//    Academic Press, 1971.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P11_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = 1.0 / ( 1.0 + exp ( x[i] ) );
  }
  return fx;
}
/******************************************************************************/

void p11_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P11_LIM returns the integration limits for problem 11.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p12_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P12_EXACT returns the estimated integral for problem 12.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P12_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = 0.77750463411224827642;

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

double *p12_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P12_FUN evaluates the integrand for problem 12.
//
//  Discussion:
//
//    The integrand has a removable singularity at x = 0.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    x / ( exp ( x ) - 1 )
//
//  Antiderivative:
//
//    The Debye function.
//
//  Approximate Integral (20 digits):
//
//    0.77750463411224827642...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software, edited by John R Rice,
//    Academic Press, 1971.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P12_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] == 0.0 )
    {
      fx[i] = 1.0;
    }
    else
    {
      fx[i] = x[i] / ( exp ( x[i] ) - 1.0 );
    }
  }
  return fx;
}
/******************************************************************************/

void p12_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P12_LIM returns the integration limits for problem 12.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p13_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P13_EXACT returns the estimated integral for problem 13.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P13_EXACT, the estimated value of the integral.
//
{
  double a;
  double b;
  double exact;
 
  p13_lim ( &a, &b );

  exact = r8_si ( b ) - r8_si ( a );

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

double *p13_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P13_FUN evaluates the integrand for problem 13.
//
//  Interval:
//
//    0 <= x <= 10
//
//  Integrand:
//
//    sin ( x ) / x
//
//  Approximate Integral (20 digits):
//
//    1.6583475942188740493...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P13_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] == 0.0 )
    {
      fx[i] = 1.0;
    }
    else
    {
      fx[i] = sin ( x[i] ) / x[i];
    }
  }
  return fx;
}
/******************************************************************************/

void p13_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P13_LIM returns the integration limits for problem 13.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 10.0;

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

double p14_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P14_EXACT returns the estimated integral for problem 14.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P14_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = 0.500000211166;

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

double *p14_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P14_FUN evaluates the integrand for problem 14.
//
//  Discussion:
//
//    For X's that aren't actually very big, the function becomes very 
//    small.  Some compilers may product code that fails in these cases.
//    An attempt has been made to return a value of 0 when the computed
//    value of F(X) would be extremely small.
//
//  Interval:
//
//    0 <= x <= 10
//
//  Integrand:
//
//    sqrt ( 50 ) * exp ( - 50 * pi * x * x )
//
//  Exact Integral:
//
//    0.5 * erf ( 50 * sqrt ( 2 * pi ) )
//
//  Approximate Integral:
//
//    0.500000211166...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software, edited by John R Rice,
//    Academic Press, 1971.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P14_FUN[N], the integrand values.
//
{
  double *fx;
  int i;
  const double r8_pi = 3.141592653589793;
  static double x_max = 0.0;

  if ( x_max == 0.0 )
  {
    x_max = sqrt ( log ( r8_max ( r8_epsilon ( ), 1.0E-10 ) ) 
      / ( - 50.0 * r8_pi ) );
  }

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {

    if ( x_max < r8_abs ( x[i] ) )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = sqrt ( 50.0 ) * exp ( - 50.0 * r8_pi * x[i] * x[i] );
    }
  }
  return fx;
}
/******************************************************************************/

void p14_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P14_LIM returns the integration limits for problem 14.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 10.0;

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

double p15_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P15_EXACT returns the exact integral for problem 15.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P15_EXACT, the value of the integral.
//
{
  double exact;

  exact = 1.0;

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

double *p15_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P15_FUN evaluates the integrand for problem 15.
//
//  Interval:
//
//    0 <= x <= 10
//
//  Integrand:
//
//    25 * exp ( - 25 * x )
//
//  Antiderivative:
//
//    - exp ( - 25 * x )
//
//  Exact Integral:
//
//    1 - exp ( - 250 )
//
//  Approximate Integral:
//
//    1.00000000...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software, edited by John R Rice,
//    Academic Press, 1971.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P15_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = 25.0 * exp ( - 25.0 * x[i] );
  }
  return fx;
}
/******************************************************************************/

void p15_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P15_LIM returns the integration limits for problem 15.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 10.0;

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

double p16_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P16_EXACT returns the exact integral for problem 16.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P16_EXACT, the value of the integral.
//
{
  double exact;

  exact = 0.49936338107645674464;

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

double *p16_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P16_FUN evaluates the integrand for problem 16.
//
//  Interval:
//
//    0 <= x <= 10
//
//  Integrand:
//
//    50.0 / ( pi * ( 2500.0 * x * x + 1.0 ) )
//
//  Antiderivative:
//
//    ( 1 / pi ) * arctan ( 50 * x )
//
//  Approximate Integral (20 digits):
//
//    0.49936338107645674464...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software, edited by John R Rice,
//    Academic Press, 1971.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P16_FUN[N], the integrand values.
//
{
  double *fx;
  int i;
  const double r8_pi = 3.141592653589793;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = 50.0 / r8_pi / ( 2500.0 * x[i] * x[i] + 1.0 );
  }
  return fx;
}
/******************************************************************************/

void p16_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P16_LIM returns the integration limits for problem 16.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p17_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P17_EXACT returns the estimated integral for problem 17.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P17_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = 0.5;

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

double *p17_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P17_FUN evaluates the integrand for problem 17.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    ( sin ( 50 * pi * x ) )^2
//
//  Antiderivative:
//
//    1/2 x - sin ( 100 * pi * x ) / ( 200 * pi )
//
//  Approximate Integral:
//
//    0.5
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P17_FUN[N], the integrand values.
//
{
  double *fx;
  int i;
  const double r8_pi = 3.141592653589793;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = pow ( sin ( 50.0 * r8_pi * x[i] ), 2 );
  }
  return fx;
}
/******************************************************************************/

void p17_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P17_LIM returns the integration limits for problem 17.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p18_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P18_EXACT returns the estimated integral for problem 18.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P18_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = 0.17055734950243820437;

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

double *p18_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P18_FUN evaluates the integrand for problem 18.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    x / ( exp ( x ) + 1 )
//
//  Approximate Integral (20 digits):
//
//    0.17055734950243820437...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Hermann Engels,
//    Numerical Quadrature and Cubature,
//    Academic Press, 1980.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P18_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = x[i] / ( exp ( x[i] ) + 1.0 );
  }
  return fx;
}
/******************************************************************************/

void p18_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P18_LIM returns the integration limits for problem 18.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p19_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P19_EXACT returns the exact integral for problem 19.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P19_EXACT, the value of the integral.
//
{
  double exact;

  exact = - 1.0;

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

double *p19_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P19_FUN evaluates the integrand for problem 19.
//
//  Discussion:
//
//    The integrand is singular at x = 0.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    ln ( x )
//
//  Antiderivative:
//
//    x * ln ( x ) - x
//
//  Exact Integral:
//
//    -1
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software, edited by John R Rice,
//    Academic Press, 1971.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P19_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] <= 1.0E-15 )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = log ( x[i] );
    }
  }
  return fx;
}
/******************************************************************************/

void p19_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P19_LIM returns the integration limits for problem 19.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p20_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P20_EXACT returns the estimated integral for problem 20.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P20_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = 1.5643964440690497731;

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

double *p20_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P20_FUN evaluates the integrand for problem 20.
//
//  Interval:
//
//    -1 <= x <= 1
//
//  Integrand:
//
//    1 / ( x^2 + 1.005 )
//
//  Antiderivative:
//
//    ( 1 / sqrt ( 1.005 ) ) * arctan ( x / sqrt ( 1.005 ) )
//
//  Approximate Integral (20 digits):
//
//    1.5643964440690497731...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software, edited by John R Rice,
//    Academic Press, 1971.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P20_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = 1.0 / ( pow ( x[i], 2 ) + 1.005 );
  }
  return fx;
}
/******************************************************************************/

void p20_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P20_LIM returns the integration limits for problem 20.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = -1.0;
  *b = 1.0;

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

double p21_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P21_EXACT returns the estimated integral for problem 21.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P21_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = 0.21080273631018169851;

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

double *p21_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P21_FUN evaluates the integrand for problem 21.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//       ( sech (   10.0 * ( x - 0.2 ) ) )^2
//     + ( sech (  100.0 * ( x - 0.4 ) ) )^4
//     + ( sech ( 1000.0 * ( x - 0.6 ) ) )^6
//
//  Exact Integral:
//
//    ( 1 + tanh ( 8 ) * tanh ( 2 ) ) / 10.0 + 2 / 150 + 2 / 1875
//
//  Approximate Integral (20 digits):
//
//    0.21080273631018169851...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software, edited by John R Rice,
//    Academic Press, 1971.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P21_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = 
      pow ( r8_sech (   10.0 * ( x[i] - 0.2 ) ), 2 )
    + pow ( r8_sech (  100.0 * ( x[i] - 0.4 ) ), 4 )
    + pow ( r8_sech ( 1000.0 * ( x[i] - 0.6 ) ), 6 );
  }
  return fx;
}
/******************************************************************************/

void p21_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P21_LIM returns the integration limits for problem 21.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p22_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P22_EXACT returns the estimated integral for problem 22.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P22_EXACT, the estimated value of the integral.
//
{
  double exact;
  const double r8_pi = 3.141592653589793;

  exact = 0.125 * log ( 9.0 ) + r8_pi / sqrt ( 48.0 );

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

double *p22_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P22_FUN evaluates the integrand for problem 22.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    1 / ( x^4 + x^2 + 1 )
//
//  Exact integral:
//
//    ln ( 9 ) / 8 + pi / sqrt ( 48 )
//
//  Approximate Integral (20 digits):
//
//    0.72810291322558188550...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  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 N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P22_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = 1.0 / ( pow ( x[i], 4 ) + pow ( x[i], 2 ) + 1.0 );
  }
  return fx;
}
/******************************************************************************/

void p22_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P22_LIM returns the integration limits for problem 22.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p23_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P23_EXACT returns the estimated integral for problem 23.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P23_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = 0.62471325642771360429;

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

double *p23_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P23_FUN evaluates the integrand for problem 23.
//
//  Discussion:
//
//    The integrand has a singularity at X = 0.
//    The integrand is discontinuous at X = 0.
//    The integrand is arbitrarily oscillatory as X decreases to 0.
//    The integrand becomes unbounded as X decreases to 0.
//
//    Integral ( 0 < X < 1 ) ( 1 / X ) sin ( 1 / X ) dX
//    = Integral ( 1 < X < Infinity ) ( 1 / X ) * sin ( X ) dX.
//
//  Interval:
//
//    0 <= X <= 1
//
//  Integrand:
//
//    ( 1 / x ) sin ( 1 / x )
//
//  Approximate Integral (20 digits):
//
//    0.62471325642771360429...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  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 N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P23_FUN[N], the integrand values.
//
{
  double *fx;
  int i;
  
  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] == 0.0 )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = ( 1.0 / x[i] ) * sin ( 1.0 / x[i] );
    }
  }
  return fx;
}
/******************************************************************************/

void p23_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P23_LIM returns the integration limits for problem 23.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p24_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P24_EXACT returns the estimated integral for problem 24.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P24_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = - 0.0067547455;

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

double *p24_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P24_FUN evaluates the integrand for problem 24.
//
//  Discussion:
//
//    The integrand is continuous, but nowhere differentiable.
//
//  Interval:
//
//    0 <= X <= 0.5
//
//  Integrand:
//
//    ( 1 / pi ) * sum ( 1 <= I < Infinity ) 2^(-I) * cos ( 7^I * pi * X )
//
//  Approximate Integral:
//
//    - 0.0067547455
//
//  Antiderivative:
//
//    ( 1 / pi^2 ) * sum ( 1 <= I < Infinity ) 14^(-I) * sin ( 7^I * pi * X )
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Philip Davis, Philip Rabinowitz,
//    Methods of Numerical Integration,
//    Second Edition,
//    Dover, 2007,
//    ISBN: 0486453391,
//    LC: QA299.3.D28.
//
//    Herbert Salzer, Norman Levine,
//    Table of a Weierstrass Continuous Nondifferentiable Function,
//    Mathematics of Computation,
//    Volume 15, pages 120 - 130, 1961.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P24_FUN[N], the integrand values.
//
{
  double *fx;
  int i;
  int j;
  static int n_term = 40;
  const double r8_pi = 3.141592653589793;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = 0.0;
    for ( j = 1; j <= n_term; j++ )
    {
      fx[i] = fx[i] + cos ( pow ( 7.0, j ) * r8_pi * x[i] ) / pow ( 2.0, j );
    }
    fx[i] = fx[i] / r8_pi;
  }
  return fx;
}
/******************************************************************************/

void p24_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P24_LIM returns the integration limits for problem 24.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 0.5;

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

double p25_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P25_EXACT returns the estimated integral for problem 25.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P25_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = 0.3 * log ( 0.3 ) + 0.7 * log ( 0.7 ) - 1.0;

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

double *p25_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P25_FUN evaluates the integrand for problem 25.
//
//  Interval:
//
//    0 <= X <= 1.
//
//  Integrand:
//
//    ln ( abs ( x - 0.7 ) )
//
//  Exact Integral:
//
//    0.3 * ln ( 0.3 ) + 0.7 * ln ( 0.7 ) - 1
//
//  Approximate Integral (20 digits):
//
//    -1.6108643020548934630
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Kendall Atkinson,
//    An Introduction to Numerical Analysis,
//    Prentice Hall, 1984, page 303.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P25_FUN[N], the integrand values.
//
{
  double *fx;
  int i;
  
  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] == 0.7 )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = log ( r8_abs ( x[i] - 0.7 ) );
    }
  }
  return fx;
}
/******************************************************************************/

void p25_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P25_LIM returns the integration limits for problem 25.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p26_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P26_EXACT returns the exact integral for problem 26.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P26_EXACT, the value of the integral.
//
{
  double exact;

  exact = 7.9549265210128452745;

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

double *p26_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P26_FUN evaluates the integrand for problem 26.
//
//  Interval:
//
//    0 <= x <= 2 pi
//
//  Integrand:
//
//    exp ( cos ( x ) )
//
//  Exact Integral:
//
//    2 * Pi * I0(1)
//
//  Approximate Integral (20 digits):
//
//    7.9549265210128452745...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Kendall Atkinson,
//    An Introduction to Numerical Analysis,
//    Prentice Hall, 1984, page 262.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P26_FUN[N], the integrand values.
//
{
  double *fx;
  int i;
  
  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = exp ( cos ( x[i] ) );
  }
  return fx;
}
/******************************************************************************/

void p26_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P26_LIM returns the integration limits for problem 26.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  const double r8_pi = 3.141592653589793;

  *a = 0.0;
  *b = 2.0 * r8_pi;

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

double p27_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P27_EXACT returns the exact integral for problem 27.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P27_EXACT, the value of the integral.
//
{
  double exact;

  exact = 5.0 - 6.0 * log ( 2.0 );

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

double *p27_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P27_FUN evaluates the integrand for problem 27.
//
//  Interval:
//
//    0 <= X <= 1
//
//  Integrand:
//
//    1 / ( X^(1/2) + X^(1/3) )
//
//  Exact Integral:
//
//    5 - 6 * ln ( 2 )
//
//  Approximate Integral (20 digits):
//
//    0.84111691664032814350...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  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 N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P27_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] == 0.0 )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = 1.0 / ( sqrt ( x[i] ) + pow ( x[i], 1.0 / 3.0 ) );
    }
  }
  return fx;
}
/******************************************************************************/

void p27_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P27_LIM returns the integration limits for problem 27.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p28_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P28_EXACT returns the exact integral for problem 28.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P28_EXACT, the value of the integral.
//
{
  double exact;
  const double r8_pi = 3.141592653589793;

  exact = ( 50.0 / 2501.0 ) * ( 1.0 - exp ( - 2.0 * r8_pi ) );

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

double *p28_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P28_FUN evaluates the integrand for problem 28.
//
//  Interval:
//
//    0 <= X <= 2 PI
//
//  Integrand:
//
//    exp ( - X ) * sin ( 50 * X )
//
//  Exact Integral:
//
//    50 / ( 2501 ) * ( 1 - exp ( - 2 * PI ) )
//
//  Approximate Integral (20 digits):
//
//    0.019954669277654778312...
//
//  Reference:
//
//    Kendall Atkinson,
//    An Introduction to Numerical Analysis,
//    Prentice Hall, 1984, page 303.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P28_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = exp ( - x[i] ) * sin ( 50.0 * x[i] );
  }
  return fx;
}
/******************************************************************************/

void p28_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P28_LIM returns the integration limits for problem 28.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  const double r8_pi = 3.141592653589793;

  *a = 0.0;
  *b = 2.0 * r8_pi;

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

double p29_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P29_EXACT returns the exact integral for problem 29.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P29_EXACT, the value of the integral.
//
{
  double exact;

  exact = 1.0 - log ( 2.0 );

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

double *p29_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P29_FUN evaluates the integrand for problem 29.
//
//  Interval:
//
//    0 <= X <= 1
//
//  Integrand:
//
//    F ( X ) = 1 / ( X + 2 )   for 0 < X < E - 2
//            = 0               otherwise
//
//  Exact Integral:
//
//    1 - ln ( 2 )
//
//  Approximate Integral (20 digits):
//
//    0.30685281944005469058...
//
//  Reference:
//
//    Philip Davis, Philip Rabinowitz,
//    Methods of Numerical Integration,
//    Second Edition,
//    Dover, 2007,
//    ISBN: 0486453391,
//    LC: QA299.3.D28.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P29_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( 0.0 <= x[i] && x[i] <= exp ( 1.0 ) - 2.0 )
    {
      fx[i] = 1.0 / ( x[i] + 2.0 );
    }
    else
    {
      fx[i] = 0.0;
    }
  }
  return fx;
}
/******************************************************************************/

void p29_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P29_LIM returns the integration limits for problem 29.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p30_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P30_EXACT returns the exact integral for problem 30.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P30_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = -4.5275696251606720278;

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

double *p30_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P30_FUN evaluates the integrand for problem 30.
//
//  Interval:
//
//    2 <= x <= 7
//
//  Integrand:
//
//          cos (       x )
//    + 5 * cos ( 1.6 * x )
//    - 2 * cos ( 2.0 * x )
//    + 5 * cos ( 4.5 * x )
//    + 7 * cos ( 9.0 * x )
//
//  Antiderivative:
//
//          sin (       x )
//    + 5 * sin ( 1.6 * x ) / 1.6
//    - 2 * sin ( 2.0 * x ) / 2.0
//    + 5 * sin ( 4.5 * x ) / 4.5
//    + 7 * sin ( 9.0 * x ) / 9.0
//
//  Exact Integral:
//
//    -4.5275696251606720278
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Dianne OLeary,
//    Scientific Computing with Case Studies,
//    SIAM, 2008,
//    ISBN13: 978-0-898716-66-5,
//    LC: QA401.O44.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P30_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = 
              cos (       x[i] ) 
      + 5.0 * cos ( 1.6 * x[i] ) 
      - 2.0 * cos ( 2.0 * x[i] ) 
      + 5.0 * cos ( 4.5 * x[i] ) 
      + 7.0 * cos ( 9.0 * x[i] );
  }
  return fx;
}
/******************************************************************************/

void p30_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P30_LIM returns the integration limits for problem 30.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 2.0;
  *b = 7.0;

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

double p31_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P31_EXACT returns the exact integral for problem 31.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P31_EXACT, the value of the integral.
//
{
  double exact;

  exact = 2.0 * atan ( 4.0 );

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

double *p31_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P31_FUN evaluates the integrand for problem 31.
//
//  Discussion:
//
//    A simple Newton-Cotes quadrature rule, in which the order of the
//    rule is increased, but the interval is not subdivided, diverges
//    for this integrand.
//
//    This is Runge's function, a standard example of the perils of
//    using high order polynomial interpolation at equally spaced nodes.
//    Since this is exactly what is really going on in a Newton Cotes
//    rule, it is little wonder that the result is so poor.
//
//  Interval:
//
//    -4 <= x <= 4
//
//  Integrand:
//
//    1 / ( 1 + x^2 )
//
//  Antiderivative:
//
//    arctan ( x )
//
//  Exact Integral:
//
//    2 * arctan ( 4 )
//
//  Approximate Integral (20 digits):
//
//    2.6516353273360649301...
//
//  Reference:
//
//    Kendall Atkinson,
//    An Introduction to Numerical Analysis,
//    Prentice Hall, 1984, page 266.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P31_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = 1.0 / ( 1.0 + x[i] * x[i] );
  }
  return fx;
}
/******************************************************************************/

void p31_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P31_LIM returns the integration limits for problem 31.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = - 4.0;
  *b =   4.0;

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

double p32_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P32_EXACT returns the exact integral for problem 32.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P32_EXACT, the value of the integral.
//
{
  double exact;
  const double r8_pi = 3.141592653589793;

  exact = - 0.5 * ( exp ( r8_pi ) + 1.0 );

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

double *p32_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P32_FUN evaluates the integrand for problem 32.
//
//  Interval:
//
//    0 <= X <= PI
//
//  Integrand:
//
//    exp ( X ) * cos ( X )
//
//  Antiderivative:
//
//    0.5 * exp ( X ) * ( sin ( X ) + cos ( X ) )
//
//  Exact Integral:
//
//    - 0.5 * ( exp ( PI ) + 1 )
//
//  Approximate Integral (20 digits):
//
//    -12.070346316389634503...
//
//  Reference:
//
//    Kendall Atkinson,
//    An Introduction to Numerical Analysis,
//    Prentice Hall, 1984, page 254, 277, 297.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P32_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = exp ( x[i] ) * cos ( x[i] );
  }
  return fx;
}
/******************************************************************************/

void p32_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P32_LIM returns the integration limits for problem 32.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  const double r8_pi = 3.141592653589793;

  *a = 0.0;
  *b = r8_pi;

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

double p33_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P33_EXACT returns the exact integral for problem 33.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P33_EXACT, the value of the integral.
//
{
  double exact;
  const double r8_pi = 3.141592653589793;

  exact = 0.5 * sqrt ( r8_pi );

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

double *p33_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P33_FUN evaluates the integrand for problem 33.
//
//  Discussion:
//
//    The integrand is singular at both endpoints of the interval.
//
//  Interval:
//
//    0 <= X <= 1
//
//  Integrand:
//
//    sqrt ( - ln ( X ) )
//
//  Exact Integral:
//
//    sqrt ( pi ) / 2
//
//  Approximate Integral (20 digits):
//
//    0.88622692545275801365...
//
//  Reference:
//
//    Kendall Atkinson,
//    An Introduction to Numerical Analysis,
//    Prentice Hall, 1984, page 307.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P33_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] <= 0.0 )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = sqrt ( - log ( x[i] ) );
    }
  }
  return fx;
}
/******************************************************************************/

void p33_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P33_LIM returns the integration limits for problem 33.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p34_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P34_EXACT returns the exact integral for problem 34.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P34_EXACT, the value of the integral.
//
{
  double exact;

  exact = 1627879.0 / 1500.0;

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

double *p34_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P34_FUN evaluates the integrand for problem 34.
//
//  Interval:
//
//    0 <= X <= 1
//
//  Integrand:
//
//    ( 10 * X - 1 ) * ( 10 * X - 1.1 ) * ( 10 * X - 1.2 ) * ( 10 * X - 1.3 )
//
//  Exact Integral:
//
//    1627879 / 1500
//
//  Approximate Integral (20 digits):
//
//    1085.2526666666666666...
//
//  Reference:
//
//    Hermann Engels,
//    Numerical Quadrature and Cubature,
//    Academic Press, 1980.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P34_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = ( 10.0 * x[i] - 1.0 ) * ( 10.0 * x[i] - 1.1 )
      * ( 10.0 * x[i] - 1.2 ) * ( 10.0 * x[i] - 1.3 );
  }
  return fx;
}
/******************************************************************************/

void p34_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P34_LIM returns the integration limits for problem 34.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p35_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P35_EXACT returns the exact integral for problem 35.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P35_EXACT, the value of the integral.
//
{
  double exact;

  exact = 26.0;

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

double *p35_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P35_FUN evaluates the integrand for problem 35.
//
//  Interval:
//
//    -9 <= X <= 100
//
//  Integrand:
//
//    1 / sqrt ( abs ( X ) )
//
//  Exact Integral:
//
//    26
//
//  Reference:
//
//    Hermann Engels,
//    Numerical Quadrature and Cubature,
//    Academic Press, 1980.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P35_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] == 0.0 )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = 1.0 / sqrt ( r8_abs ( x[i] ) );
    }
  }
  return fx;
}
/******************************************************************************/

void p35_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P35_LIM returns the integration limits for problem 35.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = -9.0;
  *b = 100.0;

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

double p36_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P36_EXACT returns the exact integral for problem 36.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P36_EXACT, the value of the integral.
//
{
  double alpha;
  double exact;

  alpha = p36_param_get ( );

  exact = 1.0 / pow ( alpha + 1.0, 2 );

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

double *p36_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P36_FUN evaluates the integrand for problem 36.
//
//  Discussion:
//
//    The problem has a parameter ALPHA that can be set by calling
//    P36_PARAM_SET.  It had a default value of -0.9.
//
//    The integrand has an endpoint singularity at X=0.
//
//    Suggested values of ALPHA include -0.9 through 2.6.
//
//  Interval:
//
//    0 <= X <= 1
//
//  Integrand:
//
//    X^alpha * ln ( 1 / X )
//
//  Exact Integral:
//
//    1 / ( alpha + 1 )^2
//
//  Reference:
//
//    Robert Piessens, Elise de Doncker-Kapenga, 
//    Christian Ueberhuber, David Kahaner,
//    QUADPACK: A Subroutine Package for Automatic Integration,
//    Springer, 1983, page 83.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P36_FUN[N], the integrand values.
//
{
  double alpha;
  double *fx;
  int i;

  alpha = p36_param_get ( );

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] <= 0.0 )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = pow ( x[i], alpha ) * log ( 1.0 / x[i] );
    }
  }
  return fx;
}
/******************************************************************************/

void p36_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P36_LIM returns the integration limits for problem 36.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

void p36_param ( char *action, char *name, double *value )

/******************************************************************************/
//
//  Purpose:
//
//    P36_PARAM gets or sets the parameter values for problem 36.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *ACTION.
//    "GET" to get the value.
//    "SET" to set the value.
//
//    Input, char *NAME, the name of the parameter.
//    "ALPHA" is the only option.
//
//    Input/output, double *VALUE.
//    If the action is "GET", then VALUE returns the current parameter value.
//    If ACTION is "SET", then the parameter value is set to VALUE.
//
{
  static double alpha = - 0.9;

  if ( s_eqi ( action, "get" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      *value = alpha;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P36_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else if ( s_eqi ( action, "set" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      alpha = *value;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P36_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P36_PARAM - Fatal error!\n" );
    fprintf ( stderr, "  Unrecognized action.\n" );
    exit ( 1 );
  }

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

double p36_param_get ( )

/******************************************************************************/
//
//  Purpose:
//
//    P36_PARAM_GET returns the parameter values for problem 36.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    29 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double ALPHA, the current value of the parameter.
//
{
  double alpha;

  p36_param ( "get", "alpha", &alpha );

  return alpha;
}
/******************************************************************************/

void p36_param_set ( double alpha )

/******************************************************************************/
//
//  Purpose:
//
//    P36_PARAM_SET sets the parameter values for problem 36.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ALPHA, the new value of the parameter.
//
{
  p36_param ( "set", "alpha", &alpha );

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

double p37_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P37_EXACT returns the exact integral for problem 37.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P37_EXACT, the value of the integral.
//
{
  double alpha;
  double exact;
  const double r8_pi = 3.141592653589793;

  alpha = p37_param_get ( );

  exact = atan ( ( 4.0 - r8_pi ) * pow ( 4.0, alpha - 1.0 ) ) 
        + atan (         r8_pi   * pow ( 4.0, alpha - 1.0 ) );

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

double *p37_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P37_FUN evaluates the integrand for problem 37.
//
//  Discussion:
//
//    The problem has a parameter ALPHA that can be set by calling
//    P37_PARAM_SET.  It had a default value of 5.0.
//
//    The integrand has a peak of height 4^ALPHA at X = PI/4.
//
//    Suggested values of ALPHA include 0 through 20.
//
//  Interval:
//
//    0 <= X <= 1
//
//  Integrand:
//
//    4^(-ALPHA) / ( (X-PI/4)^2 + 16^(-ALPHA) )
//
//  Exact Integral:
//
//    atan ( ( 4 - PI ) * 4^(ALPHA-1) ) + atan ( PI * 4^(ALPHA-1) )
//
//  Reference:
//
//    Robert Piessens, Elise de Doncker-Kapenga, 
//    Christian Ueberhuber, David Kahaner,
//    QUADPACK: A Subroutine Package for Automatic Integration,
//    Springer, 1983, page 83.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P37_FUN[N], the integrand values.
//
{
  double alpha;
  double *fx;
  int i;
  const double r8_pi = 3.141592653589793;

  alpha = p37_param_get ( );

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = pow ( 4.0, -alpha ) 
      / ( pow ( x[i] - 0.25 * r8_pi, 2 ) + pow ( 16.0, - alpha ) );
  }
  return fx;
}
/******************************************************************************/

void p37_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P37_LIM returns the integration limits for problem 37.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

void p37_param ( char *action, char *name, double *value )

/******************************************************************************/
//
//  Purpose:
//
//    P37_PARAM gets or sets the parameter values for problem 37.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *ACTION.
//    "GET" to get the value.
//    "SET" to set the value.
//
//    Input, char *NAME, the name of the parameter.
//    "ALPHA" is the only option.
//
//    Input/output, double *VALUE.
//    If the action is "GET", then VALUE returns the current parameter value.
//    If ACTION is "SET", then the parameter value is set to VALUE.
//
{
  static double alpha = 5.0;

  if ( s_eqi ( action, "get" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      *value = alpha;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P37_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else if ( s_eqi ( action, "set" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      alpha = *value;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P37_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P37_PARAM - Fatal error!\n" );
    fprintf ( stderr, "  Unrecognized action.\n" );
    exit ( 1 );
  }

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

double p37_param_get ( )

/******************************************************************************/
//
//  Purpose:
//
//    P37_PARAM_GET returns the parameter values for problem 37.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double ALPHA, the current value of the parameter.
//
{
  double alpha;

  p37_param ( "get", "alpha", &alpha );

  return alpha;
}
/******************************************************************************/

void p37_param_set ( double alpha )

/******************************************************************************/
//
//  Purpose:
//
//    P37_PARAM_SET sets the parameter values for problem 37.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ALPHA, the new value of the parameter.
//
{
  p37_param ( "set", "alpha", &alpha );

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

double p38_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P38_EXACT returns the exact integral for problem 38.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P38_EXACT, the value of the integral.
//
{
  double alpha;
  double exact;
  const double r8_pi = 3.141592653589793;
  double x;

  alpha = p38_param_get ( );

  x = pow ( 2.0, alpha );

  exact = r8_pi * r8_besj0 ( x );

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

double *p38_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P38_FUN evaluates the integrand for problem 38.
//
//  Discussion:
//
//    The problem has a parameter ALPHA that can be set by calling
//    P38_PARAM_SET.
//
//    The integrand oscillates more strongly as ALPHA is increased.
//
//    The suggested range for ALPHA is 0 to 10.
//
//  Interval:
//
//    0 <= X <= PI
//
//  Integrand:
//
//    cos ( 2^ALPHA * sin ( x ) )
//
//  Exact Integral:
//
//    pi * J0 ( 2^ALPHA )
//
//    where J0 ( x ) is the J Bessel function of order 0.
//
//  Reference:
//
//    Robert Piessens, Elise de Doncker-Kapenga, 
//    Christian Ueberhuber, David Kahaner,
//    QUADPACK: A Subroutine Package for Automatic Integration,
//    Springer, 1983, page 83.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P38_FUN[N], the integrand values.
//
{
  double alpha;
  double *fx;
  int i;

  alpha = p38_param_get ( );

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = cos ( pow ( 2.0, alpha ) * sin ( x[i] ) );
  }
  return fx;
}
/******************************************************************************/

void p38_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P38_LIM returns the integration limits for problem 38.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  const double r8_pi = 3.141592653589793;

  *a = 0.0;
  *b = r8_pi;

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

void p38_param ( char *action, char *name, double *value )

/******************************************************************************/
//
//  Purpose:
//
//    P38_PARAM gets or sets the parameter values for problem 38.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *ACTION.
//    "GET" to get the value.
//    "SET" to set the value.
//
//    Input, char *NAME, the name of the parameter.
//    "ALPHA" is the only option.
//
//    Input/output, double *VALUE.
//    If the action is "GET", then VALUE returns the current parameter value.
//    If ACTION is "SET", then the parameter value is set to VALUE.
//
{
  static double alpha = 3.0;

  if ( s_eqi ( action, "get" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      *value = alpha;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P38_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else if ( s_eqi ( action, "set" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      alpha = *value;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P38_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P38_PARAM - Fatal error!\n" );
    fprintf ( stderr, "  Unrecognized action.\n" );
    exit ( 1 );
  }

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

double p38_param_get ( )

/******************************************************************************/
//
//  Purpose:
//
//    P38_PARAM_GET returns the parameter values for problem 38.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double ALPHA, the current value of the parameter.
//
{
  double alpha;

  p38_param ( "get", "alpha", &alpha );

  return alpha;
}
/******************************************************************************/

void p38_param_set ( double alpha )

/******************************************************************************/
//
//  Purpose:
//
//    P38_PARAM_SET sets the parameter values for problem 38.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ALPHA, the new value of the parameter.
//
{
  p38_param ( "set", "alpha", &alpha );

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

double p39_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P39_EXACT returns the exact integral for problem 39.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P39_EXACT, the value of the integral.
//
{
  double alpha;
  double exact;

  alpha = p39_param_get ( );

  exact = ( pow ( 2.0 / 3.0, alpha + 1.0 ) 
          + pow ( 1.0 / 3.0, alpha + 1.0 ) ) / ( alpha + 1.0 );

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

double *p39_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P39_FUN evaluates the integrand for problem 39.
//
//  Discussion:
//
//    The problem has a parameter ALPHA that can be set by calling
//    P39_PARAM_SET.
//
//    The integrand has a singularity at an internal point ( x = 1/3 )
//    with small binary period.
//
//    The suggested range for ALPHA is -0.8 through 2.1.
//
//  Interval:
//
//    0 <= X <= 1
//
//  Integrand:
//
//    ( abs ( x - 1/3 ) )^alpha
//
//  Exact Integral:
//
//    ( (2/3)^(alpha+1) + (1/3)^(alpha+1) ) / ( alpha + 1 )
//
//  Reference:
//
//    Robert Piessens, Elise de Doncker-Kapenga, 
//    Christian Ueberhuber, David Kahaner,
//    QUADPACK: A Subroutine Package for Automatic Integration,
//    Springer, 1983, page 83.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P39_FUN[N], the integrand values.
//
{
  double alpha;
  double *fx;
  int i;

  alpha = p39_param_get ( );

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] - 1.0 / 3.0 == 0.0 )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = pow ( r8_abs ( x[i] - 1.0 / 3.0 ), alpha );
    }
  }
  return fx;
}
/******************************************************************************/

void p39_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P39_LIM returns the integration limits for problem 39.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

void p39_param ( char *action, char *name, double *value )

/******************************************************************************/
//
//  Purpose:
//
//    P39_PARAM gets or sets the parameter values for problem 39.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *ACTION.
//    "GET" to get the value.
//    "SET" to set the value.
//
//    Input, char *NAME, the name of the parameter.
//    "ALPHA" is the only option.
//
//    Input/output, double *VALUE.
//    If the action is "GET", then VALUE returns the current parameter value.
//    If ACTION is "SET", then the parameter value is set to VALUE.
//
{
  static double alpha = - 0.5;

  if ( s_eqi ( action, "get" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      *value = alpha;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P39_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else if ( s_eqi ( action, "set" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      alpha = *value;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P39_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P39_PARAM - Fatal error!\n" );
    fprintf ( stderr, "  Unrecognized action.\n" );
    exit ( 1 );
  }
  return;
}
/******************************************************************************/

double p39_param_get ( )

/******************************************************************************/
//
//  Purpose:
//
//    P39_PARAM_GET returns the parameter values for problem 39.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double ALPHA, the current value of the parameter.
//
{
  double alpha;

  p39_param ( "get", "alpha", &alpha );

  return alpha;
}
/******************************************************************************/

void p39_param_set ( double alpha )

/******************************************************************************/
//
//  Purpose:
//
//    P39_PARAM_SET sets the parameter values for problem 39.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ALPHA, the new value of the parameter.
//
{
  p39_param ( "set", "alpha", &alpha );

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

double p40_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P40_EXACT returns the exact integral for problem 40.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P40_EXACT, the value of the integral.
//
{
  double alpha;
  double exact;
  const double r8_pi = 3.141592653589793;

  alpha = p40_param_get ( );

  exact = ( pow ( 1.0 - 0.25 * r8_pi, alpha + 1.0 ) 
          + pow (     + 0.25 * r8_pi, alpha + 1.0 ) ) 
          / ( alpha + 1.0 );

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

double *p40_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P40_FUN evaluates the integrand for problem 40.
//
//  Discussion:
//
//    The problem has a parameter ALPHA that can be set by calling
//    P40_PARAM_SET.
//
//    The integrand has a singularity at an internal point ( x = pi/4 ).
//
//    The suggested range for ALPHA is -0.8 through 2.1.
//
//  Interval:
//
//    0 <= X <= 1
//
//  Integrand:
//
//    ( abs ( x - pi/4 ) )^alpha
//
//  Exact Integral:
//
//    ( (1-pi/4)^(alpha+1) + (pi/4)^(alpha+1) ) / ( alpha + 1 )
//
//  Reference:
//
//    Robert Piessens, Elise de Doncker-Kapenga, 
//    Christian Ueberhuber, David Kahaner,
//    QUADPACK: A Subroutine Package for Automatic Integration,
//    Springer, 1983, page 83.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P40_FUN[N], the integrand values.
//
{
  double alpha;
  double *fx;
  int i;
  const double r8_pi = 3.141592653589793;

  alpha = p40_param_get ( );

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = pow ( r8_abs ( x[i] - 0.25 * r8_pi ), alpha );
  }
  return fx;
}
/******************************************************************************/

void p40_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P40_LIM returns the integration limits for problem 40.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

void p40_param ( char *action, char *name, double *value )

/******************************************************************************/
//
//  Purpose:
//
//    P40_PARAM gets or sets the parameter values for problem 40.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *ACTION.
//    "GET" to get the value.
//    "SET" to set the value.
//
//    Input, char *NAME, the name of the parameter.
//    "ALPHA" is the only option.
//
//    Input/output, double *VALUE.
//    If the action is "GET", then VALUE returns the current parameter value.
//    If ACTION is "SET", then the parameter value is set to VALUE.
//
{
  static double alpha = - 0.5;

  if ( s_eqi ( action, "get" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      *value = alpha;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P40_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else if ( s_eqi ( action, "set" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      alpha = *value;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P40_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P40_PARAM - Fatal error!\n" );
    fprintf ( stderr, "  Unrecognized action.\n" );
    exit ( 1 );
  }

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

double p40_param_get ( )

/******************************************************************************/
//
//  Purpose:
//
//    P40_PARAM_GET returns the parameter values for problem 40.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double ALPHA, the current value of the parameter.
//
{
  double alpha;

  p40_param ( "get", "alpha", &alpha );

  return alpha;
}
/******************************************************************************/

void p40_param_set ( double alpha )

/******************************************************************************/
//
//  Purpose:
//
//    P40_PARAM_SET sets the parameter values for problem 40.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ALPHA, the new value of the parameter.
//
{
  p40_param ( "set", "alpha", &alpha );

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

double p41_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P41_EXACT returns the exact integral for problem 41.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P41_EXACT, the value of the integral.
//
{
  double alpha;
  double exact;
  const double r8_pi = 3.141592653589793;

  alpha = p41_param_get ( );

  exact = r8_pi / sqrt ( pow ( 1.0 + pow ( 2.0, - alpha ), 2 ) - 1.0 );

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

double *p41_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P41_FUN evaluates the integrand for problem 41.
//
//  Discussion:
//
//    The problem has a parameter ALPHA that can be set by calling
//    P41_PARAM_SET.
//
//    The integrand has a singularity at both endpoints, whose
//    severity increases with ALPHA.
//
//    The suggested range for ALPHA is 1 through 20.
//
//  Interval:
//
//    -1 <= X <= 1
//
//  Integrand:
//
//    1 / ( sqrt ( 1 - x^2 ) * ( x + 1 + 2^(-alpha) ) )
//
//  Exact Integral:
//
//    pi / sqrt ( ( 1 + 2^(-alpha) ) - 1 )
//
//  Reference:
//
//    Robert Piessens, Elise de Doncker-Kapenga, 
//    Christian Ueberhuber, David Kahaner,
//    QUADPACK: A Subroutine Package for Automatic Integration,
//    Springer, 1983, page 83.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P41_FUN[N], the integrand values.
//
{
  double alpha;
  double *fx;
  int i;

  alpha = p41_param_get ( );

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( 1.0 - x[i] * x[i] == 0.0 || x[i] + 1.0 + pow ( 0.5, alpha ) == 0.0 )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = 1.0 / 
        ( sqrt ( 1.0 - x[i] * x[i] ) * ( x[i] + 1.0 + pow ( 0.5, alpha ) ) );
    }
  }
  return fx;
}
/******************************************************************************/

void p41_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P41_LIM returns the integration limits for problem 41.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = -1.0;
  *b = 1.0;

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

void p41_param ( char *action, char *name, double *value )

/******************************************************************************/
//
//  Purpose:
//
//    P41_PARAM gets or sets the parameter values for problem 41.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *ACTION.
//    "GET" to get the value.
//    "SET" to set the value.
//
//    Input, char *NAME, the name of the parameter.
//    "ALPHA" is the only option.
//
//    Input/output, double *VALUE.
//    If the action is "GET", then VALUE returns the current parameter value.
//    If ACTION is "SET", then the parameter value is set to VALUE.
//
{
  static double alpha = 3.0;

  if ( s_eqi ( action, "get" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      *value = alpha;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P41_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else if ( s_eqi ( action, "set" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      alpha = *value;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P41_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P41_PARAM - Fatal error!\n" );
    fprintf ( stderr, "  Unrecognized action.\n" );
    exit ( 1 );
  }

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

double p41_param_get ( )

/******************************************************************************/
//
//  Purpose:
//
//    P41_PARAM_GET returns the parameter values for problem 41.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double ALPHA, the current value of the parameter.
//
{
  double alpha;

  p41_param ( "get", "alpha", &alpha );

  return alpha;
}
/******************************************************************************/

void p41_param_set ( double alpha )

/******************************************************************************/
//
//  Purpose:
//
//    P41_PARAM_SET sets the parameter values for problem 41.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ALPHA, the new value of the parameter.
//
{
  p41_param ( "set", "alpha", &alpha );

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

double p42_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P42_EXACT returns the exact integral for problem 42.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P42_EXACT, the value of the integral.
//
{
  double alpha;
  double exact;

  alpha = p42_param_get ( );

  exact = pow ( 2.0, alpha - 2.0 ) * pow ( r8_gamma ( alpha / 2.0 ), 2 )
    / r8_gamma ( alpha );

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

double *p42_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P42_FUN evaluates the integrand for problem 42.
//
//  Discussion:
//
//    The problem has a parameter ALPHA that can be set by calling
//    P42_PARAM_SET.
//
//    The integrand has a singularity at X = 0 if ALPHA < 1.
//
//    The suggested range for ALPHA is 0.1 through 2.
//
//  Interval:
//
//    0 <= X <= pi/2
//
//  Integrand:
//
//    ( sin(x) )^( alpha - 1 )
//
//  Exact Integral:
//
//    2^( alpha - 2 ) * ( Gamma(alpha/2) )^2 / Gamma(alpha)
//
//  Reference:
//
//    Robert Piessens, Elise de Doncker-Kapenga, 
//    Christian Ueberhuber, David Kahaner,
//    QUADPACK: A Subroutine Package for Automatic Integration,
//    Springer, 1983, page 83.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P42_FUN[N], the integrand values.
//
{
  double alpha;
  double base;
  double *fx;
  int i;

  alpha = p42_param_get ( );

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    base = sin ( x[i] );

    if ( base == 0.0 )
    {
      if ( 1.0 < alpha )
      {
        fx[i] = 0.0;
      }
      else if ( alpha == 1.0 )
      {
        fx[i] = 1.0;
      }
      else
      {
        fx[i] = 0.0;
      }
    }
    else
    {
      fx[i] = pow ( base, alpha - 1.0 );
    }
  }
  return fx;
}
/******************************************************************************/

void p42_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P42_LIM returns the integration limits for problem 42.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  const double r8_pi = 3.141592653589793;

  *a = 0.0;
  *b = r8_pi / 2.0;

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

void p42_param ( char *action, char *name, double *value )

/******************************************************************************/
//
//  Purpose:
//
//    P42_PARAM gets or sets the parameter values for problem 42.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *ACTION.
//    "GET" to get the value.
//    "SET" to set the value.
//
//    Input, char *NAME, the name of the parameter.
//    "ALPHA" is the only option.
//
//    Input/output, double *VALUE.
//    If the action is "GET", then VALUE returns the current parameter value.
//    If ACTION is "SET", then the parameter value is set to VALUE.
//
{
  static double alpha = 0.3;

  if ( s_eqi ( action, "get" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      *value = alpha;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P42_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else if ( s_eqi ( action, "set" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      alpha = *value;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P42_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P42_PARAM - Fatal error!\n" );
    fprintf ( stderr, "  Unrecognized action.\n" );
    exit ( 1 );
  }

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

double p42_param_get ( )

/******************************************************************************/
//
//  Purpose:
//
//    P42_PARAM_GET returns the parameter values for problem 42.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double ALPHA, the current value of the parameter.
//
{
  double alpha;

  p42_param ( "get", "alpha", &alpha );

  return alpha;
}
/******************************************************************************/

void p42_param_set ( double alpha )

/******************************************************************************/
//
//  Purpose:
//
//    P42_PARAM_SET sets the parameter values for problem 42.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ALPHA, the new value of the parameter.
//
{
  p42_param ( "set", "alpha", &alpha );

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

double p43_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P43_EXACT returns the exact integral for problem 43.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P43_EXACT, the value of the integral.
//
{
  double alpha;
  double exact;

  alpha = p43_param_get ( );

  exact = r8_gamma ( alpha );

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

double *p43_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P43_FUN evaluates the integrand for problem 43.
//
//  Discussion:
//
//    The problem has a parameter ALPHA that can be set by calling
//    P43_PARAM_SET.
//
//    The suggested parameter range is 0.1 <= ALPHA <= 2.0.
//
//    The integrand has an algebraic endpoint singularity at X = 1
//    times a singular factor.
//
//  Interval:
//
//    0 <= X <= 1
//
//  Integrand:
//
//    ( ln ( 1 / x ) )^( alpha - 1 )
//
//  Exact Integral:
//
//    Gamma(alpha)
//
//  Reference:
//
//    Robert Piessens, Elise de Doncker-Kapenga, 
//    Christian Ueberhuber, David Kahaner,
//    QUADPACK: A Subroutine Package for Automatic Integration,
//    Springer, 1983, page 84.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P43_FUN[N], the integrand values.
//
{
  double alpha;
  double *fx;
  int i;

  alpha = p43_param_get ( );

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] <= 0.0 )
    {
      fx[i] = 0.0;
    }
    else if ( x[i] == 0.0 )
    {
      if ( alpha - 1.0 < 0.0 )
      {
        fx[i] = 0.0;
      }
      else if ( alpha - 1.0 == 0.0 )
      {
        fx[i] = 1.0;
      }
      else
      {
        fx[i] = 0.0;
      }
    }
    else if ( x[i] == 1.0 )
    {
      if ( alpha - 1.0 < 0.0 )
      {
        fx[i] = 0.0;
      }
      else if ( alpha - 1.0 == 0.0 )
      {
        fx[i] = 1.0;
      }
      else
      {
        fx[i] = 0.0;
      }
    }
    else
    {
      fx[i] = pow ( log ( 1.0 / x[i] ), alpha - 1.0 );
    }
  }
  return fx;
}
/******************************************************************************/

void p43_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P43_LIM returns the integration limits for problem 43.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

void p43_param ( char *action, char *name, double *value )

/******************************************************************************/
//
//  Purpose:
//
//    P43_PARAM gets or sets the parameter values for problem 43.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *ACTION.
//    "GET" to get the value.
//    "SET" to set the value.
//
//    Input, char *NAME, the name of the parameter.
//    "ALPHA" is the only option.
//
//    Input/output, double *VALUE.
//    If the action is "GET", then VALUE returns the current parameter value.
//    If ACTION is "SET", then the parameter value is set to VALUE.
//
{
  static double alpha = 0.3;

  if ( s_eqi ( action, "get" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      *value = alpha;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P43_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else if ( s_eqi ( action, "set" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      alpha = *value;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P43_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P43_PARAM - Fatal error!\n" );
    fprintf ( stderr, "  Unrecognized action.\n" );
    exit ( 1 );
  }

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

double p43_param_get ( )

/******************************************************************************/
//
//  Purpose:
//
//    P43_PARAM_GET returns the parameter values for problem 43.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double ALPHA, the current value of the parameter.
//
{
  double alpha;

  p43_param ( "get", "alpha", &alpha );

  return alpha;
}
/******************************************************************************/

void p43_param_set ( double alpha )

/******************************************************************************/
//
//  Purpose:
//
//    P43_PARAM_SET sets the parameter values for problem 43.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ALPHA, the new value of the parameter.
//
{
  p43_param ( "set", "alpha", &alpha );

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

double p44_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P44_EXACT returns the exact integral for problem 44.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P44_EXACT, the value of the integral.
//
{
  double alpha;
  double exact;
  double t;

  alpha = p44_param_get ( );

  t = pow ( 2.0, alpha );

  exact = ( 20.0 * sin ( t ) - t * cos ( t ) 
    + t * exp ( - 20.0 ) ) / ( 400.0 + t * t );

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

double *p44_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P44_FUN evaluates the integrand for problem 44.
//
//  Discussion:
//
//    The problem has a parameter ALPHA that can be set by calling
//    P44_PARAM_SET.
//
//    The suggested parameter range is 0.0 <= ALPHA <= 9.0.
//
//    As ALPHA increases, the integrand becomes more oscillatory.
//
//  Interval:
//
//    0 <= X <= 1
//
//  Integrand:
//
//    exp ( 20 * ( x - 1 ) ) * sin ( 2^alpha * x )
//
//  Exact Integral:
//
//    ( 20 sin ( 2^alpha ) - 2^alpha cos ( 2^alpha )
//    + 2^alpha exp ( -20 ) ) / ( 400 + 4^alpha )
//
//  Reference:
//
//    Robert Piessens, Elise de Doncker-Kapenga, 
//    Christian Ueberhuber, David Kahaner,
//    QUADPACK: A Subroutine Package for Automatic Integration,
//    Springer, 1983, page 84.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P44_FUN[N], the integrand values.
//
{
  double alpha;
  double *fx;
  int i;

  alpha = p44_param_get ( );

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = exp ( 20.0 * ( x[i] - 1.0 ) ) * sin ( pow ( 2.0, alpha ) * x[i] );
  }
  return fx;
}
/******************************************************************************/

void p44_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P44_LIM returns the integration limits for problem 44.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

void p44_param ( char *action, char *name, double *value )

/******************************************************************************/
//
//  Purpose:
//
//    P44_PARAM gets or sets the parameter values for problem 44.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *ACTION.
//    "GET" to get the value.
//    "SET" to set the value.
//
//    Input, char *NAME, the name of the parameter.
//    "ALPHA" is the only option.
//
//    Input/output, double *VALUE.
//    If the action is "GET", then VALUE returns the current parameter value.
//    If ACTION is "SET", then the parameter value is set to VALUE.
//
{
  static double alpha = 2.0;

  if ( s_eqi ( action, "get" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      *value = alpha;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P44_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else if ( s_eqi ( action, "set" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      alpha = *value;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P44_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P44_PARAM - Fatal error!\n" );
    fprintf ( stderr, "  Unrecognized action.\n" );
    exit ( 1 );
  }

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

double p44_param_get ( )

/******************************************************************************/
//
//  Purpose:
//
//    P44_PARAM_GET returns the parameter values for problem 44.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double ALPHA, the current value of the parameter.
//
{
  double alpha;

  p44_param ( "get", "alpha", &alpha );

  return alpha;
}
/******************************************************************************/

void p44_param_set ( double alpha )

/******************************************************************************/
//
//  Purpose:
//
//    P44_PARAM_SET sets the parameter values for problem 44.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ALPHA, the new value of the parameter.
//
{
  p44_param ( "set", "alpha", &alpha );

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

double p45_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P45_EXACT returns the exact integral for problem 45.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P45_EXACT, the value of the integral.
//
{
  double alpha;
  double exact;
  const double r8_pi = 3.141592653589793;
  double t;

  alpha = p45_param_get ( );

  t = pow ( 2.0, alpha - 1.0 );

  exact = r8_pi * cos ( t ) * r8_besj0 ( t );

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

double *p45_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P45_FUN evaluates the integrand for problem 45.
//
//  Discussion:
//
//    The problem has a parameter ALPHA that can be set by calling
//    P45_PARAM_SET.
//
//    The suggested parameter range is 0.0 <= ALPHA <= 8.0.
//
//    The function is singular at 0 and 1.
//
//  Interval:
//
//    0 <= X <= 1
//
//  Integrand:
//
//    cos ( 2^alpha * x ) / sqrt ( x * ( 1 - x ) )
//
//  Exact Integral:
//
//    pi * cos ( 2^(alpha-1) ) * J0 ( 2^(alpha-1) )
//
//  Reference:
//
//    Robert Piessens, Elise de Doncker-Kapenga, 
//    Christian Ueberhuber, David Kahaner,
//    QUADPACK: A Subroutine Package for Automatic Integration,
//    Springer, 1983, page 84.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P45_FUN[N], the integrand values.
//
{
  double alpha;
  double *fx;
  int i;

  alpha = p45_param_get ( );

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] == 0.0 )
    {
      fx[i] = 0.0;
    }
    else if ( x[i] == 1.0 )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = cos ( pow ( 2.0, alpha ) * x[i] ) 
        / sqrt ( x[i] * ( 1.0 - x[i] ) );
    }
  }
  return fx;
}
/******************************************************************************/

void p45_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P45_LIM returns the integration limits for problem 45.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

void p45_param ( char *action, char *name, double *value )

/******************************************************************************/
//
//  Purpose:
//
//    P45_PARAM gets or sets the parameter values for problem 45.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *ACTION.
//    "GET" to get the value.
//    "SET" to set the value.
//
//    Input, char *NAME, the name of the parameter.
//    "ALPHA" is the only option.
//
//    Input/output, double *VALUE.
//    If the action is "GET", then VALUE returns the current parameter value.
//    If ACTION is "SET", then the parameter value is set to VALUE.
//
{
  static double alpha = 2.0;

  if ( s_eqi ( action, "get" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      *value = alpha;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P45_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else if ( s_eqi ( action, "set" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      alpha = *value;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P45_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P45_PARAM - Fatal error!\n" );
    fprintf ( stderr, "  Unrecognized action.\n" );
    exit ( 1 );
  }

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

double p45_param_get ( )

/******************************************************************************/
//
//  Purpose:
//
//    P45_PARAM_GET returns the parameter values for problem 45.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double ALPHA, the current value of the parameter.
//
{
  double alpha;

  p45_param ( "get", "alpha", &alpha );

  return alpha;
}
/******************************************************************************/

void p45_param_set ( double alpha )

/******************************************************************************/
//
//  Purpose:
//
//    P45_PARAM_SET sets the parameter values for problem 45.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ALPHA, the new value of the parameter.
//
{
  p45_param ( "set", "alpha", &alpha );

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

double p46_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P46_EXACT returns the exact integral for problem 46.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    26 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P46_EXACT, the value of the integral.
//
{
  double exact;

  exact = 6.0690909595647754101;

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

double *p46_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P46_FUN evaluates the integrand for problem 46.
//
//  Discussion:
//
//    The problem has a parameter ALPHA that can be set by calling
//    P46_PARAM_SET.
//
//    The integrand is the radius of an ellipse as a function of angle.
//
//    The integral represents the arc length of the ellipse.
//
//    The suggested parameter range is 0.0 <= ALPHA < 1.0.  ALPHA is
//    the eccentricity of the ellipse.
//
//  Interval:
//
//    0 <= theta <= 2 pi
//
//  Integrand:
//
//    r(theta) = ( 1 - alpha^2 ) / ( 1 - alpha * cos ( theta ) )
//
//  Exact Integral:
//
//    When alpha = sin ( pi / 12 ), then
//
//      6.0690909595647754101
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    26 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Richard Crandall,
//    Projects in Scientific Computing,
//    Springer, 2000, page 47.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P46_FUN[N], the integrand values.
//
{
  double alpha;
  double *fx;
  int i;

  alpha = p46_param_get ( );

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = ( 1.0 - alpha * alpha ) / ( 1.0 - alpha * cos ( x[i] ) );
  }
  return fx;
}
/******************************************************************************/

void p46_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P46_LIM returns the integration limits for problem 46.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    20 March 2002
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  const double r8_pi = 3.141592653589793;

  *a = 0.0;
  *b = 2.0 * r8_pi;

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

void p46_param ( char *action, char *name, double *value )

/******************************************************************************/
//
//  Purpose:
//
//    P46_PARAM gets or sets the parameter values for problem 46.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *ACTION.
//    "GET" to get the value.
//    "SET" to set the value.
//
//    Input, char *NAME, the name of the parameter.
//    "ALPHA" is the only option.
//
//    Input/output, double *VALUE.
//    If the action is "GET", then VALUE returns the current parameter value.
//    If ACTION is "SET", then the parameter value is set to VALUE.
//
{
  static double alpha = 0.0;
  const double r8_pi = 3.141592653589793;
  static int set = 0;

  if ( !set )
  {
    alpha = sin ( r8_pi / 12.0 );
    set = 1;
  }

  if ( s_eqi ( action, "get" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      *value = alpha;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P46_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else if ( s_eqi ( action, "set" ) )
  {
    if ( s_eqi ( name, "alpha" ) )
    {
      alpha = *value;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P46_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P46_PARAM - Fatal error!\n" );
    fprintf ( stderr, "  Unrecognized action.\n" );
    exit ( 1 );
  }

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

double p46_param_get ( )

/******************************************************************************/
//
//  Purpose:
//
//    P46_PARAM_GET returns the parameter values for problem 46.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    20 March 2002
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double ALPHA, the current value of the parameter.
//
{
  double alpha;

  p46_param ( "get", "alpha", &alpha );

  return alpha;
}
/******************************************************************************/

void p46_param_set ( double alpha )

/******************************************************************************/
//
//  Purpose:
//
//    P46_PARAM_SET sets the parameter values for problem 46.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    20 March 2002
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ALPHA, the new value of the parameter.
//
{
  p46_param ( "set", "alpha", &alpha );

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

double p47_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P47_EXACT returns the exact integral for problem 47.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P47_EXACT, the value of the integral.
//
{
  double exact;

  exact = - 4.0 / 9.0;

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

double *p47_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P47_FUN evaluates the integrand for problem 47.
//
//  Discussion:
//
//    The function is singular at the left endpoint.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    sqrt ( x ) * ln ( x )
//
//  Exact Integral:
//
//    -4/9 = -0.4444...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Robert Piessens, Elise de Doncker-Kapenga, 
//    Christian Ueberhuber, David Kahaner,
//    QUADPACK: A Subroutine Package for Automatic Integration,
//    Springer, 1983, page 101.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P47_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] == 0.0 )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = sqrt ( x[i] ) * log ( x[i] );
    }
  }
  return fx;
}
/******************************************************************************/

void p47_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P47_LIM returns the integration limits for problem 47.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p48_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P48_EXACT returns the exact integral for problem 48.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P48_EXACT, the value of the integral.
//
{
  double exact;

  exact = -4.0;

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

double *p48_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P48_FUN evaluates the integrand for problem 48.
//
//  Discussion:
//
//    The function is singular at the left endpoint.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    ln ( x ) / sqrt ( x )
//
//  Exact Integral:
//
//    -4
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Robert Piessens, Elise de Doncker-Kapenga, 
//    Christian Ueberhuber, David Kahaner,
//    QUADPACK: A Subroutine Package for Automatic Integration,
//    Springer, 1983, page 103.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P48_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] == 0.0 )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = log ( x[i] ) / sqrt ( x[i] );
    }
  }
  return fx;
}
/******************************************************************************/

void p48_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P48_LIM returns the integration limits for problem 48.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p49_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P49_EXACT returns the exact integral for problem 49.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P49_EXACT, the value of the integral.
//
{
  double exact;

  exact = 61.0 * log ( 2.0 ) + 77.0 * log ( 7.0 ) / 4.0 - 27.0;

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

double *p49_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P49_FUN evaluates the integrand for problem 49.
//
//  Discussion:
//
//    The function is singular at two internal points, 1 and sqrt(2).
//
//  Interval:
//
//    0 <= x <= 3
//
//  Integrand:
//
//    x^3 * log ( abs ( ( x^2 - 1 ) * ( x^2 - 2 ) ) )
//
//  Exact Integral:
//
//    61 log ( 2 ) + (77/4) log ( 7 ) - 27 = 52.7408...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Robert Piessens, Elise de Doncker-Kapenga, 
//    Christian Ueberhuber, David Kahaner,
//    QUADPACK: A Subroutine Package for Automatic Integration,
//    Springer, 1983, page 104.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P49_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( ( x[i] * x[i] - 1.0 ) * ( x[i] * x[i] - 2.0 ) == 0.0 )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = pow ( x[i], 3 ) * log ( r8_abs ( ( x[i] * x[i] - 1.0 ) 
        * ( x[i] * x[i] - 2.0 ) ) );
    }
  }
  return fx;
}
/******************************************************************************/

void p49_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P49_LIM returns the integration limits for problem 49.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 3.0;

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

double p50_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P50_EXACT returns the exact integral for problem 50.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P50_EXACT, the value of the integral.
//
{
  double exact;
  const double r8_pi = 3.141592653589793;
  double t;

  t = 10.0 * r8_pi;

  exact = ( - euler_constant ( ) - log ( t ) + r8_ci ( t ) ) / t;

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

double *p50_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P50_FUN evaluates the integrand for problem 50.
//
//  Discussion:
//
//    The function has a removable singularity at x = 0.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    log ( x ) * sin ( 10 * pi * x )
//
//  Exact Integral:
//
//    ( - gamma - log ( 10 * pi ) + Ci ( 10 * pi ) ) / 10 * pi = -0.1281316...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Robert Piessens, Elise de Doncker-Kapenga, 
//    Christian Ueberhuber, David Kahaner,
//    QUADPACK: A Subroutine Package for Automatic Integration,
//    Springer, 1983, page 106.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P50_FUN[N], the integrand values.
//
{
  double *fx;
  int i;
  const double r8_pi = 3.141592653589793;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] == 0.0 )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = log ( x[i] ) * sin ( 10.0 * r8_pi * x[i] );
    }
  }
  return fx;
}
/******************************************************************************/

void p50_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P50_LIM returns the integration limits for problem 50.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p51_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P51_EXACT returns the exact integral for problem 51.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P51_EXACT, the value of the integral.
//
{
  double exact;
  const double r8_pi = 3.141592653589793;

  exact = - ( r8_ci ( 1.0 ) * sin ( 1.0 ) + 
    ( 0.5 * r8_pi  - r8_si ( 1.0 ) ) * cos ( 1.0 ) ) / r8_pi;

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

double *p51_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P51_FUN evaluates the integrand for problem 51.
//
//  Discussion:
//
//    The function has a removable singularity at x = 0.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    ln ( x ) / ( 1 + ( ln(x) )^2 )^2
//
//  Exact Integral:
//
//    - ( ci(1) * sin(1) + ( pi/2 - si(1) ) * cos(1) ) / pi = - 0.1892752...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Robert Piessens, Elise de Doncker-Kapenga, 
//    Christian Ueberhuber, David Kahaner,
//    QUADPACK: A Subroutine Package for Automatic Integration,
//    Springer, 1983, page 108.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P51_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] == 0.0 )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = log ( x[i] ) / pow ( 1.0 + pow ( log ( x[i] ), 2 ), 2 );
    }
  }
  return fx;
}
/******************************************************************************/

void p51_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P51_LIM returns the integration limits for problem 51.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p52_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P52_EXACT returns the exact integral for problem 52.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P52_EXACT, the value of the integral.
//
{
  double exact;

  exact = log ( 125.0 / 631.0 ) / 18.0;

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

double *p52_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P52_FUN evaluates the integrand for problem 52.
//
//  Discussion:
//
//    The function has a singularity at x = 0.
//
//  Interval:
//
//    -1 <= x <= 5
//
//  Integrand:
//
//    1 / ( x * ( 5 * x^3 + 6 ) )
//
//  Exact Integral:
//
//    ln ( 125 / 631 ) / 18 = -0.08994401...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Robert Piessens, Elise de Doncker-Kapenga, 
//    Christian Ueberhuber, David Kahaner,
//    QUADPACK: A Subroutine Package for Automatic Integration,
//    Springer, 1983, page 109.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P52_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] == 0.0 )
    {
      fx[i] = 0.0;
    }
    else
    {
      fx[i] = 1.0 / ( x[i] * ( 5.0 * pow ( x[i], 3 ) + 6.0 ) );
    }
  }
  return fx;
}
/******************************************************************************/

void p52_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P52_LIM returns the integration limits for problem 52.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = -1.0;
  *b = 5.0;

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

double p53_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P53_EXACT returns the exact integral for problem 53.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    26 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P53_EXACT, the value of the integral.
//
{
  double exact;
  const double r8_pi = 3.141592653589793;

  exact = 0.5 * r8_pi - atan ( 1.0 / sqrt ( 2.0 ) ) + log ( 3.0 ) / 2.0;

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

double *p53_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P53_FUN evaluates the integrand for problem 53.
//
//  Discussion:
//
//    The integrand is singular at x = -1 + sqrt ( 3 ) = 0.732...
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    1 / sqrt ( abs ( x^2 + 2 * x - 2 ) )
//
//  Exact Integral:
//
//    pi / 2 - arctan ( 1 / sqrt ( 2 ) ) + ln ( 3 ) / 2 = 1.504622...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Robert Piessens, Elise de Doncker-Kapenga, 
//    Christian Ueberhuber, David Kahaner,
//    QUADPACK: A Subroutine Package for Automatic Integration,
//    Springer, 1983, page 110.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P53_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = 1.0 / sqrt ( r8_abs ( x[i] * x[i] + 2.0 * x[i] - 2.0 ) );
  }
  return fx;
}
/******************************************************************************/

void p53_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P53_LIM returns the integration limits for problem 53.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p54_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P54_EXACT returns the exact integral for problem 54.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P54_EXACT, the value of the integral.
//
{
  double exact;

  exact = 2.0 / sqrt ( 3.0 );

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

double *p54_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P54_FUN evaluates the integrand for problem 54.
//
//  Discussion:
//
//    The reference claims that this integrand is more closely approximated
//    by the trapezoid rule than by Gauss-Legendre quadrature.
//
//    Points  Trapezoid  Gauss-Legendre
//     4      1.91667    2.53883
//    12      2.1594     2.25809
//
//    However, the stated results hardly give one confindence in
//    the convergence of the trapezoid results, and I am unable to
//    confirm them, because my results for 4 points give good results
//    (about 1.14) for BOTH Trapezoid and Gauss-Legendre//
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    2 / ( 2 + sin ( 10 * PI * X ) )
//
//  Exact Integral:
//
//    2 / sqrt ( 3 ) = 1.1547...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Prem Kythe, Pratap Puri,
//    Computational Methods for Linear Integral Equations,
//    Birkhaeuser, 2002.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P54_FUN[N], the integrand values.
//
{
  double *fx;
  int i;
  const double r8_pi = 3.141592653589793;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = 2.0 / ( 2.0 + sin ( 10.0 * r8_pi * x[i] ) );
  }
  return fx;
}
/******************************************************************************/

void p54_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P54_LIM returns the integration limits for problem 54.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double p55_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P55_EXACT returns the exact integral for problem 55.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 April 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P55_EXACT, the value of the integral.
//
{
  double a;
  double b;
  double c;
  double exact;
  const double r8_pi = 3.141592653589793;
  double x0;

  p55_lim ( &a, &b );
  p55_param ( "get", "c", &c );
  p55_param ( "get", "x0", &x0 );

  exact = sqrt ( r8_pi ) * 
    ( r8_erf ( c * ( b - x0 ) ) - r8_erf ( c * ( a - x0 ) ) ) 
    / ( 2.0 * c );

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

double *p55_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P55_FUN evaluates the integrand for problem 55.
//
//  Interval:
//
//    a = 0 <= x <= 1 = b
//
//  Integrand:
//
//    exp ( - c^2 * ( x - x0 )^2 )
//
//  Exact Integral:
//
//    sqrt ( pi ) 
//    * ( erf ( c * ( b - x0 ) ) - erf ( c * ( a - x0 ) ) )
//    / ( 2 * c )
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    20 April 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P55_FUN[N], the integrand values.
//
{
  double c;
  double *fx;
  int i;
  double x0;

  p55_param ( "get", "c", &c );
  p55_param ( "get", "x0", &x0 );

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = exp ( - c * c * pow ( x[i] - x0, 2 ) );
  }
  return fx;
}
/******************************************************************************/

void p55_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P55_LIM returns the integration limits for problem 55.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

void p55_param ( char *action, char *name, double *value )

/******************************************************************************/
//
//  Purpose:
//
//    P55_PARAM sets or gets real scalar parameters for problem 55.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *ACTION,
//    "get" to get a parameter.
//    "set" to set a parameter.
//
//    Input, char *NAME, the name of the variable.
//    "c" is the coefficient.
//    "x0" is the base point.
//
//    Input/output, double *VALUE.
//    * If ACTION = "set", then VALUE is an input quantity, and is the
//      new value to be assigned to NAME.
//    * If ACTION = "get", then VALUE is an output quantity, and is the
//      current value of NAME.
//
{
  static double c = 3.0;
  static double x0 = 0.75;

  if ( s_eqi ( action, "get" ) )
  {
    if ( s_eqi ( name, "c" ) )
    {
      *value = c;
    }
    else if ( s_eqi ( name, "x0" ) )
    {
      *value = x0;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P55_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else if ( s_eqi ( action, "set" ) )
  {
    if ( s_eqi ( name, "c" ) )
    {
      c = *value;
    }
    else if ( s_eqi ( name, "x0" ) )
    {
      x0 = *value;
    }
    else
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P55_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Unrecognized name.\n" );
      exit ( 1 );
    }
  }
  else
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P55_PARAM - Fatal error!\n" );
    fprintf ( stderr, "  Unrecognized action.\n" );
    exit ( 1 );
  }
  return;
}
/******************************************************************************/

double p56_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P56_EXACT returns the estimated integral for problem 56.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P56_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = 1.9922524079504000171;

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

double *p56_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P56_FUN evaluates the integrand for problem 56.
//
//  Interval:
//
//    -1 <= x <= 1
//
//  Integrand:
//
//    1 / ( x^6 + 0.9 )
//
//  Approximate Integral (20 digits):
//
//    1.9922524079504000171...
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software,
//    edited by John R Rice,
//    Academic Press, 1971.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P56_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = 1.0 / ( pow ( x[i], 6 ) + 0.9 );
  }
  return fx;
}
/******************************************************************************/

void p56_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P56_LIM returns the integration limits for problem 56.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = - 1.0;
  *b = 1.0;

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

double p57_exact ( )

/******************************************************************************/
//
//  Purpose:
//
//    P57_EXACT returns the exact integral for problem 57.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double P57_EXACT, the estimated value of the integral.
//
{
  double exact;

  exact = 0.4;

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

double *p57_fun ( int n, double x[] )

/******************************************************************************/
//
//  Purpose:
//
//    P57_FUN evaluates the integrand for problem 57.
//
//  Interval:
//
//    0 <= x <= 1
//
//  Integrand:
//
//    x^(3/2)
//
//  Antiderivative:
//
//    (2/5) * x^(5/2)
//
//  Exact Integral:
//
//    0.4
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    David Kahaner,
//    Comparison of Numerical Quadrature Formulas,
//    in Mathematical Software, edited by John R Rice,
//    Academic Press, 1971.
//
//  Parameters:
//
//    Input, int N, the number of evaluation points.
//
//    Input, double X[N], the evaluation points.
//
//    Output, double P57_FUN[N], the integrand values.
//
{
  double *fx;
  int i;

  fx = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    fx[i] = sqrt ( pow ( x[i], 3 ) );
  }
  return fx;
}
/******************************************************************************/

void p57_lim ( double *a, double *b )

/******************************************************************************/
//
//  Purpose:
//
//    P57_LIM returns the integration limits for problem 57.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 December 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double *A, *B, the limits of integration.
//
{
  *a = 0.0;
  *b = 1.0;

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

double r8_abs ( double x )

/******************************************************************************/
/*
  Purpose:

    R8_ABS returns the absolute value of an R8.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    07 May 2006

  Author:

    John Burkardt

  Parameters:

    Input, double X, the quantity whose absolute value is desired.

    Output, double R8_ABS, the absolute value of X.
*/
{
  double value;

  if ( 0.0 <= x )
  {
    value = + x;
  }
  else
  {
    value = - x;
  }
  return value;
}
/******************************************************************************/

double r8_aint ( double x )

/******************************************************************************/
/*
  Purpose:

    R8_AINT truncates an R8 argument to an integer.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    1 September 2011

  Author:

    John Burkardt.

  Parameters:

    Input, double X, the argument.

    Output, double R8_AINT, the truncated version of X.
*/
{
  double value;

  if ( x < 0.0E+00 )
  {
    value = - ( double ) ( ( int ) ( r8_abs ( x ) ) );
  }
  else
  {
    value =   ( double ) ( ( int ) ( r8_abs ( x ) ) );
  }

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

void r8_b0mp ( double x, double *ampl, double *theta )

/******************************************************************************/
/*
  Purpose:

    R8_B0MP evaluates the modulus and phase for the Bessel J0 and Y0 functions.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 January 2012

  Author:

    Original FORTRAN77 version by Wayne Fullerton.
    C version by John Burkardt.

  Reference:

    Wayne Fullerton,
    Portable Special Function Routines,
    in Portability of Numerical Software,
    edited by Wayne Cowell,
    Lecture Notes in Computer Science, Volume 57,
    Springer 1977,
    ISBN: 978-3-540-08446-4,
    LC: QA297.W65.

  Parameters:

    Input, double X, the argument.

    Output, double *AMPL, *THETA, the modulus and phase.
*/
{
  static double bm0cs[37] = {
    +0.9211656246827742712573767730182E-01,
    -0.1050590997271905102480716371755E-02,
    +0.1470159840768759754056392850952E-04,
    -0.5058557606038554223347929327702E-06,
    +0.2787254538632444176630356137881E-07,
    -0.2062363611780914802618841018973E-08,
    +0.1870214313138879675138172596261E-09,
    -0.1969330971135636200241730777825E-10,
    +0.2325973793999275444012508818052E-11,
    -0.3009520344938250272851224734482E-12,
    +0.4194521333850669181471206768646E-13,
    -0.6219449312188445825973267429564E-14,
    +0.9718260411336068469601765885269E-15,
    -0.1588478585701075207366635966937E-15,
    +0.2700072193671308890086217324458E-16,
    -0.4750092365234008992477504786773E-17,
    +0.8615128162604370873191703746560E-18,
    -0.1605608686956144815745602703359E-18,
    +0.3066513987314482975188539801599E-19,
    -0.5987764223193956430696505617066E-20,
    +0.1192971253748248306489069841066E-20,
    -0.2420969142044805489484682581333E-21,
    +0.4996751760510616453371002879999E-22,
    -0.1047493639351158510095040511999E-22,
    +0.2227786843797468101048183466666E-23,
    -0.4801813239398162862370542933333E-24,
    +0.1047962723470959956476996266666E-24,
    -0.2313858165678615325101260800000E-25,
    +0.5164823088462674211635199999999E-26,
    -0.1164691191850065389525401599999E-26,
    +0.2651788486043319282958336000000E-27,
    -0.6092559503825728497691306666666E-28,
    +0.1411804686144259308038826666666E-28,
    -0.3298094961231737245750613333333E-29,
    +0.7763931143074065031714133333333E-30,
    -0.1841031343661458478421333333333E-30,
    +0.4395880138594310737100799999999E-31 };
  static double bm02cs[40] = {
    +0.9500415145228381369330861335560E-01,
    -0.3801864682365670991748081566851E-03,
    +0.2258339301031481192951829927224E-05,
    -0.3895725802372228764730621412605E-07,
    +0.1246886416512081697930990529725E-08,
    -0.6065949022102503779803835058387E-10,
    +0.4008461651421746991015275971045E-11,
    -0.3350998183398094218467298794574E-12,
    +0.3377119716517417367063264341996E-13,
    -0.3964585901635012700569356295823E-14,
    +0.5286111503883857217387939744735E-15,
    -0.7852519083450852313654640243493E-16,
    +0.1280300573386682201011634073449E-16,
    -0.2263996296391429776287099244884E-17,
    +0.4300496929656790388646410290477E-18,
    -0.8705749805132587079747535451455E-19,
    +0.1865862713962095141181442772050E-19,
    -0.4210482486093065457345086972301E-20,
    +0.9956676964228400991581627417842E-21,
    -0.2457357442805313359605921478547E-21,
    +0.6307692160762031568087353707059E-22,
    -0.1678773691440740142693331172388E-22,
    +0.4620259064673904433770878136087E-23,
    -0.1311782266860308732237693402496E-23,
    +0.3834087564116302827747922440276E-24,
    -0.1151459324077741271072613293576E-24,
    +0.3547210007523338523076971345213E-25,
    -0.1119218385815004646264355942176E-25,
    +0.3611879427629837831698404994257E-26,
    -0.1190687765913333150092641762463E-26,
    +0.4005094059403968131802476449536E-27,
    -0.1373169422452212390595193916017E-27,
    +0.4794199088742531585996491526437E-28,
    -0.1702965627624109584006994476452E-28,
    +0.6149512428936330071503575161324E-29,
    -0.2255766896581828349944300237242E-29,
    +0.8399707509294299486061658353200E-30,
    -0.3172997595562602355567423936152E-30,
    +0.1215205298881298554583333026514E-30,
    -0.4715852749754438693013210568045E-31 };
  static double bt02cs[39] = {
    -0.24548295213424597462050467249324,
    +0.12544121039084615780785331778299E-02,
    -0.31253950414871522854973446709571E-04,
    +0.14709778249940831164453426969314E-05,
    -0.99543488937950033643468850351158E-07,
    +0.85493166733203041247578711397751E-08,
    -0.86989759526554334557985512179192E-09,
    +0.10052099533559791084540101082153E-09,
    -0.12828230601708892903483623685544E-10,
    +0.17731700781805131705655750451023E-11,
    -0.26174574569485577488636284180925E-12,
    +0.40828351389972059621966481221103E-13,
    -0.66751668239742720054606749554261E-14,
    +0.11365761393071629448392469549951E-14,
    -0.20051189620647160250559266412117E-15,
    +0.36497978794766269635720591464106E-16,
    -0.68309637564582303169355843788800E-17,
    +0.13107583145670756620057104267946E-17,
    -0.25723363101850607778757130649599E-18,
    +0.51521657441863959925267780949333E-19,
    -0.10513017563758802637940741461333E-19,
    +0.21820381991194813847301084501333E-20,
    -0.46004701210362160577225905493333E-21,
    +0.98407006925466818520953651199999E-22,
    -0.21334038035728375844735986346666E-22,
    +0.46831036423973365296066286933333E-23,
    -0.10400213691985747236513382399999E-23,
    +0.23349105677301510051777740800000E-24,
    -0.52956825323318615788049749333333E-25,
    +0.12126341952959756829196287999999E-25,
    -0.28018897082289428760275626666666E-26,
    +0.65292678987012873342593706666666E-27,
    -0.15337980061873346427835733333333E-27,
    +0.36305884306364536682359466666666E-28,
    -0.86560755713629122479172266666666E-29,
    +0.20779909972536284571238399999999E-29,
    -0.50211170221417221674325333333333E-30,
    +0.12208360279441714184191999999999E-30,
    -0.29860056267039913454250666666666E-31 };
  static double bth0cs[44] = {
    -0.24901780862128936717709793789967,
    +0.48550299609623749241048615535485E-03,
    -0.54511837345017204950656273563505E-05,
    +0.13558673059405964054377445929903E-06,
    -0.55691398902227626227583218414920E-08,
    +0.32609031824994335304004205719468E-09,
    -0.24918807862461341125237903877993E-10,
    +0.23449377420882520554352413564891E-11,
    -0.26096534444310387762177574766136E-12,
    +0.33353140420097395105869955014923E-13,
    -0.47890000440572684646750770557409E-14,
    +0.75956178436192215972642568545248E-15,
    -0.13131556016891440382773397487633E-15,
    +0.24483618345240857495426820738355E-16,
    -0.48805729810618777683256761918331E-17,
    +0.10327285029786316149223756361204E-17,
    -0.23057633815057217157004744527025E-18,
    +0.54044443001892693993017108483765E-19,
    -0.13240695194366572724155032882385E-19,
    +0.33780795621371970203424792124722E-20,
    -0.89457629157111779003026926292299E-21,
    +0.24519906889219317090899908651405E-21,
    -0.69388422876866318680139933157657E-22,
    +0.20228278714890138392946303337791E-22,
    -0.60628500002335483105794195371764E-23,
    +0.18649748964037635381823788396270E-23,
    -0.58783732384849894560245036530867E-24,
    +0.18958591447999563485531179503513E-24,
    -0.62481979372258858959291620728565E-25,
    +0.21017901684551024686638633529074E-25,
    -0.72084300935209253690813933992446E-26,
    +0.25181363892474240867156405976746E-26,
    -0.89518042258785778806143945953643E-27,
    +0.32357237479762298533256235868587E-27,
    -0.11883010519855353657047144113796E-27,
    +0.44306286907358104820579231941731E-28,
    -0.16761009648834829495792010135681E-28,
    +0.64292946921207466972532393966088E-29,
    -0.24992261166978652421207213682763E-29,
    +0.98399794299521955672828260355318E-30,
    -0.39220375242408016397989131626158E-30,
    +0.15818107030056522138590618845692E-30,
    -0.64525506144890715944344098365426E-31,
    +0.26611111369199356137177018346367E-31 };
  double eta;
  static int nbm0 = 0;
  static int nbm02 = 0;
  static int nbt02 = 0;
  static int nbth0 = 0;
  const double pi4 = 0.785398163397448309615660845819876;
  double z;

  if ( nbm0 == 0 )
  {
    eta = 0.1 * r8_mach ( 3 );
    nbm0 = r8_inits ( bm0cs, 37, eta );
    nbt02 = r8_inits ( bt02cs, 39, eta );
    nbm02 = r8_inits ( bm02cs, 40, eta );
    nbth0 = r8_inits ( bth0cs, 44, eta );
  }

  if ( x < 4.0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8_B0MP - Fatal error!\n" );
    fprintf ( stderr, "  X < 4.\n" );
    exit ( 1 );
  }
  else if ( x <= 8.0 )
  {
    z = ( 128.0 / x / x - 5.0 ) / 3.0;
    *ampl = ( 0.75 + r8_csevl ( z, bm0cs, nbm0 ) ) / sqrt ( x );
    *theta = x - pi4 + r8_csevl ( z, bt02cs, nbt02 ) / x;
  }
  else
  {
    z = 128.0 / x / x - 1.0;
    *ampl = ( 0.75 + r8_csevl ( z, bm02cs, nbm02) ) / sqrt ( x );
    *theta = x - pi4 + r8_csevl ( z, bth0cs, nbth0 ) / x;
  }
  return;
}
/******************************************************************************/

double r8_besj0 ( double x )

/******************************************************************************/
/*
  Purpose:

    R8_BESJ0 evaluates the Bessel function J of order 0 of an R8 argument.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 January 2012

  Author:

    Original FORTRAN77 version by Wayne Fullerton.
    C version by John Burkardt.

  Reference:

    Wayne Fullerton,
    Portable Special Function Routines,
    in Portability of Numerical Software,
    edited by Wayne Cowell,
    Lecture Notes in Computer Science, Volume 57,
    Springer 1977,
    ISBN: 978-3-540-08446-4,
    LC: QA297.W65.

  Parameters:

    Input, double X, the argument.

    Output, double R8_BESJ0, the Bessel function J of order 0 of X.
*/
{
  double ampl;
  static double bj0cs[19] = {
    +0.10025416196893913701073127264074,
    -0.66522300776440513177678757831124,
    +0.24898370349828131370460468726680,
    -0.33252723170035769653884341503854E-01,
    +0.23114179304694015462904924117729E-02,
    -0.99112774199508092339048519336549E-04,
    +0.28916708643998808884733903747078E-05,
    -0.61210858663032635057818407481516E-07,
    +0.98386507938567841324768748636415E-09,
    -0.12423551597301765145515897006836E-10,
    +0.12654336302559045797915827210363E-12,
    -0.10619456495287244546914817512959E-14,
    +0.74706210758024567437098915584000E-17,
    -0.44697032274412780547627007999999E-19,
    +0.23024281584337436200523093333333E-21,
    -0.10319144794166698148522666666666E-23,
    +0.40608178274873322700800000000000E-26,
    -0.14143836005240913919999999999999E-28,
    +0.43910905496698880000000000000000E-31 };
  static int ntj0 = 0;
  double theta;
  double value;
  static double xsml = 0.0;
  double y;

  if ( ntj0 == 0 )
  {
    ntj0 = r8_inits ( bj0cs, 19, 0.1 * r8_mach ( 3 ) );
    xsml = sqrt ( 4.0 * r8_mach ( 3 ) );
  }

  y = r8_abs ( x );

  if ( y <= xsml )
  {
    value = 1.0;
  }
  else if ( y <= 4.0 )
  {
    value = r8_csevl ( 0.125 * y * y - 1.0, bj0cs, ntj0 );
  }
  else
  {
    r8_b0mp ( y, &ampl, &theta );
    value = ampl * cos ( theta );
  }
  return value;
}
/******************************************************************************/

double r8_ci ( double x )

/******************************************************************************/
/*
  Purpose:

    R8_CI evaluates the cosine integral Ci of an R8 argument.

  Discussion:

    The cosine integral is defined by

      CI(X) = - integral ( X <= T < Infinity ) ( cos ( T ) ) / T  dT

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 January 2012

  Author:

    Original FORTRAN77 version by Wayne Fullerton.
    C version by John Burkardt.

  Reference:

    Wayne Fullerton,
    Portable Special Function Routines,
    in Portability of Numerical Software,
    edited by Wayne Cowell,
    Lecture Notes in Computer Science, Volume 57,
    Springer 1977,
    ISBN: 978-3-540-08446-4,
    LC: QA297.W65.

  Parameters:

    Input, double X, the argument.

    Output, double R8_CI, the cosine integral Ci evaluated at X.
*/
{
  static double cics[19] = {
    -0.34004281856055363156281076633129873,
    -1.03302166401177456807159271040163751,
     0.19388222659917082876715874606081709,
    -0.01918260436019865893946346270175301,
     0.00110789252584784967184098099266118,
    -0.00004157234558247208803840231814601,
     0.00000109278524300228715295578966285,
    -0.00000002123285954183465219601280329,
     0.00000000031733482164348544865129873,
    -0.00000000000376141547987683699381798,
     0.00000000000003622653488483964336956,
    -0.00000000000000028911528493651852433,
     0.00000000000000000194327860676494420,
    -0.00000000000000000001115183182650184,
     0.00000000000000000000005527858887706,
    -0.00000000000000000000000023907013943,
     0.00000000000000000000000000091001612,
    -0.00000000000000000000000000000307233,
     0.00000000000000000000000000000000926 };
  double f;
  double g;
  static int nci = 0;
  double sinx;
  double value;
  static double xsml = 0.0;
  double y;

  if ( nci == 0 )
  {
    nci = r8_inits ( cics, 19, 0.1 * r8_mach ( 3 ) );
    xsml = sqrt ( r8_mach ( 3 ) );
  }

  if ( x <= 0.0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8_CI - Fatal error!\n" );
    fprintf ( stderr, "  X <= 0.0.\n" );
    exit ( 1 );
  }
  else if ( x <= xsml )
  {
    y = - 1.0;
    value = log ( x ) - 0.5 + r8_csevl ( y, cics, nci );
  }
  else if ( x <= 4.0 )
  {
    y = ( x * x - 8.0 ) * 0.125;
    value = log ( x ) - 0.5 + r8_csevl ( y, cics, nci );
  }
  else
  {
    r8_sifg ( x, &f, &g );
    sinx = sin ( x );
    value = f * sinx - g * cos ( x );
  }
  return value;
}
/******************************************************************************/

 double r8_csevl ( double x, double a[], int n )

/******************************************************************************/
/*
  Purpose:

    R8_CSEVL evaluates a Chebyshev series.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 January 2012

  Author:

    C version by John Burkardt.

  Reference:

    Roger Broucke,
    Algorithm 446:
    Ten Subroutines for the Manipulation of Chebyshev Series,
    Communications of the ACM,
    Volume 16, Number 4, April 1973, pages 254-256.

  Parameters:

    Input, double X, the evaluation point.

    Input, double CS[N], the Chebyshev coefficients.

    Input, int N, the number of Chebyshev coefficients.

    Output, double R8_CSEVL, the Chebyshev series evaluated at X.
*/
{
  double b0;
  double b1;
  double b2;
  int i;
  double twox;
  double value;

  if ( n < 1 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8_CSEVL - Fatal error!\n" );
    fprintf ( stderr, "  Number of terms <= 0.\n" );
    exit ( 1 );
  }

  if ( 1000 < n )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8_CSEVL - Fatal error!\n" );
    fprintf ( stderr, "  Number of terms greater than 1000.\n" );
    exit ( 1 );
 }

  if ( x < -1.1 || 1.1 < x )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8_CSEVL - Fatal error!\n" );
    fprintf ( stderr, "  X outside (-1,+1).\n" );
    exit ( 1 );
  }

  twox = 2.0 * x;
  b1 = 0.0;
  b0 = 0.0;

  for ( i = n - 1; 0 <= i; i-- )
  {
    b2 = b1;
    b1 = b0;
    b0 = twox * b1 - b2 + a[i];
  }

  value = 0.5 * ( b0 - b2 );

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

double r8_epsilon ( )

/******************************************************************************/
/*
  Purpose:

    R8_EPSILON returns the R8 round off unit.

  Discussion:

    R8_EPSILON is a number R which is a power of 2 with the property that,
    to the precision of the computer's arithmetic,
      1 < 1 + R
    but
      1 = ( 1 + R / 2 )

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    01 September 2012

  Author:

    John Burkardt

  Parameters:

    Output, double R8_EPSILON, the R8 round-off unit.
*/
{
  const double value = 2.220446049250313E-016;

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

double r8_erf ( double x )

/******************************************************************************/
/*
  Purpose:

    R8_ERF evaluates the error function of an R8 argument.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 January 2012

  Author:

    Original FORTRAN77 version by Wayne Fullerton.
    C version by John Burkardt.

  Reference:

    Wayne Fullerton,
    Portable Special Function Routines,
    in Portability of Numerical Software,
    edited by Wayne Cowell,
    Lecture Notes in Computer Science, Volume 57,
    Springer 1977,
    ISBN: 978-3-540-08446-4,
    LC: QA297.W65.

  Parameters:

    Input, double X, the argument.

    Output, double R8_ERF, the error function of X.
*/
{
  static double erfcs[21] = {
    -0.49046121234691808039984544033376E-01,
    -0.14226120510371364237824741899631,
    +0.10035582187599795575754676712933E-01,
    -0.57687646997674847650827025509167E-03,
    +0.27419931252196061034422160791471E-04,
    -0.11043175507344507604135381295905E-05,
    +0.38488755420345036949961311498174E-07,
    -0.11808582533875466969631751801581E-08,
    +0.32334215826050909646402930953354E-10,
    -0.79910159470045487581607374708595E-12,
    +0.17990725113961455611967245486634E-13,
    -0.37186354878186926382316828209493E-15,
    +0.71035990037142529711689908394666E-17,
    -0.12612455119155225832495424853333E-18,
    +0.20916406941769294369170500266666E-20,
    -0.32539731029314072982364160000000E-22,
    +0.47668672097976748332373333333333E-24,
    -0.65980120782851343155199999999999E-26,
    +0.86550114699637626197333333333333E-28,
    -0.10788925177498064213333333333333E-29,
    +0.12811883993017002666666666666666E-31 };
  static int nterf = 0;
  static double sqeps = 0.0;
  static double sqrtpi = 1.77245385090551602729816748334115;
  double value;
  static double xbig = 0.0;
  double y;

  if ( nterf == 0 )
  {
    nterf = r8_inits ( erfcs, 21, 0.1 * r8_mach ( 3 ) );
    xbig = sqrt ( - log ( sqrtpi * r8_mach ( 3 ) ) );
    sqeps = sqrt ( 2.0 * r8_mach ( 3 ) );
  }

  y = r8_abs ( x );

  if ( y <= sqeps )
  {
    value = 2.0 * x / sqrtpi;
  }
  else if ( y <= 1.0 )
  {
    value = x * ( 1.0 + r8_csevl ( 2.0 * x * x - 1.0, erfcs, nterf ) );
  }
  else if ( y <= xbig )
  {
    value = 1.0 - r8_erfc ( y );
    if ( x < 0.0 )
    {
      value = - value;
    }
  }
  else
  {
    value = 1.0;
    if ( x < 0.0 )
    {
      value = - value;
    }
  }
  return value;
}
/******************************************************************************/

double r8_erfc ( double x )

/******************************************************************************/
/*
  Purpose:

    R8_ERFC evaluates the co-error function of an R8 argument.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 January 2012

  Author:

    Original FORTRAN77 version by Wayne Fullerton.
    C version by John Burkardt.

  Reference:

    Wayne Fullerton,
    Portable Special Function Routines,
    in Portability of Numerical Software,
    edited by Wayne Cowell,
    Lecture Notes in Computer Science, Volume 57,
    Springer 1977,
    ISBN: 978-3-540-08446-4,
    LC: QA297.W65.

  Parameters:

    Input, double X, the argument.

    Output, double R8_ERFC, the co-error function of X.
*/
{
  static double erc2cs[49] = {
    -0.6960134660230950112739150826197E-01,
    -0.4110133936262089348982212084666E-01,
    +0.3914495866689626881561143705244E-02,
    -0.4906395650548979161280935450774E-03,
    +0.7157479001377036380760894141825E-04,
    -0.1153071634131232833808232847912E-04,
    +0.1994670590201997635052314867709E-05,
    -0.3642666471599222873936118430711E-06,
    +0.6944372610005012589931277214633E-07,
    -0.1371220902104366019534605141210E-07,
    +0.2788389661007137131963860348087E-08,
    -0.5814164724331161551864791050316E-09,
    +0.1238920491752753181180168817950E-09,
    -0.2690639145306743432390424937889E-10,
    +0.5942614350847910982444709683840E-11,
    -0.1332386735758119579287754420570E-11,
    +0.3028046806177132017173697243304E-12,
    -0.6966648814941032588795867588954E-13,
    +0.1620854541053922969812893227628E-13,
    -0.3809934465250491999876913057729E-14,
    +0.9040487815978831149368971012975E-15,
    -0.2164006195089607347809812047003E-15,
    +0.5222102233995854984607980244172E-16,
    -0.1269729602364555336372415527780E-16,
    +0.3109145504276197583836227412951E-17,
    -0.7663762920320385524009566714811E-18,
    +0.1900819251362745202536929733290E-18,
    -0.4742207279069039545225655999965E-19,
    +0.1189649200076528382880683078451E-19,
    -0.3000035590325780256845271313066E-20,
    +0.7602993453043246173019385277098E-21,
    -0.1935909447606872881569811049130E-21,
    +0.4951399124773337881000042386773E-22,
    -0.1271807481336371879608621989888E-22,
    +0.3280049600469513043315841652053E-23,
    -0.8492320176822896568924792422399E-24,
    +0.2206917892807560223519879987199E-24,
    -0.5755617245696528498312819507199E-25,
    +0.1506191533639234250354144051199E-25,
    -0.3954502959018796953104285695999E-26,
    +0.1041529704151500979984645051733E-26,
    -0.2751487795278765079450178901333E-27,
    +0.7290058205497557408997703680000E-28,
    -0.1936939645915947804077501098666E-28,
    +0.5160357112051487298370054826666E-29,
    -0.1378419322193094099389644800000E-29,
    +0.3691326793107069042251093333333E-30,
    -0.9909389590624365420653226666666E-31,
    +0.2666491705195388413323946666666E-31 };
  static double erfccs[59] = {
    +0.715179310202924774503697709496E-01,
    -0.265324343376067157558893386681E-01,
    +0.171115397792085588332699194606E-02,
    -0.163751663458517884163746404749E-03,
    +0.198712935005520364995974806758E-04,
    -0.284371241276655508750175183152E-05,
    +0.460616130896313036969379968464E-06,
    -0.822775302587920842057766536366E-07,
    +0.159214187277090112989358340826E-07,
    -0.329507136225284321486631665072E-08,
    +0.722343976040055546581261153890E-09,
    -0.166485581339872959344695966886E-09,
    +0.401039258823766482077671768814E-10,
    -0.100481621442573113272170176283E-10,
    +0.260827591330033380859341009439E-11,
    -0.699111056040402486557697812476E-12,
    +0.192949233326170708624205749803E-12,
    -0.547013118875433106490125085271E-13,
    +0.158966330976269744839084032762E-13,
    -0.472689398019755483920369584290E-14,
    +0.143587337678498478672873997840E-14,
    -0.444951056181735839417250062829E-15,
    +0.140481088476823343737305537466E-15,
    -0.451381838776421089625963281623E-16,
    +0.147452154104513307787018713262E-16,
    -0.489262140694577615436841552532E-17,
    +0.164761214141064673895301522827E-17,
    -0.562681717632940809299928521323E-18,
    +0.194744338223207851429197867821E-18,
    -0.682630564294842072956664144723E-19,
    +0.242198888729864924018301125438E-19,
    -0.869341413350307042563800861857E-20,
    +0.315518034622808557122363401262E-20,
    -0.115737232404960874261239486742E-20,
    +0.428894716160565394623737097442E-21,
    -0.160503074205761685005737770964E-21,
    +0.606329875745380264495069923027E-22,
    -0.231140425169795849098840801367E-22,
    +0.888877854066188552554702955697E-23,
    -0.344726057665137652230718495566E-23,
    +0.134786546020696506827582774181E-23,
    -0.531179407112502173645873201807E-24,
    +0.210934105861978316828954734537E-24,
    -0.843836558792378911598133256738E-25,
    +0.339998252494520890627359576337E-25,
    -0.137945238807324209002238377110E-25,
    +0.563449031183325261513392634811E-26,
    -0.231649043447706544823427752700E-26,
    +0.958446284460181015263158381226E-27,
    -0.399072288033010972624224850193E-27,
    +0.167212922594447736017228709669E-27,
    -0.704599152276601385638803782587E-28,
    +0.297976840286420635412357989444E-28,
    -0.126252246646061929722422632994E-28,
    +0.539543870454248793985299653154E-29,
    -0.238099288253145918675346190062E-29,
    +0.109905283010276157359726683750E-29,
    -0.486771374164496572732518677435E-30,
    +0.152587726411035756763200828211E-30 };
  static double erfcs[21] = {
    -0.49046121234691808039984544033376E-01,
    -0.14226120510371364237824741899631,
    +0.10035582187599795575754676712933E-01,
    -0.57687646997674847650827025509167E-03,
    +0.27419931252196061034422160791471E-04,
    -0.11043175507344507604135381295905E-05,
    +0.38488755420345036949961311498174E-07,
    -0.11808582533875466969631751801581E-08,
    +0.32334215826050909646402930953354E-10,
    -0.79910159470045487581607374708595E-12,
    +0.17990725113961455611967245486634E-13,
    -0.37186354878186926382316828209493E-15,
    +0.71035990037142529711689908394666E-17,
    -0.12612455119155225832495424853333E-18,
    +0.20916406941769294369170500266666E-20,
    -0.32539731029314072982364160000000E-22,
    +0.47668672097976748332373333333333E-24,
    -0.65980120782851343155199999999999E-26,
    +0.86550114699637626197333333333333E-28,
    -0.10788925177498064213333333333333E-29,
    +0.12811883993017002666666666666666E-31 };
  double eta;
  static int nterc2 = 0;
  static int nterf = 0;
  static int nterfc = 0;
  static double sqeps = 0.0;
  static double sqrtpi = 1.77245385090551602729816748334115;
  double value;
  static double xmax = 0.0;
  static double xsml = 0.0;
  double y;

  if ( nterf == 0 )
  {
    eta = 0.1 * r8_mach ( 3 );
    nterf = r8_inits ( erfcs, 21, eta );
    nterfc = r8_inits ( erfccs, 59, eta );
    nterc2 = r8_inits ( erc2cs, 49, eta );

    xsml = - sqrt ( - log ( sqrtpi * r8_mach ( 3 ) ) );
    xmax = sqrt (- log ( sqrtpi * r8_mach ( 1 ) ) );
    xmax = xmax - 0.5 * log ( xmax ) / xmax - 0.01;
    sqeps = sqrt ( 2.0 * r8_mach ( 3 ) );
  }

  if ( x <= xsml )
  {
    value = 2.0;
    return value;
  }

  if ( xmax < x )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8_ERFC - Warning!\n" );
    fprintf ( stderr, "  X so big that ERFC underflows.\n" );
    value = 0.0;
    return value;
  }

  y = r8_abs ( x );

  if ( y < sqeps )
  {
    value = 1.0 - 2.0 * x / sqrtpi;
    return value;
  }
  else if ( y <= 1.0 )
  {
    value = 1.0 - x * ( 1.0 
      + r8_csevl ( 2.0 * x * x - 1.0, erfcs, nterf ) );
    return value;
  }

  y = y * y;

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

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

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

void r8_gaml ( double *xmin, double *xmax )

/******************************************************************************/
/*
  Purpose:

    R8_GAML evaluates bounds for an R8 argument of the gamma function.

  Discussion:

    This function calculates the minimum and maximum legal bounds 
    for X in the evaluation of GAMMA ( X ).

    XMIN and XMAX are not the only bounds, but they are the only 
    non-trivial ones to calculate.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 January 2012

  Author:

    Original FORTRAN77 version by Wayne Fullerton.
    C version by John Burkardt.

  Reference:

    Wayne Fullerton,
    Portable Special Function Routines,
    in Portability of Numerical Software,
    edited by Wayne Cowell,
    Lecture Notes in Computer Science, Volume 57,
    Springer 1977,
    ISBN: 978-3-540-08446-4,
    LC: QA297.W65.

  Parameters:

    Output, double *XMIN, *XMAX, the bounds.
*/
{
  double alnbig;
  double alnsml;
  int i;
  int j;
  double xln;
  double xold;

  alnsml = log ( r8_mach ( 1 ) );
  *xmin = - alnsml;

  for ( i = 1; i <= 10; i++ )
  {
    xold = *xmin;
    xln = log ( *xmin );
    *xmin = *xmin - *xmin * ( ( *xmin + 0.5 ) * xln - *xmin 
      - 0.2258 + alnsml ) / ( *xmin * xln + 0.5 );

    if ( r8_abs ( *xmin - xold ) < 0.005 )
    {
      *xmin = - *xmin + 0.01;

      alnbig = log ( r8_mach ( 2 ) );
      *xmax = alnbig;

      for ( j = 1; j <= 10; j++ )
      {
        xold = *xmax;
        xln = log ( *xmax );
        *xmax = *xmax - *xmax * ( ( *xmax - 0.5 ) * xln - *xmax 
          + 0.9189 - alnbig ) / ( *xmax * xln - 0.5 );

        if ( r8_abs ( *xmax - xold ) < 0.005 )
        {
          *xmax = *xmax - 0.01;
          *xmin = r8_max ( *xmin, - *xmax + 1.0 );
          return;
        }
      }
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "R8_GAML - Fatal error!\n" );
      fprintf ( stderr, "  Unable to find XMAX.\n" );
      exit ( 1 );
    }
  }
  fprintf ( stderr, "\n" );
  fprintf ( stderr, "R8_GAML - Fatal error!\n" );
  fprintf ( stderr, "  Unable to find XMIN.\n" );
  exit ( 1 );
}
/******************************************************************************/

double r8_gamma ( double x )

/******************************************************************************/
/*
  Purpose:

    R8_GAMMA evaluates the gamma function of an R8 argument.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 January 2012

  Author:

    Original FORTRAN77 version by Wayne Fullerton.
    C version by John Burkardt.

  Reference:

    Wayne Fullerton,
    Portable Special Function Routines,
    in Portability of Numerical Software,
    edited by Wayne Cowell,
    Lecture Notes in Computer Science, Volume 57,
    Springer 1977,
    ISBN: 978-3-540-08446-4,
    LC: QA297.W65.

  Parameters:

    Input, double X, the argument.

    Output, double R8_GAMMA, the gamma function of X.
*/
{
  static double dxrel = 0.0;
  static double gcs[42] = {
    +0.8571195590989331421920062399942E-02,
    +0.4415381324841006757191315771652E-02,
    +0.5685043681599363378632664588789E-01,
    -0.4219835396418560501012500186624E-02,
    +0.1326808181212460220584006796352E-02,
    -0.1893024529798880432523947023886E-03,
    +0.3606925327441245256578082217225E-04,
    -0.6056761904460864218485548290365E-05,
    +0.1055829546302283344731823509093E-05,
    -0.1811967365542384048291855891166E-06,
    +0.3117724964715322277790254593169E-07,
    -0.5354219639019687140874081024347E-08,
    +0.9193275519859588946887786825940E-09,
    -0.1577941280288339761767423273953E-09,
    +0.2707980622934954543266540433089E-10,
    -0.4646818653825730144081661058933E-11,
    +0.7973350192007419656460767175359E-12,
    -0.1368078209830916025799499172309E-12,
    +0.2347319486563800657233471771688E-13,
    -0.4027432614949066932766570534699E-14,
    +0.6910051747372100912138336975257E-15,
    -0.1185584500221992907052387126192E-15,
    +0.2034148542496373955201026051932E-16,
    -0.3490054341717405849274012949108E-17,
    +0.5987993856485305567135051066026E-18,
    -0.1027378057872228074490069778431E-18,
    +0.1762702816060529824942759660748E-19,
    -0.3024320653735306260958772112042E-20,
    +0.5188914660218397839717833550506E-21,
    -0.8902770842456576692449251601066E-22,
    +0.1527474068493342602274596891306E-22,
    -0.2620731256187362900257328332799E-23,
    +0.4496464047830538670331046570666E-24,
    -0.7714712731336877911703901525333E-25,
    +0.1323635453126044036486572714666E-25,
    -0.2270999412942928816702313813333E-26,
    +0.3896418998003991449320816639999E-27,
    -0.6685198115125953327792127999999E-28,
    +0.1146998663140024384347613866666E-28,
    -0.1967938586345134677295103999999E-29,
    +0.3376448816585338090334890666666E-30,
    -0.5793070335782135784625493333333E-31 };
  int i;
  int n;
  static int ngcs = 0;
  const double r8_pi = 3.14159265358979323846264338327950;
  double sinpiy;
  static double sq2pil = 0.91893853320467274178032973640562;
  double value;
  static double xmax = 0.0;
  static double xmin = 0.0;
  static double xsml = 0.0;
  double y;

  if ( ngcs == 0 )
  {
    ngcs = r8_inits ( gcs, 42, 0.1 * r8_mach ( 3 ) );
    r8_gaml ( &xmin, &xmax );
    xsml = exp ( r8_max ( log ( r8_mach ( 1 ) ),
      - log ( r8_mach ( 2 ) ) ) + 0.01 );
    dxrel = sqrt ( r8_mach ( 4 ) );
  }

  y = r8_abs ( x );

  if ( y <= 10.0 )
  {
    n = ( int ) ( x );
    if ( x < 0.0 )
    {
      n = n - 1;
    }
    y = x - ( double ) ( n );
    n = n - 1;
    value = 0.9375 + r8_csevl ( 2.0 * y - 1.0, gcs, ngcs );

    if ( n == 0 )
    {
      return value;
    }
    else if ( n < 0 )
    {
      n = - n;

      if ( x == 0.0 )
      {
        fprintf ( stderr, "\n" );
        fprintf ( stderr, "R8_GAMMA - Fatal error!\n" );
        fprintf ( stderr, "  X is 0.\n" );
        exit ( 1 );
      }

      if ( x < 0.0 && x + ( double ) ( n - 2 ) == 0.0 )
      {
        fprintf ( stderr, "\n" );
        fprintf ( stderr, "R8_GAMMA - Fatal error!\n" );
        fprintf ( stderr, "  X is a negative int.\n" );
        exit ( 1 );
      }

      if ( x < - 0.5 && r8_abs ( ( x - r8_aint ( x - 0.5 ) ) / x ) < dxrel )
      {
        fprintf ( stderr, "\n" );
        fprintf ( stderr, "R8_GAMMA - Warning!\n" );
        fprintf ( stderr, "  X too near a negative int,\n" );
        fprintf ( stderr, "  answer is half precision.\n" );
      }

      if ( y < xsml )
      {
        fprintf ( stderr, "\n" );
        fprintf ( stderr, "R8_GAMMA - Fatal error!\n" );
        fprintf ( stderr, "  X is so close to zero that Gamma overflows.\n" );
        exit ( 1 );
      }

      for ( i = 1; i <= n; i++ )
      {
        value = value / ( x + ( double ) ( i - 1 ) );
      }

    }
    else if ( n == 0 )
    {
    }
    else
    {
      for ( i = 1; i <= n; i++ )
      {
        value = ( y + ( double ) ( i ) ) * value;
      }
    }
  }
  else
  {
    if ( xmax < x )
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "R8_GAMMA - Fatal error!\n" );
      fprintf ( stderr, "  X so big that Gamma overflows.\n" );
      exit ( 1 );
    }
/*
  Underflow.
*/
    if ( x < xmin )
    {
      value = 0.0;
      return value;
    }

    value = exp ( ( y - 0.5 ) * log ( y ) - y + sq2pil + r8_lgmc ( y ) );

    if ( 0.0 < x )
    {
      return value;
    }

    if ( r8_abs ( ( x - r8_aint ( x - 0.5 ) ) / x ) < dxrel )
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "R8_GAMMA - Warning!\n" );
      fprintf ( stderr, "  X too near a negative int,\n" );
      fprintf ( stderr, "  answer is half precision.\n" );
    }

    sinpiy = sin ( r8_pi * y );

    if ( sinpiy == 0.0 )
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "R8_GAMMA - Fatal error!\n" );
      fprintf ( stderr, "  X is a negative int.\n" );
      exit ( 1 );
    }
    value = - r8_pi / ( y * sinpiy * value );
  }
  return value;
}
/******************************************************************************/

double r8_huge ( )

/******************************************************************************/
/*
  Purpose:

    R8_HUGE returns a "huge" R8.

  Discussion:

    The value returned by this function is NOT required to be the
    maximum representable R8.  This value varies from machine to machine,
    from compiler to compiler, and may cause problems when being printed.
    We simply want a "very large" but non-infinite number.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    06 October 2007

  Author:

    John Burkardt

  Parameters:

    Output, double R8_HUGE, a "huge" R8 value.
*/
{
  double value;

  value = 1.0E+30;

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

int r8_inits ( double dos[], int nos, double eta )

/******************************************************************************/
/*
  Purpose:

    R8_INITS initializes a Chebyshev series.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 January 2012

  Author:

    C version by John Burkardt.

  Reference:

    Roger Broucke,
    Algorithm 446:
    Ten Subroutines for the Manipulation of Chebyshev Series,
    Communications of the ACM,
    Volume 16, Number 4, April 1973, pages 254-256.

  Parameters:

    Input, double DOS[NOS], the Chebyshev coefficients.

    Input, int NOS, the number of coefficients.

    Input, double ETA, the desired accuracy.

    Output, int R8_INITS, the number of terms of the series needed
    to ensure the requested accuracy.
*/
{
  double err;
  int i;
  int value;

  if ( nos < 1 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8_INITS - Fatal error!\n" );
    fprintf ( stderr, "  Number of coefficients < 1.\n" );
    exit ( 1 );
  }

  err = 0.0;

  for ( i = nos - 1; 0 <= i; i-- )
  {
    err = err + r8_abs ( dos[i] );
    if ( eta < err )
    {
      value = i + 1;
      return value;
    }
  }

  value = i;
  fprintf ( stderr, "\n" );
  fprintf ( stderr, "R8_INITS - Warning!\n" );
  fprintf ( stderr, "  ETA may be too small.\n" );

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

double r8_lgmc ( double x )

/******************************************************************************/
/*
  Purpose:

    R8_LGMC evaluates the log gamma correction factor for an R8 argument.

  Discussion:

    For 10 <= X, compute the log gamma correction factor so that

      log ( gamma ( x ) ) = log ( sqrt ( 2 * pi ) ) 
                          + ( x - 0.5 ) * log ( x ) - x 
                          + r8_lgmc ( x )

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 January 2012

  Author:

    Original FORTRAN77 version by Wayne Fullerton.
    C version by John Burkardt.

  Reference:

    Wayne Fullerton,
    Portable Special Function Routines,
    in Portability of Numerical Software,
    edited by Wayne Cowell,
    Lecture Notes in Computer Science, Volume 57,
    Springer 1977,
    ISBN: 978-3-540-08446-4,
    LC: QA297.W65.

  Parameters:

    Input, double X, the argument.

    Output, double R8_LGMC, the correction factor.
*/
{
  static double algmcs[15] = {
    +0.1666389480451863247205729650822,
    -0.1384948176067563840732986059135E-04,
    +0.9810825646924729426157171547487E-08,
    -0.1809129475572494194263306266719E-10,
    +0.6221098041892605227126015543416E-13,
    -0.3399615005417721944303330599666E-15,
    +0.2683181998482698748957538846666E-17,
    -0.2868042435334643284144622399999E-19,
    +0.3962837061046434803679306666666E-21,
    -0.6831888753985766870111999999999E-23,
    +0.1429227355942498147573333333333E-24,
    -0.3547598158101070547199999999999E-26,
    +0.1025680058010470912000000000000E-27,
    -0.3401102254316748799999999999999E-29,
    +0.1276642195630062933333333333333E-30 };
  static int nalgm = 0;
  double value;
  static double xbig = 0.0;
  static double xmax = 0.0;

  if ( nalgm == 0 )
  {
    nalgm = r8_inits ( algmcs, 15, r8_mach ( 3 ) );
    xbig = 1.0 / sqrt ( r8_mach ( 3 ) );
    xmax = exp ( r8_min ( log ( r8_mach ( 2 ) / 12.0 ), 
      - log ( 12.0 * r8_mach ( 1 ) ) ) );
  }

  if ( x < 10.0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8_LGMC - Fatal error!\n" );
    fprintf ( stderr, "  X must be at least 10.\n" );
    exit ( 1 );
  }
  else if ( x < xbig )
  {
    value = r8_csevl ( 2.0 * ( 10.0 / x ) 
      * ( 10.0 / x ) - 1.0, algmcs, nalgm ) / x;
  }
  else if ( x < xmax )
  {
    value = 1.0 / ( 12.0 * x );
  }
  else
  {
    value = 0.0;
  }
  return value;
}
/******************************************************************************/

double r8_mach ( int i )

/******************************************************************************/
/*
  Purpose:

    R8_MACH returns double precision real machine constants.

  Discussion:

    Assuming that the internal representation of a double precision real
    number is in base B, with T the number of base-B digits in the mantissa,
    and EMIN the smallest possible exponent and EMAX the largest possible 
    exponent, then

      R8_MACH(1) = B^(EMIN-1), the smallest positive magnitude.
      R8_MACH(2) = B^EMAX*(1-B^(-T)), the largest magnitude.
      R8_MACH(3) = B^(-T), the smallest relative spacing.
      R8_MACH(4) = B^(1-T), the largest relative spacing.
      R8_MACH(5) = log10(B).

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    24 April 2007

  Author:

    Original FORTRAN77 version by Phyllis Fox, Andrew Hall, Norman Schryer.
    C version by John Burkardt.

  Reference:

    Phyllis Fox, Andrew Hall, Norman Schryer,
    Algorithm 528:
    Framework for a Portable Library,
    ACM Transactions on Mathematical Software,
    Volume 4, Number 2, June 1978, page 176-188.

  Parameters:

    Input, int I, chooses the parameter to be returned.
    1 <= I <= 5.

    Output, double R8_MACH, the value of the chosen parameter.
*/
{
  double value;

  if ( i < 1 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8_MACH - Fatal error!\n" );
    fprintf ( stderr, "  The input argument I is out of bounds.\n" );
    fprintf ( stderr, "  Legal values satisfy 1 <= I <= 5.\n" );
    fprintf ( stderr, "  I = %d\n", i );
    value = 0.0;
    exit ( 1 );
  }
  else if ( i == 1 )
  {
    value = 4.450147717014403E-308;
  }
  else if ( i == 2 )
  {
    value = 8.988465674311579E+307;
  }
  else if ( i == 3 )
  {
    value = 1.110223024625157E-016;
  }
  else if ( i == 4 )
  {
    value = 2.220446049250313E-016;
  }
  else if ( i == 5 )
  {
    value = 0.301029995663981E+000;
  }
  else if ( 5 < i )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8_MACH - Fatal error!\n" );
    fprintf ( stderr, "  The input argument I is out of bounds.\n" );
    fprintf ( stderr, "  Legal values satisfy 1 <= I <= 5.\n" );
    fprintf ( stderr, "  I = %d\n", i );
    value = 0.0;
    exit ( 1 );
  }

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

double r8_max ( double x, double y )

/******************************************************************************/
/*
  Purpose:

    R8_MAX returns the maximum of two R8's.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    07 May 2006

  Author:

    John Burkardt

  Parameters:

    Input, double X, Y, the quantities to compare.

    Output, double R8_MAX, the maximum of X and Y.
*/
{
  double value;

  if ( y < x )
  {
    value = x;
  }
  else
  {
    value = y;
  }
  return value;
}
/******************************************************************************/

double r8_min ( double x, double y )

/******************************************************************************/
/*
  Purpose:

    R8_MIN returns the minimum of two R8's.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    07 May 2006

  Author:

    John Burkardt

  Parameters:

    Input, double X, Y, the quantities to compare.

    Output, double R8_MIN, the minimum of X and Y.
*/
{
  double value;

  if ( y < x )
  {
    value = y;
  }
  else
  {
    value = x;
  }
  return value;
}
/******************************************************************************/

double r8_sech ( double x )

/******************************************************************************/
/*
  Purpose:

    R8_SECH evaluates the hyperbolic secant, while avoiding COSH overflow.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    21 December 2011

  Author:

    John Burkardt

  Parameters:

    Input, double X, the argument of the function.

    Output, double R8_SECH, the value of the function.
*/
{
  double log_huge = 80.0;
  double value;

  if ( log_huge < r8_abs ( x ) )
  {
    value = 0.0;
  }
  else
  {
    value = 1.0 / cosh ( x );
  }
  return value;
}
/******************************************************************************/

double r8_si ( double x )

/******************************************************************************/
/*
  Purpose:

    R8_SI evaluates the sine integral Si of an R8 argument.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 January 2012

  Author:

    Original FORTRAN77 version by Wayne Fullerton.
    C version by John Burkardt.

  Reference:

    Wayne Fullerton,
    Portable Special Function Routines,
    in Portability of Numerical Software,
    edited by Wayne Cowell,
    Lecture Notes in Computer Science, Volume 57,
    Springer 1977,
    ISBN: 978-3-540-08446-4,
    LC: QA297.W65.

  Parameters:

    Input, double X, the argument.

    Output, double R8_SI, the sine integral Si evaluated at X.
*/
{
  double absx;
  double cosx;
  double f;
  double g;
  static int nsi = 0;
  const double pi2 = 1.57079632679489661923132169163975;
  static double sics[18] = {
    -0.1315646598184841928904275173000457,
    -0.2776578526973601892048287660157299,
     0.0354414054866659179749135464710086,
    -0.0025631631447933977658752788361530,
     0.0001162365390497009281264921482985,
    -0.0000035904327241606042670004347148,
     0.0000000802342123705710162308652976,
    -0.0000000013562997692540250649931846,
     0.0000000000179440721599736775567759,
    -0.0000000000001908387343087145490737,
     0.0000000000000016669989586824330853,
    -0.0000000000000000121730988368503042,
     0.0000000000000000000754181866993865,
    -0.0000000000000000000004014178842446,
     0.0000000000000000000000018553690716,
    -0.0000000000000000000000000075166966,
     0.0000000000000000000000000000269113,
    -0.0000000000000000000000000000000858 };
  double value;
  static double xsml = 0.0;

  if ( nsi == 0 )
  {
    nsi = r8_inits ( sics, 18, 0.1 * r8_mach ( 3 ) );
    xsml = sqrt ( r8_mach ( 3 ) );
  }

  absx = r8_abs ( x );

  if ( absx < xsml )
  {
    value = x;
  }
  else if ( absx <= 4.0 )
  {
    value = x * ( 0.75 + r8_csevl ( ( x * x - 8.0 ) * 0.125, sics, nsi ) );
  }
  else
  {
    r8_sifg ( absx, &f, &g );
    cosx = cos ( absx );
    value = pi2 - f * cosx - g * sin ( x );
    if ( x < 0.0 )
    {
      value = - value;
    }
  }
  return value;
}
/******************************************************************************/

void r8_sifg ( double x, double *f, double *g )

/******************************************************************************/
/*
  Purpose:

    R8_SIFG is a utility routine.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 January 2012

  Author:

    Original FORTRAN77 version by Wayne Fullerton.
    C version by John Burkardt.

  Reference:

    Wayne Fullerton,
    Portable Special Function Routines,
    in Portability of Numerical Software,
    edited by Wayne Cowell,
    Lecture Notes in Computer Science, Volume 57,
    Springer 1977,
    ISBN: 978-3-540-08446-4,
    LC: QA297.W65.

  Parameters:

    Input, double X, the argument.

    Output, double *F, *G.
*/
{
  static double f1cs[43] = {
-0.1191081969051363610348201965828918,
-0.0247823144996236247590074150823133,
 0.0011910281453357821268120363054457,
-0.0000927027714388561748308600360706,
 0.0000093373141568270996868204582766,
-0.0000011058287820557143938979426306,
 0.0000001464772071460162169336550799,
-0.0000000210694496287689532601227548,
 0.0000000032293492366848236382857374,
-0.0000000005206529617529375828014986,
 0.0000000000874878884570278750268316,
-0.0000000000152176187056123668294574,
 0.0000000000027257192405419573900583,
-0.0000000000005007053075968556290255,
 0.0000000000000940240902726068511779,
-0.0000000000000180014444791803678336,
 0.0000000000000035062621432741785826,
-0.0000000000000006935282926769149709,
 0.0000000000000001390925136454216568,
-0.0000000000000000282486885074170585,
 0.0000000000000000058031305693579081,
-0.0000000000000000012046901573375820,
 0.0000000000000000002525052443655940,
-0.0000000000000000000533980268805594,
 0.0000000000000000000113855786274122,
-0.0000000000000000000024462861505259,
 0.0000000000000000000005293659320439,
-0.0000000000000000000001153184940277,
 0.0000000000000000000000252786568318,
-0.0000000000000000000000055738645378,
 0.0000000000000000000000012358245621,
-0.0000000000000000000000002754350842,
 0.0000000000000000000000000616906808,
-0.0000000000000000000000000138817443,
 0.0000000000000000000000000031375329,
-0.0000000000000000000000000007121249,
 0.0000000000000000000000000001622778,
-0.0000000000000000000000000000371206,
 0.0000000000000000000000000000085221,
-0.0000000000000000000000000000019633,
 0.0000000000000000000000000000004538,
-0.0000000000000000000000000000001052,
 0.0000000000000000000000000000000245 };
  static double f2cs[99] = {
-0.03484092538970132330836049733745577,
-0.01668422056779596873246786312278676,
 0.00067529012412377385045207859239727,
-0.00005350666225447013628785577557429,
 0.00000626934217790075267050759431626,
-0.00000095266388019916680677790414293,
 0.00000017456292242509880425504427666,
-0.00000003687954030653093307097646628,
 0.00000000872026777051395264075816938,
-0.00000000226019703919738748530423167,
 0.00000000063246249765250612520444877,
-0.00000000018889118884717869240911480,
 0.00000000005967746729997813372620472,
-0.00000000001980443117372239011196007,
 0.00000000000686413954772103383713264,
-0.00000000000247310193070199106074890,
 0.00000000000092263594549941404196042,
-0.00000000000035523634999261784497297,
 0.00000000000014076049625351591461820,
-0.00000000000005726228499747652794311,
 0.00000000000002386537545413171810106,
-0.00000000000001017141890764597142232,
 0.00000000000000442594531078364424968,
-0.00000000000000196344933049189761979,
 0.00000000000000088688748314810461024,
-0.00000000000000040743345027311546948,
 0.00000000000000019016837215675339859,
-0.00000000000000009009707297478042442,
 0.00000000000000004329211274095668667,
-0.00000000000000002108144465322479526,
 0.00000000000000001039637907026452274,
-0.00000000000000000518891007948931936,
 0.00000000000000000261955324869899371,
-0.00000000000000000133690399951301570,
 0.00000000000000000068941057702931664,
-0.00000000000000000035905362610437250,
 0.00000000000000000018878077255791706,
-0.00000000000000000010016125265594380,
 0.00000000000000000005360725691578228,
-0.00000000000000000002893198974944827,
 0.00000000000000000001574065100202625,
-0.00000000000000000000863027106431206,
 0.00000000000000000000476715602862288,
-0.00000000000000000000265222739998504,
 0.00000000000000000000148582865063866,
-0.00000000000000000000083797235923135,
 0.00000000000000000000047565916422711,
-0.00000000000000000000027169073353112,
 0.00000000000000000000015612738881686,
-0.00000000000000000000009024555078347,
 0.00000000000000000000005246097049119,
-0.00000000000000000000003066450818697,
 0.00000000000000000000001801996250957,
-0.00000000000000000000001064443050752,
 0.00000000000000000000000631942158881,
-0.00000000000000000000000377013812246,
 0.00000000000000000000000225997542918,
-0.00000000000000000000000136100844814,
 0.00000000000000000000000082333232003,
-0.00000000000000000000000050025986091,
 0.00000000000000000000000030526245684,
-0.00000000000000000000000018705164021,
 0.00000000000000000000000011508404393,
-0.00000000000000000000000007108714611,
 0.00000000000000000000000004408065533,
-0.00000000000000000000000002743760867,
 0.00000000000000000000000001714144851,
-0.00000000000000000000000001074768860,
 0.00000000000000000000000000676259777,
-0.00000000000000000000000000426981348,
 0.00000000000000000000000000270500637,
-0.00000000000000000000000000171933331,
 0.00000000000000000000000000109636138,
-0.00000000000000000000000000070132573,
 0.00000000000000000000000000045001784,
-0.00000000000000000000000000028963835,
 0.00000000000000000000000000018697009,
-0.00000000000000000000000000012104646,
 0.00000000000000000000000000007859065,
-0.00000000000000000000000000005116867,
 0.00000000000000000000000000003340627,
-0.00000000000000000000000000002186851,
 0.00000000000000000000000000001435340,
-0.00000000000000000000000000000944523,
 0.00000000000000000000000000000623117,
-0.00000000000000000000000000000412101,
 0.00000000000000000000000000000273208,
-0.00000000000000000000000000000181558,
 0.00000000000000000000000000000120934,
-0.00000000000000000000000000000080737,
 0.00000000000000000000000000000054022,
-0.00000000000000000000000000000036227,
 0.00000000000000000000000000000024348,
-0.00000000000000000000000000000016401,
 0.00000000000000000000000000000011074,
-0.00000000000000000000000000000007497,
 0.00000000000000000000000000000005091,
-0.00000000000000000000000000000003470,
 0.00000000000000000000000000000002377 };
  static double g1cs[44] = {
-0.3040578798253495954499726682091083,
-0.0566890984597120587731339156118269,
 0.0039046158173275643919984071554082,
-0.0003746075959202260618619339867489,
 0.0000435431556559843679552220840065,
-0.0000057417294453025046561970723475,
 0.0000008282552104502629741937616492,
-0.0000001278245892594642727883913223,
 0.0000000207978352948687884439257529,
-0.0000000035313205921990798042032682,
 0.0000000006210824236308951068631449,
-0.0000000001125215474446292649336987,
 0.0000000000209088917684421605267019,
-0.0000000000039715831737681727689158,
 0.0000000000007690431314272089939005,
-0.0000000000001514696742731613519826,
 0.0000000000000302892146552359684119,
-0.0000000000000061399703834708825400,
 0.0000000000000012600605829510933553,
-0.0000000000000002615029250939483683,
 0.0000000000000000548278844891796821,
-0.0000000000000000116038182129526571,
 0.0000000000000000024771654107129795,
-0.0000000000000000005330672753223389,
 0.0000000000000000001155666075598465,
-0.0000000000000000000252280547744957,
 0.0000000000000000000055429038550786,
-0.0000000000000000000012252208421297,
 0.0000000000000000000002723664318684,
-0.0000000000000000000000608707831422,
 0.0000000000000000000000136724874476,
-0.0000000000000000000000030856626806,
 0.0000000000000000000000006995212319,
-0.0000000000000000000000001592587569,
 0.0000000000000000000000000364051056,
-0.0000000000000000000000000083539465,
 0.0000000000000000000000000019240303,
-0.0000000000000000000000000004446816,
 0.0000000000000000000000000001031182,
-0.0000000000000000000000000000239887,
 0.0000000000000000000000000000055976,
-0.0000000000000000000000000000013100,
 0.0000000000000000000000000000003074,
-0.0000000000000000000000000000000723 };
  static double g2cs[44] = {
-0.1211802894731646263541834046858267,
-0.0316761386394950286701407923505610,
 0.0013383199778862680163819429492182,
-0.0000895511011392252425531905069518,
 0.0000079155562961718213115249467924,
-0.0000008438793322241520181418982080,
 0.0000001029980425677530146647227274,
-0.0000000139295750605183835795834444,
 0.0000000020422703959875980400677594,
-0.0000000003196534694206427035434752,
 0.0000000000528147832657267698615312,
-0.0000000000091339554672671033735289,
 0.0000000000016426251238967760444819,
-0.0000000000003055897039322660002410,
 0.0000000000000585655825785779717892,
-0.0000000000000115229197730940120563,
 0.0000000000000023209469119988537310,
-0.0000000000000004774355834177535025,
 0.0000000000000001000996765800180573,
-0.0000000000000000213533778082256704,
 0.0000000000000000046277190777367671,
-0.0000000000000000010175807410227657,
 0.0000000000000000002267657399884672,
-0.0000000000000000000511630776076426,
 0.0000000000000000000116767014913108,
-0.0000000000000000000026935427672470,
 0.0000000000000000000006275665841146,
-0.0000000000000000000001475880557531,
 0.0000000000000000000000350145314739,
-0.0000000000000000000000083757732152,
 0.0000000000000000000000020191815152,
-0.0000000000000000000000004903567705,
 0.0000000000000000000000001199123348,
-0.0000000000000000000000000295170610,
 0.0000000000000000000000000073113112,
-0.0000000000000000000000000018217843,
 0.0000000000000000000000000004565148,
-0.0000000000000000000000000001150151,
 0.0000000000000000000000000000291267,
-0.0000000000000000000000000000074125,
 0.0000000000000000000000000000018953,
-0.0000000000000000000000000000004868,
 0.0000000000000000000000000000001256,
-0.0000000000000000000000000000000325 };
  static double g3cs[56] = {
-0.0280574367809472928402815264335299,
-0.0137271597162236975409100508089556,
 0.0002894032638760296027448941273751,
-0.0000114129239391197145908743622517,
 0.0000006813965590726242997720207302,
-0.0000000547952289604652363669058052,
 0.0000000055207429918212529109406521,
-0.0000000006641464199322920022491428,
 0.0000000000922373663487041108564960,
-0.0000000000144299088886682862611718,
 0.0000000000024963904892030710248705,
-0.0000000000004708240675875244722971,
 0.0000000000000957217659216759988140,
-0.0000000000000207889966095809030537,
 0.0000000000000047875099970877431627,
-0.0000000000000011619070583377173759,
 0.0000000000000002956508969267836974,
-0.0000000000000000785294988256492025,
 0.0000000000000000216922264368256612,
-0.0000000000000000062113515831676342,
 0.0000000000000000018384568838450977,
-0.0000000000000000005610887482137276,
 0.0000000000000000001761862805280062,
-0.0000000000000000000568111050541451,
 0.0000000000000000000187786279582313,
-0.0000000000000000000063531694151124,
 0.0000000000000000000021968802368238,
-0.0000000000000000000007754666550395,
 0.0000000000000000000002791018356581,
-0.0000000000000000000001023178525247,
 0.0000000000000000000000381693403919,
-0.0000000000000000000000144767895606,
 0.0000000000000000000000055779512634,
-0.0000000000000000000000021817239071,
 0.0000000000000000000000008656646309,
-0.0000000000000000000000003482157895,
 0.0000000000000000000000001419188130,
-0.0000000000000000000000000585714314,
 0.0000000000000000000000000244660482,
-0.0000000000000000000000000103387099,
 0.0000000000000000000000000044177299,
-0.0000000000000000000000000019080079,
 0.0000000000000000000000000008326038,
-0.0000000000000000000000000003669553,
 0.0000000000000000000000000001632875,
-0.0000000000000000000000000000733357,
 0.0000000000000000000000000000332327,
-0.0000000000000000000000000000151906,
 0.0000000000000000000000000000070020,
-0.0000000000000000000000000000032539,
 0.0000000000000000000000000000015240,
-0.0000000000000000000000000000007193,
 0.0000000000000000000000000000003420,
-0.0000000000000000000000000000001638,
 0.0000000000000000000000000000000790,
-0.0000000000000000000000000000000383 };
  static int nf1 = 0;
  static int nf2 = 0;
  static int ng1 = 0;
  static int ng2 = 0;
  static int ng3 = 0;
  double tol;
  static double xbig = 0.0;
  static double xbnd = 0.0;
  static double xbndg = 0.0;
  static double xmaxf = 0.0;
  static double xmaxg = 0.0;

  if ( nf1 == 0 )
  {
    tol = 0.1 * r8_mach ( 3 );
    nf1 = r8_inits ( f1cs, 43, tol );
    nf2 = r8_inits ( f2cs, 99, tol );
    ng1 = r8_inits ( g1cs, 44, tol );
    ng2 = r8_inits ( g2cs, 44, tol );
    ng3 = r8_inits ( g3cs, 56, tol );
    xbig = sqrt ( 1.0 / r8_mach ( 3 ) );
    xmaxf = exp ( r8_min ( - log ( r8_mach ( 1 ) ), 
      log ( r8_mach ( 2 ) ) ) - 0.01 );
    xmaxg = 1.0 / sqrt ( r8_mach ( 1 ) );
    xbnd = sqrt ( 50.0 );
    xbndg = sqrt ( 200.0 );
  }

  if ( x < 4.0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8_SIFG - Fatal error!\n" );
    fprintf ( stderr, "  Approximation invalid for X < 4.\n" );
    exit ( 1 );
  }
  else if ( x <= xbnd )
  {
    *f = ( 1.0 + r8_csevl ( ( 1.0 / x / x - 0.04125 )
      / 0.02125, f1cs, nf1 ) ) / x;
    *g = ( 1.0 + r8_csevl ( ( 1.0 / x / x - 0.04125 )
      / 0.02125, g1cs, ng1 ) ) / x / x;
  }
  else if ( x <= xbig )
  {
    *f = ( 1.0 + r8_csevl ( 100. / x / x - 1.0, f2cs, nf2 ) ) / x;
    if ( x <= xbndg )
    {
      *g = ( 1.0 + r8_csevl ( ( 10000.0 / x / x - 125.0 ) 
        / 75.0, g2cs, ng2 ) ) / x / x;
    }
    else
    {
      *g = ( 1.0 + r8_csevl ( 400.0 / x / x - 1.0, g3cs, ng3 ) ) / x / x;
    }
  }
  else
  {
    if ( x < xmaxf )
    {
      *f = 1.0 / x;
    }
    else
    {
      *f = 0.0;
    }
    if ( x < xmaxg )
    {
      *g = 1.0 / x / x;
    }
    else
    {
      *g = 0.0;
    }
  }
  return;
}
/******************************************************************************/

double r8vec_dot_product ( int n, double a1[], double a2[] )

/******************************************************************************/
/*
  Purpose:

    R8VEC_DOT_PRODUCT computes the dot product of a pair of R8VEC's.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    26 July 2007

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of entries in the vectors.

    Input, double A1[N], A2[N], the two vectors to be considered.

    Output, double R8VEC_DOT_PRODUCT, the dot product of the vectors.
*/
{
  int i;
  double value;

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

double *r8vec_linspace_new ( int n, double a_first, double a_last )

/******************************************************************************/
/*
  Purpose:

    R8VEC_LINSPACE_NEW creates a vector of linearly spaced values.

  Discussion:

    An R8VEC is a vector of R8's.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    29 March 2011

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of entries in the vector.

    Input, double A_FIRST, A_LAST, the first and last entries.

    Output, double R8VEC_LINSPACE_NEW[N], a vector of linearly spaced data.
*/
{
  double *a;
  int i;

  a = ( double * ) malloc ( n * sizeof ( double ) );

  if ( n == 1 )
  {
    a[0] = ( a_first + a_last ) / 2.0;
  }
  else
  {
    for ( i = 0; i < n; i++ )
    {
      a[i] = ( ( double ) ( n - 1 - i ) * a_first 
             + ( double ) (         i ) * a_last ) 
             / ( double ) ( n - 1     );
    }
  }
  return a;
}
/******************************************************************************/

double r8vec_sum ( int n, double a[] )

/******************************************************************************/
/*
  Purpose:

    R8VEC_SUM returns the sum of an R8VEC.

  Discussion:

    An R8VEC is a vector of R8's.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    26 August 2008

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of entries in the vector.

    Input, double A[N], the vector.

    Output, double R8VEC_SUM, the sum of the vector.
*/
{
  int i;
  double value;

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

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

double *r8vec_uniform_01_new ( int n, int *seed )

/******************************************************************************/
/*
  Purpose:

    R8VEC_UNIFORM_01_NEW returns a unit pseudorandom R8VEC.

  Discussion:

    This routine implements the recursion

      seed = 16807 * seed mod ( 2^31 - 1 )
      unif = seed / ( 2^31 - 1 )

    The integer arithmetic never requires more than 32 bits,
    including a sign bit.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    19 August 2004

  Author:

    John Burkardt

  Reference:

    Paul Bratley, Bennett Fox, Linus Schrage,
    A Guide to Simulation,
    Second Edition,
    Springer, 1987,
    ISBN: 0387964673,
    LC: QA76.9.C65.B73.

    Bennett Fox,
    Algorithm 647:
    Implementation and Relative Efficiency of Quasirandom
    Sequence Generators,
    ACM Transactions on Mathematical Software,
    Volume 12, Number 4, December 1986, pages 362-376.

    Pierre L'Ecuyer,
    Random Number Generation,
    in Handbook of Simulation,
    edited by Jerry Banks,
    Wiley, 1998,
    ISBN: 0471134031,
    LC: T57.62.H37.

    Peter Lewis, Allen Goodman, James Miller,
    A Pseudo-Random Number Generator for the System/360,
    IBM Systems Journal,
    Volume 8, Number 2, 1969, pages 136-143.

  Parameters:

    Input, int N, the number of entries in the vector.

    Input/output, int *SEED, a seed for the random number generator.

    Output, double R8VEC_UNIFORM_01_NEW[N], the vector of pseudorandom values.
*/
{
  int i;
  int i4_huge = 2147483647;
  int k;
  double *r;

  if ( *seed == 0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8VEC_UNIFORM_01_NEW - Fatal error!\n" );
    fprintf ( stderr, "  Input value of SEED = 0.\n" );
    exit ( 1 );
  }

  r = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    k = *seed / 127773;

    *seed = 16807 * ( *seed - k * 127773 ) - k * 2836;

    if ( *seed < 0 )
    {
      *seed = *seed + i4_huge;
    }

    r[i] = ( double ) ( *seed ) * 4.656612875E-10;
  }

  return r;
}
/******************************************************************************/

int s_eqi ( char *s1, char *s2 )

/******************************************************************************/
/*
  Purpose:

    S_EQI reports whether two strings are equal, ignoring case.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    03 June 2008

  Author:

    John Burkardt

  Parameters:

    Input, char *S1, char *S2, pointers to two strings.

    Output, int S_EQI, is true if the strings are equal.
*/
{
  int i;
  int nchar;
  int nchar1;
  int nchar2;

  nchar1 = strlen ( s1 );
  nchar2 = strlen ( s2 );
  if ( nchar1 < nchar2 )
  {
    nchar = nchar1;
  }
  else
  {
    nchar = nchar2;
  }
/*
  The strings are not equal if they differ over their common length.
*/
  for ( i = 0; i < nchar; i++ )
  {

    if ( ch_cap ( s1[i] ) != ch_cap ( s2[i] ) )
    {
      return 0;
    }
  }
/*
  The strings are not equal if the longer one includes nonblanks
  in the tail.
*/
  if ( nchar < nchar1 )
  {
    for ( i = nchar; i < nchar1; i++ )
    {
      if ( s1[i] != ' ' )
      {
        return 0;
      }
    }
  }
  else if ( nchar < nchar2 )
  {
    for ( i = nchar; i < nchar2; i++ )
    {
      if ( s2[i] != ' ' )
      {
        return 0;
      }
    }
  }

  return 1;
}
/******************************************************************************/

void timestamp ( void )

/******************************************************************************/
/*
  Purpose:

    TIMESTAMP prints the current YMDHMS date as a time stamp.

  Example:

    31 May 2001 09:45:54 AM

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    24 September 2003

  Author:

    John Burkardt

  Parameters:

    None
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  printf ( "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}