[ceph.git] / ceph / src / boost / libs / math / test / gauss_kronrod_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 Gauss Kronrod_quadrature_test

#include <complex>
#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/math/tools/test_value.hpp>
#include <boost/test/included/unit_test.hpp>
#include <boost/test/tools/floating_point_comparison.hpp>
#include <boost/math/quadrature/gauss_kronrod.hpp>
#include <boost/math/special_functions/sinc.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/multiprecision/debug_adaptor.hpp>

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

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

#ifdef _MSC_VER
#pragma warning(disable:4127)  // Conditional expression is constant
#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_kronrod;
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;
using boost::multiprecision::cpp_dec_float_50;
using boost::multiprecision::debug_adaptor;
using boost::multiprecision::number;

//
// 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<15>(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<17>(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<21>(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<31>(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<41>(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<51>(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<61>(unsigned id)
{
   switch (id)
   {
   case test_ca_error_id:
      return 5e-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_kronrod on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
    Real tol = boost::math::tools::epsilon<Real>() * 10;
    Real error;
    auto f = [](const Real& x)->Real
    {
       return 5*x + 7;
    };
    Real L1;
    Real Q = gauss_kronrod<Real, Points>::integrate(f, (Real) 0, (Real) 1, 0, 0, &error, &L1);
    BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol);
    BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol);

    Q = gauss_kronrod<Real, Points>::integrate(f, (Real) 1, (Real) 0, 0, 0, &error, &L1);
    BOOST_CHECK_CLOSE_FRACTION(Q, -9.5, tol);
    BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol);

    Q = gauss_kronrod<Real, Points>::integrate(f, (Real) 0, (Real) 0, 0, 0, &error, &L1);
    BOOST_CHECK_CLOSE(Q, Real(0), tol);
}

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

    auto f = [](const Real& x)->Real { return 5*x*x + 7*x + 12; };
    Real L1;
    Real Q = gauss_kronrod<Real, Points>::integrate(f, 0, 1, 0, 0, &error, &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;
    Real error;

    auto f1 = [](const Real& x)->Real { return atan(x)/(x*(x*x + 1)) ; };
    Real Q = gauss_kronrod<Real, Points>::integrate(f1, 0, 1, 0, 0, &error, &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_kronrod<Real, Points>::integrate(f2, 0 , 1, 0, 0, &error, &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_kronrod<Real, Points>::integrate(f5, 0, 1, 0);
    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)->Real { return t*boost::math::log1p(t); };
    Q = gauss_kronrod<Real, Points>::integrate(f1, 0 , 1, 0);
    Q_expected = half<Real>()*half<Real>();
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);


    // Example 2:
    auto f2 = [](const Real& t)->Real { return t*t*atan(t); };
    Q = gauss_kronrod<Real, Points>::integrate(f2, 0 , 1, 0);
    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)->Real { return exp(t)*cos(t); };
    Q = gauss_kronrod<Real, Points>::integrate(f3, 0, half_pi<Real>(), 0);
    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_kronrod<Real, Points>::integrate(f4, 0 , 1, 0);
    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)->Real { return sqrt(t)*log(t); };
    Q = gauss_kronrod<Real, Points>::integrate(f5, 0 , 1, 0);
    Q_expected = -4/ (Real) 9;
    BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol);

    // Example 6:
    auto f6 = [](const Real& t)->Real { return sqrt(1 - t*t); };
    Q = gauss_kronrod<Real, Points>::integrate(f6, 0 , 1, 0);
    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;
    Real error;

    auto f1 = [](const Real& t)->Real { return 1/(1+t*t);};
    Q = gauss_kronrod<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), boost::math::tools::max_value<Real>(), 0, 0, &error, &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 Gauss Kronrod 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;
    Real error;

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

    auto f4 = [](const Real& t)->Real { return 1/(1+t*t); };
    Q = gauss_kronrod<Real, Points>::integrate(f4, 1, boost::math::tools::max_value<Real>(), 0, 0, &error, &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 Gauss Kronrod 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)->Real { return 1/(1+t*t);};
    Q = gauss_kronrod<Real, Points>::integrate(f1, -boost::math::tools::max_value<Real>(), Real(0), 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]
    boost::math::quadrature::gauss_kronrod<Real, 61> integrator;
    Complex Q = integrator.integrate(lw, (Real) 0, pi<Real>());
    BOOST_CHECK_CLOSE_FRACTION(Q.real(), BOOST_MATH_TEST_VALUE(Real, 1.0900765344857908463017778267816696498710210863535777805644), tol);
    BOOST_CHECK_CLOSE_FRACTION(Q.imag(), BOOST_MATH_TEST_VALUE(Real, 0.5301397207748388014268602135741217419287056313827031782979), tol);
}

BOOST_AUTO_TEST_CASE(gauss_quadrature_test)
{
#ifdef TEST1
    std::cout << "Testing 15 point approximation:\n";
    test_linear<double, 15>();
    test_quadratic<double, 15>();
    test_ca<double, 15>();
    test_three_quadrature_schemes_examples<double, 15>();
    test_integration_over_real_line<double, 15>();
    test_right_limit_infinite<double, 15>();
    test_left_limit_infinite<double, 15>();

    //  test one case where we do not have pre-computed constants:
    std::cout << "Testing 17 point approximation:\n";
    test_linear<double, 17>();
    test_quadratic<double, 17>();
    test_ca<double, 17>();
    test_three_quadrature_schemes_examples<double, 17>();
    test_integration_over_real_line<double, 17>();
    test_right_limit_infinite<double, 17>();
    test_left_limit_infinite<double, 17>();
    test_complex_lambert_w<std::complex<double>>();
    test_complex_lambert_w<std::complex<long double>>();
#endif
#ifdef TEST1A
    std::cout << "Testing 21 point approximation:\n";
    test_linear<cpp_bin_float_quad, 21>();
    test_quadratic<cpp_bin_float_quad, 21>();
    test_ca<cpp_bin_float_quad, 21>();
    test_three_quadrature_schemes_examples<cpp_bin_float_quad, 21>();
    test_integration_over_real_line<cpp_bin_float_quad, 21>();
    test_right_limit_infinite<cpp_bin_float_quad, 21>();
    test_left_limit_infinite<cpp_bin_float_quad, 21>();

    std::cout << "Testing 31 point approximation:\n";
    test_linear<cpp_bin_float_quad, 31>();
    test_quadratic<cpp_bin_float_quad, 31>();
    test_ca<cpp_bin_float_quad, 31>();
    test_three_quadrature_schemes_examples<cpp_bin_float_quad, 31>();
    test_integration_over_real_line<cpp_bin_float_quad, 31>();
    test_right_limit_infinite<cpp_bin_float_quad, 31>();
    test_left_limit_infinite<cpp_bin_float_quad, 31>();
#endif
#ifdef TEST2
    std::cout << "Testing 41 point approximation:\n";
    test_linear<cpp_bin_float_quad, 41>();
    test_quadratic<cpp_bin_float_quad, 41>();
    test_ca<cpp_bin_float_quad, 41>();
    test_three_quadrature_schemes_examples<cpp_bin_float_quad, 41>();
    test_integration_over_real_line<cpp_bin_float_quad, 41>();
    test_right_limit_infinite<cpp_bin_float_quad, 41>();
    test_left_limit_infinite<cpp_bin_float_quad, 41>();

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