[ceph.git] / ceph / src / boost / libs / math / test / gauss_quadrature_test.cpp

// Copyright Nick Thompson, 2017
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)

#define BOOST_TEST_MODULE tanh_sinh_quadrature_test

#include <complex>
//#include <boost/multiprecision/mpc.hpp>
#include <boost/config.hpp>
#include <boost/detail/workaround.hpp>

#if !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_TRAILING_RESULT_TYPES) && !defined(BOOST_NO_SFINAE_EXPR)

#include <boost/math/concepts/real_concept.hpp>
#include <boost/test/included/unit_test.hpp>
#include <boost/test/tools/floating_point_comparison.hpp>
#include <boost/math/quadrature/gauss.hpp>
#include <boost/math/special_functions/sinc.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/multiprecision/cpp_complex.hpp>

#ifdef BOOST_HAS_FLOAT128
#include <boost/multiprecision/complex128.hpp>
#endif

#ifdef _MSC_VER
#pragma warning(disable:4127)  // Conditional expression is constant
#endif

#if !defined(TEST1) && !defined(TEST2) && !defined(TEST3)
#  define TEST1
#  define TEST2
#  define TEST3
#endif

using std::expm1;
using std::atan;
using std::tan;
using std::log;
using std::log1p;
using std::asinh;
using std::atanh;
using std::sqrt;
using std::isnormal;
using std::abs;
using std::sinh;
using std::tanh;
using std::cosh;
using std::pow;
using std::exp;
using std::sin;
using std::cos;
using std::string;
using boost::math::quadrature::gauss;
using boost::math::constants::pi;
using boost::math::constants::half_pi;
using boost::math::constants::two_div_pi;
using boost::math::constants::two_pi;
using boost::math::constants::half;
using boost::math::constants::third;
using boost::math::constants::half;
using boost::math::constants::third;
using boost::math::constants::catalan;
using boost::math::constants::ln_two;
using boost::math::constants::root_two;
using boost::math::constants::root_two_pi;
using boost::math::constants::root_pi;
using boost::multiprecision::cpp_bin_float_quad;

//
// Error rates depend only on the number of points in the approximation, not the type being tested,
// define all our expected errors here:
//

enum
{
   test_ca_error_id,
   test_ca_error_id_2,
   test_three_quad_error_id,
   test_three_quad_error_id_2,
   test_integration_over_real_line_error_id,
   test_right_limit_infinite_error_id,
   test_left_limit_infinite_error_id
};

template <unsigned Points>
double expected_error(unsigned)
{
   return 0; // placeholder, all tests will fail
}

template <>
double expected_error<7>(unsigned id)
{
   switch (id)
   {
   case test_ca_error_id:
      return 1e-7;
   case test_ca_error_id_2:
      return 2e-5;
   case test_three_quad_error_id:
      return 1e-8;
   case test_three_quad_error_id_2:
      return 3.5e-3;
   case test_integration_over_real_line_error_id:
      return 6e-3;
   case test_right_limit_infinite_error_id:
   case test_left_limit_infinite_error_id:
      return 1e-5;
   }
   return 0;  // placeholder, all tests will fail
}

template <>
double expected_error<9>(unsigned id)
{
   switch (id)
   {
   case test_ca_error_id:
      return 1e-7;
   case test_ca_error_id_2:
      return 2e-5;
   case test_three_quad_error_id:
      return 1e-8;
   case test_three_quad_error_id_2:
      return 3.5e-3;
   case test_integration_over_real_line_error_id:
      return 6e-3;
   case test_right_limit_infinite_error_id:
   case test_left_limit_infinite_error_id:
      return 1e-5;
   }
   return 0;  // placeholder, all tests will fail
}

template <>
double expected_error<10>(unsigned id)
{
   switch (id)
   {
   case test_ca_error_id:
      return 1e-12;
   case test_ca_error_id_2:
      return 3e-6;
   case test_three_quad_error_id:
      return 2e-13;
   case test_three_quad_error_id_2:
      return 2e-3;
   case test_integration_over_real_line_error_id:
      return 6e-3;  // doesn't get any better with more points!
   case test_right_limit_infinite_error_id:
   case test_left_limit_infinite_error_id:
      return 5e-8;
   }
   return 0;  // placeholder, all tests will fail
}

template <>
double expected_error<15>(unsigned id)
{
   switch (id)
   {
   case test_ca_error_id:
      return 6e-20;
   case test_ca_error_id_2:
      return 3e-7;
   case test_three_quad_error_id:
      return 1e-19;
   case test_three_quad_error_id_2:
      return 6e-4;
   case test_integration_over_real_line_error_id:
      return 6e-3;  // doesn't get any better with more points!
   case test_right_limit_infinite_error_id:
   case test_left_limit_infinite_error_id:
      return 5e-11;
   }
   return 0;  // placeholder, all tests will fail
}

