diff --git a/include/boost/math/concepts/real_concept.hpp b/include/boost/math/concepts/real_concept.hpp index 3a759955a..d43a3c50e 100644 --- a/include/boost/math/concepts/real_concept.hpp +++ b/include/boost/math/concepts/real_concept.hpp @@ -31,6 +31,8 @@ #include #include #include +#include +#include #if defined(__SGI_STL_PORT) # include #endif @@ -255,6 +257,24 @@ inline real_concept sqrt(real_concept a) inline real_concept tanh(real_concept a) { return std::tanh(a.value()); } +// +// C++11 ism's +// Note that these must not actually call the std:: versions as that precludes using this +// header to test in C++03 mode, call the Boost versions instead: +// +inline boost::math::concepts::real_concept asinh(boost::math::concepts::real_concept a) +{ + return boost::math::asinh(a.value(), boost::math::policies::make_policy(boost::math::policies::overflow_error())); +} +inline boost::math::concepts::real_concept acosh(boost::math::concepts::real_concept a) +{ + return boost::math::acosh(a.value(), boost::math::policies::make_policy(boost::math::policies::overflow_error())); +} +inline boost::math::concepts::real_concept atanh(boost::math::concepts::real_concept a) +{ + return boost::math::atanh(a.value(), boost::math::policies::make_policy(boost::math::policies::overflow_error())); +} + // // Conversion and truncation routines: // diff --git a/include/boost/math/concepts/std_real_concept.hpp b/include/boost/math/concepts/std_real_concept.hpp index d0eb1d89e..d855d2d55 100644 --- a/include/boost/math/concepts/std_real_concept.hpp +++ b/include/boost/math/concepts/std_real_concept.hpp @@ -226,12 +226,17 @@ inline boost::math::concepts::std_real_concept sqrt(boost::math::concepts::std_r { return std::sqrt(a.value()); } inline boost::math::concepts::std_real_concept tanh(boost::math::concepts::std_real_concept a) { return std::tanh(a.value()); } -inline bool isfinite(boost::math::concepts::std_real_concept a) -{ return std::isfinite(a.value()); } +// +// C++11 ism's +// Note that these must not actually call the std:: versions as that precludes using this +// header to test in C++03 mode, call the Boost versions instead: +// inline boost::math::concepts::std_real_concept asinh(boost::math::concepts::std_real_concept a) -{ return std::asinh(a.value()); } +{ return boost::math::asinh(a.value(), boost::math::policies::make_policy(boost::math::policies::overflow_error())); } +inline boost::math::concepts::std_real_concept acosh(boost::math::concepts::std_real_concept a) +{ return boost::math::acosh(a.value(), boost::math::policies::make_policy(boost::math::policies::overflow_error())); } inline boost::math::concepts::std_real_concept atanh(boost::math::concepts::std_real_concept a) -{ return std::atanh(a.value()); } +{ return boost::math::atanh(a.value(), boost::math::policies::make_policy(boost::math::policies::overflow_error())); } } // namespace std @@ -405,4 +410,12 @@ using concepts::llround; } // namespace math } // namespace boost +// +// These must go at the end, as they include stuff that won't compile until +// after std_real_concept has been defined: +// +#include +#include +#include + #endif // BOOST_MATH_STD_REAL_CONCEPT_HPP diff --git a/include/boost/math/quadrature/detail/exp_sinh_detail.hpp b/include/boost/math/quadrature/detail/exp_sinh_detail.hpp index 510572a5e..13c305da9 100644 --- a/include/boost/math/quadrature/detail/exp_sinh_detail.hpp +++ b/include/boost/math/quadrature/detail/exp_sinh_detail.hpp @@ -68,7 +68,7 @@ exp_sinh_detail::exp_sinh_detail(Real tol, size_t max_refinements) for (size_t i = 0; i < max_refinements; ++i) { Real h = (Real) 1/ (Real) (1<::sinh_sinh_detail(Real tol, size_t max_refinement for (size_t i = 0; i < max_refinements; ++i) { Real h = (Real) 1/ (Real) (1<::integrate(const F f, Real* error, Real* L1) "The function you are trying to integrate does not go to zero at infinity, and instead evaluates to %1%", y_max, Policy()); } - Real y_min = f(std::numeric_limits::lowest()); + Real y_min = f(-boost::math::tools::max_value()); if(abs(y_min) > boost::math::tools::epsilon() || !(boost::math::isfinite)(y_min)) { return policies::raise_domain_error(function, @@ -151,7 +151,7 @@ Real sinh_sinh_detail::integrate(const F f, Real* error, Real* L1) Real absum = 0; Real abterm1 = 1; - Real eps = std::numeric_limits::epsilon()*L1_I1; + Real eps = boost::math::tools::epsilon()*L1_I1; for(size_t j = 0; j < m_weights[i].size(); ++j) { Real x = m_abscissas[i][j]; diff --git a/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp b/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp index 062671483..5650beebf 100644 --- a/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp +++ b/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp @@ -12,6 +12,7 @@ #include #include #include +#include namespace boost{ namespace math{ namespace quadrature { namespace detail{ diff --git a/include/boost/math/quadrature/tanh_sinh.hpp b/include/boost/math/quadrature/tanh_sinh.hpp index 36b6eca74..a691980d8 100644 --- a/include/boost/math/quadrature/tanh_sinh.hpp +++ b/include/boost/math/quadrature/tanh_sinh.hpp @@ -59,57 +59,60 @@ Real tanh_sinh::integrate(const F f, Real a, Real b, Real* error, static const char* function = "tanh_sinh<%1%>::integrate"; - if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b)) + if (!(boost::math::isnan)(a) && !(boost::math::isnan)(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 diff = (b-a)*half(); - auto u = [&](Real z) { return f(avg + diff*z); }; - Real Q = diff*m_imp->integrate(u, error, L1, function); - if(L1) - { - *L1 *= diff; - } - return Q; + // Infinite limits: + if ((a <= -tools::max_value()) && (b >= tools::max_value())) + { + auto u = [&](Real t) { auto t_sq = t*t; auto inv = 1 / (1 - t_sq); return f(t*inv)*(1 + t_sq)*inv*inv; }; + return m_imp->integrate(u, error, L1, function); + } + + // Right limit is infinite: + if ((boost::math::isfinite)(a) && (b >= tools::max_value())) + { + auto u = [&](Real t) { auto z = 1 / (t + 1); auto arg = 2 * z + a - 1; return f(arg)*z*z; }; + Real Q = 2 * m_imp->integrate(u, error, L1, function); + if (L1) + { + *L1 *= 2; + } + + return Q; + } + + if ((boost::math::isfinite)(b) && (a <= -tools::max_value())) + { + auto u = [&](Real t) { return f(b - t); }; + auto v = [&](Real t) { auto z = 1 / (t + 1); auto arg = 2 * z - 1; return u(arg)*z*z; }; + + Real Q = 2 * m_imp->integrate(v, error, L1, function); + if (L1) + { + *L1 *= 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 diff = (b - a)*half(); + auto u = [&](Real z) { return f(avg + diff*z); }; + Real Q = diff*m_imp->integrate(u, error, L1, function); + + if (L1) + { + *L1 *= diff; + } + return Q; + } } - - // Infinite limits: - if ((a <= -tools::max_value()) && (b >= tools::max_value())) - { - auto u = [&](Real t) { auto t_sq = t*t; auto inv = 1/(1 - t_sq); return f(t*inv)*(1+t_sq)*inv*inv; }; - return m_imp->integrate(u, error, L1, function); - } - - // Right limit is infinite: - if ((boost::math::isfinite)(a) && (b >= tools::max_value())) - { - auto u = [&](Real t) { auto z = 1/(t+1); auto arg = 2*z + a - 1; return f(arg)*z*z; }; - Real Q = 2*m_imp->integrate(u, error, L1, function); - if(L1) - { - *L1 *= 2; - } - - return Q; - } - - if ((boost::math::isfinite)(b) && (a <= -tools::max_value())) - { - auto u = [&](Real t) { return f(b-t);}; - auto v = [&](Real t) { auto z = 1/(t+1); auto arg = 2*z - 1; return u(arg)*z*z; }; - - Real Q = 2*m_imp->integrate(v, error, L1, function); - if (L1) - { - *L1 *= 2; - } - return Q; - } - return policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy()); } diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 2ca9e6ce4..1ec2b3c40 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -769,23 +769,33 @@ run test_uniform.cpp pch ../../test/build//boost_unit_test_framework ; run test_weibull.cpp ../../test/build//boost_unit_test_framework ; run test_zeta.cpp ../../test/build//boost_unit_test_framework test_instances//test_instances pch_light ; run tanh_sinh_quadrature_test.cpp ../../test/build//boost_unit_test_framework - : : : TEST1 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] : + : : : TEST1 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr ] : tanh_sinh_quadrature_test_1 ; run tanh_sinh_quadrature_test.cpp ../../test/build//boost_unit_test_framework - : : : TEST2 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] : + : : : TEST2 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr ] : tanh_sinh_quadrature_test_2 ; run tanh_sinh_quadrature_test.cpp ../../test/build//boost_unit_test_framework - : : : TEST3 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] : + : : : TEST3 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr ] : tanh_sinh_quadrature_test_3 ; run tanh_sinh_quadrature_test.cpp ../../test/build//boost_unit_test_framework - : : : release TEST4 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] : + : : : release TEST4 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr ] : tanh_sinh_quadrature_test_4 ; run tanh_sinh_quadrature_test.cpp ../../test/build//boost_unit_test_framework - : : : release TEST5 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] : + : : : release TEST5 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr ] : tanh_sinh_quadrature_test_5 ; run tanh_sinh_quadrature_test.cpp ../../test/build//boost_unit_test_framework - : : : TEST6 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] : + : : : TEST6 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr ] : tanh_sinh_quadrature_test_6 ; +run sinh_sinh_quadrature_test.cpp ../../test/build//boost_unit_test_framework + : : : [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr ] ; +run exp_sinh_quadrature_test.cpp ../../test/build//boost_unit_test_framework + : : : TEST1 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr ] : exp_sinh_quadrature_test_1 ; +run exp_sinh_quadrature_test.cpp ../../test/build//boost_unit_test_framework + : : : release TEST2 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr ] : exp_sinh_quadrature_test_2 ; +run exp_sinh_quadrature_test.cpp ../../test/build//boost_unit_test_framework + : : : TEST3 [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr ] : exp_sinh_quadrature_test_3 ; + + run test_policy.cpp ../../test/build//boost_unit_test_framework ; run test_policy_2.cpp ../../test/build//boost_unit_test_framework ; @@ -886,9 +896,9 @@ run compile_test/dist_triangular_incl_test.cpp compile_test_main ; run compile_test/dist_uniform_incl_test.cpp compile_test_main ; run compile_test/dist_weibull_incl_test.cpp compile_test_main ; run compile_test/distribution_concept_check.cpp ; -run compile_test/exp_sinh_incl_test.cpp compile_test_main ; -run compile_test/sinh_sinh_incl_test.cpp compile_test_main ; -run compile_test/tanh_sinh_incl_test.cpp compile_test_main ; +run compile_test/exp_sinh_incl_test.cpp compile_test_main : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr ] ; +run compile_test/sinh_sinh_incl_test.cpp compile_test_main : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr ] ; +run compile_test/tanh_sinh_incl_test.cpp compile_test_main : : : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr ] ; run compile_test/sf_beta_incl_test.cpp compile_test_main ; run compile_test/sf_bernoulli_incl_test.cpp compile_test_main ; run compile_test/sf_bessel_incl_test.cpp compile_test_main ; @@ -969,11 +979,11 @@ compile compile_test/tools_test_inc_test.cpp ; compile compile_test/tools_toms748_inc_test.cpp ; compile compile_test/cubic_spline_concept_test.cpp : [ requires cxx11_smart_ptr cxx11_defaulted_functions ] ; compile compile_test/barycentric_rational_concept_test.cpp : [ requires cxx11_smart_ptr cxx11_defaulted_functions ] ; -compile compile_test/sf_legendre_stieltjes_concept_test.cpp : [ requires cxx11_auto_declarations cxx11_defaulted_functions cxx11_lambdas ] ; +compile compile_test/sf_legendre_stieltjes_concept_test.cpp : [ requires cxx11_auto_declarations cxx11_defaulted_functions cxx11_lambdas cxx11_auto_declarations ] ; compile compile_test/trapezoidal_concept_test.cpp ; -compile compile_test/cubic_spline_concept_test.cpp : [ requires cxx11_smart_ptr ] ; -compile compile_test/barycentric_rational_concept_test.cpp : [ requires cxx11_smart_ptr ] ; -compile compile_test/sf_legendre_stieltjes_concept_test.cpp : [ requires cxx11_auto_declarations ] ; +compile compile_test/exp_sinh_concept_test.cpp : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr ] ; +compile compile_test/sinh_sinh_concept_test.cpp : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr ] ; +compile compile_test/tanh_sinh_concept_test.cpp : [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr ] ; run octonion_test.cpp ../../test/build//boost_unit_test_framework ; diff --git a/test/exp_sinh_quadrature_test.cpp b/test/exp_sinh_quadrature_test.cpp index 08d7ad5e6..d2bef59b8 100644 --- a/test/exp_sinh_quadrature_test.cpp +++ b/test/exp_sinh_quadrature_test.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -47,11 +48,17 @@ using boost::math::constants::root_two_pi; using boost::math::constants::root_pi; using boost::math::quadrature::exp_sinh; +#if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) +# define TEST1 +# define TEST2 +# define TEST3 +#endif + 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 = sqrt(std::numeric_limits::epsilon()); + Real tol = sqrt(boost::math::tools::epsilon()); Real Q; Real Q_expected; Real error; @@ -60,19 +67,19 @@ void test_right_limit_infinite() // Example 12 const auto f2 = [](Real t) { return exp(-t)/sqrt(t); }; - Q = integrator.integrate(f2, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f2, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = root_pi(); BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); // The integrand is strictly positive, so it coincides with the value of the integral: BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); auto f3 = [](Real t) { Real z = exp(-t); if (z == 0) { return z; } return z*cos(t); }; - Q = integrator.integrate(f3, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f3, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = half(); BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); auto f4 = [](Real t) { return 1/(1+t*t); }; - Q = integrator.integrate(f4, 1, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f4, 1, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = pi()/4; BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); @@ -82,7 +89,7 @@ template void test_left_limit_infinite() { std::cout << "Testing left limit infinite for 1/(1+t^2) on type " << boost::typeindex::type_id().pretty_name() << "\n"; - Real tol = sqrt(std::numeric_limits::epsilon()); + Real tol = sqrt(boost::math::tools::epsilon()); Real Q; Real Q_expected; Real error; @@ -91,7 +98,7 @@ void test_left_limit_infinite() // Example 11: auto f1 = [](Real t) { return 1/(1+t*t);}; - Q = integrator.integrate(f1, -std::numeric_limits::infinity(), 0, &error, &L1); + Q = integrator.integrate(f1, std::numeric_limits::has_infinity ? -std::numeric_limits::infinity() : -boost::math::tools::max_value(), 0, &error, &L1); Q_expected = half_pi(); BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); @@ -107,7 +114,7 @@ void test_nr_examples() using std::exp; using std::sqrt; std::cout << "Testing type " << boost::typeindex::type_id().pretty_name() << "\n"; - Real tol = sqrt(std::numeric_limits::epsilon()); + Real tol = sqrt(boost::math::tools::epsilon()); std::cout << std::setprecision(std::numeric_limits::digits10); Real Q; Real Q_expected; @@ -116,13 +123,13 @@ void test_nr_examples() exp_sinh integrator(tol, 12); auto f0 = [](Real) { return (Real) 0; }; - Q = integrator.integrate(f0, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f0, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = 0; BOOST_CHECK_CLOSE_FRACTION(Q, 0.0f, tol); BOOST_CHECK_CLOSE_FRACTION(L1, 0.0f, tol); auto f = [](Real x) { return 1/(1+x*x); }; - Q = integrator.integrate(f, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = half_pi(); BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); @@ -141,33 +148,33 @@ void test_nr_examples() return sin(x*half())*z2; }; - Q = integrator.integrate(f1, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f1, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = sqrt(pi()*(sqrt((Real) 5) - 2)); // The integrand is oscillatory; the accuracy is low. BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); auto f2 = [](Real x) { return pow(x, -(Real) 2/(Real) 7)*exp(-x*x); }; - Q = integrator.integrate(f2, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f2, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = half()*boost::math::tgamma((Real) 5/ (Real) 14); BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); auto f3 = [](Real x) { return (Real) 1/ (sqrt(x)*(1+x)); }; - Q = integrator.integrate(f3, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f3, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = pi(); - BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10*std::numeric_limits::epsilon()); - BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, 10*std::numeric_limits::epsilon()); + BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10*boost::math::tools::epsilon()); + BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, 10*boost::math::tools::epsilon()); auto f4 = [](Real t) { return exp(-t*t*half()); }; - Q = integrator.integrate(f4, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f4, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = root_two_pi()/2; BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); auto f5 = [](Real t) { return 1/cosh(t);}; - Q = integrator.integrate(f5, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f5, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = half_pi(); BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); @@ -183,7 +190,7 @@ void test_crc() using std::sqrt; using std::log; std::cout << "Testing integral from CRC handbook on type " << boost::typeindex::type_id().pretty_name() << "\n"; - Real tol = sqrt(std::numeric_limits::epsilon()); + Real tol = sqrt(boost::math::tools::epsilon()); std::cout << std::setprecision(std::numeric_limits::digits10); Real Q; Real Q_expected; @@ -192,7 +199,7 @@ void test_crc() exp_sinh integrator(tol, 14); auto f0 = [](Real x) { return log(x)*exp(-x); }; - Q = integrator.integrate(f0, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f0, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = -boost::math::constants::euler(); BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); @@ -205,7 +212,7 @@ void test_crc() return pow(t, (Real) 12 - 1)*x; }; - Q = integrator.integrate(f1, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f1, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = boost::math::tgamma(12.0f); BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); @@ -219,7 +226,7 @@ void test_crc() } return x*cosh(5*t); }; - Q = integrator.integrate(f2, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f2, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = boost::math::cyl_bessel_k(5, 12); BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); // Laplace transform of cos(at) @@ -235,7 +242,7 @@ void test_crc() }; // For high oscillation frequency, the quadrature sum is ill-conditioned. - Q = integrator.integrate(f3, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f3, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = s/(a*a+s*s); // Since the integrand is oscillatory, we increase the tolerance: BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 100*tol); @@ -250,12 +257,12 @@ void test_crc() return boost::math::cyl_bessel_j(0, t)*x; }; - Q = integrator.integrate(f4, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f4, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = 1/sqrt(1+s*s); BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); auto f6 = [](Real t) { return exp(-t*t)*log(t);}; - Q = integrator.integrate(f6, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f6, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = -boost::math::constants::root_pi()*(boost::math::constants::euler() + 2*ln_two())/4; BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); } @@ -263,23 +270,34 @@ void test_crc() BOOST_AUTO_TEST_CASE(exp_sinh_quadrature_test) { - test_left_limit_infinite(); +#ifdef TEST1 + + test_left_limit_infinite(); test_left_limit_infinite(); test_left_limit_infinite(); - test_left_limit_infinite(); - test_right_limit_infinite(); test_right_limit_infinite(); test_right_limit_infinite(); - test_right_limit_infinite(); - test_nr_examples(); test_nr_examples(); test_nr_examples(); - //test_nr_examples(); - test_crc(); test_crc(); test_crc(); - //test_crc(); + +#endif +#ifdef TEST2 + + test_left_limit_infinite(); + test_right_limit_infinite(); + test_nr_examples(); + test_crc(); +#endif + +#ifdef TEST3 + test_left_limit_infinite(); + test_right_limit_infinite(); + test_nr_examples(); + test_crc(); +#endif } diff --git a/test/sinh_sinh_quadrature_test.cpp b/test/sinh_sinh_quadrature_test.cpp index 4754edc31..524e3e8d4 100644 --- a/test/sinh_sinh_quadrature_test.cpp +++ b/test/sinh_sinh_quadrature_test.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -56,7 +57,7 @@ template void test_nr_examples() { std::cout << "Testing type " << boost::typeindex::type_id().pretty_name() << "\n"; - Real tol = sqrt(std::numeric_limits::epsilon()); + Real tol = sqrt(boost::math::tools::epsilon()); std::cout << std::setprecision(std::numeric_limits::digits10); Real Q; Real Q_expected; @@ -102,7 +103,7 @@ template void test_crc() { std::cout << "Testing CRC formulas on type " << boost::typeindex::type_id().pretty_name() << "\n"; - Real tol = sqrt(std::numeric_limits::epsilon()); + Real tol = sqrt(boost::math::tools::epsilon()); std::cout << std::setprecision(std::numeric_limits::digits10); Real Q; Real Q_expected; @@ -141,10 +142,12 @@ BOOST_AUTO_TEST_CASE(sinh_sinh_quadrature_test) test_nr_examples(); test_nr_examples(); test_nr_examples(); - //test_nr_examples(); + test_nr_examples(); + test_nr_examples(); test_crc(); test_crc(); test_crc(); - //test_crc(); + test_crc(); + test_crc(); } diff --git a/test/tanh_sinh_quadrature_test.cpp b/test/tanh_sinh_quadrature_test.cpp index a84114a02..61b09ba27 100644 --- a/test/tanh_sinh_quadrature_test.cpp +++ b/test/tanh_sinh_quadrature_test.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -24,6 +25,10 @@ #include #include +#ifdef _MSC_VER +#pragma warning(disable:4127) // Conditional expression is constant +#endif + #if !defined(TEST1) && !defined(TEST2) && !defined(TEST3) && !defined(TEST4) && !defined(TEST5) && !defined(TEST6) # define TEST1 # define TEST2 @@ -157,7 +162,7 @@ template void test_detail() { std::cout << "Testing tanh_sinh_detail on type " << boost::typeindex::type_id().pretty_name() << "\n"; - Real tol = sqrt(std::numeric_limits::epsilon()); + Real tol = sqrt(boost::math::tools::epsilon()); tanh_sinh_detail > integrator(tol, 20); auto f = [](Real x) { return x*x; }; Real err; @@ -172,7 +177,7 @@ template void test_linear() { std::cout << "Testing linear functions are integrated properly by tanh_sinh on type " << boost::typeindex::type_id().pretty_name() << "\n"; - Real tol = 100*std::numeric_limits::epsilon(); + Real tol = 100*boost::math::tools::epsilon(); tanh_sinh integrator(tol, 20); auto f = [](Real x) { return 5*x + 7; }; Real error; @@ -187,7 +192,7 @@ 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 = 100*std::numeric_limits::epsilon(); + Real tol = 100*boost::math::tools::epsilon(); tanh_sinh integrator(tol, 20); auto f = [](Real x) { return 5*x*x + 7*x + 12; }; Real error; @@ -202,7 +207,7 @@ template void test_singular() { std::cout << "Testing integration of endpoint singularities on type " << boost::typeindex::type_id().pretty_name() << "\n"; - Real tol = 100*std::numeric_limits::epsilon(); + Real tol = 100*boost::math::tools::epsilon(); Real error; Real L1; tanh_sinh integrator(tol, 20); @@ -221,7 +226,7 @@ template void test_ca() { std::cout << "Testing integration of C(a) on type " << boost::typeindex::type_id().pretty_name() << "\n"; - Real tol = sqrt(std::numeric_limits::epsilon()); + Real tol = sqrt(boost::math::tools::epsilon()); Real error; Real L1; @@ -272,13 +277,13 @@ 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 = sqrt(std::numeric_limits::epsilon()); + Real tol = sqrt(boost::math::tools::epsilon()); Real Q; Real Q_expected; tanh_sinh integrator(tol, 20); // Example 1: - auto f1 = [](Real t) { return t*log1p(t); }; + auto f1 = [](Real t) { return t*boost::math::log1p(t); }; Q = integrator.integrate(f1, (Real) 0 , (Real) 1); Q_expected = half()*half(); BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); @@ -293,7 +298,7 @@ void test_three_quadrature_schemes_examples() // Example 3: auto f3 = [](Real t) { return exp(t)*cos(t); }; Q = integrator.integrate(f3, (Real) 0, half_pi()); - Q_expected = expm1(half_pi())*half(); + Q_expected = boost::math::expm1(half_pi())*half(); BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); // Example 4: @@ -320,7 +325,7 @@ 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 = sqrt(std::numeric_limits::epsilon()); + Real tol = sqrt(boost::math::tools::epsilon()); Real Q; Real Q_expected; Real error; @@ -328,13 +333,13 @@ void test_integration_over_real_line() tanh_sinh integrator(tol, 20); auto f1 = [](Real t) { return 1/(1+t*t);}; - Q = integrator.integrate(f1, -std::numeric_limits::infinity(), std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f1, std::numeric_limits::has_infinity ? -std::numeric_limits::infinity() : -boost::math::tools::max_value(), std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = pi(); BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); auto f2 = [](Real t) { return exp(-t*t*half()); }; - Q = integrator.integrate(f2, -std::numeric_limits::infinity(), std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f2, std::numeric_limits::has_infinity ? -std::numeric_limits::infinity() : -boost::math::tools::max_value(), std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = root_two_pi(); BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); @@ -348,7 +353,7 @@ void test_integration_over_real_line() //BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); auto f4 = [](Real t) { return 1/cosh(t);}; - Q = integrator.integrate(f4, -std::numeric_limits::infinity(), std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f4, std::numeric_limits::has_infinity ? -std::numeric_limits::infinity() : -boost::math::tools::max_value(), std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = pi(); BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); BOOST_CHECK_CLOSE_FRACTION(L1, Q_expected, tol); @@ -359,7 +364,7 @@ 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 = sqrt(std::numeric_limits::epsilon()); + Real tol = sqrt(boost::math::tools::epsilon()); Real Q; Real Q_expected; Real error; @@ -368,23 +373,23 @@ void test_right_limit_infinite() // Example 11: auto f1 = [](Real t) { return 1/(1+t*t);}; - Q = integrator.integrate(f1, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f1, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = half_pi(); BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); // Example 12 auto f2 = [](Real t) { return exp(-t)/sqrt(t); }; - Q = integrator.integrate(f2, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f2, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = root_pi(); BOOST_CHECK_CLOSE(Q, Q_expected, 1000*tol); auto f3 = [](Real t) { return exp(-t)*cos(t); }; - Q = integrator.integrate(f3, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f3, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = half(); BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); auto f4 = [](Real t) { return 1/(1+t*t); }; - Q = integrator.integrate(f4, 1, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f4, 1, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = pi()/4; BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); } @@ -393,14 +398,14 @@ 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 = sqrt(std::numeric_limits::epsilon()); + Real tol = sqrt(boost::math::tools::epsilon()); Real Q; Real Q_expected; tanh_sinh integrator; // Example 11: auto f1 = [](Real t) { return 1/(1+t*t);}; - Q = integrator.integrate(f1, -std::numeric_limits::infinity(), 0); + Q = integrator.integrate(f1, std::numeric_limits::has_infinity ? -std::numeric_limits::infinity() : -boost::math::tools::max_value(), 0); Q_expected = half_pi(); BOOST_CHECK_CLOSE(Q, Q_expected, 100*tol); } @@ -432,7 +437,7 @@ template void test_nr_examples() { std::cout << "Testing singular integrals from NR 4.5.4 on type " << boost::typeindex::type_id().pretty_name() << "\n"; - Real tol = sqrt(std::numeric_limits::epsilon()); + Real tol = sqrt(boost::math::tools::epsilon()); Real Q; Real Q_expected; Real error; @@ -440,12 +445,12 @@ void test_nr_examples() tanh_sinh integrator(tol, 15); auto f1 = [](Real x) { return sin(x*half())*pow(x, -3*half())*exp(-x); }; - Q = integrator.integrate(f1, 0, std::numeric_limits::infinity(), &error, &L1); + Q = integrator.integrate(f1, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value(), &error, &L1); Q_expected = sqrt(pi()*(sqrt((Real) 5) - 2)); BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, 10*tol); auto f2 = [](Real x) { return pow(x, -(Real) 2/(Real) 7)*exp(-x*x); }; - Q = integrator.integrate(f2, 0, std::numeric_limits::infinity()); + Q = integrator.integrate(f2, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value()); Q_expected = half()*boost::math::tgamma((Real) 5/ (Real) 14); BOOST_CHECK_CLOSE_FRACTION(Q, Q_expected, tol); @@ -456,7 +461,7 @@ template void test_early_termination() { std::cout << "Testing Clenshaw & Curtis's example of integrand which fools termination schemes on type " << boost::typeindex::type_id().pretty_name() << "\n"; - Real tol = sqrt(std::numeric_limits::epsilon()); + Real tol = sqrt(boost::math::tools::epsilon()); Real Q; Real Q_expected; Real error; @@ -475,7 +480,7 @@ template void test_crc() { std::cout << "Testing CRC formulas on type " << boost::typeindex::type_id().pretty_name() << "\n"; - Real tol = sqrt(std::numeric_limits::epsilon()); + Real tol = sqrt(boost::math::tools::epsilon()); Real Q; Real Q_expected; Real error; @@ -501,7 +506,7 @@ void test_sf() using std::sqrt; // Test some special functions that we already know how to evaluate: std::cout << "Testing special functions on type " << boost::typeindex::type_id().pretty_name() << "\n"; - Real tol = sqrt(std::numeric_limits::epsilon()); + Real tol = sqrt(boost::math::tools::epsilon()); tanh_sinh integrator; // incomplete beta: if (std::numeric_limits::digits < 120) // Otherwise too slow @@ -510,15 +515,15 @@ void test_sf() } Real x = 2, y = 3, z = 0.5, p = 0.25; - // This one has much larger (10^-9) error than exp_sinh quadrature: - //BOOST_CHECK_CLOSE_FRACTION(integrator.integrate([&](Real t) { return 1 / (sqrt(t + x) * (t + y)); }, 0, std::numeric_limits::infinity()) / 2, boost::math::ellint_rc(x, y), tol); + // This one has a higher than expected error rate, even though exp_sinh quadrature is just fine: + //BOOST_CHECK_CLOSE_FRACTION(integrator.integrate([&](Real t) { return 1 / (sqrt(t + x) * (t + y)); }, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value()) / 2, boost::math::ellint_rc(x, y), tol); - BOOST_CHECK_CLOSE_FRACTION((Real(3) / 2) * integrator.integrate([&](Real t) { return 1 / (sqrt((t + x) * (t + y) * (t + z)) * (t + p)); }, 0, std::numeric_limits::infinity()), boost::math::ellint_rj(x, y, z, p), tol); + BOOST_CHECK_CLOSE_FRACTION((Real(3) / 2) * integrator.integrate([&](Real t) { return 1 / (sqrt((t + x) * (t + y) * (t + z)) * (t + p)); }, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value()), boost::math::ellint_rj(x, y, z, p), tol); z = 5.5; - BOOST_CHECK_CLOSE_FRACTION(integrator.integrate([&](Real t) { using std::pow; using std::exp; return t > 10000 ? 0 : pow(t, z - 1) * exp(-t); }, 0, std::numeric_limits::infinity()), + BOOST_CHECK_CLOSE_FRACTION(integrator.integrate([&](Real t) { using std::pow; using std::exp; return t > 10000 ? 0 : pow(t, z - 1) * exp(-t); }, 0, std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value()), boost::math::tgamma(z), tol); - BOOST_CHECK_CLOSE_FRACTION(integrator.integrate([](Real t) { using std::exp; return exp(-t*t); }, -std::numeric_limits::infinity(), std::numeric_limits::infinity()), + BOOST_CHECK_CLOSE_FRACTION(integrator.integrate([](Real t) { using std::exp; return exp(-t*t); }, std::numeric_limits::has_infinity ? -std::numeric_limits::infinity() : -boost::math::tools::max_value(), std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : boost::math::tools::max_value()), boost::math::constants::root_pi(), tol); } @@ -603,8 +608,20 @@ BOOST_AUTO_TEST_CASE(tanh_sinh_quadrature_test) #endif #ifdef TEST6 - //test_linear(); - //test_quadratic(); + test_detail(); + test_right_limit_infinite(); + test_left_limit_infinite(); + test_linear(); + test_quadratic(); + test_singular(); + test_ca(); + test_three_quadrature_schemes_examples(); + test_horrible(); + test_integration_over_real_line(); + test_nr_examples(); + test_early_termination(); + test_crc(); + test_sf(); #endif }