# include <cstdlib>
# include <iostream>
# include <iomanip>
# include <cmath>
# include <ctime>
using namespace std;
# include "hermite_test_int.hpp"
//****************************************************************************80
void hermite_compute ( int order, double xtab[], double weight[] )
//****************************************************************************80
//
// 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:
//
// 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:
//
// 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 ( false )
{
for ( i = 0; i < order/2; i++ )
{
temp = xtab[i];
xtab[i] = xtab[order-1-i];
xtab[order-1-i] = temp;
}
}
return;
}
//****************************************************************************80
double hermite_integral ( int n )
//****************************************************************************80
//
// 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;
}
//****************************************************************************80
void hermite_recur ( double *p2, double *dp2, double *p1, double x, int order )
//****************************************************************************80
//
// 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;
}
//****************************************************************************80
void hermite_root ( double *x, int order, double *dp2, double *p1 )
//****************************************************************************80
//
// Purpose:
//
// HERMITE_ROOT improves an approximate root of a Hermite polynomial.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 09 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 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;
}
//****************************************************************************80
int i4_factorial2 ( int n )
//****************************************************************************80
//
// 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;
}
//****************************************************************************80
double p00_exact ( int problem )
//****************************************************************************80
//
// 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:
//
// 31 Julyl 2010
//
// 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
{
cout << "\n";
cout << "P00_EXACT - Fatal error!\n";
cout << " Illegal problem number = " << problem << "\n";
exit ( 1 );
}
return exact;
}
//****************************************************************************80
void p00_fun ( int problem, int option, int n, double x[], double f[] )
//****************************************************************************80
//
// Purpose:
//
// P00_FUN evaluates the integrand for any problem.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 31 July 2010
//
// 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
{
cout << "\n";
cout << "P00_FUN - Fatal error!\n";
cout << " Illegal problem number = " << problem << "\n";
exit ( 1 );
}
return;
}
//****************************************************************************80
double p00_gauss_hermite ( int problem, int order )
//****************************************************************************80
//
// Purpose:
//
// P00_GAUSS_HERMITE applies a Gauss-Hermite quadrature rule.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 26 May 2009
//
// 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 i;
int option;
double result;
double *weight;
double *xtab;
f_vec = new double[order];
weight = new double[order];
xtab = new double[order];
hermite_compute ( order, xtab, weight );
option = 1;
p00_fun ( problem, option, order, xtab, f_vec );
result = r8vec_dot_product ( order, weight, f_vec );
delete [] f_vec;
delete [] weight;
delete [] xtab;
return result;
}
//****************************************************************************80
double p00_monte_carlo ( int problem, int order )
//****************************************************************************80
//
// 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:
//
// 17 May 2010
//
// 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 = new double[order];
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;
delete [] f_vec;
delete [] x_vec;
return result;
}
//****************************************************************************80
int p00_problem_num ( )
//****************************************************************************80
//
// Purpose:
//
// P00_PROBLEM_NUM returns the number of test integration problems.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 31 July 2010
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Output, int P00_PROBLEM_NUM, the number of test problems.
//
{
int problem_num;
problem_num = 8;
return problem_num;
}
//****************************************************************************80
string p00_title ( int problem )
//****************************************************************************80
//
// Purpose:
//
// P00_TITLE returns the title for any problem.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 02 February 2010
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Input, int PROBLEM, the index of the problem.
//
// Output, string P00_TITLE, the title of the problem.
//
{
string 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
{
cout << "\n";
cout << "P00_TITLE - Fatal error!\n";
cout << " Illegal problem number = " << problem << "\n";
exit ( 1 );
}
return title;
}
//****************************************************************************80
double p00_turing ( int problem, double h, double tol, int *n )
//****************************************************************************80
//
// 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:
//
// 26 May 2009
//
// 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 i;
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 )
{
cout << "\n";
cout << "P00_TURING - Warning!\n";
cout << " Number of steps exceeded N_TOO_MANY = " << n_too_many << "\n";
break;
}
}
return result;
}
//****************************************************************************80
double p01_exact ( )
//****************************************************************************80
//
// Purpose:
//
// P01_EXACT returns the exact integral for problem 1.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 28 July 2007
//
// 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;
}
//****************************************************************************80
void p01_fun ( int option, int n, double x[], double f[] )
//****************************************************************************80
//
// 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:
//
// 26 May 2009
//
// 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;
}
//****************************************************************************80
string p01_title ( )
//****************************************************************************80
//
// Purpose:
//
// P01_TITLE returns the title for problem 1.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 26 May 2009
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Output, string P01_TITLE, the title of the problem.
//
{
string title;
title = "exp(-x*x) * cos(2*omega*x)";
return title;
}
//****************************************************************************80
double p02_exact ( )
//****************************************************************************80
//
// Purpose:
//
// P02_EXACT returns the exact integral for problem 2.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 31 July 2007
//
// 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;
}
//****************************************************************************80
void p02_fun ( int option, int n, double x[], double f[] )
//****************************************************************************80
//
// 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:
//
// 26 May 2009
//
// 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;
}
//****************************************************************************80
string p02_title ( )
//****************************************************************************80
//
// Purpose:
//
// P02_TITLE returns the title for problem 2.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 26 May 2009
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Output, string P02_TITLE, the title of the problem.
//
{
string title;
title = "exp(-x*x)";
return title;
}
//****************************************************************************80
double p03_exact ( )
//****************************************************************************80
//
// Purpose:
//
// P03_EXACT returns the exact integral for problem 3.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 28 July 2007
//
// 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;
}
//****************************************************************************80
void p03_fun ( int option, int n, double x[], double f[] )
//****************************************************************************80
//
// 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:
//
// 26 May 2009
//
// 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;
}
//****************************************************************************80
string p03_title ( )
//****************************************************************************80
//
// Purpose:
//
// P03_TITLE returns the title for problem 3.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 26 May 2009
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Output, string P03_TITLE, the title of the problem.
//
{
string title;
title = "exp(-px) / ( 1 + exp(-qx) )";
return title;
}
//****************************************************************************80
double p04_exact ( )
//****************************************************************************80
//
// Purpose:
//
// P04_EXACT returns the estimated integral for problem 4.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 28 July 2007
//
// 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;
}
//****************************************************************************80
void p04_fun ( int option, int n, double x[], double f[] )
//****************************************************************************80
//
// 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:
//
// 26 May 2009
//
// 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;
}
//****************************************************************************80
string p04_title ( )
//****************************************************************************80
//
// Purpose:
//
// P04_TITLE returns the title for problem 4.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 26 May 2009
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Output, string P04_TITLE, the title of the problem.
//
{
string title;
title = "sin(x^2)";
return title;
}
//****************************************************************************80
double p05_exact ( )
//****************************************************************************80
//
// Purpose:
//
// P05_EXACT returns the estimated integral for problem 5.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 06 October 2006
//
// 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;
}
//****************************************************************************80
void p05_fun ( int option, int n, double x[], double f[] )
//****************************************************************************80
//
// 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:
//
// 26 May 2009
//
// 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;
}
//****************************************************************************80
string p05_title ( )
//****************************************************************************80
//
// Purpose:
//
// P05_TITLE returns the title for problem 5.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 26 May 2009
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Output, string P05_TITLE, the title of the problem.
//
{
string title;
title = "1/( (1+x^2) sqrt(4+3x^2) )";
return title;
}
//****************************************************************************80
double p06_exact ( )
//****************************************************************************80
//
// Purpose:
//
// P06_EXACT returns the exact integral for problem 6.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 26 May 2009
//
// 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;
}
//****************************************************************************80
void p06_fun ( int option, int n, double x[], double f[] )
//****************************************************************************80
//
// 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:
//
// 26 May 2009
//
// 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;
}
//****************************************************************************80
void p06_param ( char action, char name, int *value )
//****************************************************************************80
//
// 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:
//
// 26 May 2009
//
// 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 )
{
cerr << "\n";
cerr << "P06_PARAM - Fatal error!\n";
cerr << " Parameter M must be greater than -1.\n";
exit ( 1 );
}
m = *value;
}
else if ( action == 'G' || action == 'g' )
{
*value = m;
}
else
{
cerr << "\n";
cerr << "P06_PARAM - Fatal error!\n";
cerr << " Unrecognized value of ACTION = \"" << action << "\".\n";
exit ( 1 );
}
return;
}
//****************************************************************************80
string p06_title ( )
//****************************************************************************80
//
// Purpose:
//
// P06_TITLE returns the title for problem 6.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 26 May 2009
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Output, string P06_TITLE, the title of the problem.
//
{
string title;
title = "x^m exp(-x*x)";
return title;
}
//****************************************************************************80
double p07_exact ( )
//****************************************************************************80
//
// Purpose:
//
// P07_EXACT returns the exact integral for problem 7.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 02 February 2010
//
// 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;
}
//****************************************************************************80
void p07_fun ( int option, int n, double x[], double f[] )
//****************************************************************************80
//
// 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:
//
// 02 February 2010
//
// 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;
}
//****************************************************************************80
string p07_title ( )
//****************************************************************************80
//
// Purpose:
//
// P07_TITLE returns the title for problem 7.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 02 February 2010
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Output, string P07_TITLE, the title of the problem.
//
{
string title;
title = "x^2 cos ( x ) exp(-x*x)";
return title;
}
//****************************************************************************80
double p08_exact ( )
//****************************************************************************80
//
// 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:
//
// 31 July 2010
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Output, double P08_EXACT, the value of the integral.
//
{
double exact;
exact = 3.0088235661136433510;
return exact;
}
//****************************************************************************80
void p08_fun ( int option, int n, double x[], double f[] )
//****************************************************************************80
//
// 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:
//
// 31 July 2010
//
// 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;
}
//****************************************************************************80
string p08_title ( )
//****************************************************************************80
//
// Purpose:
//
// P08_TITLE returns the title for problem 8.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 31 July 2010
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Output, string P08_TITLE, the title of the problem.
//
{
string title;
title = "sqrt(1+x*x/2) * exp(-x*x/2)";
return title;
}
//****************************************************************************80
double r8_abs ( double x )
//****************************************************************************80
//
// Purpose:
//
// R8_ABS returns the absolute value of an R8.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 14 November 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;
}
//****************************************************************************80
double r8_epsilon ( )
//****************************************************************************80
//
// Purpose:
//
// R8_EPSILON returns the R8 roundoff unit.
//
// Discussion:
//
// The roundoff unit 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;
}
//****************************************************************************80
double r8_gamma ( double x )
//****************************************************************************80
//
// 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 };
bool 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 = false;
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 = true;
}
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;
}
//****************************************************************************80
double r8_huge ( )
//****************************************************************************80
//
// Purpose:
//
// R8_HUGE returns a "huge" R8.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 31 August 2004
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Output, double R8_HUGE, a "huge" R8 value.
//
{
double value;
value = 0.1797693134862E+309;
return value;
}
//****************************************************************************80
double r8vec_dot_product ( int n, double a1[], double a2[] )
//****************************************************************************80
//
// 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:
//
// 03 July 2005
//
// 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;
}
//****************************************************************************80
double *r8vec_normal_01_new ( int n, int *seed )
//****************************************************************************80
//
// 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 October 2004
//
// 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, real 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, real 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 = new double[n];
//
// 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;
delete [] 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;
delete [] 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;
delete [] r;
}
return x;
}
//****************************************************************************80
void r8vec_reverse ( int n, double a[] )
//****************************************************************************80
//
// Purpose:
//
// R8VEC_REVERSE reverses the elements of an R8VEC.
//
// Discussion:
//
// An R8VEC is a vector of double precision values.
//
// 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:
//
// 30 April 2006
//
// 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 = 0; i < n / 2; i++ )
{
temp = a[i];
a[i] = a[n-1-i];
a[n-1-i] = temp;
}
return;
}
//****************************************************************************80
double r8vec_sum ( int n, double a[] )
//****************************************************************************80
//
// 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:
//
// 15 October 2004
//
// 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;
}
//****************************************************************************80
double *r8vec_uniform_01_new ( int n, int *seed )
//****************************************************************************80
//
// Purpose:
//
// R8VEC_UNIFORM_01_NEW returns a new unit pseudorandom R8VEC.
//
// Discussion:
//
// This routine implements the recursion
//
// seed = ( 16807 * seed ) mod ( 2^31 - 1 )
// u = 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 )
{
cerr << "\n";
cerr << "R8VEC_UNIFORM_01_NEW - Fatal error!\n";
cerr << " Input value of SEED = 0.\n";
exit ( 1 );
}
r = new double[n];
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;
}
//****************************************************************************80
void timestamp ( )
//****************************************************************************80
//
// 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;
size_t len;
time_t now;
now = time ( NULL );
tm = localtime ( &now );
len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
cout << time_buffer << "\n";
return;
# undef TIME_SIZE
}