template <>
double expected_error<20>(unsigned id)
{
   switch (id)
   {
   case test_ca_error_id:
      return 1e-26;
   case test_ca_error_id_2:
      return 1e-7;
   case test_three_quad_error_id:
      return 3e-27;
   case test_three_quad_error_id_2:
      return 3e-4;
   case test_integration_over_real_line_error_id:
      return 5e-5;  // doesn't get any better with more points!
   case test_right_limit_infinite_error_id:
   case test_left_limit_infinite_error_id:
      return 1e-15;
   }
   return 0;  // placeholder, all tests will fail
}

template <>
double expected_error<25>(unsigned id)
{
   switch (id)
   {
   case test_ca_error_id:
      return 5e-33;
   case test_ca_error_id_2:
      return 1e-8;
   case test_three_quad_error_id:
      return 1e-32;
   case test_three_quad_error_id_2:
      return 3e-4;
   case test_integration_over_real_line_error_id:
      return 1e-14;
   case test_right_limit_infinite_error_id:
   case test_left_limit_infinite_error_id:
      return 3e-19;
   }
   return 0;  // placeholder, all tests will fail
}

template <>
double expected_error<30>(unsigned id)
{
   switch (id)
   {
   case test_ca_error_id:
      return 2e-34;
   case test_ca_error_id_2:
      return 5e-9;
   case test_three_quad_error_id:
      return 4e-34;
   case test_three_quad_error_id_2:
      return 1e-4;
   case test_integration_over_real_line_error_id:
      return 1e-16;
   case test_right_limit_infinite_error_id:
   case test_left_limit_infinite_error_id:
      return 3e-23;
   }
   return 0;  // placeholder, all tests will fail
}


template<class Real, unsigned Points>
void test_linear()
{
    std::cout << "Testing linear functions are integrated properly by gauss on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
    Real tol = boost::math::tools::epsilon<Real>() * 10;
    auto f = [](const Real& x)
    {
       return 5*x + 7;
    };
    Real L1;
    Real Q = gauss<Real, Points>::integrate(f, (Real) 0, (Real) 1, &L1);
    BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol);
    BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol);
    Q = gauss<Real, Points>::integrate(f, (Real) 0, (Real) 0, &L1);
    BOOST_CHECK_CLOSE(Q, 0, tol);
    Q = gauss<Real, Points>::integrate(f, (Real) 1, (Real) 0, &L1);
    BOOST_CHECK_CLOSE_FRACTION(Q, -9.5, tol);
}

template<class Real, unsigned Points>
void test_quadratic()
{
    std::cout << "Testing quadratic functions are integrated properly by Gaussian quadrature on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
    Real tol = boost::math::tools::epsilon<Real>() * 10;

    auto f = [](const Real& x) { return 5*x*x + 7*x + 12; };
    Real L1;
    Real Q = gauss<Real, Points>::integrate(f, 0, 1, &L1);
    BOOST_CHECK_CLOSE_FRACTION(Q, (Real) 17 + half<Real>()*third<Real>(), tol);
    BOOST_CHECK_CLOSE_FRACTION(L1, (Real) 17 + half<Real>()*third<Real>(), tol);
}

// Examples taken from
//http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
template<class Real, unsigned Points>
void test_ca()
{
    std::cout << "Testing integration of C(a) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
    Real tol = expected_error<Points>(test_ca_error_id);
    Real L1;

    auto f1 = [](const Real& x) { return atan(x)/(x*(x*x + 1)) ; };
    Real Q = gauss<Real, Points>::integrate(f1, 0, 1, &L1);
    Real Q_expected = pi<Real>()*ln_two<Real>()/8 + catalan<Real>()*half<Real>();
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
    BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);

    auto f2 = [](Real x)->Real { Real t0 = x*x + 1; Real t1 = sqrt(t0); return atan(t1)/(t0*t1); };
    Q = gauss<Real, Points>::integrate(f2, 0 , 1, &L1);
    Q_expected = pi<Real>()/4 - pi<Real>()/root_two<Real>() + 3*atan(root_two<Real>())/root_two<Real>();
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
    BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);

    tol = expected_error<Points>(test_ca_error_id_2);
    auto f5 = [](Real t)->Real { return t*t*log(t)/((t*t - 1)*(t*t*t*t + 1)); };
    Q = gauss<Real, Points>::integrate(f5, 0 , 1);
    Q_expected = pi<Real>()*pi<Real>()*(2 - root_two<Real>())/32;
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
}

