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

# include "laguerre_integrands.h"

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

void laguerre_compute ( int order, double xtab[], double weight[], 
  double alpha )

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

    LAGUERRE_COMPUTE computes a Gauss-Laguerre quadrature rule.

  Discussion:

    In the simplest case, ALPHA is 0, and we are approximating the
    integral from 0 to +oo of EXP(-X) * F(X).  When this is so,
    it is easy to modify the rule to approximate the integral from
    A to +oo as well.

    If ALPHA is nonzero, then there is no simple way to extend the
    rule to approximate the integral from A to +oo.  The simplest
    procedures would be to approximate the integral from 0 to A.

    The integration interval is [ A, +oo ) or [ 0, +oo ).

    The weight function is w(x) = exp ( -x ) or exp ( -x ) * x^alpha.


    If the integral to approximate is:

        Integral ( A <= X < +oo ) EXP ( - X ) * F(X) dX
      or
        Integral ( 0 <= X < +oo ) EXP ( - X ) * X^ALPHA * F(X) dX

    then the quadrature rule is:

      EXP ( - A ) * Sum ( 1 <= I <= ORDER ) WEIGHT(I) * F ( A+XTAB(I) )
    or
      sum ( 1 <= I <= ORDER ) WEIGHT(I) * F ( XTAB(I) )

    If the integral to approximate is:

        Integral ( A <= X < +oo ) F(X) dX
      or
        Integral ( 0 <= X < +oo ) X^ALPHA * F(X) dX

    then the quadrature rule is:

      EXP ( - A ) * Sum ( 1 <= I <= ORDER ) 
        WEIGHT(I) * EXP(A+XTAB(I)) * F ( A+XTAB(I) )
    or
      sum ( 1 <= I <= ORDER ) WEIGHT(I) * EXP(XTAB(I)) * F ( XTAB(I) )

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    02 May 2006

  Author:

    Original FORTRAN77 version by Arthur Stroud, Don Secrest.
    C version by John Burkardt.

  Reference:

    Arthur Stroud, Don Secrest,
    Gaussian Quadrature Formulas,
    Prentice Hall, 1966,
    LC: QA299.4G3S7.

  Parameters:

    Input, int ORDER, the order of the quadrature rule to be computed.
    ORDER must be at least 1.

    Output, double XTAB[ORDER], the Gauss-Laguerre abscissas.

    Output, double WEIGHT[ORDER], the Gauss-Laguerre weights.

    Input, double ALPHA, the exponent of the X factor.
    Set ALPHA = 0.0 for the simplest rule.
    ALPHA must be nonnegative.
