From 7d2002db80c6e7cdffc69fa15d2c251a20077ad9 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Thu, 31 Aug 2017 19:42:26 +0100 Subject: [PATCH] Quadrature: add gauss and gauss-kronrod quadrature. --- include/boost/math/quadrature/gauss.hpp | 183 +++++++++ .../boost/math/quadrature/gauss_kronrod.hpp | 281 +++++++++++++ ...adaptive_gauss_kronrod_quadrature_test.cpp | 259 ++++++++++++ test/gauss_kronrod_quadrature_test.cpp | 377 ++++++++++++++++++ test/gauss_quadrature_test.cpp | 372 +++++++++++++++++ 5 files changed, 1472 insertions(+) create mode 100644 include/boost/math/quadrature/gauss.hpp create mode 100644 include/boost/math/quadrature/gauss_kronrod.hpp create mode 100644 test/adaptive_gauss_kronrod_quadrature_test.cpp create mode 100644 test/gauss_kronrod_quadrature_test.cpp create mode 100644 test/gauss_quadrature_test.cpp diff --git a/include/boost/math/quadrature/gauss.hpp b/include/boost/math/quadrature/gauss.hpp new file mode 100644 index 000000000..5c75844cf --- /dev/null +++ b/include/boost/math/quadrature/gauss.hpp @@ -0,0 +1,183 @@ +// Copyright John Maddock 2015. +// 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) + +#ifndef BOOST_MATH_QUADRATURE_GAUSS_HPP +#define BOOST_MATH_QUADRATURE_GAUSS_HPP + +#ifdef _MSC_VER +#pragma once +#endif + +#include +#include + +namespace boost { namespace math{ namespace quadrature{ namespace detail{ + +template +class gauss_detail +{ + static std::vector calculate_weights() + { + std::vector result(abscissa().size(), 0); + for (unsigned i = 0; i < abscissa().size(); ++i) + { + Real x = abscissa()[i]; + Real p = boost::math::legendre_p_prime(N, x); + result[i] = 2 / ((1 - x * x) * p * p); + } + return result; + } +public: + static const std::vector& abscissa() + { + static std::vector data = boost::math::legendre_p_zeros(N); + return data; + } + static const std::vector& weights() + { + static std::vector data = calculate_weights(); + return data; + } +}; + +template <> +class gauss_detail +{ +public: + static constexpr std::array const & abscissa() + { + static constexpr std::array data = { + 0.000000000000000000000000000000000e+00, + 4.058451513773971669066064120769615e-01, + 7.415311855993944398638647732807884e-01, + 9.491079123427585245261896840478513e-01, + }; + return data; + } + static constexpr std::array const & weights() + { + static constexpr std::array data = { + 4.179591836734693877551020408163265e-01, + 3.818300505051189449503697754889751e-01, + 2.797053914892766679014677714237796e-01, + 1.294849661688696932706114326790820e-01, + }; + return data; + } +}; + +} + +template > +class gauss : public detail::gauss_detail +{ +public: + typedef Real value_type; + + template + static value_type integrate(F f, Real* pL1 = nullptr) + { + using std::fabs; + unsigned non_zero_start = 1; + value_type result = 0; + if (N & 1) + result = f(value_type(0)) * weights()[0]; + else + non_zero_start = 0; + value_type L1 = fabs(result); + for (unsigned i = non_zero_start; i < abscissa().size(); ++i) + { + value_type fp = f(abscissa()[i]); + value_type fm = f(-abscissa()[i]); + result += (fp + fm) * weights()[i]; + L1 += (fabs(fp) + fabs(fm)) * weights()[i]; + } + if (pL1) + *pL1 = L1; + return result; + } + template + static value_type integrate(F f, Real a, Real b, Real* pL1 = nullptr) + { + static const char* function = "boost::math::quadrature::gauss<%1%>::integrate(f, %1%, %1%)"; + if (!(boost::math::isnan)(a) && !(boost::math::isnan)(b)) + { + // Infinite limits: + if ((a <= -tools::max_value()) && (b >= tools::max_value())) + { + auto u = [&](const Real& t)->Real + { + Real t_sq = t*t; + Real inv = 1 / (1 - t_sq); + return f(t*inv)*(1 + t_sq)*inv*inv; + }; + return integrate(u, pL1); + } + + // Right limit is infinite: + if ((boost::math::isfinite)(a) && (b >= tools::max_value())) + { + auto u = [&](const Real& t)->Real + { + Real z = 1 / (t + 1); + Real arg = 2 * z + a - 1; + return f(arg)*z*z; + }; + Real Q = 2 * integrate(u, pL1); + if (pL1) + { + *pL1 *= 2; + } + return Q; + } + + if ((boost::math::isfinite)(b) && (a <= -tools::max_value())) + { + auto v = [&](const Real& t)->Real + { + Real z = 1 / (t + 1); + Real arg = 2 * z - 1; + return f(b - arg) * z * z; + }; + Real Q = 2 * integrate(v, pL1); + if (pL1) + { + *pL1 *= 2; + } + return Q; + } + + if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b)) + { + if (b <= a) + { + return policies::raise_domain_error(function, "Arguments to integrate are in wrong order; integration over [a,b] must have b > a.", a, Policy()); + } + Real avg = (a + b)*half(); + Real scale = (b - a)*half(); + + auto u = [&](Real z)->Real + { + return f(avg + scale*z); + }; + Real Q = scale*integrate(u, pL1); + + if (pL1) + { + *pL1 *= scale; + } + return Q; + } + } + return policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy()); + } +}; + +} // namespace quadrature +} // namespace math +} // namespace boost + +#endif // BOOST_MATH_QUADRATURE_GAUSS_HPP + diff --git a/include/boost/math/quadrature/gauss_kronrod.hpp b/include/boost/math/quadrature/gauss_kronrod.hpp new file mode 100644 index 000000000..08811d801 --- /dev/null +++ b/include/boost/math/quadrature/gauss_kronrod.hpp @@ -0,0 +1,281 @@ +// Copyright John Maddock 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) + +#ifndef BOOST_MATH_QUADRATURE_GAUSS_KRONROD_HPP +#define BOOST_MATH_QUADRATURE_GAUSS_KRONROD_HPP + +#ifdef _MSC_VER +#pragma once +#pragma warning(push) +#pragma warning(disable: 4127) +#endif + +#include +#include +#include +#include + +namespace boost { namespace math{ namespace quadrature{ namespace detail{ + +template +class gauss_kronrod_detail +{ + static legendre_stieltjes const& get_legendre_stieltjes() + { + static const legendre_stieltjes data((N - 1) / 2 + 1); + return data; + } + static std::vector calculate_abscissa() + { + static std::vector result = boost::math::legendre_p_zeros((N - 1) / 2); + const legendre_stieltjes E = get_legendre_stieltjes(); + std::vector ls_zeros = E.zeros(); + result.insert(result.end(), ls_zeros.begin(), ls_zeros.end()); + std::sort(result.begin(), result.end()); + return result; + } + static std::vector calculate_weights() + { + std::vector result(abscissa().size(), 0); + unsigned gauss_order = (N - 1) / 2; + unsigned gauss_start = gauss_order & 1 ? 0 : 1; + const legendre_stieltjes& E = get_legendre_stieltjes(); + + for (unsigned i = gauss_start; i < abscissa().size(); i += 2) + { + Real x = abscissa()[i]; + Real p = boost::math::legendre_p_prime(gauss_order, x); + Real gauss_weight = 2 / ((1 - x * x) * p * p); + result[i] = gauss_weight + static_cast(2) / (static_cast(gauss_order + 1) * legendre_p_prime(gauss_order, x) * E(x)); + } + for (unsigned i = gauss_start ? 0 : 1; i < abscissa().size(); i += 2) + { + Real x = abscissa()[i]; + result[i] = static_cast(2) / (static_cast(gauss_order + 1) * legendre_p(gauss_order, x) * E.prime(x)); + } + return result; + } +public: + static const std::vector& abscissa() + { + static std::vector data = calculate_abscissa(); + return data; + } + static const std::vector& weights() + { + static std::vector data = calculate_weights(); + return data; + } +}; + +template <> +class gauss_kronrod_detail +{ +public: + static constexpr std::array const & abscissa() + { + static constexpr std::array data = { + 0.000000000000000000000000000000000e+00, + 2.077849550078984676006894037732449e-01, + 4.058451513773971669066064120769615e-01, + 5.860872354676911302941448382587296e-01, + 7.415311855993944398638647732807884e-01, + 8.648644233597690727897127886409262e-01, + 9.491079123427585245261896840478513e-01, + 9.914553711208126392068546975263285e-01, + }; + return data; + } + static constexpr std::array const & weights() + { + static constexpr std::array data = { + 2.094821410847278280129991748917143e-01, + 2.044329400752988924141619992346491e-01, + 1.903505780647854099132564024210137e-01, + 1.690047266392679028265834265985503e-01, + 1.406532597155259187451895905102379e-01, + 1.047900103222501838398763225415180e-01, + 6.309209262997855329070066318920429e-02, + 2.293532201052922496373200805896959e-02, + }; + return data; + } +}; + +} + +template > +class gauss_kronrod : public detail::gauss_kronrod_detail +{ +public: + typedef Real value_type; +private: + template + static value_type integrate_non_adaptive_m1_1(F f, Real* error = nullptr, Real* pL1 = nullptr) + { + using std::fabs; + unsigned gauss_start = 2; + unsigned kronrod_start = 1; + unsigned gauss_order = (N - 1) / 2; + value_type kronrod_result = 0; + value_type gauss_result = 0; + value_type fp, fm; + if (gauss_order & 1) + { + fp = f(value_type(0)); + kronrod_result = fp * weights()[0]; + gauss_result += fp * gauss::weights()[0]; + } + else + { + fp = f(value_type(0)); + kronrod_result = fp * weights()[0]; + gauss_start = 1; + kronrod_start = 2; + } + value_type L1 = fabs(kronrod_result); + for (unsigned i = gauss_start; i < abscissa().size(); i += 2) + { + fp = f(abscissa()[i]); + fm = f(-abscissa()[i]); + kronrod_result += (fp + fm) * weights()[i]; + L1 += (fabs(fp) + fabs(fm)) * weights()[i]; + gauss_result += (fp + fm) * gauss::weights()[i / 2]; + } + for (unsigned i = kronrod_start; i < abscissa().size(); i += 2) + { + fp = f(abscissa()[i]); + fm = f(-abscissa()[i]); + kronrod_result += (fp + fm) * weights()[i]; + L1 += (fabs(fp) + fabs(fm)) * weights()[i]; + } + if (pL1) + *pL1 = L1; + if (error) + *error = fabs(kronrod_result - gauss_result); + return kronrod_result; + } + + template + struct recursive_info + { + F f; + Real tol; + }; + + template + static value_type recursive_adaptive_integrate(const recursive_info* info, Real a, Real b, unsigned max_levels, Real abs_tol, Real* error, Real* L1) + { + Real error_local; + Real mean = (b + a) / 2; + Real scale = (b - a) / 2; + auto ff = [&](const Real& x) + { + return info->f(scale * x + mean); + }; + Real estimate = scale * integrate_non_adaptive_m1_1(ff, &error_local, L1); + + Real abs_tol1 = fabs(estimate * info->tol); + if (abs_tol == 0) + abs_tol = abs_tol1; + + if (max_levels && (abs_tol1 < error_local) && (abs_tol < error_local)) + { + Real mid = (a + b) / 2; + Real L1_local; + estimate = recursive_adaptive_integrate(info, a, mid, max_levels - 1, abs_tol / 2, error, L1); + estimate += recursive_adaptive_integrate(info, mid, b, max_levels - 1, abs_tol / 2, &error_local, &L1_local); + if (error) + *error += error_local; + if (L1) + *L1 += L1_local; + return estimate; + } + if(L1) + *L1 *= scale; + if (error) + *error = error_local; + return estimate; + } + +public: + template + static value_type integrate(F f, Real a, Real b, unsigned max_depth = 15, Real tol = tools::root_epsilon(), Real* error = nullptr, Real* pL1 = nullptr) + { + static const char* function = "boost::math::quadrature::gauss_kronrod<%1%>::integrate(f, %1%, %1%)"; + if (!(boost::math::isnan)(a) && !(boost::math::isnan)(b)) + { + // Infinite limits: + if ((a <= -tools::max_value()) && (b >= tools::max_value())) + { + auto u = [&](const Real& t)->Real + { + Real t_sq = t*t; + Real inv = 1 / (1 - t_sq); + return f(t*inv)*(1 + t_sq)*inv*inv; + }; + recursive_info info = { u, tol }; + return recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1); + } + + // Right limit is infinite: + if ((boost::math::isfinite)(a) && (b >= tools::max_value())) + { + auto u = [&](const Real& t)->Real + { + Real z = 1 / (t + 1); + Real arg = 2 * z + a - 1; + return f(arg)*z*z; + }; + recursive_info info = { u, tol }; + Real Q = 2 * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1); + if (pL1) + { + *pL1 *= 2; + } + return Q; + } + + if ((boost::math::isfinite)(b) && (a <= -tools::max_value())) + { + auto v = [&](const Real& t)->Real + { + Real z = 1 / (t + 1); + Real arg = 2 * z - 1; + return f(b - arg) * z * z; + }; + recursive_info info = { v, tol }; + Real Q = 2 * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1); + if (pL1) + { + *pL1 *= 2; + } + return Q; + } + + if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b)) + { + if (b <= a) + { + return policies::raise_domain_error(function, "Arguments to integrate are in wrong order; integration over [a,b] must have b > a.", a, Policy()); + } + recursive_info info = { f, tol }; + return recursive_adaptive_integrate(&info, a, b, max_depth, Real(0), error, pL1); + } + } + return policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy()); + } +}; + +} // namespace quadrature +} // namespace math +} // namespace boost + +#ifdef _MSC_VER +#pragma warning(pop) +#endif + +#endif // BOOST_MATH_QUADRATURE_GAUSS_KRONROD_HPP + diff --git a/test/adaptive_gauss_kronrod_quadrature_test.cpp b/test/adaptive_gauss_kronrod_quadrature_test.cpp new file mode 100644 index 000000000..6ac2577d0 --- /dev/null +++ b/test/adaptive_gauss_kronrod_quadrature_test.cpp @@ -0,0 +1,259 @@ +// 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 +#include + +#if !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_TRAILING_RESULT_TYPES) && !defined(BOOST_NO_SFINAE_EXPR) + +#include +#include +#include +#include +#include +#include + +#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; + +template +Real get_termination_condition() +{ + return boost::math::tools::epsilon() * 1000; +} + + +template +void test_linear() +{ + std::cout << "Testing linear functions are integrated properly by gauss_kronrod on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = boost::math::tools::epsilon() * 10; + Real error; + auto f = [](const Real& x) + { + return 5*x + 7; + }; + Real L1; + Real Q = gauss_kronrod::integrate(f, (Real) 0, (Real) 1, 15, get_termination_condition(), &error, &L1); + BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol); + BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol); +} + +template +void test_quadratic() +{ + std::cout << "Testing quadratic functions are integrated properly by tanh_sinh on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = boost::math::tools::epsilon() * 10; + Real error; + + auto f = [](const Real& x) { return 5*x*x + 7*x + 12; }; + Real L1; + Real Q = gauss_kronrod::integrate(f, 0, 1, 15, get_termination_condition(), &error, &L1); + BOOST_CHECK_CLOSE_FRACTION(Q, (Real) 17 + half()*third(), tol); + BOOST_CHECK_CLOSE_FRACTION(L1, (Real) 17 + half()*third(), tol); +} + +// Examples taken from +//http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf +template +void test_ca() +{ + std::cout << "Testing integration of C(a) on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = boost::math::tools::epsilon() * 10; + Real L1; + Real error; + + auto f1 = [](const Real& x) { return atan(x)/(x*(x*x + 1)) ; }; + Real Q = gauss_kronrod::integrate(f1, 0, 1, 15, get_termination_condition(), &error, &L1); + Real Q_expected = pi()*ln_two()/8 + catalan()*half(); + 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::integrate(f2, 0 , 1, 15, get_termination_condition(), &error, &L1); + Q_expected = pi()/4 - pi()/root_two() + 3*atan(root_two())/root_two(); + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); + BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); + + auto f5 = [](Real t)->Real { return t*t*log(t)/((t*t - 1)*(t*t*t*t + 1)); }; + Q = gauss_kronrod::integrate(f5, 0, 1, 25); + Q_expected = pi()*pi()*(2 - root_two())/32; + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 100 * tol); +} + +template +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().pretty_name() << "\n"; + Real tol = boost::math::tools::epsilon() * 10; + Real Q; + Real Q_expected; + + // Example 1: + auto f1 = [](const Real& t) { return t*boost::math::log1p(t); }; + Q = gauss_kronrod::integrate(f1, 0 , 1); + Q_expected = half()*half(); + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); + + + // Example 2: + auto f2 = [](const Real& t) { return t*t*atan(t); }; + Q = gauss_kronrod::integrate(f2, 0, 1); + Q_expected = (pi() -2 + 2*ln_two())/12; + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 2 * tol); + + // Example 3: + auto f3 = [](const Real& t) { return exp(t)*cos(t); }; + Q = gauss_kronrod::integrate(f3, 0, half_pi()); + Q_expected = boost::math::expm1(half_pi())*half(); + 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::integrate(f4, 0, 1); + Q_expected = 5*pi()*pi()/96; + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); +} + + +template +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().pretty_name() << "\n"; + Real tol = boost::math::tools::epsilon() * 10; + Real Q; + Real Q_expected; + Real L1; + Real error; + + auto f1 = [](const Real& t) { return 1/(1+t*t);}; + Q = gauss_kronrod::integrate(f1, -boost::math::tools::max_value(), boost::math::tools::max_value(), 15, get_termination_condition(), &error, &L1); + Q_expected = pi(); + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); + BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); + + auto f4 = [](const Real& t) { return 1/cosh(t);}; + Q = gauss_kronrod::integrate(f4, -boost::math::tools::max_value(), boost::math::tools::max_value(), 15, get_termination_condition(), &error, &L1); + Q_expected = pi(); + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); + BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); + +} + +template +void test_right_limit_infinite() +{ + std::cout << "Testing right limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = boost::math::tools::epsilon() * 10; + Real Q; + Real Q_expected; + Real L1; + Real error; + + // Example 11: + auto f1 = [](const Real& t) { return 1/(1+t*t);}; + Q = gauss_kronrod::integrate(f1, 0, boost::math::tools::max_value(), 15, get_termination_condition(), &error, &L1); + Q_expected = half_pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + auto f4 = [](const Real& t) { return 1/(1+t*t); }; + Q = gauss_kronrod::integrate(f4, 1, boost::math::tools::max_value(), 15, get_termination_condition(), &error, &L1); + Q_expected = pi()/4; + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); +} + +template +void test_left_limit_infinite() +{ + std::cout << "Testing left limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = boost::math::tools::epsilon() * 10; + Real Q; + Real Q_expected; + + // Example 11: + auto f1 = [](const Real& t) { return 1/(1+t*t);}; + Q = gauss_kronrod::integrate(f1, -boost::math::tools::max_value(), 0); + Q_expected = half_pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); +} + +BOOST_AUTO_TEST_CASE(gauss_quadrature_test) +{ + test_linear(); + test_quadratic(); + test_ca(); + test_three_quadrature_schemes_examples(); + test_integration_over_real_line(); + test_right_limit_infinite(); + test_left_limit_infinite(); + + test_linear(); + test_quadratic(); + test_ca(); + test_three_quadrature_schemes_examples(); + test_integration_over_real_line(); + test_right_limit_infinite(); + test_left_limit_infinite(); + + test_linear(); + test_quadratic(); + test_ca(); + test_three_quadrature_schemes_examples(); + test_integration_over_real_line(); + test_right_limit_infinite(); + test_left_limit_infinite(); + + test_linear(); + test_quadratic(); + test_ca(); + test_three_quadrature_schemes_examples(); + test_integration_over_real_line(); + test_right_limit_infinite(); + test_left_limit_infinite(); +} + +#else + +int main() { return 0; } + +#endif diff --git a/test/gauss_kronrod_quadrature_test.cpp b/test/gauss_kronrod_quadrature_test.cpp new file mode 100644 index 000000000..d5df066e4 --- /dev/null +++ b/test/gauss_kronrod_quadrature_test.cpp @@ -0,0 +1,377 @@ +// 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 +#include + +#if !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_TRAILING_RESULT_TYPES) && !defined(BOOST_NO_SFINAE_EXPR) + +#include +#include +#include +#include +#include +#include + +#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; + +// +// 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 +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<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 +void test_linear() +{ + std::cout << "Testing linear functions are integrated properly by gauss_kronrod on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = boost::math::tools::epsilon() * 10; + Real error; + auto f = [](const Real& x) + { + return 5*x + 7; + }; + Real L1; + Real Q = gauss_kronrod::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); +} + +template +void test_quadratic() +{ + std::cout << "Testing quadratic functions are integrated properly by tanh_sinh on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = boost::math::tools::epsilon() * 10; + Real error; + + auto f = [](const Real& x) { return 5*x*x + 7*x + 12; }; + Real L1; + Real Q = gauss_kronrod::integrate(f, 0, 1, 0, 0, &error, &L1); + BOOST_CHECK_CLOSE_FRACTION(Q, (Real) 17 + half()*third(), tol); + BOOST_CHECK_CLOSE_FRACTION(L1, (Real) 17 + half()*third(), tol); +} + +// Examples taken from +//http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf +template +void test_ca() +{ + std::cout << "Testing integration of C(a) on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = expected_error(test_ca_error_id); + Real L1; + Real error; + + auto f1 = [](const Real& x) { return atan(x)/(x*(x*x + 1)) ; }; + Real Q = gauss_kronrod::integrate(f1, 0, 1, 0, 0, &error, &L1); + Real Q_expected = pi()*ln_two()/8 + catalan()*half(); + 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::integrate(f2, 0 , 1, 0, 0, &error, &L1); + Q_expected = pi()/4 - pi()/root_two() + 3*atan(root_two())/root_two(); + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); + BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); + + tol = expected_error(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::integrate(f5, 0, 1, 0); + Q_expected = pi()*pi()*(2 - root_two())/32; + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); +} + +template +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().pretty_name() << "\n"; + Real tol = expected_error(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_kronrod::integrate(f1, 0 , 1, 0); + Q_expected = half()*half(); + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); + + + // Example 2: + auto f2 = [](const Real& t) { return t*t*atan(t); }; + Q = gauss_kronrod::integrate(f2, 0 , 1, 0); + Q_expected = (pi() -2 + 2*ln_two())/12; + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 2 * tol); + + // Example 3: + auto f3 = [](const Real& t) { return exp(t)*cos(t); }; + Q = gauss_kronrod::integrate(f3, 0, half_pi(), 0); + Q_expected = boost::math::expm1(half_pi())*half(); + 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::integrate(f4, 0 , 1, 0); + Q_expected = 5*pi()*pi()/96; + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); + + tol = expected_error(test_three_quad_error_id_2); + // Example 5: + auto f5 = [](const Real& t) { return sqrt(t)*log(t); }; + Q = gauss_kronrod::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) { return sqrt(1 - t*t); }; + Q = gauss_kronrod::integrate(f6, 0 , 1, 0); + Q_expected = pi()/4; + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); +} + + +template +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().pretty_name() << "\n"; + Real tol = expected_error(test_integration_over_real_line_error_id); + Real Q; + Real Q_expected; + Real L1; + Real error; + + auto f1 = [](const Real& t) { return 1/(1+t*t);}; + Q = gauss_kronrod::integrate(f1, -boost::math::tools::max_value(), boost::math::tools::max_value(), 0, 0, &error, &L1); + Q_expected = pi(); + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); + BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); + + auto f4 = [](const Real& t) { return 1/cosh(t);}; + Q = gauss_kronrod::integrate(f4, -boost::math::tools::max_value(), boost::math::tools::max_value(), 0, 0, &error, &L1); + Q_expected = pi(); + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); + BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); + +} + +template +void test_right_limit_infinite() +{ + std::cout << "Testing right limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = expected_error(test_right_limit_infinite_error_id); + Real Q; + Real Q_expected; + Real L1; + Real error; + + // Example 11: + auto f1 = [](const Real& t) { return 1/(1+t*t);}; + Q = gauss_kronrod::integrate(f1, 0, boost::math::tools::max_value(), 0, 0, &error, &L1); + Q_expected = half_pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + auto f4 = [](const Real& t) { return 1/(1+t*t); }; + Q = gauss_kronrod::integrate(f4, 1, boost::math::tools::max_value(), 0, 0, &error, &L1); + Q_expected = pi()/4; + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); +} + +template +void test_left_limit_infinite() +{ + std::cout << "Testing left limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = expected_error(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_kronrod::integrate(f1, -boost::math::tools::max_value(), Real(0), 0); + Q_expected = half_pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); +} + +BOOST_AUTO_TEST_CASE(gauss_quadrature_test) +{ + test_linear(); + test_quadratic(); + test_ca(); + test_three_quadrature_schemes_examples(); + test_integration_over_real_line(); + test_right_limit_infinite(); + test_left_limit_infinite(); + + test_linear(); + test_quadratic(); + test_ca(); + test_three_quadrature_schemes_examples(); + test_integration_over_real_line(); + test_right_limit_infinite(); + test_left_limit_infinite(); + + test_linear(); + test_quadratic(); + test_ca(); + test_three_quadrature_schemes_examples(); + test_integration_over_real_line(); + test_right_limit_infinite(); + test_left_limit_infinite(); + + test_linear(); + test_quadratic(); + test_ca(); + test_three_quadrature_schemes_examples(); + test_integration_over_real_line(); + test_right_limit_infinite(); + test_left_limit_infinite(); +} + +#else + +int main() { return 0; } + +#endif diff --git a/test/gauss_quadrature_test.cpp b/test/gauss_quadrature_test.cpp new file mode 100644 index 000000000..edff0dec0 --- /dev/null +++ b/test/gauss_quadrature_test.cpp @@ -0,0 +1,372 @@ +// 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 +#include + +#if !defined(BOOST_NO_CXX11_DECLTYPE) && !defined(BOOST_NO_CXX11_TRAILING_RESULT_TYPES) && !defined(BOOST_NO_SFINAE_EXPR) + +#include +#include +#include +#include +#include +#include + +#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; +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 +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<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 +void test_linear() +{ + std::cout << "Testing linear functions are integrated properly by gauss on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = boost::math::tools::epsilon() * 10; + auto f = [](const Real& x) + { + return 5*x + 7; + }; + Real L1; + Real Q = gauss::integrate(f, (Real) 0, (Real) 1, &L1); + BOOST_CHECK_CLOSE_FRACTION(Q, 9.5, tol); + BOOST_CHECK_CLOSE_FRACTION(L1, 9.5, tol); +} + +template +void test_quadratic() +{ + std::cout << "Testing quadratic functions are integrated properly by tanh_sinh on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = boost::math::tools::epsilon() * 10; + + auto f = [](const Real& x) { return 5*x*x + 7*x + 12; }; + Real L1; + Real Q = gauss::integrate(f, 0, 1, &L1); + BOOST_CHECK_CLOSE_FRACTION(Q, (Real) 17 + half()*third(), tol); + BOOST_CHECK_CLOSE_FRACTION(L1, (Real) 17 + half()*third(), tol); +} + +// Examples taken from +//http://crd-legacy.lbl.gov/~dhbailey/dhbpapers/quadrature.pdf +template +void test_ca() +{ + std::cout << "Testing integration of C(a) on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = expected_error(test_ca_error_id); + Real L1; + + auto f1 = [](const Real& x) { return atan(x)/(x*(x*x + 1)) ; }; + Real Q = gauss::integrate(f1, 0, 1, &L1); + Real Q_expected = pi()*ln_two()/8 + catalan()*half(); + 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::integrate(f2, 0 , 1, &L1); + Q_expected = pi()/4 - pi()/root_two() + 3*atan(root_two())/root_two(); + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); + BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); + + tol = expected_error(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::integrate(f5, 0 , 1); + Q_expected = pi()*pi()*(2 - root_two())/32; + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); +} + +template +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().pretty_name() << "\n"; + Real tol = expected_error(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::integrate(f1, 0 , 1); + Q_expected = half()*half(); + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); + + + // Example 2: + auto f2 = [](const Real& t) { return t*t*atan(t); }; + Q = gauss::integrate(f2, 0 , 1); + Q_expected = (pi() -2 + 2*ln_two())/12; + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 2 * tol); + + // Example 3: + auto f3 = [](const Real& t) { return exp(t)*cos(t); }; + Q = gauss::integrate(f3, 0, half_pi()); + Q_expected = boost::math::expm1(half_pi())*half(); + 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::integrate(f4, 0 , 1); + Q_expected = 5*pi()*pi()/96; + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); + + tol = expected_error(test_three_quad_error_id_2); + // Example 5: + auto f5 = [](const Real& t) { return sqrt(t)*log(t); }; + Q = gauss::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::integrate(f6, 0 , 1); + Q_expected = pi()/4; + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); +} + + +template +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().pretty_name() << "\n"; + Real tol = expected_error(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::integrate(f1, -boost::math::tools::max_value(), boost::math::tools::max_value(), &L1); + Q_expected = pi(); + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); + BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); + + auto f4 = [](const Real& t) { return 1/cosh(t);}; + Q = gauss::integrate(f4, -boost::math::tools::max_value(), boost::math::tools::max_value(), &L1); + Q_expected = pi(); + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); + BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); + +} + +template +void test_right_limit_infinite() +{ + std::cout << "Testing right limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = expected_error(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::integrate(f1, 0, boost::math::tools::max_value(), &L1); + Q_expected = half_pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); + + auto f4 = [](const Real& t) { return 1/(1+t*t); }; + Q = gauss::integrate(f4, 1, boost::math::tools::max_value(), &L1); + Q_expected = pi()/4; + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); +} + +template +void test_left_limit_infinite() +{ + std::cout << "Testing left limit infinite for tanh_sinh in 'A Comparison of Three High Precision Quadrature Schemes' on type " << boost::typeindex::type_id().pretty_name() << "\n"; + Real tol = expected_error(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::integrate(f1, -boost::math::tools::max_value(), Real(0)); + Q_expected = half_pi(); + BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); +} + +BOOST_AUTO_TEST_CASE(gauss_quadrature_test) +{ + test_linear(); + test_quadratic(); + test_ca(); + test_three_quadrature_schemes_examples(); + test_integration_over_real_line(); + test_right_limit_infinite(); + test_left_limit_infinite(); + + test_linear(); + test_quadratic(); + test_ca(); + test_three_quadrature_schemes_examples(); + test_integration_over_real_line(); + test_right_limit_infinite(); + test_left_limit_infinite(); + + test_linear(); + test_quadratic(); + test_ca(); + test_three_quadrature_schemes_examples(); + test_integration_over_real_line(); + test_right_limit_infinite(); + test_left_limit_infinite(); + + test_linear(); + test_quadratic(); + test_ca(); + test_three_quadrature_schemes_examples(); + test_integration_over_real_line(); + test_right_limit_infinite(); + test_left_limit_infinite(); +} + +#else + +int main() { return 0; } + +#endif