template<class Real, unsigned Points>
void test_three_quadrature_schemes_examples()
{
    std::cout << "Testing integral in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
    Real tol = expected_error<Points>(test_three_quad_error_id);
    Real Q;
    Real Q_expected;

    // Example 1:
    auto f1 = [](const Real& t) { return t*boost::math::log1p(t); };
    Q = gauss<Real, Points>::integrate(f1, 0 , 1);
    Q_expected = half<Real>()*half<Real>();
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);


    // Example 2:
    auto f2 = [](const Real& t) { return t*t*atan(t); };
    Q = gauss<Real, Points>::integrate(f2, 0 , 1);
    Q_expected = (pi<Real>() -2 + 2*ln_two<Real>())/12;
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 2 * tol);

    // Example 3:
    auto f3 = [](const Real& t) { return exp(t)*cos(t); };
    Q = gauss<Real, Points>::integrate(f3, 0, half_pi<Real>());
    Q_expected = boost::math::expm1(half_pi<Real>())*half<Real>();
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);

    // Example 4:
    auto f4 = [](Real x)->Real { Real t0 = sqrt(x*x + 2); return atan(t0)/(t0*(x*x+1)); };
    Q = gauss<Real, Points>::integrate(f4, 0 , 1);
    Q_expected = 5*pi<Real>()*pi<Real>()/96;
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);

    tol = expected_error<Points>(test_three_quad_error_id_2);
    // Example 5:
    auto f5 = [](const Real& t) { return sqrt(t)*log(t); };
    Q = gauss<Real, Points>::integrate(f5, 0 , 1);
    Q_expected = -4/ (Real) 9;
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);

    // Example 6:
    auto f6 = [](const Real& t) { return sqrt(1 - t*t); };
    Q = gauss<Real, Points>::integrate(f6, 0 , 1);
    Q_expected = pi<Real>()/4;
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
}


template<class Real, unsigned Points>
void test_integration_over_real_line()
{
    std::cout << "Testing integrals over entire real line in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
    Real tol = expected_error<Points>(test_integration_over_real_line_error_id);
    Real Q;
    Real Q_expected;
    Real L1;

    auto f1 = [](const Real& t) { return 1/(1+t*t);};
    Q = gauss<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), boost::math::tools::max_value<Real>(), &L1);
    Q_expected = pi<Real>();
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
    BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
}

template<class Real, unsigned Points>
void test_right_limit_infinite()
{
    std::cout << "Testing right limit infinite for Gaussian quadrature in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
    Real tol = expected_error<Points>(test_right_limit_infinite_error_id);
    Real Q;
    Real Q_expected;
    Real L1;

    // Example 11:
    auto f1 = [](const Real& t) { return 1/(1+t*t);};
    Q = gauss<Real, Points>::integrate(f1, 0, boost::math::tools::max_value<Real>(), &L1);
    Q_expected = half_pi<Real>();
    BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);

    auto f4 = [](const Real& t) { return 1/(1+t*t); };
    Q = gauss<Real, Points>::integrate(f4, 1, boost::math::tools::max_value<Real>(), &L1);
    Q_expected = pi<Real>()/4;
    BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
}

template<class Real, unsigned Points>
void test_left_limit_infinite()
{
    std::cout << "Testing left limit infinite for Gaussian quadrature in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
    Real tol = expected_error<Points>(test_left_limit_infinite_error_id);
    Real Q;
    Real Q_expected;

    // Example 11:
    auto f1 = [](const Real& t) { return 1/(1+t*t);};
    Q = gauss<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), Real(0));
    Q_expected = half_pi<Real>();
    BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
}

template<class Complex>
void test_complex_lambert_w()
{
    std::cout << "Testing that complex-valued integrands are integrated correctly by Gaussian quadrature on type " << boost::typeindex::type_id<Complex>().pretty_name() << "\n";
    typedef typename Complex::value_type Real;
    Real tol = 10e-9;
    using boost::math::constants::pi;
    Complex z{2, 3};
    auto lw = [&z](Real v)->Complex {
      using std::cos;
      using std::sin;
      using std::exp;
      Real sinv = sin(v);
      Real cosv = cos(v);

      Real cotv = cosv/sinv;
      Real cscv = 1/sinv;
      Real t = (1-v*cotv)*(1-v*cotv) + v*v;
      Real x = v*cscv*exp(-v*cotv);
      Complex den = z + x;
      Complex num = t*(z/pi<Real>());
      Complex res = num/den;
      return res;
    };

    //N[ProductLog[2+3*I], 150]
    Complex Q = gauss<Real, 30>::integrate(lw, (Real) 0, pi<Real>());
    BOOST_CHECK_CLOSE_FRACTION(Q.real(), boost::lexical_cast<Real>("1.09007653448579084630177782678166964987102108635357778056449870727913321296238687023915522935120701763447787503167111962008709116746523970476893277703"), tol);
    BOOST_CHECK_CLOSE_FRACTION(Q.imag(), boost::lexical_cast<Real>("0.530139720774838801426860213574121741928705631382703178297940568794784362495390544411799468140433404536019992695815009036975117285537382995180319280835"), tol);
}

