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

# include "hermite_test_int.h"

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

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

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

    HERMITE_COMPUTE computes a Gauss-Hermite quadrature rule.

  Discussion:

    The abscissas are the zeros of the N-th order Hermite polynomial.

    The integration interval is ( -oo, +oo ).

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

    The integral to approximate:

      Integral ( -oo < X < +oo ) exp ( - X*X ) * F(X) dX

    The quadrature rule:

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

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  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 formula to be computed.

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

    Output, double WEIGHT[ORDER], the Gauss-Hermite weights.
*/
{
  double cc;
  double dp2;
  int i;
  double p1;
  double s;
  double temp;
  double x;

  cc = 1.7724538509 * r8_gamma ( ( double ) ( order ) ) 
    / pow ( 2.0, order - 1 );

  s = pow ( 2.0 * ( double ) ( order ) + 1.0, 1.0 / 6.0 );

  for ( i = 0; i < ( order + 1 ) / 2; i++ )
  {
    if ( i == 0 )
    {
      x = s * s * s - 1.85575 / s;
    }
    else if ( i == 1 )
    {
      x = x - 1.14 * pow ( ( double ) ( order ), 0.426 ) / x;
    }
    else if ( i == 2 )
    {
      x = 1.86 * x - 0.86 * xtab[0];
    }
    else if ( i == 3 )
    {
      x = 1.91 * x - 0.91 * xtab[1];
    }
    else
    {
      x = 2.0 * x - xtab[i-2];
    }

    hermite_root ( &x, order, &dp2, &p1 );

    xtab[i] = x;
    weight[i] = ( cc / dp2 ) / p1;

    xtab[order-i-1] = -x;
    weight[order-i-1] = weight[i];
  }
/*
  Reverse the order of the XTAB values.
*/
  r8vec_reverse ( order, xtab );

  if ( 0 )
  {
  for ( i = 0; i < order/2; i++ )
  {
    temp            = xtab[i];
    xtab[i]         = xtab[order-1-i];
    xtab[order-1-i] = temp;
  }
  }
  return;
}
/******************************************************************************/

double hermite_integral ( int n )

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

    HERMITE_INTEGRAL returns the value of a Hermite polynomial integral.

  Discussion:

    H(n) = Integral ( -oo < x < +oo ) x^n exp(-x*x) dx

    H(n) is 0 for n odd.

    H(n) = (n-1)!! * sqrt(pi) / 2^(n/2) for n even.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    25 July 2007

  Author:

    John Burkardt

  Parameters:

    Input, int N, the order of the integral.  
    0 <= N.

    Output, double VALUE, the value of the integral.