*/
{
  double *b;
  double *c;
  double cc;
  double dp2;
  int i;
  double p1;
  double prod;
  double r1;
  double r2;
  double ratio;
  double x;

  b = ( double * ) malloc ( order * sizeof ( double ) );
  c = ( double * ) malloc ( order * sizeof ( double ) );
/*
  Set the recursion coefficients.
*/
  for ( i = 0; i < order; i++ )
  {
    b[i] = ( alpha + ( double ) ( 2 * i + 1 ) );
  }

  for ( i = 0; i < order; i++ )
  {
    c[i] = ( double ) ( i ) * ( alpha + ( double ) ( i ) );
  }
  prod = 1.0;
  for ( i = 1; i < order; i++ )
  {
    prod = prod * c[i];
  }
  cc = tgamma ( alpha + 1.0 ) * prod;

  for ( i = 0; i < order; i++ )
  {
/*
  Compute an estimate for the root.
*/
    if ( i == 0 )
    {
      x = ( 1.0 + alpha ) * ( 3.0+ 0.92 * alpha ) / 
        ( 1.0 + 2.4 * ( double ) ( order ) + 1.8 * alpha );
    }
    else if ( i == 1 )
    {
      x = x + ( 15.0 + 6.25 * alpha ) / 
        ( 1.0 + 0.9 * alpha + 2.5 * ( double ) ( order ) );
    }
    else
    {
      r1 = ( 1.0 + 2.55 * ( double ) ( i - 1 ) ) 
        / ( 1.9 * ( double ) ( i - 1 ) );

      r2 = 1.26 * ( double ) ( i - 1 ) * alpha / 
        ( 1.0 + 3.5 * ( double ) ( i - 1 ) );

      ratio = ( r1 + r2 ) / ( 1.0 + 0.3 * alpha );

      x = x + ratio * ( x - xtab[i-2] );
    }
/*
  Use iteration to find the root.
*/
    laguerre_root ( &x, order, alpha, &dp2, &p1, b, c );
/*
  Set the abscissa and weight.
*/
    xtab[i] = x;
    weight[i] = ( cc / dp2 ) / p1;
  }

  free ( b );
  free ( c );

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

void laguerre_recur ( double *p2, double *dp2, double *p1, double x, 
  int order, double alpha, double b[], double c[] )

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

    LAGUERRE_RECUR finds the value and derivative of a Laguerre polynomial.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    03 May 2006

  Author:

    Original FORTRAN77 version by Arthur Stroud, Don Secrest.
    C version by John Burkardt.

  Reference:

    Arthur Stroud, Don Secrest,
    Gaussian Quadrature Formulas,
    Prentice Hall, 1966,
    LC: QA299.4G3S7.

  Parameters:

    Output, double *P2, the value of L(ORDER)(X).

    Output, double *DP2, the value of L'(ORDER)(X).

    Output, double *P1, the value of L(ORDER-1)(X).

    Input, double X, the point at which polynomials are evaluated.

    Input, int ORDER, the order of the polynomial to be computed.

    Input, double ALPHA, the exponent of the X factor in the
    integrand.

    Input, double B[ORDER], C[ORDER], the recursion coefficients.
*/
{
  double dp0;
  double dp1;
  int i;
  double p0;

  *p1 = 1.0;
  dp1 = 0.0;

  *p2 = x - alpha - 1.0;
  *dp2 = 1.0;

  for ( i = 1; i < order; i++ )
  {
    p0 = *p1;
    dp0 = dp1;

    *p1 = *p2;
    dp1 = *dp2;

    *p2  = ( x - b[i] ) * ( *p1 ) - c[i] * p0;
    *dp2 = ( x - b[i] ) * dp1 + ( *p1 ) - c[i] * dp0;
  }

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

void laguerre_root ( double *x, int order, double alpha, double *dp2, 
  double *p1, double b[], double c[] )

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

    LAGUERRE_ROOT improves an approximate root of a Laguerre polynomial.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    03 May 2006

  Author:

    Original FORTRAN77 version by Arthur Stroud, Don Secrest.
    C version by John Burkardt.

  Reference:

    Arthur Stroud, Don Secrest,
    Gaussian Quadrature Formulas,
    Prentice Hall, 1966,
    LC: QA299.4G3S7.

  Parameters:

    Input/output, double *X, the approximate root, which
    should be improved on output.

    Input, int ORDER, the order of the polynomial to be computed.

    Input, double ALPHA, the exponent of the X factor.

    Output, double *DP2, the value of L'(ORDER)(X).

    Output, double *P1, the value of L(ORDER-1)(X).

    Input, double B[ORDER], C[ORDER], the recursion coefficients.
*/
{
  double d;
  double eps;
  double p2;
  int step;
  int step_max = 10;

  eps = DBL_EPSILON;

  for ( step = 1; step <= step_max; step++ )
  {
    laguerre_recur ( &p2, dp2, p1, *x, order, alpha, b, c );

    d = p2 / ( *dp2 );
    *x = *x - d;

    if ( fabs ( d ) <= eps * ( fabs ( *x ) + 1.0 ) )
    {
      break;
    }
  }

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

void legendre_compute ( int order, double xtab[], double weight[] )

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

    LEGENDRE_COMPUTE computes a Gauss-Legendre quadrature rule.

  Discussion:

    The integration interval is [ -1, 1 ].

    The weight function is w(x) = 1.0.

    The integral to approximate:

      Integral ( -1 <= X <= 1 ) F(X) dX

    The quadrature rule:

      Sum ( 1 <= I <= ORDER ) WEIGHT(I) * F ( XTAB(I) )

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 April 2006

  Author:

    Original FORTRAN77 version by Philip Davis, Philip Rabinowitz.
    C version by John Burkardt.

  Reference:

    Philip Davis, Philip Rabinowitz,
    Methods of Numerical Integration,
    Second Edition,
    Dover, 2007,
    ISBN: 0486453391,
    LC: QA299.3.D28.

  Parameters:

    Input, int ORDER, the order of the rule.
    ORDER must be greater than 0.

    Output, double XTAB[ORDER], the abscissas of the rule.

    Output, double WEIGHT[ORDER], the weights of the rule.
    The weights are positive, symmetric, and should sum to 2.
*/
{
  double d1;
  double d2pn;
  double d3pn;
  double d4pn;
  double dp;
  double dpn;
  double e1;
  double fx;
  double h;
  int i;
  int iback;
  int k;
  int m;
  int mp1mi;
  int ncopy;
  int nmove;
  double p;
  double pk;
  double pkm1;
  double pkp1;
  const double r8_pi = 3.1415926535897932385;
  double t;
  double u;
  double v;
  double x0;
  double xtemp;

  if ( order < 1 )
  {
    printf ( "\n" );
    printf ( "LEGENDRE_COMPUTE - Fatal error!\n" );
    printf ( "  Illegal value of ORDER = %d\n", order );
    exit ( 1 );
  }

  e1 = ( double ) ( order * ( order + 1 ) );

  m = ( order + 1 ) / 2;

  for ( i = 1; i <= m; i++ )
  {
    mp1mi = m + 1 - i;

    t = ( double ) ( 4 * i - 1 ) * r8_pi / ( double ) ( 4 * order + 2 );

    x0 = cos ( t ) * ( 1.0 - ( 1.0 - 1.0 / ( double ) ( order ) ) 
      / ( double ) ( 8 * order * order ) );

    pkm1 = 1.0;
    pk = x0;

    for ( k = 2; k <= order; k++ )
    {
      pkp1 = 2.0 * x0 * pk - pkm1 - ( x0 * pk - pkm1 ) / ( double ) ( k );
      pkm1 = pk;
      pk = pkp1;
    }

    d1 = ( double ) ( order ) * ( pkm1 - x0 * pk );

    dpn = d1 / ( 1.0 - x0 * x0 );

    d2pn = ( 2.0 * x0 * dpn - e1 * pk ) / ( 1.0 - x0 * x0 );

    d3pn = ( 4.0 * x0 * d2pn + ( 2.0 - e1 ) * dpn ) / ( 1.0 - x0 * x0 );

    d4pn = ( 6.0 * x0 * d3pn + ( 6.0 - e1 ) * d2pn ) / ( 1.0 - x0 * x0 );

    u = pk / dpn;
    v = d2pn / dpn;
/*
  Initial approximation H:
*/
    h = -u * ( 1.0 + 0.5 * u * ( v + u * ( v * v - d3pn / ( 3.0 * dpn ) ) ) );
/*
  Refine H using one step of Newton's method:
*/
    p = pk + h * ( dpn + 0.5 * h * ( d2pn + h / 3.0 
      * ( d3pn + 0.25 * h * d4pn ) ) );

    dp = dpn + h * ( d2pn + 0.5 * h * ( d3pn + h * d4pn / 3.0 ) );

    h = h - p / dp;

    xtemp = x0 + h;

    xtab[mp1mi-1] = xtemp;

    fx = d1 - h * e1 * ( pk + 0.5 * h * ( dpn + h / 3.0 
      * ( d2pn + 0.25 * h * ( d3pn + 0.2 * h * d4pn ) ) ) );

    weight[mp1mi-1] = 2.0 * ( 1.0 - xtemp * xtemp ) / ( fx * fx );
  }

  if ( ( order % 2 ) == 1 )
  {
    xtab[0] = 0.0;
  }
/*
  Shift the data up.
*/
  nmove = ( order + 1 ) / 2;
  ncopy = order - nmove;

  for ( i = 1; i <= nmove; i++ )
  {
    iback = order + 1 - i;
    xtab[iback-1] = xtab[iback-ncopy-1];
    weight[iback-1] = weight[iback-ncopy-1];
  }
/*
  Reflect values for the negative abscissas.
*/
  for ( i = 1; i <= order - nmove; i++ )
  {
    xtab[i-1] = - xtab[order-i];
    weight[i-1] = weight[order-i];
  }

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

double p00_alpha ( int problem )

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

    P00_ALPHA returns the value of ALPHA for any problem.

  Discussion:

    ALPHA is the lower, finite limit of integration in the integral.

    The typical or default value is 0.0.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 December 2011

  Author:

    John Burkardt

  Parameters:

    Input, int PROBLEM, the index of the problem.

    Output, double P00_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  if ( problem == 1 )
  {
    alpha = p01_alpha ( );
  }
  else if ( problem == 2 )
  {
    alpha = p02_alpha ( );
  }
  else if ( problem == 3 )
  {
    alpha = p03_alpha ( );
  }
  else if ( problem == 4 )
  {
    alpha = p04_alpha ( );
  }
  else if ( problem == 5 )
  {
    alpha = p05_alpha ( );
  }
  else if ( problem == 6 )
  {
    alpha = p06_alpha ( );
  }
  else if ( problem == 7 )
  {
    alpha = p07_alpha ( );
  }
  else if ( problem == 8 )
  {
    alpha = p08_alpha ( );
  }
  else if ( problem == 9 )
  {
    alpha = p09_alpha ( );
  }
  else if ( problem == 10 )
  {
    alpha = p10_alpha ( );
  }
  else if ( problem == 11 )
  {
    alpha = p11_alpha ( );
  }
  else if ( problem == 12 )
  {
    alpha = p12_alpha ( );
  }
  else if ( problem == 13 )
  {
    alpha = p13_alpha ( );
  }
  else if ( problem == 14 )
  {
    alpha = p14_alpha ( );
  }
  else if ( problem == 15 )
  {
    alpha = p15_alpha ( );
  }
  else if ( problem == 16 )
  {
    alpha = p16_alpha ( );
  }
  else if ( problem == 17 )
  {
    alpha = p17_alpha ( );
  }
  else if ( problem == 18 )
  {
    alpha = p18_alpha ( );
  }
  else if ( problem == 19 )
  {
    alpha = p19_alpha ( );
  }
  else if ( problem == 20 )
  {
    alpha = p20_alpha ( );
  }
  else
  {
    printf ( "\n" );
    printf ( "P00_ALPHA - Fatal error!\n" );
    printf ( "  Illegal problem number = %d\n", problem );
    exit ( 1 );
  }

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

double p00_exact ( int problem )

/******************************************************************************/
/*
  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 index (PROBLEM) rather than by name.

    In most cases, the "exact" value of the integral is not given;
    instead a "respectable" approximation is available.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 December 2011

  Author:

    John Burkardt

  Parameters:

    Input, int PROBLEM, the index of the problem.

    Output, double P00_EXACT, the exact value of the integral.
*/
{
  double exact;

  if ( problem == 1 )
  {
    exact = p01_exact ( );
  }
  else if ( problem == 2 )
  {
    exact = p02_exact ( );
  }
  else if ( problem == 3 )
  {
    exact = p03_exact ( );
  }
  else if ( problem == 4 )
  {
    exact = p04_exact ( );
  }
  else if ( problem == 5 )
  {
    exact = p05_exact ( );
  }
  else if ( problem == 6 )
  {
    exact = p06_exact ( );
  }
  else if ( problem == 7 )
  {
    exact = p07_exact ( );
  }
  else if ( problem == 8 )
  {
    exact = p08_exact ( );
  }
  else if ( problem == 9 )
  {
    exact = p09_exact ( );
  }
  else if ( problem == 10 )
  {
    exact = p10_exact ( );
  }
  else if ( problem == 11 )
  {
    exact = p11_exact ( );
  }
  else if ( problem == 12 )
  {
    exact = p12_exact ( );
  }
  else if ( problem == 13 )
  {
    exact = p13_exact ( );
  }
  else if ( problem == 14 )
  {
    exact = p14_exact ( );
  }
  else if ( problem == 15 )
  {
    exact = p15_exact ( );
  }
  else if ( problem == 16 )
  {
    exact = p16_exact ( );
  }
  else if ( problem == 17 )
  {
    exact = p17_exact ( );
  }
  else if ( problem == 18 )
  {
    exact = p18_exact ( );
  }
  else if ( problem == 19 )
  {
    exact = p19_exact ( );
  }
  else if ( problem == 20 )
  {
    exact = p20_exact ( );
  }
  else
  {
    printf ( "\n" );
    printf ( "P00_EXACT - Fatal error!\n" );
    printf ( "  Illegal problem number = %d\n", problem );
    exit ( 1 );
  }

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

double p00_exp_transform ( int problem, int order )

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

    P00_EXP_TRANSFORM applies an exponential transform and Gauss-Legendre rule.

  Discussion:

    To approximate:

      Integral ( alpha <= x < +oo ) f(x) dx

    Transform:

      u = exp ( -x )
      du = - exp ( -x ) dx

      x = - log ( u )
      dx = - du / u

      x = alpha    => u = exp ( -alpha )
      x = Infinity => u = 0

    Transformed integral:

      Integral ( 0 <= u <= exp ( -alpha ) ) f ( -log(u) ) du / u

    We apply a Gauss-Legendre rule here, but we could easily use any rule
    that avoids evaluation at U = 0.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Reference:

    Arthur Stroud, Don Secrest,
    Gaussian Quadrature Formulas,
    Prentice Hall, 1966,
    LC: QA299.4G3S7.

  Parameters:

    Input, int PROBLEM, the index of the problem.

    Input, int ORDER, the order of the Gauss-Legendre rule 
    to apply.

    Output, double P00_EXP_TRANSFORM, the approximate integral.
*/
{
  double alpha;
  double *fu;
  int i;
  double result;
  double *u;
  double *u_log;
  double *weight;

  u = ( double * ) malloc ( order * sizeof ( double ) );
  u_log = ( double * ) malloc ( order * sizeof ( double ) );
  weight = ( double * ) malloc ( order * sizeof ( double ) );

  alpha = p00_alpha ( problem );
/*
  Get the abscissas and weights for Gauss-Legendre quadrature.
*/
  legendre_compute ( order, u, weight );
/*
  Modify the weights from [-1,1] to [0,exp(-alpha)].
*/
  for ( i = 0; i < order; i++ )
  {
    weight[i] = exp ( -alpha ) * weight[i] / 2.0;
  }
/*
  Linear transform of abscissas from [-1,1] to [0,exp(-alpha)].
*/
  for ( i = 0; i < order; i++ )
  {
    u[i] = ( ( 1.0 + u[i] ) * exp ( - alpha )
           + ( 1.0 - u[i] ) * 0.0 )
           / ( 2.0        );
  }
/*
  Define U_LOG = - log ( U )
*/
  for ( i = 0; i < order; i++ )
  {
    u_log[i] = - log ( u[i] );
  }
/*
  Evaluate F ( -LOG(U) ).
*/
  fu = p00_fun ( problem, order, u_log );
/*
  The integrand is F ( -LOG(U) ) / U
*/
  for ( i = 0; i < order; i++ )
  {
    fu[i] = fu[i] / u[i];
  }
/*
  Sum.
*/
  result = r8vec_dot_product ( order, weight, fu );

  free ( fu );
  free ( u );
  free ( u_log );
  free ( weight );

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

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

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

    P00_FUN evaluates the integrand for any problem.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 December 2011

  Author:

    John Burkardt

  Parameters:

    Input, int PROBLEM, the index of the problem.

    Input, int N, the number of points.

    Input, double X[N], the evaluation points.

    Output, double P00_FUN[N], the function values.
*/
{
  double *f;

  if ( problem == 1 )
  {
    f = p01_fun ( n, x );
  }
  else if ( problem == 2 )
  {
    f = p02_fun ( n, x );
  }
  else if ( problem == 3 )
  {
    f = p03_fun ( n, x );
  }
  else if ( problem == 4 )
  {
    f = p04_fun ( n, x );
  }
  else if ( problem == 5 )
  {
    f = p05_fun ( n, x );
  }
  else if ( problem == 6 )
  {
    f = p06_fun ( n, x );
  }
  else if ( problem == 7 )
  {
    f = p07_fun ( n, x );
  }
  else if ( problem == 8 )
  {
    f = p08_fun ( n, x );
  }
  else if ( problem == 9 )
  {
    f = p09_fun ( n, x );
  }
  else if ( problem == 10 )
  {
    f = p10_fun ( n, x );
  }
  else if ( problem == 11 )
  {
    f = p11_fun ( n, x );
  }
  else if ( problem == 12 )
  {
    f = p12_fun ( n, x );
  }
  else if ( problem == 13 )
  {
    f = p13_fun ( n, x );
  }
  else if ( problem == 14 )
  {
    f = p14_fun ( n, x );
  }
  else if ( problem == 15 )
  {
    f = p15_fun ( n, x );
  }
  else if ( problem == 16 )
  {
    f = p16_fun ( n, x );
  }
  else if ( problem == 17 )
  {
    f = p17_fun ( n, x );
  }
  else if ( problem == 18 )
  {
    f = p18_fun ( n, x );
  }
  else if ( problem == 19 )
  {
    f = p19_fun ( n, x );
  }
  else if ( problem == 20 )
  {
    f = p20_fun ( n, x );
  }
  else
  {
    printf ( "\n" );
    printf ( "P00_FUN - Fatal error!\n" );
    printf ( "  Illegal problem number = %d\n", problem );
    exit ( 1 );
  }

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

double p00_gauss_laguerre ( int problem, int order )

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

    P00_GAUSS_LAGUERRE applies a Gauss-Laguerre rule.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Input, int PROBLEM, the index of the problem.

    Input, int ORDER, the order of the Gauss-Laguerre rule 
    to apply.

    Output, double P00_GAUSS_LAGUERRE, the approximate integral.
*/
{
  double alpha;
  double alpha2;
  double *fx;
  int i;
  double result;
  double *weight;
  double *xtab;

  weight = ( double * ) malloc ( order * sizeof ( double ) );
  xtab = ( double * ) malloc ( order * sizeof ( double ) );

  alpha = p00_alpha ( problem );

  alpha2 = 0.0;
  laguerre_compute ( order, xtab, weight, alpha2 );

  for ( i = 0; i < order; i++ )
  {
    xtab[i] = xtab[i] + alpha;
  }

   fx = p00_fun ( problem, order, xtab );
/*
  The Gauss-Laguerre rule assumes a weight of EXP(-X).

  We need to multiply each F(X) by EXP(X) to implicitly 
  adjust for this weight.
*/
  for ( i = 0; i < order; i++ )
  {
    fx[i] = fx[i] * exp ( xtab[i] );
  }

  result = exp ( -alpha ) * r8vec_dot_product ( order, weight, fx );

  free ( fx );
  free ( weight );
  free ( xtab );

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

int p00_problem_num ( )

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

    P00_PROBLEM_NUM returns the number of test integration problems.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 December 2011

  Author:

    John Burkardt

  Parameters:

    Output, int P00_PROBLEM_NUM, the number of test problems.
*/
{
  int problem_num;

  problem_num = 20;

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

double p00_rat_transform ( int problem, int order )

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

    P00_RAT_TRANSFORM applies a rational transform and Gauss-Legendre rule.

  Discussion:

    To approximate:

      Integral ( alpha <= x < +oo ) f(x) dx

    Transform:

      u = 1 / ( 1 + x )
      du = - dx / ( 1 + x )^2

      x = ( 1 - u ) / u
      dx = - du / u^2

      x = alpha    => u = 1 / ( 1 + alpha )
      x = Infinity => u = 0

    Transformed integral:

      Integral ( 0 < u <= 1 / ( 1 + alpha ) ) f ( ( 1 - u ) / u ) du / u^2

    We apply a Gauss-Legendre rule here, but we could easily use any rule
    that avoids evaluation at U = 0.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    30 July 2007

  Author:

    John Burkardt

  Reference:

    Arthur Stroud, Don Secrest,
    Gaussian Quadrature Formulas,
    Prentice Hall, 1966,
    LC: QA299.4G3S7.

  Parameters:

    Input, int PROBLEM, the index of the problem.

    Input, int ORDER, the order of the Gauss-Legendre rule 
    to apply.

    Output, double P00_RAT_TRANSFORM, the approximate integral.
*/
{
  double alpha;
  double *fu;
  int i;
  double result;
  double *u;
  double *u_rat;
  double *weight;

  u = ( double * ) malloc ( order * sizeof ( double ) );
  u_rat = ( double * ) malloc ( order * sizeof ( double ) );
  weight = ( double * ) malloc ( order * sizeof ( double ) );

  alpha = p00_alpha ( problem );
/*
  Get the abscissas and weights for Gauss-Legendre quadrature.
*/
  legendre_compute ( order, u, weight );
/*
  Modify the weights from [-1,1] to [0,1/(1+alpha)].
*/
  for ( i = 0; i < order; i++ )
  {
    weight[i] = weight[i] / 2.0 / ( 1.0 + alpha );
  }
/*
  Linear transform of abscissas from [-1,1] to [0,exp(-alpha)].
*/
  for ( i = 0; i < order; i++ )
  {
    u[i] = ( ( 1.0 + u[i] ) / ( 1.0 + alpha )
           + ( 1.0 - u[i] ) * 0.0 )
           / ( 2.0        );
  }
/*
  Define U_RAT = ( 1 - U ) / U.
*/
  for ( i = 0; i < order; i++ )
  {
    u_rat[i] = ( 1.0 - u[i] ) / u[i];
  }
/*
  Evaluate F ( ( 1 - U ) / U ).
*/
  fu = p00_fun ( problem, order, u_rat );
/*
  The integrand is F ( ( 1 - U ) / U ) / U^2
*/
  for ( i = 0; i < order; i++ )
  {
    fu[i] = fu[i] / u[i] / u[i];
  }
/*
  Sum.
*/
  result = r8vec_dot_product ( order, weight, fu );

  free ( fu );
  free ( u );
  free ( u_rat );
  free ( weight );

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

char *p00_title ( int problem )

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

    P00_TITLE returns the title for any problem.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 December 2011

  Author:

    John Burkardt

  Parameters:

    Input, int PROBLEM, the index of the problem.

    Output, char *P00_TITLE, the title of the problem.
*/
{
  char *title;

  if ( problem == 1 )
  {
    title = p01_title ( );
  }
  else if ( problem == 2 )
  {
    title = p02_title ( );
  }
  else if ( problem == 3 )
  {
    title = p03_title ( );
  }
  else if ( problem == 4 )
  {
    title = p04_title ( );
  }
  else if ( problem == 5 )
  {
    title = p05_title ( );
  }
  else if ( problem == 6 )
  {
    title = p06_title ( );
  }
  else if ( problem == 7 )
  {
    title = p07_title ( );
  }
  else if ( problem == 8 )
  {
    title = p08_title ( );
  }
  else if ( problem == 9 )
  {
    title = p09_title ( );
  }
  else if ( problem == 10 )
  {
    title = p10_title ( );
  }
  else if ( problem == 11 )
  {
    title = p11_title ( );
  }
  else if ( problem == 12 )
  {
    title = p12_title ( );
  }
  else if ( problem == 13 )
  {
    title = p13_title ( );
  }
  else if ( problem == 14 )
  {
    title = p14_title ( );
  }
  else if ( problem == 15 )
  {
    title = p15_title ( );
  }
  else if ( problem == 16 )
  {
    title = p16_title ( );
  }
  else if ( problem == 17 )
  {
    title = p17_title ( );
  }
  else if ( problem == 18 )
  {
    title = p18_title ( );
  }
  else if ( problem == 19 )
  {
    title = p19_title ( );
  }
  else if ( problem == 20 )
  {
    title = p20_title ( );
  }
  else
  {
    printf ( "\n" );
    printf ( "P00_TITLE - Fatal error!\n" );
    printf ( "  Illegal problem number = %d\n", problem );
    exit ( 1 );
  }

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

double p01_alpha ( )

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

    P01_ALPHA returns ALPHA for problem 1.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P01_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 2.0;

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

double p01_exact ( )

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

    P01_EXACT returns the exact integral for problem 1.

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P01_EXACT, the value of the integral.
*/
{
  double exact;

  exact = 0.19524754198276439152;

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

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

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

    P01_FUN evaluates the integrand for problem 1.

  Discussion:

    D&R gives "exact" value as 0.19524753.
    Mathematica returns        0.19524754198276439152...
    D&R gives Laguerre(16) as  0.16623627...

  Integral:

    exp ( -2 ) Integral ( 2 <= x < +oo ) / ( x * log(x)^2 ) dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Reference:

    Philip Davis, Philip Rabinowitz,
    Methods of Numerical Integration,
    Second Edition,
    Dover, 2007,
    ISBN: 0486453391,
    LC: QA299.3.D28.

  Parameters:

    Input, int N, the number of points.

    Input, double X[N], the evaluation points.

    Output, double P01_FUN[N], the function values.
*/
{
  double *f;
  int i;

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

  for ( i = 0; i < n; i++ )
  {
    f[i] = exp ( -2.0 ) / ( x[i] * pow ( log ( x[i] ), 2 ) );
  }

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

char *p01_title ( )

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

    P01_TITLE returns the title for problem 1.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, char *P01_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "1 / ( x * log ( x )^2 )" );

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

double p02_alpha ( )

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

    P02_ALPHA returns ALPHA for problem 2.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P02_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 2.0;

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

double p02_exact ( )

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

    P02_EXACT returns the exact integral for problem 2.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P02_EXACT, the value of the integral.
*/
{
  double exact;

  exact = 0.32510848278991335198;

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

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

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

    P02_FUN evaluates the integrand for problem 2.

  Discussion:

    D&R gives "exact" value as 0.32510855.
    Mathematica returns        0.32510848278991335198...
    D&R gives Laguerre(16) as  0.19142399...

  Integral:

    exp ( -2 ) Integral ( 2 <= x < +oo ) / ( x * log(x)^(3/2) ) dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Reference:

    Philip Davis, Philip Rabinowitz,
    Methods of Numerical Integration,
    Second Edition,
    Dover, 2007,
    ISBN: 0486453391,
    LC: QA299.3.D28.

  Parameters:

    Input, int N, the number of points.

    Input, double X[N], the evaluation points.

    Output, double P02_FUN[N], the function values.
*/
{
  double *f;
  int i;

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

  for ( i = 0; i < n; i++ )
  {
    f[i] = exp ( -2.0 ) / ( x[i] * pow ( log ( x[i] ), 1.5 ) );
  }

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

char *p02_title ( )

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

    P02_TITLE returns the title for problem 2.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, char *P02_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "1 / ( x * log ( x )^(3/2) )" );

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

double p03_alpha ( )

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

    P03_ALPHA returns ALPHA for problem 3.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P03_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 2.0;

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

double p03_exact ( )

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

    P03_EXACT returns the exact integral for problem 3.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P03_EXACT, the value of the integral.
*/
{
  double exact;

  exact = 13.628;

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

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

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

    P03_FUN evaluates the integrand for problem 3.

  Discussion:

    D&R gives "exact" value as 13.628...
    Mathematica returns        13.440045415012575106...
    D&R gives Laguerre(16) as   0.44996932...

  Integral:

    exp ( -2 ) Integral ( 2 <= x < +oo ) / ( x^1.01 ) dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Reference:

    Philip Davis, Philip Rabinowitz,
    Methods of Numerical Integration,
    Second Edition,
    Dover, 2007,
    ISBN: 0486453391,
    LC: QA299.3.D28.

  Parameters:

    Input, int N, the number of points.

    Input, double X[N], the evaluation points.

    Output, double P03_FUN[N], the function values.
*/
{
  double *f;
  int i;

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

  for ( i = 0; i < n; i++ )
  {
    f[i] = exp ( -2.0 ) * 1.0 / pow ( x[i], 1.01 );
  }

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

char *p03_title ( )

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

    P03_TITLE returns the title for problem 3.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, char *P03_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "1 / ( x^1.01 )" );

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

double p04_alpha ( )

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

    P04_ALPHA returns ALPHA for problem 4.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P04_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 2.0;

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

double p04_exact ( )

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

    P04_EXACT returns the estimated integral for problem 4.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P04_EXACT, the estimated value of the integral.
*/
{
  double exact;

  exact = -0.0046848541335080643181;

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

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

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

    P04_FUN evaluates the integrand for problem 4.

  Discussion:

    D&R gives "exact" value as -0.0046984...
    Mathematica returns        -0.0046848541335080643181...
    D&R gives Laguerre(16) as  -0.039258696...

  Integral:

    exp ( -2 ) Integral ( 2 <= x < +oo ) sin ( x ) / x dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Reference:

    Philip Davis, Philip Rabinowitz,
    Methods of Numerical Integration,
    Second Edition,
    Dover, 2007,
    ISBN: 0486453391,
    LC: QA299.3.D28.

  Parameters:

    Input, int N, the number of points.

    Input, double X[N], the evaluation points.

    Output, double P04_FUN[N], the function values.
*/
{
  double *f;
  int i;

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

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

char *p04_title ( )

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

    P04_TITLE returns the title for problem 4.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, char *P04_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "Sine integral" );

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

double p05_alpha ( )

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

    P05_ALPHA returns ALPHA for problem 5.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P05_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 2.0;

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

double p05_exact ( )

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

    P05_EXACT returns the estimated integral for problem 5.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    06 October 2006

  Author:

    John Burkardt

  Parameters:

    Output, double P05_EXACT, the estimated value of the integral.
*/
{
  double exact;

  exact = 0.0015897286158592328774;

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

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

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

    P05_FUN evaluates the integrand for problem 5.

  Discussion:

    D&R gives "exact" value as  0.00158973...
    Mathematica returns         0.0015897286158592328774...
    D&R gives Laguerre(16) as  -0.067859545...

  Integral:

    exp ( -2 ) Integral ( 2 <= x < +oo ) cos ( pi * x^2 / 2 ) dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Reference:

    Philip Davis, Philip Rabinowitz,
    Methods of Numerical Integration,
    Second Edition,
    Dover, 2007,
    ISBN: 0486453391,
    LC: QA299.3.D28.

  Parameters:

    Input, int N, the number of points.

    Input, double X[N], the evaluation points.

    Output, double P05_FUN[N], the function values.
*/
{
  double *f;
  int i;
  const double r8_pi = 3.1415926535897932385;

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

  for ( i = 0; i < n; i++ )
  {
    f[i] = exp ( -2.0 ) * cos ( 0.5 * r8_pi * x[i] * x[i] );
  }

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

char *p05_title ( )

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

    P05_TITLE returns the title for problem 5.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, char *P05_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "Fresnel integral" );

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

double p06_alpha ( )

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

    P06_ALPHA returns ALPHA for problem 6.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P06_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 2.0;

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

double p06_exact ( )

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

    P06_EXACT returns the exact integral for problem 6.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P06_EXACT, the estimated value of the integral.
*/
{
  double exact;

  exact = 0.00056103711148387120640;

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

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

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

    P06_FUN evaluates the integrand for problem 6.

  Discussion:

    D&R gives "exact" value as 0.0005610371...
    Mathematica returns        0.00056103711148387120640...
    D&R gives Laguerre(16) as  0.00056100775...

  Integral:

    exp ( -2 ) Integral ( 2 <= x < +oo ) exp ( -x^2 ) dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    30 July 2007

  Author:

    John Burkardt

  Reference:

    Philip Davis, Philip Rabinowitz,
    Methods of Numerical Integration,
    Second Edition,
    Dover, 2007,
    ISBN: 0486453391,
    LC: QA299.3.D28.

  Parameters:

    Input, int N, the number of points.

    Input, double X[N], the evaluation points.

    Output, double P06_FUN[N], the function values.
*/
{
  double exponent_min = -80.0;
  double *f;
  int i;

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

  for ( i = 0; i < n; i++ )
  {
    if ( - x[i] * x[i] < exponent_min )
    {
      f[i] = 0.0;
    }
    else
    {
      f[i] = exp ( -2.0 ) * exp ( - x[i] * x[i] );
    }
  }

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

char *p06_title ( )

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

    P06_TITLE returns the title for problem 6.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, char *P06_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "Complementary error function" );

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

double p07_alpha ( )

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

    P07_ALPHA returns ALPHA for problem 7.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P07_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 2.0;

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

double p07_exact ( )

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

    P07_EXACT returns the exact integral for problem 7.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P07_EXACT, the value of the integral.
*/
{
  double exact;

  exact = 0.16266891;

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

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

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

    P07_FUN evaluates the integrand for problem 7.

  Discussion:

    D&R gives "exact" value as 0.16266891...
    Mathematica does not return a value.
    D&R gives Laguerre(16) as  0.097083064...

  Integral:

    exp ( -2 ) Integral ( 2 <= x < +oo ) sin ( x - 1 ) 
      / sqrt ( x * ( x - 2 ) ) dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Reference:

    Philip Davis, Philip Rabinowitz,
    Methods of Numerical Integration,
    Second Edition,
    Dover, 2007,
    ISBN: 0486453391,
    LC: QA299.3.D28.

  Parameters:

    Input, int N, the number of points.

    Input, double X[N], the evaluation points.

    Output, double P07_FUN[N], the function values.
*/
{
  double *f;
  int i;

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

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

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

char *p07_title ( )

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

    P07_TITLE returns the title for problem 7.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, char *P07_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "Bessel function" );

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

double p08_alpha ( )

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

    P08_ALPHA returns ALPHA for problem 8.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 May 2014

  Author:

    John Burkardt

  Parameters:

    Output, double P08_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 0.0;

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

double p08_exact ( )

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

    P08_EXACT returns the estimated integral for problem 8.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    18 May 2014

  Author:

    John Burkardt

  Parameters:

    Output, double P08_EXACT, the estimated value of the integral.
*/
{
  double exact;
  const double r8_pi = 3.1415926535897932385;

  exact = r8_pi * r8_pi / 6.0;

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

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

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

    P08_FUN evaluates the integrand for problem 8.

  Integral:

    Integral ( 0 <= x < +oo ) x / ( exp ( x ) - 1 ) dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    30 July 2007

  Author:

    John Burkardt

  Reference:

    Philip Davis, Philip Rabinowitz,
    Methods of Numerical Integration,
    Second Edition,
    Dover, 2007,
    ISBN: 0486453391,
    LC: QA299.3.D28.

  Parameters:

    Input, int N, the number of points.

    Input, double X[N], the evaluation points.

    Output, double P08_FUN[N], the function values.
*/
{
  double exponent_max = 80.0;
  double *f;
  int i;

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

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

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

char *p08_title ( )

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

    P08_TITLE returns the title for problem 8.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, char *P08_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "Debye function" );

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

double p09_alpha ( )

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

    P09_ALPHA returns ALPHA for problem 9.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P09_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 0.0;

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

double p09_exact ( )

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

    P09_EXACT returns the estimated integral for problem 9.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P09_EXACT, the estimated value of the integral.
*/
{
  double exact;

  exact = 24.0;

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

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

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

    P09_FUN evaluates the integrand for problem 9.

  Discussion:

    The integral is the definition of the Gamma function for
    Z = 5, with exact value (Z-1)! = 24.

  Integral:

    Integral ( 0 <= x < +oo ) x^4 exp ( -x ) dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    30 July 2007

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of points.

    Input, double X[N], the evaluation points.

    Output, double P09_FUN[N], the function values.
*/
{
  double exponent_min = -80.0;
  double *f;
  int i;

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

  for ( i = 0; i < n; i++ )
  {
    if ( -x[i] < exponent_min )
    {
      f[i] = 0.0;
    }
    else
    {
      f[i] = pow ( x[i], 4 ) * exp ( -x[i] );
    }
  }

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

char *p09_title ( )

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

    P09_TITLE returns the title for problem 9.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    28 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, char *P09_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "Gamma(Z=5) function" );

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

double p10_alpha ( )

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

    P10_ALPHA returns ALPHA for problem 10.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P10_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 0.0;

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

double p10_exact ( )

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

    P10_EXACT returns the estimated integral for problem 10.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double EXACT, the estimated value of the integral.
*/
{
  double exact;
  const double r8_pi = 3.1415926535897932385;

  exact = r8_pi / 2.0;

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

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

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

    P10_FUN evaluates the integrand for problem 10.

  Discussion:

    S&S gives exact value as pi/2 = 1.5707963267948966192...
    S&S gives Laguerre(16) as       1.5537377347...
    S&S gives EXP_TRANSFORM(16) as  1.4293043007...

  Integral:

    Integral ( 0 <= x < +oo ) 1/(1+x*x) dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Reference:

    Arthur Stroud, Don Secrest,
    Gaussian Quadrature Formulas,
    Prentice Hall, 1966,
    LC: QA299.4G3S7.

  Parameters:

    Input, int N, the number of points.

    Input, double X[N], the evaluation points.

    Output, double P10_FUN[N], the function values.
*/
{
  double *f;
  int i;

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

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

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

char *p10_title ( )

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

    P10_TITLE returns the title for problem 10.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, char *P10_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "1 / ( 1 + x*x )" );

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

double p11_alpha ( )

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

    P11_ALPHA returns ALPHA for problem 11.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P11_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 0.0;

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

double p11_exact ( )

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

    P11_EXACT returns the estimated integral for problem 11.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double EXACT, the estimated value of the integral.
*/
{
  double exact;
  const double r8_pi = 3.1415926535897932385;

  exact = r8_pi;

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

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

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

    P11_FUN evaluates the integrand for problem 11.

  Discussion:

    S&S gives exact value as pi =  3.1415926535897932385...
    S&S gives Laguerre(16) as      2.6652685196...
    S&S gives EXP_TRANSFORM(16) as 2.3629036166... 

  Integral:

    Integral ( 0 <= x < +oo ) 1/((1+x)*sqrt(x)) dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Reference:

    Arthur Stroud, Don Secrest,
    Gaussian Quadrature Formulas,
    Prentice Hall, 1966,
    LC: QA299.4G3S7.

  Parameters:

    Input, int N, the number of points.

    Input, double X[N], the evaluation points.

    Output, double P11_FUN[N], the function values.
*/
{
  double *f;
  int i;

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

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

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

char *p11_title ( )

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

    P11_TITLE returns the title for problem 11.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, char *P11_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "1 / ( (1+x) * sqrt(x) )" );

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

double p12_alpha ( )

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

    P12_ALPHA returns ALPHA for problem 12.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P12_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 0.0;

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

double p12_exact ( )

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

    P12_EXACT returns the estimated integral for problem 12.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double EXACT, the estimated value of the integral.
*/
{
  double exact;

  exact = 0.5;

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

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

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

    P12_FUN evaluates the integrand for problem 12.

  Discussion:

    S&S gives exact value as pi =  0.5
    S&S gives Laguerre(16) as      0.5000000000...
    S&S gives EXP_TRANSFORM(16) as 0.5019065783... 

  Integral:

    Integral ( 0 <= x < +oo ) exp ( -x ) * cos ( x ) dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    30 July 2007

  Author:

    John Burkardt

  Reference:

    Arthur Stroud, Don Secrest,
    Gaussian Quadrature Formulas,
    Prentice Hall, 1966,
    LC: QA299.4G3S7.

  Parameters:

    Input, int N, the number of points.

    Input, double X[N], the evaluation points.

    Output, double P12_FUN[N], the function values.
*/
{
  double exponent_min = -80.0;
  double *f;
  int i;

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

  for ( i = 0; i < n; i++ )
  {
    if ( - x[i] < exponent_min )
    {
      f[i] = 0.0;
    }
    else
    {
      f[i] = exp ( -x[i] ) * cos ( x[i] );
    }
  }

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

char *p12_title ( )

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

    P12_TITLE returns the title for problem 12.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, char *P12_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "exp ( - x ) * cos ( x )" );

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

double p13_alpha ( )

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

    P13_ALPHA returns ALPHA for problem 13.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P13_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 0.0;

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

double p13_exact ( )

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

    P13_EXACT returns the estimated integral for problem 13.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double EXACT, the estimated value of the integral.
*/
{
  double exact;
  const double r8_pi = 3.1415926535897932385;

  exact = r8_pi / 2.0;

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

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

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

    P13_FUN evaluates the integrand for problem 13.

  Discussion:

    S&S gives exact value as pi/2 = 1.5707963267948966192...
    S&S gives Laguerre(16) as       1.4399523793...
    S&S gives EXP_TRANSFORM(16) as  1.3045186595...

  Integral:

    Integral ( 0 <= x < +oo ) sin ( x ) / x dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Reference:

    Arthur Stroud, Don Secrest,
    Gaussian Quadrature Formulas,
    Prentice Hall, 1966,
    LC: QA299.4G3S7.

  Parameters:

    Input, int N, the number of points.

    Input, double X[N], the evaluation points.

    Output, double P13_FUN[N], the function values.
*/
{
  double *f;
  int i;

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

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

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

char *p13_title ( )

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

    P13_TITLE returns the title for problem 13.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, char *P13_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "sin(x) / x" );

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

double p14_alpha ( )

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

    P14_ALPHA returns ALPHA for problem 14.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P14_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 0.0;

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

double p14_exact ( )

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

    P14_EXACT returns the estimated integral for problem 14.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, double EXACT, the estimated value of the integral.
*/
{
  double exact;

  exact = 1.0634618101722400407;

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

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

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

    P14_FUN evaluates the integrand for problem 14.

  Discussion:

    S&S gives "exact" value as     1.0634618101...
    Mathematica returns            1.0634618101722400407...
    S&S gives Laguerre(16) as      1.0634713425...
    S&S gives EXP_TRANSFORM(16) as 1.0634618101...

  Integral:

    Integral ( 0 <= x < +oo ) sin ( exp ( - x ) + exp ( - 4 x ) ) dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Reference:

    Arthur Stroud, Don Secrest,
    Gaussian Quadrature Formulas,
    Prentice Hall, 1966,
    LC: QA299.4G3S7.

  Parameters:

    Input, int N, the number of points.

    Input, double X[N], the evaluation points.

    Output, double P14_FUN[N], the function values.
*/
{
  double exponent_min = -80.0;
  double *f;
  int i;

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

  for ( i = 0; i < n; i++ )
  {
    if ( - x[i] < exponent_min )
    {
      f[i] = 0.0;
    }
    else if ( -4.0 * x[i] < exponent_min )
    {
      f[i] = sin ( exp ( -x[i] ) );
    }
    else
    {
      f[i] = sin ( exp ( -x[i] ) + exp ( -4.0 * x[i] ) );
    }
  }

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

char *p14_title ( )

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

    P14_TITLE returns the title for problem 14.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 July 2007

  Author:

    John Burkardt

  Parameters:

    Output, char *P14_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "sin ( exp(-x) + exp(-4x) )" );

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

double p15_alpha ( )

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

    P15_ALPHA returns ALPHA for problem 15.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 August 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P15_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 0.0;

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

double p15_exact ( )

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

    P15_EXACT returns the estimated integral for problem 15.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 August 2007

  Author:

    John Burkardt

  Parameters:

    Output, double EXACT, the estimated value of the integral.
*/
{
  double exact;
  const double r8_pi = 3.1415926535897932385;

  exact = - r8_pi * log ( 10.0 ) / 20.0;

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

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

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

    P15_FUN evaluates the integrand for problem 15.

  Integral:

    Integral ( 0 <= x < +oo ) log(x) / (1+100*x*x) dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 August 2007

  Author:

    John Burkardt

  Reference:

    Robert Piessens, Elise deDoncker-Kapenga, 
    Christian Ueberhuber, David Kahaner,
    QUADPACK: A Subroutine Package for Automatic Integration,
    Springer, 1983,
    ISBN: 3540125531,
    LC: QA299.3.Q36.

  Parameters:

    Input, int N, the number of points.

    Input, double X[N], the evaluation points.

    Output, double P15_FUN[N], the function values.
*/
{
  double *f;
  int i;

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

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

char *p15_title ( )

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

    P15_TITLE returns the title for problem 15.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 August 2007

  Author:

    John Burkardt

  Parameters:

    Output, char *P15_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "log(x) / ( 1 + 100 x^2 )" );

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

double p16_alpha ( )

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

    P16_ALPHA returns ALPHA for problem 16.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 August 2007

  Author:

    John Burkardt

  Parameters:

    Output, double P16_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 0.0;

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

double p16_exact ( )

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

    P16_EXACT returns the estimated integral for problem 16.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 August 2007

  Author:

    John Burkardt

  Parameters:

    Output, double EXACT, the estimated value of the integral.
*/
{
  double exact;

  exact = 1.0;

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

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

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

    P16_FUN evaluates the integrand for problem 16.

  Integral:

    Integral ( 0 <= x < +oo ) cos ( pi * x / 2 ) / sqrt ( x ) dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 August 2007

  Author:

    John Burkardt

  Reference:

    Robert Piessens, Elise deDoncker-Kapenga, 
    Christian Ueberhuber, David Kahaner,
    QUADPACK: A Subroutine Package for Automatic Integration,
    Springer, 1983,
    ISBN: 3540125531,
    LC: QA299.3.Q36.

  Parameters:

    Input, int N, the number of points.

    Input, double X[N], the evaluation points.

    Output, double P16_FUN[N], the function values.
*/
{
  double *f;
  int i;
  const double r8_pi = 3.1415926535897932385;

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

  for ( i = 0; i < n; i++ )
  {
    if ( x[i] == 0.0 )
    {
      f[i] = HUGE_VAL;
    }
    else
    {
      f[i] = cos ( r8_pi * x[i] / 2.0 ) / sqrt ( x[i] );
    }
  }
  return f;
}
/******************************************************************************/

char *p16_title ( )

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

    P16_TITLE returns the title for problem 16.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    29 August 2007

  Author:

    John Burkardt

  Parameters:

    Output, char *P16_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "cos ( pi x / 2 ) / sqrt ( x )" );

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

double p17_alpha ( )

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

    P17_ALPHA returns ALPHA for problem 17.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 December 2011

  Author:

    John Burkardt

  Parameters:

    Output, double P17_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 0.0;

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

double p17_exact ( )

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

    P17_EXACT returns the exact integral for problem 17.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 December 2011

  Author:

    John Burkardt

  Parameters:

    Output, double P17_EXACT, the value of the integral.
*/
{
  const double beta = 2.0;
  double exact;
  const double r8_pi = 3.1415926535897932385;

  exact = sqrt ( r8_pi ) * cos ( 0.5 * atan ( pow ( 2.0, beta ) ) ) 
    / sqrt ( sqrt ( 1.0 + pow ( 0.25, beta ) ) );

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

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

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

    P17_FUN evaluates the integrand for problem 17.

  Integral:

    Integral ( 0 <= x < +oo) exp ( - x / 2^beta ) * cos ( x ) / sqrt ( x ) dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 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 84.

  Parameters:

    Input, int N, the number of evaluation points.

    Input, double X[N], the evaluation points.

    Output, double P17_FUN[N], the integrand values.
*/
{
  static double beta = 2.0;
  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] = exp ( - x[i] / pow ( 2.0, beta ) ) * cos ( x[i] ) 
        / sqrt ( x[i] );
    }
  }
  return fx;
}
/******************************************************************************/

char *p17_title ( )

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

    P17_TITLE returns the title for problem 17.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 December 2011

  Author:

    John Burkardt

  Parameters:

    Output, char *P17_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "exp ( - x / 2^beta ) * cos ( x ) / sqrt ( x )" );

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

double p18_alpha ( )

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

    P18_ALPHA returns ALPHA for problem 18.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 December 2011

  Author:

    John Burkardt

  Parameters:

    Output, double P18_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 0.0;

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

double p18_exact ( )

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

    P18_EXACT returns the exact integral for problem 18.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 December 2011

  Author:

    John Burkardt

  Parameters:

    Output, double P18_EXACT, the value of the integral.
*/
{
  static double beta = 1.0;
  double exact;
 
  exact = pow ( 2.0, 3.0 * beta + 1.0 );

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

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

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

    P18_FUN evaluates the integrand for problem 18.

  Integral:

    Integral ( 0 <= x < +oo ) x^2 * exp ( - x / 2^beta ) dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 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 84.

  Parameters:

    Input, int N, the number of evaluation points.

    Input, double X[N], the evaluation points.

    Output, double P18_FUN[N], the integrand values.
*/
{
  static double beta = 1.0;
  double *fx;
  int i;

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

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

char *p18_title ( )

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

    P18_TITLE returns the title for problem 18.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 December 2011

  Author:

    John Burkardt

  Parameters:

    Output, char *P18_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "x^2 * exp ( - x / 2^beta )" );

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

double p19_alpha ( )

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

    P19_ALPHA returns ALPHA for problem 19.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 December 2011

  Author:

    John Burkardt

  Parameters:

    Output, double P19_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 0.0;

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

double p19_exact ( )

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

    P19_EXACT returns the exact integral for problem 19.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 December 2011

  Author:

    John Burkardt

  Parameters:

    Output, double P19_EXACT, the value of the integral.
*/
{
  const double beta = 0.5;
  double exact;
  const double r8_pi = 3.1415926535897932385;

  if ( beta == 1.0 )
  {
    exact = 1.0 / 10.0;
  }
  else
  {
    exact = ( 1.0 - beta ) * r8_pi 
      / ( pow ( 10.0, beta ) * sin ( r8_pi * beta ) );
  }
  return exact;
}
/******************************************************************************/

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

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

    P19_FUN evaluates the integrand for problem 19.

  Integral:

    Integral ( 0 <= x < +oo ) x^(alpha-1) / ( 1 + 10 x )^2 dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 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 84.

  Parameters:

    Input, int N, the number of evaluation points.

    Input, double X[N], the evaluation points.

    Output, double P61_FUN[N], the integrand values.
*/
{
  static double beta = 0.5;
  double *fx;
  int i;

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

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

char *p19_title ( )

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

    P19_TITLE returns the title for problem 19.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 December 2011

  Author:

    John Burkardt

  Parameters:

    Output, char *P19_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "x^(beta-1) / ( 1 + 10 x )^2" );

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

double p20_alpha ( )

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

    P20_ALPHA returns ALPHA for problem 20.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 December 2011

  Author:

    John Burkardt

  Parameters:

    Output, double P20_ALPHA, the value of ALPHA.
*/
{
  double alpha;

  alpha = 0.0;

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

double p20_exact ( )

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

    P20_EXACT returns the exact integral for problem 20.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 December 2011

  Author:

    John Burkardt

  Parameters:

    Output, double P20_EXACT, the value of the integral.
*/
{
  static double beta = 1.0;
  double exact;

  exact = 
    ( 
      log ( 1.5 ) / pow ( 2.0, beta ) 
      - 1.0 / pow ( 2.0, beta + 1.0 ) * 
      log ( ( 16.0 + pow ( 0.25, beta ) ) / ( 1.0 + pow ( 0.25, beta ) ) ) 
      - atan ( pow ( 2.0, beta + 2.0 ) ) - atan ( pow ( 2.0, beta ) ) 
    ) / ( 1.0 + pow ( 0.25, beta ) );

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

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

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

    P20_FUN evaluates the integrand for problem 20.

  Integral:

    Integral ( 0 <= x < +oo ) 
      1 / ( 2^beta * ( ( x - 1 )^2 + (1/4)^beta ) * ( x - 2 ) ) dx

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 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 84.

  Parameters:

    Input, int N, the number of evaluation points.

    Input, double X[N], the evaluation points.

    Output, double P20_FUN[N], the integrand values.
*/
{
  static double beta = 1.0;
  double *fx;
  int i;

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

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

char *p20_title ( )

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

    P20_TITLE returns the title for problem 20.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    27 December 2011

  Author:

    John Burkardt

  Parameters:

    Output, char *P20_TITLE, the title of the problem.
*/
{
  char *title;

  title = ( char * ) malloc ( 80 * sizeof ( char ) );

  strcpy ( title, "1 / ( 2^beta * ( ( x - 1 )^2 + (1/4)^beta ) * ( x - 2 ) )" );

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

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

void timestamp ( )

/******************************************************************************/
/*
  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 MIT 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 );

  fprintf ( stdout, "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}