BOOST_AUTO_TEST_CASE(gauss_quadrature_test)
{
  
#ifdef TEST1
    test_linear<double, 7>();
    test_quadratic<double, 7>();
    test_ca<double, 7>();
    test_three_quadrature_schemes_examples<double, 7>();
    test_integration_over_real_line<double, 7>();
    test_right_limit_infinite<double, 7>();
    test_left_limit_infinite<double, 7>();

    test_linear<double, 9>();
    test_quadratic<double, 9>();
    test_ca<double, 9>();
    test_three_quadrature_schemes_examples<double, 9>();
    test_integration_over_real_line<double, 9>();
    test_right_limit_infinite<double, 9>();
    test_left_limit_infinite<double, 9>();

    test_linear<cpp_bin_float_quad, 10>();
    test_quadratic<cpp_bin_float_quad, 10>();
    test_ca<cpp_bin_float_quad, 10>();
    test_three_quadrature_schemes_examples<cpp_bin_float_quad, 10>();
    test_integration_over_real_line<cpp_bin_float_quad, 10>();
    test_right_limit_infinite<cpp_bin_float_quad, 10>();
    test_left_limit_infinite<cpp_bin_float_quad, 10>();
#endif
#ifdef TEST2
    test_linear<cpp_bin_float_quad, 15>();
    test_quadratic<cpp_bin_float_quad, 15>();
    test_ca<cpp_bin_float_quad, 15>();
    test_three_quadrature_schemes_examples<cpp_bin_float_quad, 15>();
    test_integration_over_real_line<cpp_bin_float_quad, 15>();
    test_right_limit_infinite<cpp_bin_float_quad, 15>();
    test_left_limit_infinite<cpp_bin_float_quad, 15>();

    test_linear<cpp_bin_float_quad, 20>();
    test_quadratic<cpp_bin_float_quad, 20>();
    test_ca<cpp_bin_float_quad, 20>();
    test_three_quadrature_schemes_examples<cpp_bin_float_quad, 20>();
    test_integration_over_real_line<cpp_bin_float_quad, 20>();
    test_right_limit_infinite<cpp_bin_float_quad, 20>();
    test_left_limit_infinite<cpp_bin_float_quad, 20>();

    test_linear<cpp_bin_float_quad, 25>();
    test_quadratic<cpp_bin_float_quad, 25>();
    test_ca<cpp_bin_float_quad, 25>();
    test_three_quadrature_schemes_examples<cpp_bin_float_quad, 25>();
    test_integration_over_real_line<cpp_bin_float_quad, 25>();
    test_right_limit_infinite<cpp_bin_float_quad, 25>();
    test_left_limit_infinite<cpp_bin_float_quad, 25>();

    test_linear<cpp_bin_float_quad, 30>();
    test_quadratic<cpp_bin_float_quad, 30>();
    test_ca<cpp_bin_float_quad, 30>();
    test_three_quadrature_schemes_examples<cpp_bin_float_quad, 30>();
    test_integration_over_real_line<cpp_bin_float_quad, 30>();
    test_right_limit_infinite<cpp_bin_float_quad, 30>();
    test_left_limit_infinite<cpp_bin_float_quad, 30>();


#endif
#ifdef TEST3
    test_left_limit_infinite<cpp_bin_float_quad, 30>();
    test_complex_lambert_w<std::complex<double>>();
    test_complex_lambert_w<std::complex<long double>>();
#ifdef BOOST_HAS_FLOAT128
    test_left_limit_infinite<boost::multiprecision::float128, 30>();
    test_complex_lambert_w<boost::multiprecision::complex128>();
#endif
    test_complex_lambert_w<boost::multiprecision::cpp_complex_quad>();
#endif
}

#else

int main() { return 0; }

#endif
Commit	Line	Data
b32b8144 FG	1	// Copyright Nick Thompson, 2017
	2	// Use, modification and distribution are subject to the
	3	// Boost Software License, Version 1.0.
	4	// (See accompanying file LICENSE_1_0.txt
	5	// or copy at http://www.boost.org/LICENSE_1_0.txt)
	6
	7	#define BOOST_TEST_MODULE tanh_sinh_quadrature_test
	8
92f5a8d4 TL	9	#include <complex>
92f5a8d4 TL	10	//#include <boost/multiprecision/mpc.hpp>
b32b8144 FG	11	#include <boost/config.hpp>
	12	#include <boost/detail/workaround.hpp>
	13
	14	#if !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_TRAILING_RESULT_TYPES) && !defined(BOOST_NO_SFINAE_EXPR)
	15
	16	#include <boost/math/concepts/real_concept.hpp>
	17	#include <boost/test/included/unit_test.hpp>
92f5a8d4	18	#include <boost/test/tools/floating_point_comparison.hpp>
b32b8144 FG	19	#include <boost/math/quadrature/gauss.hpp>
	20	#include <boost/math/special_functions/sinc.hpp>
	21	#include <boost/multiprecision/cpp_bin_float.hpp>
92f5a8d4 TL	22	#include <boost/multiprecision/cpp_complex.hpp>
	23
	24	#ifdef BOOST_HAS_FLOAT128
	25	#include <boost/multiprecision/complex128.hpp>
	26	#endif
b32b8144 FG	27
	28	#ifdef _MSC_VER
	29	#pragma warning(disable:4127) // Conditional expression is constant
	30	#endif
	31