*/
{
  const double r8_pi = 3.141592653589793;
  double value;

  if ( n < 0 )
  {
    value = - r8_huge ( );
  }
  else if ( ( n % 2 ) == 1 )
  {
    value = 0.0;
  }
  else
  {
    value = ( double ) ( i4_factorial2 ( n - 1 ) ) * sqrt ( r8_pi ) 
      / pow ( 2.0, n / 2 );
  }

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

void hermite_recur ( double *p2, double *dp2, double *p1, double x, int order )

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

    HERMITE_RECUR finds the value and derivative of a Hermite polynomial.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    30 April 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 H(ORDER)(X).

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

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

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

    Input, int ORDER, the order of the polynomial to be computed.
*/
{
  int i;
  double dq0;
  double dq1;
  double dq2;
  double q0;
  double q1;
  double q2;

  q1 = 1.0;
  dq1 = 0.0;

  q2 = x;
  dq2 = 1.0;

  for ( i = 2; i <= order; i++ )
  {
    q0 = q1;
    dq0 = dq1;

    q1 = q2;
    dq1 = dq2;

    q2  = x * q1 - 0.5 * ( ( double ) ( i ) - 1.0 ) * q0;
    dq2 = x * dq1 + q1 - 0.5 * ( ( double ) ( i ) - 1.0 ) * dq0;
  }

  *p2 = q2;
  *dp2 = dq2;
  *p1 = q1;

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

void hermite_root ( double *x, int order, double *dp2, double *p1 )

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

    HERMITE_ROOT improves an approximate root of a Hermite polynomial.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  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 Hermite polynomial.

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

    Output, double *P1, the value of H(ORDER-1)(X).
*/
{
  double d;
  double eps = 1.0E-12;
  double p2;
  int step;
  int step_max = 10;

  for ( step = 1; step <= step_max; step++ )
  {
    hermite_recur ( &p2, dp2, p1, *x, order );

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

    if ( r8_abs ( d ) <= eps * ( r8_abs ( *x ) + 1.0 ) )
    {
      return;
    }
  }

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

int i4_factorial2 ( int n )

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

    I4_FACTORIAL2 computes the double factorial function.

  Formula:

    FACTORIAL2( N ) = Product ( N * (N-2) * (N-4) * ... * 2 )  (N even)
                    = Product ( N * (N-2) * (N-4) * ... * 1 )  (N odd)

  Example:

     N    N!!

     0     1
     1     1
     2     2
     3     3
     4     8
     5    15
     6    48
     7   105
     8   384
     9   945
    10  3840

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    12 May 2003

  Author:

    John Burkardt

  Parameters:

    Input, int N, the argument of the double factorial function.
    If N is less than 1, I4_FACTORIAL2 is returned as 1.

    Output, int I4_FACTORIAL2, the value of the double factorial of N.
*/
{
  int value;

  if ( n < 1 )
  {
    return 1;
  }

  value = 1;

  while ( 1 < n )
  {
    value = value * n;
    n = n - 2;
  }

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

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 GNU LGPL license. 

  Modified:

    11 September 2012

  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
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P00_EXACT - Fatal error!\n" );
    fprintf ( stderr, "  Illegal problem number = %d", problem );
    exit ( 1 );
  }

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

void p00_fun ( int problem, int option, int n, double x[], double f[] )

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

    P00_FUN evaluates the integrand for any problem.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

    Input, int PROBLEM, the index of the problem.

    Input, int OPTION:
    0, integrand is f(x).
    1, integrand is exp(-x*x) * f(x);
    2, integrand is exp(-x*x/2) * f(x);

    Input, int N, the number of points.

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

    Output, double F[N], the function values.
*/
{
  if ( problem == 1 )
  {
    p01_fun ( option, n, x, f );
  }
  else if ( problem == 2 )
  {
    p02_fun ( option, n, x, f );
  }
  else if ( problem == 3 )
  {
    p03_fun ( option, n, x, f );
  }
  else if ( problem == 4 )
  {
    p04_fun ( option, n, x, f );
  }
  else if ( problem == 5 )
  {
    p05_fun ( option, n, x, f );
  }
  else if ( problem == 6 )
  {
    p06_fun ( option, n, x, f );
  }
  else if ( problem == 7 )
  {
    p07_fun ( option, n, x, f );
  }
  else if ( problem == 8 )
  {
    p08_fun ( option, n, x, f );
  }
  else
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P00_FUN - Fatal error!\n" );
    fprintf ( stderr, "  Illegal problem number = %d", problem );
    exit ( 1 );
  }

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

double p00_gauss_hermite ( int problem, int order )

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

    P00_GAUSS_HERMITE applies a Gauss-Hermite quadrature rule.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

    Input, int PROBLEM, the index of the problem.

    Input, int ORDER, the order of the rule to apply.

    Output, double P00_GAUSS_HERMITE, the approximate integral.
*/
{
  double *f_vec;
  int option;
  double result;
  double *weight;
  double *xtab;

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

  hermite_compute ( order, xtab, weight );

  option = 1;
  p00_fun ( problem, option, order, xtab, f_vec );

  result = r8vec_dot_product ( order, weight, f_vec );

  free ( f_vec );
  free ( weight );
  free ( xtab );

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

double p00_monte_carlo ( int problem, int order )

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

    P00_MONTE_CARLO applies a Monte Carlo procedure to Hermite integrals.

  Discussion:

    We wish to estimate the integral:

      I(f) = integral ( -oo < x < +oo ) f(x) exp ( - x * x ) dx

    We do this by a Monte Carlo sampling procedure, in which 
    we select N points X(1:N) from a standard normal distribution,
    and estimate

      Q(f) = sum ( 1 <= I <= N ) f(x(i)) / sqrt ( pi )

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    11 September 2012

  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_MONTE_CARLO, the approximate integral.
*/
{
  double *f_vec;
  int option;
  const double r8_pi = 3.141592653589793;
  double result;
  int seed;
  double weight;
  double *x_vec;

  seed = 123456789;
  x_vec = r8vec_normal_01_new ( order, &seed );

  option = 2;
  f_vec = ( double * ) malloc ( order * sizeof ( double ) );

  p00_fun ( problem, option, order, x_vec, f_vec );

  weight = ( double ) ( order ) / sqrt ( r8_pi ) / sqrt ( 2.0 );

  result = r8vec_sum ( order, f_vec ) / weight;

  free ( f_vec );
  free ( x_vec );

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

int p00_problem_num ( )

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

    P00_PROBLEM_NUM returns the number of test integration problems.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

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

  problem_num = 8;

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

char *p00_title ( int problem )

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

    P00_TITLE returns the title for any problem.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  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
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P00_TITLE - Fatal error!\n" );
    fprintf ( stderr, "  Illegal problem number = %d", problem );
    exit ( 1 );
  }

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

double p00_turing ( int problem, double h, double tol, int *n )

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

    P00_TURING applies the Turing quadrature rule.

  Discussion

    We consider the approximation:

      Integral ( -oo < x < +oo ) f(x) dx

      = h * Sum ( -oo < i < +oo ) f(nh) + error term

    Given H and a tolerance TOL, we start summing at I = 0, and
    adding one more term in the positive and negative I directions,
    until the absolute value of the next two terms being added 
    is less than TOL.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Reference:

    Alan Turing,
    A Method for the Calculation of the Zeta Function,
    Proceedings of the London Mathematical Society,
    Volume 48, 1943, pages 180-197.

  Parameters:

    Input, int PROBLEM, the index of the problem.

    Input, double H, the spacing to use.

    Input, double TOL, the tolerance.  

    Output, int N, the number of pairs of steps taken.
    The actual number of function evaluations is 2*N+1.

    Output, double P00_TURING, the approximate integral.
*/
{
  double f_vec[2];
  int n_too_many = 100000;
  int option;
  int order;
  double result;
  double xtab[2];

  option = 0;
  *n = 0;

  result = 0.0;
  order = 1;
  xtab[0] = 0.0;
  p00_fun ( problem, option, order, xtab, f_vec );
  result = result + h * f_vec[0];

  for ( ; ; )
  {
    *n = *n + 1;

    xtab[0] =   ( double ) ( *n ) * h;
    xtab[1] = - ( double ) ( *n ) * h;

    order = 2;
    p00_fun ( problem, option, order, xtab, f_vec );

    result = result + h * ( f_vec[0] + f_vec[1] );
/*
  Just do a simple-minded absolute error tolerance check to start with.
*/
    if ( r8_abs ( f_vec[0] ) + r8_abs ( f_vec[1] ) <= tol )
    {
      break;
    }
/*
  Just in case things go crazy.
*/
    if ( n_too_many <= *n )
    {
      fprintf ( stdout, "\n" );
      fprintf ( stdout, "P00_TURING - Warning!\n" );
      fprintf ( stdout, "  Number of steps exceeded N_TOO_MANY = %d\n", n_too_many );
      break;
    }

  }

  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:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

    Output, double P01_EXACT, the value of the integral.
*/
{
  double exact;
  double omega = 1.0;
  const double r8_pi = 3.141592653589793;

  exact = sqrt ( r8_pi ) * exp ( - omega * omega );

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

void p01_fun ( int option, int n, double x[], double f[] )

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

    P01_FUN evaluates the integrand for problem 1.

  Discussion:

    Squire gives exact value as sqrt(pi) * exp(-w*w).

    Integral ( -oo < x < +oo ) exp(-x*x) cos(2*w*x) dx

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Reference:

    William Squire,
    Comparison of Gauss-Hermite and Midpoint Quadrature with Application
    to the Voigt Function,
    in Numerical Integration: 
    Recent Developments, Software and Applications,
    edited by Patrick Keast, Graeme Fairweather,
    Reidel, 1987, pages 337-340,
    ISBN: 9027725144,
    LC: QA299.3.N38.

  Parameters:

    Input, int OPTION:
    0, integrand is f(x).
    1, integrand is exp(-x*x) * f(x);
    2, integrand is exp(-x*x/2) * f(x);

    Input, int N, the number of points.

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

    Output, double F[N], the function values.
*/
{
  int i;
  double omega = 1.0;

  for ( i = 0; i < n; i++ )
  {
    f[i] = cos ( 2.0 * omega * x[i] );
  }

  if ( option == 0 )
  {
    for ( i = 0; i < n; i++ )
    {
      f[i] = f[i] * exp ( - x[i] * x[i] );
    }
  }
  else if ( option == 1 )
  {
  }
  else if ( option == 2 )
  {
    for ( i = 0; i < n; i++ )
    {
      f[i] = f[i] * exp ( - x[i] * x[i] / 2.0 );
    }
  }

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

char *p01_title ( )

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

    P01_TITLE returns the title for problem 1.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

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

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

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

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

double p02_exact ( )

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

    P02_EXACT returns the exact integral for problem 2.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

    Output, double P02_EXACT, the value of the integral.
*/
{
  double exact;
  const double r8_pi = 3.141592653589793;

  exact = sqrt ( r8_pi );

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

void p02_fun ( int option, int n, double x[], double f[] )

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

    P02_FUN evaluates the integrand for problem 2.

  Discussion:

    The exact value is sqrt(pi).

    Integral ( -oo < x < +oo ) exp(-x*x) dx

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

    Input, int OPTION:
    0, integrand is f(x).
    1, integrand is exp(-x*x) * f(x);
    2, integrand is exp(-x*x/2) * f(x);

    Input, int N, the number of points.

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

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

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

  if ( option == 0 )
  {
    for ( i = 0; i < n; i++ )
    {
      f[i] = f[i] * exp ( - x[i] * x[i] );
    }
  }
  else if ( option == 1 )
  {
  }
  else if ( option == 2 )
  {
    for ( i = 0; i < n; i++ )
    {
      f[i] = f[i] * exp ( - x[i] * x[i] / 2.0 );
    }
  }

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

char *p02_title ( )

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

    P02_TITLE returns the title for problem 2.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

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

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

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

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

double p03_exact ( )

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

    P03_EXACT returns the exact integral for problem 3.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

    Output, double P03_EXACT, the value of the integral.
*/
{
  double exact;
  double p = 1.0;
  double q = 3.0;
  const double r8_pi = 3.141592653589793;

  exact = r8_pi / ( q * sin ( r8_pi * p / q ) );

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

void p03_fun ( int option, int n, double x[], double f[] )

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

    P03_FUN evaluates the integrand for problem 3.

  Discussion:

    The exact value is pi / (q*sin(pi*p/q) ), assuming 0 < p < q.

    Integral ( -oo < x < +oo ) exp(-px) / ( 1 + exp ( -qx) ) dx

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

    Input, int OPTION:
    0, integrand is f(x).
    1, integrand is exp(-x*x) * f(x);
    2, integrand is exp(-x*x/2) * f(x);

    Input, int N, the number of points.

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

    Output, double F[N], the function values.
*/
{
  int i;
  double p = 1.0;
  double q = 3.0;

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

  if ( option == 0 )
  {
  }
  else if ( option == 1 )
  {
    for ( i = 0; i < n; i++ )
    {
      f[i] = f[i] * exp ( + x[i] * x[i] );
    }
  }
  else if ( option == 2 )
  {
    for ( i = 0; i < n; i++ )
    {
      f[i] = f[i] * exp ( + x[i] * x[i] / 2.0 );
    }
  }

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

char *p03_title ( )

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

    P03_TITLE returns the title for problem 3.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

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

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

  strcpy ( title, "exp(-px) / ( 1 + exp(-qx) )" );

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

double p04_exact ( )

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

    P04_EXACT returns the estimated integral for problem 4.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

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

  exact = sqrt ( r8_pi / 2.0 );

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

void p04_fun ( int option, int n, double x[], double f[] )

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

    P04_FUN evaluates the integrand for problem 4.

  Discussion:

    The exact value is sqrt ( pi / 2 )

    Integral ( -oo < x < +oo ) sin ( x^2 ) dx

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

    Input, int OPTION:
    0, integrand is f(x).
    1, integrand is exp(-x*x) * f(x);
    2, integrand is exp(-x*x/2) * f(x);

    Input, int N, the number of points.

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

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

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

  if ( option == 0 )
  {
  }
  else if ( option == 1 )
  {
    for ( i = 0; i < n; i++ )
    {
      f[i] = f[i] * exp ( + x[i] * x[i] );
    }
  }
  else if ( option == 2 )
  {
    for ( i = 0; i < n; i++ )
    {
      f[i] = f[i] * exp ( + x[i] * x[i] / 2.0 );
    }
  }
  return;
}
/******************************************************************************/

char *p04_title ( )

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

    P04_TITLE returns the title for problem 4.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

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

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

  strcpy ( title, "sin(x^2)" );

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

double p05_exact ( )

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

    P05_EXACT returns the estimated integral for problem 5.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

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

  exact = r8_pi / 3.0;

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

void p05_fun ( int option, int n, double x[], double f[] )

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

    P05_FUN evaluates the integrand for problem 5.

  Discussion:

    The exact value is pi / 3.

    Integral ( -oo < x < +oo ) dx / ( (1+x^2) sqrt(4+3x^2) )

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

    Input, int OPTION:
    0, integrand is f(x).
    1, integrand is exp(-x*x) * f(x);
    2, integrand is exp(-x*x/2) * f(x);

    Input, int N, the number of points.

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

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

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

  if ( option == 0 )
  {
  }
  else if ( option == 1 )
  {
    for ( i = 0; i < n; i++ )
    {
      f[i] = f[i] * exp ( + x[i] * x[i] );
    }
  }
  else if ( option == 2 )
  {
    for ( i = 0; i < n; i++ )
    {
      f[i] = f[i] * exp ( + x[i] * x[i] / 2.0 );
    }
  }

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

char *p05_title ( )

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

    P05_TITLE returns the title for problem 5.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

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

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

  strcpy ( title, "1/( (1+x^2) sqrt(4+3x^2) )" );

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

double p06_exact ( )

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

    P06_EXACT returns the exact integral for problem 6.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

    Output, double P06_EXACT, the value of the integral.
*/
{
  double exact;
  int m;
  const double r8_pi = 3.141592653589793;

  p06_param ( 'G', 'M', &m );

  if ( m <= -1 )
  {
    exact = - r8_huge ( );
  }
  else if ( ( m % 2 ) == 1 )
  {
    exact = 0.0;
  }
  else
  {
    exact = ( double ) ( i4_factorial2 ( m - 1 ) ) * sqrt ( r8_pi ) 
      / pow ( 2.0, m / 2 );
  }

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

void p06_fun ( int option, int n, double x[], double f[] )

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

    P06_FUN evaluates the integrand for problem 6.

  Discussion:

    The exact value is (m-1)!! * sqrt ( pi ) / sqrt ( 2^m ).

    Integral ( -oo < x < +oo ) x^m exp (-x*x) dx

    The parameter M is set by calling P06_PARAM.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

    Input, int OPTION:
    0, integrand is f(x).
    1, integrand is exp(-x*x) * f(x);
    2, integrand is exp(-x*x/2) * f(x);

    Input, int N, the number of points.

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

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

  p06_param ( 'G', 'M', &m );

  for ( i = 0; i < n; i++ )
  {
    f[i] = pow ( x[i], m );
  }

  if ( option == 0 ) 
  {
    for ( i = 0; i < n; i++ )
    {
      f[i] = f[i] * exp ( - x[i] * x[i] );
    }
  }
  else if ( option == 1 )
  {
  }
  else if ( option == 2 )
  {
    for ( i = 0; i < n; i++ )
    {
      f[i] = f[i] * exp ( - x[i] * x[i] / 2.0 );
    }
  }

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

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

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

    P06_PARAM gets or sets parameters for problem 6.

  Discussion:

    The parameter is named "M", and it represents the value of the exponent
    in the integrand function:

    Integral ( -oo < x < +oo ) x^m exp (-x*x) dx

    M must be greater than -1.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

    Input, char ACTION, the action.
    'S' to set the value,
    'G' to get the value.

    Input, char NAME, the parameter name.
    'M', the exponent.

    Input/output, int *VALUE, the parameter value.
    If ACTION = 'S', then VALUE is an input quantity, and M is set to VALUE.
    If ACTION = 'G', then VALUE is an output quantity, and VALUE is set to M.
*/
{
  static int m = 0;

  if ( action == 'S' || action == 's' )
  {
    if ( *value <= -1 )
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "P06_PARAM - Fatal error!\n" );
      fprintf ( stderr, "  Parameter M must be greater than -1.\n" );
      exit ( 1 );
    }
    m = *value;
  }
  else if ( action == 'G' || action == 'g' )
  {
    *value = m;
  }
  else
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "P06_PARAM - Fatal error!\n" );
    fprintf ( stderr, "  Unrecognized value of ACTION = '%c'.\n", action );
    exit ( 1 );
  }
  return;
}
/******************************************************************************/

char *p06_title ( )

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

    P06_TITLE returns the title for problem 6.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

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

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

  strcpy ( title, "x^m exp(-x*x)" );

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

double p07_exact ( )

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

    P07_EXACT returns the exact integral for problem 7.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

    Output, double P07_EXACT, the value of the integral.
*/
{
  double e_sqrt_sqrt = 1.2840254166877414841;
  double exact;
  const double r8_pi = 3.141592653589793;

  exact = 0.25 * sqrt ( r8_pi ) / e_sqrt_sqrt;

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

void p07_fun ( int option, int n, double x[], double f[] )

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

    P07_FUN evaluates the integrand for problem 7.

  Discussion:

    The exact value is (1/4) sqrt(pi) / sqrt(sqrt(e)).

    Integral ( -oo < x < +oo ) x^2 cos(x) e^(-x^2) dx

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

    Input, int OPTION:
    0, integrand is f(x).
    1, integrand is exp(-x*x) * f(x);
    2, integrand is exp(-x*x/2) * f(x);

    Input, int N, the number of points.

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

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

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

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

char *p07_title ( )

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

    P07_TITLE returns the title for problem 7.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

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

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

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

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

double p08_exact ( )

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

    P08_EXACT returns the exact integral for problem 8.

  Discussion:

    The 20 digit value of the answer was computed by Mathematica.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

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

  exact = 3.0088235661136433510;

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

void p08_fun ( int option, int n, double x[], double f[] )

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

    P08_FUN evaluates the integrand for problem 8.

  Discussion:

    The exact value is sqrt ( 2 pi ) * HypergeometricU ( -1/2, 0, 1 ).

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

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

    Input, int OPTION:
    0, integrand is f(x).
    1, integrand is exp(-x*x) * f(x);
    2, integrand is exp(-x*x/2) * f(x);

    Input, int N, the number of points.

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

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

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

  if ( option == 0 ) 
  {
    for ( i = 0; i < n; i++ )
    {
      f[i] = f[i] * exp ( - 0.5 * x[i] * x[i] );
    }
  }
  else if ( option == 1 )
  {
    for ( i = 0; i < n; i++ )
    {
      f[i] = f[i] * exp ( + 0.5 * x[i] * x[i] );
    }
  }
  else if ( option == 2 )
  {

  }

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

char *p08_title ( )

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

    P08_TITLE returns the title for problem 8.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt

  Parameters:

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

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

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

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

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_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_gamma ( double x )

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

    R8_GAMMA evaluates Gamma(X) for a real argument.

  Discussion:

    This routine calculates the gamma function for a real argument X.

    Computation is based on an algorithm outlined in reference 1.
    The program uses rational functions that approximate the gamma
    function to at least 20 significant decimal digits.  Coefficients
    for the approximation over the interval (1,2) are unpublished.
    Those for the approximation for 12 <= X are from reference 2.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    18 January 2008

  Author:

    Original FORTRAN77 version by William Cody, Laura Stoltz.
    C version by John Burkardt.

  Reference:

    William Cody,
    An Overview of Software Development for Special Functions,
    in Numerical Analysis Dundee, 1975,
    edited by GA Watson,
    Lecture Notes in Mathematics 506,
    Springer, 1976.

    John Hart, Ward Cheney, Charles Lawson, Hans Maehly,
    Charles Mesztenyi, John Rice, Henry Thatcher,
    Christoph Witzgall,
    Computer Approximations,
    Wiley, 1968,
    LC: QA297.C64.

  Parameters:

    Input, double X, the argument of the function.

    Output, double R8_GAMMA, the value of the function.
*/
{
/*
  Coefficients for minimax approximation over (12, INF).
*/
  double c[7] = {
   -1.910444077728E-03, 
    8.4171387781295E-04, 
   -5.952379913043012E-04, 
    7.93650793500350248E-04, 
   -2.777777777777681622553E-03, 
    8.333333333333333331554247E-02, 
    5.7083835261E-03 };
  double eps = 2.22E-16;
  double fact;
  double half = 0.5;
  int i;
  int n;
  double one = 1.0;
  double p[8] = {
  -1.71618513886549492533811E+00,
   2.47656508055759199108314E+01, 
  -3.79804256470945635097577E+02,
   6.29331155312818442661052E+02, 
   8.66966202790413211295064E+02,
  -3.14512729688483675254357E+04, 
  -3.61444134186911729807069E+04,
   6.64561438202405440627855E+04 };
  int parity;
  double q[8] = {
  -3.08402300119738975254353E+01,
   3.15350626979604161529144E+02, 
  -1.01515636749021914166146E+03,
  -3.10777167157231109440444E+03, 
   2.25381184209801510330112E+04,
   4.75584627752788110767815E+03, 
  -1.34659959864969306392456E+05,
  -1.15132259675553483497211E+05 };
  const double r8_pi = 3.141592653589793;
  double res;
  double sqrtpi = 0.9189385332046727417803297;
  double sum;
  double two = 2.0;
  double value;
  double xbig = 171.624;
  double xden;
  double xinf = 1.79E+308;
  double xminin = 2.23E-308;
  double xnum;
  double y;
  double y1;
  double ysq;
  double z;
  double zero = 0.0;;

  parity = 0;
  fact = one;
  n = 0;
  y = x;
/*
  Argument is negative.
*/
  if ( y <= zero )
  {
    y = - x;
    y1 = ( double ) ( int ) ( y );
    res = y - y1;

    if ( res != zero )
    {
      if ( y1 != ( double ) ( int ) ( y1 * half ) * two )
      {
        parity = 1;
      }

      fact = - r8_pi / sin ( r8_pi * res );
      y = y + one;
    }
    else
    {
      res = xinf;
      value = res;
      return value;
    }
  }
/*
  Argument is positive.
*/
  if ( y < eps )
  {
/*
  Argument < EPS.
*/
    if ( xminin <= y )
    {
      res = one / y;
    }
    else
    {
      res = xinf;
      value = res;
      return value;
    }
  }
  else if ( y < 12.0 )
  {
    y1 = y;
/*
  0.0 < argument < 1.0.
*/
    if ( y < one )
    {
      z = y;
      y = y + one;
    }
/*
  1.0 < argument < 12.0.
  Reduce argument if necessary.
*/
    else
    {
      n = ( int ) ( y ) - 1;
      y = y - ( double ) ( n );
      z = y - one;
    }
/*
  Evaluate approximation for 1.0 < argument < 2.0.
*/
    xnum = zero;
    xden = one;
    for ( i = 0; i < 8; i++ )
    {
      xnum = ( xnum + p[i] ) * z;
      xden = xden * z + q[i];
    }
    res = xnum / xden + one;
/*
  Adjust result for case  0.0 < argument < 1.0.
*/
    if ( y1 < y )
    {
      res = res / y1;
    }
/*
  Adjust result for case 2.0 < argument < 12.0.
*/
    else if ( y < y1 )
    {
      for ( i = 1; i <= n; i++ )
      {
        res = res * y;
        y = y + one;
      }
    }
  }
  else
  {
/*
  Evaluate for 12.0 <= argument.
*/
    if ( y <= xbig )
    {
      ysq = y * y;
      sum = c[6];
      for ( i = 0; i < 6; i++ )
      {
        sum = sum / ysq + c[i];
      }
      sum = sum / y - y + sqrtpi;
      sum = sum + ( y - half ) * log ( y );
      res = exp ( sum );
    }
    else
    {
      res = xinf;
      value = res;
      return value;
    }
  }
/*
  Final adjustments and return.
*/
  if ( parity )
  {
    res = - res;
  }

  if ( fact != one )
  {
    res = fact / res;
  }

  value = res;

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

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_normal_01_new ( int n, int *seed )

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

    R8VEC_NORMAL_01_NEW returns a unit pseudonormal R8VEC.

  Discussion:

    The standard normal probability distribution function (PDF) has
    mean 0 and standard deviation 1.

    This routine can generate a vector of values on one call.  It
    has the feature that it should provide the same results
    in the same order no matter how we break up the task.

    Before calling this routine, the user may call RANDOM_SEED
    in order to set the seed of the random number generator.

    The Box-Muller method is used, which is efficient, but
    generates an even number of values each time.  On any call
    to this routine, an even number of new values are generated.
    Depending on the situation, one value may be left over.
    In that case, it is saved for the next call.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    18 February 2012

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of values desired.  If N is negative,
    then the code will flush its internal memory; in particular,
    if there is a saved value to be used on the next call, it is
    instead discarded.  This is useful if the user has reset the
    random number seed, for instance.

    Input/output, int *SEED, a seed for the random number generator.

    Output, double R8VEC_NORMAL_01_NEW[N], a sample of the standard normal PDF.

  Local parameters:

    Local, int MADE, records the number of values that have
    been computed.  On input with negative N, this value overwrites
    the return value of N, so the user can get an accounting of
    how much work has been done.

    Local, double R[N+1], is used to store some uniform random values.
    Its dimension is N+1, but really it is only needed to be the
    smallest even number greater than or equal to N.

    Local, int SAVED, is 0 or 1 depending on whether there is a
    single saved value left over from the previous call.

    Local, int X_LO, X_HI, records the range of entries of
    X that we need to compute.  This starts off as 1:N, but is adjusted
    if we have a saved value that can be immediately stored in X(1),
    and so on.

    Local, double Y, the value saved from the previous call, if
    SAVED is 1.
*/
{
  int i;
  int m;
  static int made = 0;
  double *r;
  const double r8_pi = 3.141592653589793;
  static int saved = 0;
  double *x;
  int x_hi;
  int x_lo;
  static double y = 0.0;

  x = ( double * ) malloc ( n * sizeof ( double ) );
/*
  I'd like to allow the user to reset the internal data.
  But this won't work properly if we have a saved value Y.
  I'm making a crock option that allows the user to signal
  explicitly that any internal memory should be flushed,
  by passing in a negative value for N.
*/
  if ( n < 0 )
  {
    made = 0;
    saved = 0;
    y = 0.0;
    return NULL;
  }
  else if ( n == 0 )
  {
    return NULL;
  }
/*
  Record the range of X we need to fill in.
*/
  x_lo = 1;
  x_hi = n;
/*
  Use up the old value, if we have it.
*/
  if ( saved == 1 )
  {
    x[0] = y;
    saved = 0;
    x_lo = 2;
  }
/*
  Maybe we don't need any more values.
*/
  if ( x_hi - x_lo + 1 == 0 )
  {
  }
/*
  If we need just one new value, do that here to avoid null arrays.
*/
  else if ( x_hi - x_lo + 1 == 1 )
  {
    r = r8vec_uniform_01_new ( 2, seed );

    x[x_hi-1] = sqrt ( - 2.0 * log ( r[0] ) ) * cos ( 2.0 * r8_pi * r[1] );
    y =         sqrt ( - 2.0 * log ( r[0] ) ) * sin ( 2.0 * r8_pi * r[1] );

    saved = 1;

    made = made + 2;

    free ( r );
  }
/*
  If we require an even number of values, that's easy.
*/
  else if ( ( x_hi - x_lo + 1 ) % 2 == 0 )
  {
    m = ( x_hi - x_lo + 1 ) / 2;

    r = r8vec_uniform_01_new ( 2*m, seed );

    for ( i = 0; i <= 2*m-2; i = i + 2 )
    {
      x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * r8_pi * r[i+1] );
      x[x_lo+i  ] = sqrt ( - 2.0 * log ( r[i] ) ) * sin ( 2.0 * r8_pi * r[i+1] );
    }
    made = made + x_hi - x_lo + 1;

    free ( r );
  }
/*
  If we require an odd number of values, we generate an even number,
  and handle the last pair specially, storing one in X(N), and
  saving the other for later.
*/
  else
  {
    x_hi = x_hi - 1;

    m = ( x_hi - x_lo + 1 ) / 2 + 1;

    r = r8vec_uniform_01_new ( 2*m, seed );

    for ( i = 0; i <= 2*m-4; i = i + 2 )
    {
      x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * r8_pi * r[i+1] );
      x[x_lo+i  ] = sqrt ( - 2.0 * log ( r[i] ) ) * sin ( 2.0 * r8_pi * r[i+1] );
    }

    i = 2*m - 2;

    x[x_lo+i-1] = sqrt ( - 2.0 * log ( r[i] ) ) * cos ( 2.0 * r8_pi * r[i+1] );
    y           = sqrt ( - 2.0 * log ( r[i] ) ) * sin ( 2.0 * r8_pi * r[i+1] );

    saved = 1;

    made = made + x_hi - x_lo + 2;

    free ( r );
  }

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

void r8vec_reverse ( int n, double a[] )

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

    R8VEC_REVERSE reverses the elements of an R8VEC.

  Discussion:

    An R8VEC is a vector of R8's.

  Example:

    Input:

      N = 5, A = ( 11.0, 12.0, 13.0, 14.0, 15.0 ).

    Output:

      A = ( 15.0, 14.0, 13.0, 12.0, 11.0 ).

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    23 August 2010

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of entries in the array.

    Input/output, double A[N], the array to be reversed.
*/
{
  int i;
  double temp;

  for ( i = 1; i <= n/2; i++ )
  {
    temp   = a[i-1];
    a[i-1] = a[n-i];
    a[n-i] = temp;
  }

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

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

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 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 );

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

  return;
# undef TIME_SIZE
}