92f5a8d4	32	#if !defined(TEST1) && !defined(TEST2) && !defined(TEST3)
b32b8144 FG	33	# define TEST1
b32b8144 FG	34	# define TEST2
92f5a8d4	35	# define TEST3
b32b8144 FG	36	#endif
	37
	38	using std::expm1;
	39	using std::atan;
	40	using std::tan;
	41	using std::log;
	42	using std::log1p;
	43	using std::asinh;
	44	using std::atanh;
	45	using std::sqrt;
	46	using std::isnormal;
	47	using std::abs;
	48	using std::sinh;
	49	using std::tanh;
	50	using std::cosh;
	51	using std::pow;
	52	using std::exp;
	53	using std::sin;
	54	using std::cos;
	55	using std::string;
	56	using boost::math::quadrature::gauss;
	57	using boost::math::constants::pi;
	58	using boost::math::constants::half_pi;
	59	using boost::math::constants::two_div_pi;
	60	using boost::math::constants::two_pi;
	61	using boost::math::constants::half;
	62	using boost::math::constants::third;
	63	using boost::math::constants::half;
	64	using boost::math::constants::third;
	65	using boost::math::constants::catalan;
	66	using boost::math::constants::ln_two;
	67	using boost::math::constants::root_two;
	68	using boost::math::constants::root_two_pi;
	69	using boost::math::constants::root_pi;
	70	using boost::multiprecision::cpp_bin_float_quad;
	71
	72	//
	73	// Error rates depend only on the number of points in the approximation, not the type being tested,
	74	// define all our expected errors here:
	75	//
	76
	77	enum
	78	{
	79	test_ca_error_id,
	80	test_ca_error_id_2,
	81	test_three_quad_error_id,
	82	test_three_quad_error_id_2,
	83	test_integration_over_real_line_error_id,
	84	test_right_limit_infinite_error_id,
	85	test_left_limit_infinite_error_id
	86	};
	87
	88	template <unsigned Points>
	89	double expected_error(unsigned)
	90	{
	91	return 0; // placeholder, all tests will fail
	92	}
	93
	94	template <>
	95	double expected_error<7>(unsigned id)
	96	{
	97	switch (id)
	98	{
	99	case test_ca_error_id:
100	return 1e-7;
101	case test_ca_error_id_2:
102	return 2e-5;
103	case test_three_quad_error_id:
104	return 1e-8;
105	case test_three_quad_error_id_2:
106	return 3.5e-3;
107	case test_integration_over_real_line_error_id:
108	return 6e-3;
109	case test_right_limit_infinite_error_id:
110	case test_left_limit_infinite_error_id:
111	return 1e-5;
112	}
113	return 0; // placeholder, all tests will fail
114	}
115
116	template <>
117	double expected_error<9>(unsigned id)
118	{
119	switch (id)
120	{
121	case test_ca_error_id:
122	return 1e-7;
123	case test_ca_error_id_2:
124	return 2e-5;
125	case test_three_quad_error_id:
126	return 1e-8;
127	case test_three_quad_error_id_2:
128	return 3.5e-3;
129	case test_integration_over_real_line_error_id:
130	return 6e-3;
131	case test_right_limit_infinite_error_id:
132	case test_left_limit_infinite_error_id:
133	return 1e-5;
134	}
135	return 0; // placeholder, all tests will fail
136	}
137
138	template <>
139	double expected_error<10>(unsigned id)
140	{
141	switch (id)
142	{
143	case test_ca_error_id:
144	return 1e-12;
145	case test_ca_error_id_2:
146	return 3e-6;
147	case test_three_quad_error_id:
148	return 2e-13;
149	case test_three_quad_error_id_2:
150	return 2e-3;
151	case test_integration_over_real_line_error_id:
152	return 6e-3; // doesn't get any better with more points!
153	case test_right_limit_infinite_error_id:
154	case test_left_limit_infinite_error_id:
155	return 5e-8;
156	}
157	return 0; // placeholder, all tests will fail
158	}
159
160	template <>
161	double expected_error<15>(unsigned id)
162	{
163	switch (id)
164	{
165	case test_ca_error_id:
166	return 6e-20;
167	case test_ca_error_id_2:
168	return 3e-7;
169	case test_three_quad_error_id:
170	return 1e-19;
171	case test_three_quad_error_id_2:
172	return 6e-4;
173	case test_integration_over_real_line_error_id:
174	return 6e-3; // doesn't get any better with more points!
175	case test_right_limit_infinite_error_id:
176	case test_left_limit_infinite_error_id:
177	return 5e-11;
178	}
179	return 0; // placeholder, all tests will fail
180	}
181
182	template <>
183	double expected_error<20>(unsigned id)
184	{
185	switch (id)
186	{
187	case test_ca_error_id:
188	return 1e-26;
189	case test_ca_error_id_2:
190	return 1e-7;
191	case test_three_quad_error_id:
192	return 3e-27;
193	case test_three_quad_error_id_2:
194	return 3e-4;
195	case test_integration_over_real_line_error_id:
196	return 5e-5; // doesn't get any better with more points!
197	case test_right_limit_infinite_error_id:
198	case test_left_limit_infinite_error_id:
199	return 1e-15;
200	}
201	return 0; // placeholder, all tests will fail
202	}
203
204	template <>
205	double expected_error<25>(unsigned id)
206	{
207	switch (id)
208	{
209	case test_ca_error_id:
210	return 5e-33;
211	case test_ca_error_id_2:
212	return 1e-8;
213	case test_three_quad_error_id:
214	return 1e-32;
215	case test_three_quad_error_id_2:
216	return 3e-4;
217	case test_integration_over_real_line_error_id:
218	return 1e-14;
219	case test_right_limit_infinite_error_id:
220	case test_left_limit_infinite_error_id:
221	return 3e-19;
222	}
223	return 0; // placeholder, all tests will fail
224	}
225
226	template <>
227	double expected_error<30>(unsigned id)
228	{
229	switch (id)
230	{
231	case test_ca_error_id:
232	return 2e-34;
233	case test_ca_error_id_2:
234	return 5e-9;
235	case test_three_quad_error_id:
236	return 4e-34;
237	case test_three_quad_error_id_2:
238	return 1e-4;
239	case test_integration_over_real_line_error_id:
240	return 1e-16;
241	case test_right_limit_infinite_error_id:
242	case test_left_limit_infinite_error_id:
243	return 3e-23;
244	}
245	return 0; // placeholder, all tests will fail
246	}
247
248
249	template<class Real, unsigned Points>
250	void test_linear()
251	{
252	std::cout << "Testing linear functions are integrated properly by gauss on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
253	Real tol = boost::math::tools::epsilon<Real>() * 10;
254	auto f = [](const Real& x)
255	{
256	return 5*x + 7;
257	};
258	Real L1;
259	Real Q = gauss<Real, Points>::integrate(f, (Real) 0, (Real) 1, &L1);
260	BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol);
261	BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol);
f67539c2 TL	262	Q = gauss<Real, Points>::integrate(f, (Real) 0, (Real) 0, &L1);
	263	BOOST_CHECK_CLOSE(Q, 0, tol);
	264	Q = gauss<Real, Points>::integrate(f, (Real) 1, (Real) 0, &L1);
	265	BOOST_CHECK_CLOSE_FRACTION(Q, -9.5, tol);
b32b8144 FG	266	}
	267
	268	template<class Real, unsigned Points>
	269	void test_quadratic()
	270	{
92f5a8d4	271	std::cout << "Testing quadratic functions are integrated properly by Gaussian quadrature on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
b32b8144 FG	272	Real tol = boost::math::tools::epsilon<Real>() * 10;
	273
	274	auto f = [](const Real& x) { return 5xx + 7*x + 12; };
	275	Real L1;
	276	Real Q = gauss<Real, Points>::integrate(f, 0, 1, &L1);
	277	BOOST_CHECK_CLOSE_FRACTION(Q, (Real) 17 + half<Real>()*third<Real>(), tol);
	278	BOOST_CHECK_CLOSE_FRACTION(L1, (Real) 17 + half<Real>()*third<Real>(), tol);
	279	}
	280
	281	// Examples taken from
	282	//http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf
	283	template<class Real, unsigned Points>
	284	void test_ca()
	285	{
	286	std::cout << "Testing integration of C(a) on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
	287	Real tol = expected_error<Points>(test_ca_error_id);
	288	Real L1;
	289
	290	auto f1 = [](const Real& x) { return atan(x)/(x(xx + 1)) ; };
	291	Real Q = gauss<Real, Points>::integrate(f1, 0, 1, &L1);
	292	Real Q_expected = pi<Real>()ln_two<Real>()/8 + catalan<Real>()half<Real>();
	293	BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
	294	BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
	295
	296	auto f2 = [](Real x)->Real { Real t0 = xx + 1; Real t1 = sqrt(t0); return atan(t1)/(t0t1); };
	297	Q = gauss<Real, Points>::integrate(f2, 0 , 1, &L1);
	298	Q_expected = pi<Real>()/4 - pi<Real>()/root_two<Real>() + 3*atan(root_two<Real>())/root_two<Real>();
	299	BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
	300	BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
	301
	302	tol = expected_error<Points>(test_ca_error_id_2);
	303	auto f5 = [](Real t)->Real { return ttlog(t)/((tt - 1)(ttt*t + 1)); };
	304	Q = gauss<Real, Points>::integrate(f5, 0 , 1);
	305	Q_expected = pi<Real>()pi<Real>()(2 - root_two<Real>())/32;
	306	BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
	307	}
	308
	309	template<class Real, unsigned Points>
	310	void test_three_quadrature_schemes_examples()
	311	{
	312	std::cout << "Testing integral in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
	313	Real tol = expected_error<Points>(test_three_quad_error_id);
	314	Real Q;
	315	Real Q_expected;
	316
	317	// Example 1:
	318	auto f1 = [](const Real& t) { return t*boost::math::log1p(t); };
	319	Q = gauss<Real, Points>::integrate(f1, 0 , 1);
	320	Q_expected = half<Real>()*half<Real>();
	321	BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
	322
	323
	324	// Example 2:
	325	auto f2 = [](const Real& t) { return ttatan(t); };
	326	Q = gauss<Real, Points>::integrate(f2, 0 , 1);
	327	Q_expected = (pi<Real>() -2 + 2*ln_two<Real>())/12;
	328	BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 2 * tol);
	329
	330	// Example 3:
	331	auto f3 = [](const Real& t) { return exp(t)*cos(t); };
	332	Q = gauss<Real, Points>::integrate(f3, 0, half_pi<Real>());
	333	Q_expected = boost::math::expm1(half_pi<Real>())*half<Real>();
	334	BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
	335
336	// Example 4:
337	auto f4 = [](Real x)->Real { Real t0 = sqrt(xx + 2); return atan(t0)/(t0(x*x+1)); };
338	Q = gauss<Real, Points>::integrate(f4, 0 , 1);
339	Q_expected = 5pi<Real>()pi<Real>()/96;
340	BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
341
342	tol = expected_error<Points>(test_three_quad_error_id_2);
343	// Example 5:
344	auto f5 = [](const Real& t) { return sqrt(t)*log(t); };
345	Q = gauss<Real, Points>::integrate(f5, 0 , 1);
346	Q_expected = -4/ (Real) 9;
347	BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
348
349	// Example 6:
350	auto f6 = [](const Real& t) { return sqrt(1 - t*t); };
351	Q = gauss<Real, Points>::integrate(f6, 0 , 1);
352	Q_expected = pi<Real>()/4;
353	BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
354	}
355
356
357	template<class Real, unsigned Points>
358	void test_integration_over_real_line()
359	{
360	std::cout << "Testing integrals over entire real line in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
361	Real tol = expected_error<Points>(test_integration_over_real_line_error_id);
362	Real Q;
363	Real Q_expected;
364	Real L1;
365
366	auto f1 = [](const Real& t) { return 1/(1+t*t);};
367	Q = gauss<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), boost::math::tools::max_value<Real>(), &L1);
368	Q_expected = pi<Real>();
369	BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);
370	BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol);
371	}
372
373	template<class Real, unsigned Points>
374	void test_right_limit_infinite()
375	{
92f5a8d4	376	std::cout << "Testing right limit infinite for Gaussian quadrature in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
b32b8144 FG	377	Real tol = expected_error<Points>(test_right_limit_infinite_error_id);
	378	Real Q;
	379	Real Q_expected;
	380	Real L1;
	381
	382	// Example 11:
	383	auto f1 = [](const Real& t) { return 1/(1+t*t);};
	384	Q = gauss<Real, Points>::integrate(f1, 0, boost::math::tools::max_value<Real>(), &L1);
	385	Q_expected = half_pi<Real>();
	386	BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
	387
	388	auto f4 = [](const Real& t) { return 1/(1+t*t); };
	389	Q = gauss<Real, Points>::integrate(f4, 1, boost::math::tools::max_value<Real>(), &L1);
	390	Q_expected = pi<Real>()/4;
	391	BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
	392	}
	393
	394	template<class Real, unsigned Points>
	395	void test_left_limit_infinite()
	396	{
92f5a8d4	397	std::cout << "Testing left limit infinite for Gaussian quadrature in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
b32b8144 FG	398	Real tol = expected_error<Points>(test_left_limit_infinite_error_id);
	399	Real Q;
	400	Real Q_expected;
	401
	402	// Example 11:
	403	auto f1 = [](const Real& t) { return 1/(1+t*t);};
	404	Q = gauss<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), Real(0));
	405	Q_expected = half_pi<Real>();
	406	BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol);
	407	}
	408
92f5a8d4 TL	409	template<class Complex>
	410	void test_complex_lambert_w()
	411	{
	412	std::cout << "Testing that complex-valued integrands are integrated correctly by Gaussian quadrature on type " << boost::typeindex::type_id<Complex>().pretty_name() << "\n";
	413	typedef typename Complex::value_type Real;
	414	Real tol = 10e-9;
	415	using boost::math::constants::pi;
	416	Complex z{2, 3};
	417	auto lw = [&z](Real v)->Complex {
	418	using std::cos;
	419	using std::sin;
	420	using std::exp;
	421	Real sinv = sin(v);
	422	Real cosv = cos(v);
	423
	424	Real cotv = cosv/sinv;
	425	Real cscv = 1/sinv;
	426	Real t = (1-vcotv)(1-vcotv) + vv;
	427	Real x = vcscvexp(-v*cotv);
	428	Complex den = z + x;
	429	Complex num = t*(z/pi<Real>());
	430	Complex res = num/den;
	431	return res;
	432	};
	433
	434	//N[ProductLog[2+3*I], 150]
	435	Complex Q = gauss<Real, 30>::integrate(lw, (Real) 0, pi<Real>());
	436	BOOST_CHECK_CLOSE_FRACTION(Q.real(), boost::lexical_cast<Real>("1.09007653448579084630177782678166964987102108635357778056449870727913321296238687023915522935120701763447787503167111962008709116746523970476893277703"), tol);
	437	BOOST_CHECK_CLOSE_FRACTION(Q.imag(), boost::lexical_cast<Real>("0.530139720774838801426860213574121741928705631382703178297940568794784362495390544411799468140433404536019992695815009036975117285537382995180319280835"), tol);
	438	}
	439
b32b8144 FG	440	BOOST_AUTO_TEST_CASE(gauss_quadrature_test)
b32b8144 FG	441	{
92f5a8d4	442
b32b8144 FG	443	#ifdef TEST1
	444	test_linear<double, 7>();
	445	test_quadratic<double, 7>();
	446	test_ca<double, 7>();
	447	test_three_quadrature_schemes_examples<double, 7>();
	448	test_integration_over_real_line<double, 7>();
	449	test_right_limit_infinite<double, 7>();
	450	test_left_limit_infinite<double, 7>();
	451
	452	test_linear<double, 9>();
	453	test_quadratic<double, 9>();
	454	test_ca<double, 9>();
	455	test_three_quadrature_schemes_examples<double, 9>();
	456	test_integration_over_real_line<double, 9>();
	457	test_right_limit_infinite<double, 9>();
	458	test_left_limit_infinite<double, 9>();
	459
	460	test_linear<cpp_bin_float_quad, 10>();
	461	test_quadratic<cpp_bin_float_quad, 10>();
	462	test_ca<cpp_bin_float_quad, 10>();
	463	test_three_quadrature_schemes_examples<cpp_bin_float_quad, 10>();
	464	test_integration_over_real_line<cpp_bin_float_quad, 10>();
	465	test_right_limit_infinite<cpp_bin_float_quad, 10>();
	466	test_left_limit_infinite<cpp_bin_float_quad, 10>();
	467	#endif
	468	#ifdef TEST2
	469	test_linear<cpp_bin_float_quad, 15>();
	470	test_quadratic<cpp_bin_float_quad, 15>();
	471	test_ca<cpp_bin_float_quad, 15>();
	472	test_three_quadrature_schemes_examples<cpp_bin_float_quad, 15>();
	473	test_integration_over_real_line<cpp_bin_float_quad, 15>();
	474	test_right_limit_infinite<cpp_bin_float_quad, 15>();
	475	test_left_limit_infinite<cpp_bin_float_quad, 15>();
	476
	477	test_linear<cpp_bin_float_quad, 20>();
	478	test_quadratic<cpp_bin_float_quad, 20>();
	479	test_ca<cpp_bin_float_quad, 20>();
	480	test_three_quadrature_schemes_examples<cpp_bin_float_quad, 20>();
	481	test_integration_over_real_line<cpp_bin_float_quad, 20>();
	482	test_right_limit_infinite<cpp_bin_float_quad, 20>();
	483	test_left_limit_infinite<cpp_bin_float_quad, 20>();
	484
	485	test_linear<cpp_bin_float_quad, 25>();
	486	test_quadratic<cpp_bin_float_quad, 25>();
	487	test_ca<cpp_bin_float_quad, 25>();
	488	test_three_quadrature_schemes_examples<cpp_bin_float_quad, 25>();
	489	test_integration_over_real_line<cpp_bin_float_quad, 25>();
	490	test_right_limit_infinite<cpp_bin_float_quad, 25>();
	491	test_left_limit_infinite<cpp_bin_float_quad, 25>();
	492
	493	test_linear<cpp_bin_float_quad, 30>();
	494	test_quadratic<cpp_bin_float_quad, 30>();
	495	test_ca<cpp_bin_float_quad, 30>();
	496	test_three_quadrature_schemes_examples<cpp_bin_float_quad, 30>();
	497	test_integration_over_real_line<cpp_bin_float_quad, 30>();
	498	test_right_limit_infinite<cpp_bin_float_quad, 30>();
	499	test_left_limit_infinite<cpp_bin_float_quad, 30>();
92f5a8d4 TL	500
	501
	502	#endif
	503	#ifdef TEST3
	504	test_left_limit_infinite<cpp_bin_float_quad, 30>();
	505	test_complex_lambert_w<std::complex<double>>();
	506	test_complex_lambert_w<std::complex<long double>>();
	507	#ifdef BOOST_HAS_FLOAT128
	508	test_left_limit_infinite<boost::multiprecision::float128, 30>();
	509	test_complex_lambert_w<boost::multiprecision::complex128>();
	510	#endif
	511	test_complex_lambert_w<boost::multiprecision::cpp_complex_quad>();
b32b8144 FG	512	#endif
	513	}
	514
	515	#else
	516
	517	int main() { return 0; }
	518
	519	#endif