From 9ca40b9f016bf8f4a7211a8446492ed2aeb06e6c Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Wed, 6 Mar 2024 18:47:48 +0000 Subject: [PATCH 01/27] Improve erf/expm1/expint coverage. In the expm1 case, tighten up error handling and testing. --- include/boost/math/special_functions/erf.hpp | 2 +- .../boost/math/special_functions/expint.hpp | 110 ++-------------- .../boost/math/special_functions/expm1.hpp | 124 ++++++++++++------ test/log1p_expm1_test.hpp | 4 + test/test_expint.cpp | 2 +- test/test_expint.hpp | 22 ++++ 6 files changed, 129 insertions(+), 135 deletions(-) diff --git a/include/boost/math/special_functions/erf.hpp b/include/boost/math/special_functions/erf.hpp index 57ff605299..54964617bf 100644 --- a/include/boost/math/special_functions/erf.hpp +++ b/include/boost/math/special_functions/erf.hpp @@ -690,7 +690,7 @@ T erf_imp(T z, bool invert, const Policy& pol, const std::integral_constant(%1%)", "Expected a finite argument but got %1%", z, pol); + return policies::raise_domain_error("boost::math::erf<%1%>(%1%)", "Expected a finite argument but got %1%", z, pol); // LCOV_EXCL_LINE checked as covered. if(z < 0) { diff --git a/include/boost/math/special_functions/expint.hpp b/include/boost/math/special_functions/expint.hpp index 1475a9a88b..9b2b420552 100644 --- a/include/boost/math/special_functions/expint.hpp +++ b/include/boost/math/special_functions/expint.hpp @@ -58,6 +58,7 @@ T expint_1_rational(const T& z, const std::integral_constant&) // Maximum Deviation Found: 2.006e-18 // Expected Error Term: 2.006e-18 // Max error found at double precision: 2.760e-17 + // LCOV_EXCL_START static const T Y = 0.66373538970947265625F; static const T P[6] = { BOOST_MATH_BIG_CONSTANT(T, 53, 0.0865197248079397976498), @@ -75,6 +76,7 @@ T expint_1_rational(const T& z, const std::integral_constant&) BOOST_MATH_BIG_CONSTANT(T, 53, 0.000131049900798434683324), BOOST_MATH_BIG_CONSTANT(T, 53, -0.528611029520217142048e-6) }; + // LCOV_EXCL_STOP result = tools::evaluate_polynomial(P, z) / tools::evaluate_polynomial(Q, z); result += z - log(z) - Y; @@ -83,6 +85,7 @@ T expint_1_rational(const T& z, const std::integral_constant&) { // Maximum Deviation Found (interpolated): 1.444e-17 // Max error found at double precision: 3.119e-17 + // LCOV_EXCL_START static const T P[11] = { BOOST_MATH_BIG_CONSTANT(T, 53, -0.121013190657725568138e-18), BOOST_MATH_BIG_CONSTANT(T, 53, -0.999999999999998811143), @@ -110,6 +113,7 @@ T expint_1_rational(const T& z, const std::integral_constant&) BOOST_MATH_BIG_CONSTANT(T, 53, 1229.20784182403048905), BOOST_MATH_BIG_CONSTANT(T, 53, -0.776491285282330997549) }; + // LCOV_EXCL_STOP T recip = 1 / z; result = 1 + tools::evaluate_polynomial(P, recip) / tools::evaluate_polynomial(Q, recip); @@ -132,7 +136,7 @@ T expint_1_rational(const T& z, const std::integral_constant&) // Maximum Deviation Found: 3.807e-20 // Expected Error Term: 3.807e-20 // Max error found at long double precision: 6.249e-20 - + // LCOV_EXCL_START static const T Y = 0.66373538970947265625F; static const T P[6] = { BOOST_MATH_BIG_CONSTANT(T, 64, 0.0865197248079397956816), @@ -151,6 +155,7 @@ T expint_1_rational(const T& z, const std::integral_constant&) BOOST_MATH_BIG_CONSTANT(T, 64, -0.202872781770207871975e-5), BOOST_MATH_BIG_CONSTANT(T, 64, 0.52779248094603709945e-7) }; + // LCOV_EXCL_STOP result = tools::evaluate_polynomial(P, z) / tools::evaluate_polynomial(Q, z); result += z - log(z) - Y; @@ -159,6 +164,7 @@ T expint_1_rational(const T& z, const std::integral_constant&) { // Maximum Deviation Found (interpolated): 2.220e-20 // Max error found at long double precision: 1.346e-19 + // LCOV_EXCL_START static const T P[14] = { BOOST_MATH_BIG_CONSTANT(T, 64, -0.534401189080684443046e-23), BOOST_MATH_BIG_CONSTANT(T, 64, -0.999999999999999999905), @@ -191,6 +197,7 @@ T expint_1_rational(const T& z, const std::integral_constant&) BOOST_MATH_BIG_CONSTANT(T, 64, 73930.2995984054930821), BOOST_MATH_BIG_CONSTANT(T, 64, 2063.86994219629165937) }; + // LCOV_EXCL_STOP T recip = 1 / z; result = 1 + tools::evaluate_polynomial(P, recip) / tools::evaluate_polynomial(Q, recip); @@ -213,7 +220,7 @@ T expint_1_rational(const T& z, const std::integral_constant&) // Maximum Deviation Found: 2.477e-35 // Expected Error Term: 2.477e-35 // Max error found at long double precision: 6.810e-35 - + // LCOV_EXCL_START static const T Y = 0.66373538970947265625F; static const T P[10] = { BOOST_MATH_BIG_CONSTANT(T, 113, 0.0865197248079397956434879099175975937), @@ -239,6 +246,7 @@ T expint_1_rational(const T& z, const std::integral_constant&) BOOST_MATH_BIG_CONSTANT(T, 113, 0.369373328141051577845488477377890236e-9), BOOST_MATH_BIG_CONSTANT(T, 113, -0.274149801370933606409282434677600112e-12) }; + // LCOV_EXCL_STOP result = tools::evaluate_polynomial(P, z) / tools::evaluate_polynomial(Q, z); result += z - log(z) - Y; @@ -247,7 +255,7 @@ T expint_1_rational(const T& z, const std::integral_constant&) { // Max error in interpolated form: 5.614e-35 // Max error found at long double precision: 7.979e-35 - + // LCOV_EXCL_START static const T Y = 0.70190334320068359375F; static const T P[16] = { @@ -286,6 +294,7 @@ T expint_1_rational(const T& z, const std::integral_constant&) BOOST_MATH_BIG_CONSTANT(T, 113, 169.845369689596739824177412096477219), BOOST_MATH_BIG_CONSTANT(T, 113, 2.17607292280092201170768401876895354) }; + // LCOV_EXCL_STOP T recip = 1 / z; result = Y + tools::evaluate_polynomial(P, recip) / tools::evaluate_polynomial(Q, recip); @@ -295,7 +304,7 @@ T expint_1_rational(const T& z, const std::integral_constant&) { // Max error in interpolated form: 4.413e-35 // Max error found at long double precision: 8.928e-35 - + // LCOV_EXCL_START static const T P[19] = { BOOST_MATH_BIG_CONSTANT(T, 113, -0.559148411832951463689610809550083986e-40), BOOST_MATH_BIG_CONSTANT(T, 113, -0.999999999999999999999999999999999997), @@ -339,6 +348,7 @@ T expint_1_rational(const T& z, const std::integral_constant&) BOOST_MATH_BIG_CONSTANT(T, 113, 70242279152.8241187845178443118302693), BOOST_MATH_BIG_CONSTANT(T, 113, -37633302.9409263839042721539363416685) }; + // LCOV_EXCL_STOP T recip = 1 / z; result = 1 + tools::evaluate_polynomial(P, recip) / tools::evaluate_polynomial(Q, recip); @@ -1486,94 +1496,6 @@ T expint_i_imp(T z, const Policy& pol, const std::integral_constant& t return result; } -template -struct expint_i_initializer -{ - struct init - { - init() - { - do_init(tag()); - } - static void do_init(const std::integral_constant&){} - static void do_init(const std::integral_constant&) - { - boost::math::expint(T(5), Policy()); - boost::math::expint(T(7), Policy()); - boost::math::expint(T(18), Policy()); - boost::math::expint(T(38), Policy()); - boost::math::expint(T(45), Policy()); - } - static void do_init(const std::integral_constant&) - { - boost::math::expint(T(5), Policy()); - boost::math::expint(T(7), Policy()); - boost::math::expint(T(18), Policy()); - boost::math::expint(T(38), Policy()); - boost::math::expint(T(45), Policy()); - } - static void do_init(const std::integral_constant&) - { - boost::math::expint(T(5), Policy()); - boost::math::expint(T(7), Policy()); - boost::math::expint(T(17), Policy()); - boost::math::expint(T(25), Policy()); - boost::math::expint(T(40), Policy()); - boost::math::expint(T(50), Policy()); - boost::math::expint(T(80), Policy()); - boost::math::expint(T(200), Policy()); - boost::math::expint(T(220), Policy()); - } - void force_instantiate()const{} - }; - static const init initializer; - static void force_instantiate() - { - initializer.force_instantiate(); - } -}; - -template -const typename expint_i_initializer::init expint_i_initializer::initializer; - -template -struct expint_1_initializer -{ - struct init - { - init() - { - do_init(tag()); - } - static void do_init(const std::integral_constant&){} - static void do_init(const std::integral_constant&) - { - boost::math::expint(1, T(0.5), Policy()); - boost::math::expint(1, T(2), Policy()); - } - static void do_init(const std::integral_constant&) - { - boost::math::expint(1, T(0.5), Policy()); - boost::math::expint(1, T(2), Policy()); - } - static void do_init(const std::integral_constant&) - { - boost::math::expint(1, T(0.5), Policy()); - boost::math::expint(1, T(2), Policy()); - boost::math::expint(1, T(6), Policy()); - } - void force_instantiate()const{} - }; - static const init initializer; - static void force_instantiate() - { - initializer.force_instantiate(); - } -}; - -template -const typename expint_1_initializer::init expint_1_initializer::initializer; - template inline typename tools::promote_args::type expint_forwarder(T z, const Policy& /*pol*/, std::true_type const&) @@ -1594,8 +1516,6 @@ inline typename tools::promote_args::type precision_type::value <= 113 ? 113 : 0 > tag_type; - expint_i_initializer::force_instantiate(); - return policies::checked_narrowing_cast(detail::expint_i_imp( static_cast(z), forwarding_policy(), @@ -1631,8 +1551,6 @@ inline typename tools::promote_args::type precision_type::value <= 113 ? 113 : 0 > tag_type; - detail::expint_1_initializer::force_instantiate(); - return policies::checked_narrowing_cast(detail::expint_imp( n, static_cast(z), diff --git a/include/boost/math/special_functions/expm1.hpp b/include/boost/math/special_functions/expm1.hpp index eec6356031..f53e42ab24 100644 --- a/include/boost/math/special_functions/expm1.hpp +++ b/include/boost/math/special_functions/expm1.hpp @@ -69,37 +69,6 @@ namespace detail expm1_series& operator=(const expm1_series&) = delete; }; -template -struct expm1_initializer -{ - struct init - { - init() - { - do_init(tag()); - } - template - static void do_init(const std::integral_constant&){} - static void do_init(const std::integral_constant&) - { - expm1(T(0.5)); - } - static void do_init(const std::integral_constant&) - { - expm1(T(0.5)); - } - void force_instantiate()const{} - }; - static const init initializer; - static void force_instantiate() - { - initializer.force_instantiate(); - } -}; - -template -const typename expm1_initializer::init expm1_initializer::initializer; - // // Algorithm expm1 is part of C99, but is not yet provided by many compilers. // @@ -142,6 +111,10 @@ T expm1_imp(T x, const std::integral_constant&, const P& pol) BOOST_MATH_STD_USING T a = fabs(x); + if ((boost::math::isnan)(a)) + { + return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol); + } if(a > T(0.5L)) { if(a >= tools::log_max_value()) @@ -169,6 +142,10 @@ T expm1_imp(T x, const std::integral_constant&, const P& pol) BOOST_MATH_STD_USING T a = fabs(x); + if ((boost::math::isnan)(a)) + { + return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol); + } if(a > T(0.5L)) { if(a >= tools::log_max_value()) @@ -212,6 +189,10 @@ T expm1_imp(T x, const std::integral_constant&, const P& pol) BOOST_MATH_STD_USING T a = fabs(x); + if ((boost::math::isnan)(a)) + { + return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", a, pol); + } if(a > T(0.5L)) { if(a >= tools::log_max_value()) @@ -278,8 +259,6 @@ inline typename tools::promote_args::type expm1(T x, const Policy& /* pol */) precision_type::value <= 113 ? 113 : 0 > tag_type; - detail::expm1_initializer::force_instantiate(); - return policies::checked_narrowing_cast(detail::expm1_imp( static_cast(x), tag_type(), forwarding_policy()), "boost::math::expm1<%1%>(%1%)"); @@ -294,14 +273,85 @@ inline typename tools::promote_args::type expm1(T x, const Policy& /* pol */) #if defined(BOOST_HAS_EXPM1) && !(defined(__osf__) && defined(__DECCXX_VER)) # ifdef BOOST_MATH_USE_C99 -inline float expm1(float x, const policies::policy<>&){ return ::expm1f(x); } +template +inline float expm1(float x, const Policy&) +{ + BOOST_MATH_IF_CONSTEXPR(Policy::domain_error_type::value != boost::math::policies::ignore_error && Policy::domain_error_type::value != boost::math::policies::errno_on_error) + { + if((boost::math::isnan)(x)) + return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy()); + } + BOOST_MATH_IF_CONSTEXPR(Policy::overflow_error_type::value != boost::math::policies::ignore_error && Policy::overflow_error_type::value != boost::math::policies::errno_on_error) + { + if (x >= tools::log_max_value()) + return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, Policy()); + } + return ::expm1f(x); +} # ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS -inline long double expm1(long double x, const policies::policy<>&){ return ::expm1l(x); } +template +inline long double expm1(long double x, const Policy&) +{ + BOOST_MATH_IF_CONSTEXPR(Policy::domain_error_type::value != boost::math::policies::ignore_error && Policy::domain_error_type::value != boost::math::policies::errno_on_error) + { + if ((boost::math::isnan)(x)) + return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy()); + } + BOOST_MATH_IF_CONSTEXPR(Policy::overflow_error_type::value != boost::math::policies::ignore_error && Policy::overflow_error_type::value != boost::math::policies::errno_on_error) + { + if (x >= tools::log_max_value()) + return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, Policy()); + } + return ::expm1l(x); +} # endif # else -inline float expm1(float x, const policies::policy<>&){ return static_cast(::expm1(x)); } +template +inline float expm1(float x, const Policy&) +{ + BOOST_MATH_IF_CONSTEXPR(Policy::domain_error_type::value != boost::math::policies::ignore_error && Policy::domain_error_type::value != boost::math::policies::errno_on_error) + { + if ((boost::math::isnan)(x)) + return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy()); + } + BOOST_MATH_IF_CONSTEXPR(Policy::overflow_error_type::value != boost::math::policies::ignore_error && Policy::overflow_error_type::value != boost::math::policies::errno_on_error) + { + if (x >= tools::log_max_value()) + return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, Policy()); + } + return static_cast(::expm1(x)); +} +template +inline typename std::enable_if::type expm1(long double x, const Policy&) +{ + BOOST_MATH_IF_CONSTEXPR(Policy::domain_error_type::value != boost::math::policies::ignore_error && Policy::domain_error_type::value != boost::math::policies::errno_on_error) + { + if ((boost::math::isnan)(x)) + return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy()); + } + BOOST_MATH_IF_CONSTEXPR(Policy::overflow_error_type::value != boost::math::policies::ignore_error && Policy::overflow_error_type::value != boost::math::policies::errno_on_error) + { + if (x >= tools::log_max_value()) + return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, Policy()); + } + return ::expm1(x); +} # endif -inline double expm1(double x, const policies::policy<>&){ return ::expm1(x); } +template +inline double expm1(double x, const Policy&) +{ + BOOST_MATH_IF_CONSTEXPR(Policy::domain_error_type::value != boost::math::policies::ignore_error && Policy::domain_error_type::value != boost::math::policies::errno_on_error) + { + if ((boost::math::isnan)(x)) + return policies::raise_domain_error("boost::math::expm1<%1%>(%1%)", "expm1 requires a finite argument, but got %1%", x, Policy()); + } + BOOST_MATH_IF_CONSTEXPR(Policy::overflow_error_type::value != boost::math::policies::ignore_error && Policy::overflow_error_type::value != boost::math::policies::errno_on_error) + { + if (x >= tools::log_max_value()) + return policies::raise_overflow_error("boost::math::expm1<%1%>(%1%)", nullptr, Policy()); + } + return ::expm1(x); +} #endif template diff --git a/test/log1p_expm1_test.hpp b/test/log1p_expm1_test.hpp index 56b6cd12f3..886e26d254 100644 --- a/test/log1p_expm1_test.hpp +++ b/test/log1p_expm1_test.hpp @@ -89,5 +89,9 @@ void test(T, const char* type_name) #endif #endif } + if (std::numeric_limits::has_quiet_NaN) + { + BOOST_MATH_CHECK_THROW(boost::math::expm1(std::numeric_limits::quiet_NaN()), std::domain_error); + } } diff --git a/test/test_expint.cpp b/test/test_expint.cpp index 3f44a80915..3db40d70e4 100644 --- a/test/test_expint.cpp +++ b/test/test_expint.cpp @@ -79,7 +79,7 @@ void expected_results() "float|double|long double", // test type(s) ".*Ei.*", // test data group ".*", 6, 3); // test function - if(std::numeric_limits::digits > 100) + BOOST_IF_CONSTEXPR (std::numeric_limits::digits > 100) { add_expected_result( ".*", // compiler diff --git a/test/test_expint.hpp b/test/test_expint.hpp index 491db2fcdc..d398fa1620 100644 --- a/test/test_expint.hpp +++ b/test/test_expint.hpp @@ -190,5 +190,27 @@ void test_spots(T, const char* t) BOOST_CHECK_CLOSE(::boost::math::expint(static_cast(-0.5)), static_cast(-0.559773594776160811746795939315085235226846890316353515248293L), tolerance); BOOST_CHECK_CLOSE(::boost::math::expint(static_cast(-1)), static_cast(-0.219383934395520273677163775460121649031047293406908207577979L), tolerance); BOOST_CHECK_CLOSE(::boost::math::expint(static_cast(-50.5)), static_cast(-2.27237132932219350440719707268817831250090574830769670186618e-24L), tolerance); + // + // Extra coverage cases: + // + BOOST_CHECK_EQUAL(boost::math::expint(-boost::math::tools::max_value()), T(0)); + BOOST_CHECK_THROW(boost::math::expint(1, T(-1)), std::domain_error); + BOOST_CHECK_THROW(boost::math::expint(2, T(-1)), std::domain_error); + BOOST_CHECK_EQUAL(boost::math::expint(2, T(0)), T(1)); + BOOST_CHECK_EQUAL(boost::math::expint(3, T(0)), T(0.5)); + BOOST_IF_CONSTEXPR(std::numeric_limits::has_infinity) + { + BOOST_CHECK_EQUAL(boost::math::expint(1, T(0)), std::numeric_limits::infinity()); + BOOST_CHECK_EQUAL(boost::math::expint(T(0)), -std::numeric_limits::infinity()); + BOOST_CHECK_EQUAL(boost::math::expint(boost::math::tools::log_max_value() * 2), std::numeric_limits::infinity()); + BOOST_CHECK_EQUAL(boost::math::expint(boost::math::tools::log_max_value() + T(38)), std::numeric_limits::infinity()); + } + else + { + BOOST_CHECK_GE(boost::math::expint(1, T(0)), boost::math::tools::max_value()); + BOOST_CHECK_LE(boost::math::expint(T(0)), -boost::math::tools::max_value()); + BOOST_CHECK_GE(boost::math::expint(boost::math::tools::log_max_value() * 2), boost::math::tools::max_value()); + BOOST_CHECK_GE(boost::math::expint(boost::math::tools::log_max_value() + T(38)), boost::math::tools::max_value()); + } } From fbf553d49850bd422be64e40211740693c2c1ca2 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Wed, 6 Mar 2024 19:17:44 +0000 Subject: [PATCH 02/27] Add missing #include. --- include/boost/math/special_functions/expm1.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/include/boost/math/special_functions/expm1.hpp b/include/boost/math/special_functions/expm1.hpp index f53e42ab24..fccf97bc34 100644 --- a/include/boost/math/special_functions/expm1.hpp +++ b/include/boost/math/special_functions/expm1.hpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128) From c446d7969b91aa52b3ba931df34270a5f0c88cf9 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 8 Mar 2024 17:12:11 +0000 Subject: [PATCH 03/27] expm1 coverage. --- .../boost/math/special_functions/expm1.hpp | 4 ++ test/Jamfile.v2 | 1 + test/test_expm1_extra.cpp | 67 +++++++++++++++++++ 3 files changed, 72 insertions(+) create mode 100644 test/test_expm1_extra.cpp diff --git a/include/boost/math/special_functions/expm1.hpp b/include/boost/math/special_functions/expm1.hpp index fccf97bc34..94e0921286 100644 --- a/include/boost/math/special_functions/expm1.hpp +++ b/include/boost/math/special_functions/expm1.hpp @@ -160,6 +160,7 @@ T expm1_imp(T x, const std::integral_constant&, const P& pol) if(a < tools::epsilon()) return x; + // LCOV_EXCL_START static const float Y = 0.10281276702880859375e1f; static const T n[] = { BOOST_MATH_BIG_CONSTANT(T, 64, -0.281276702880859375e-1), @@ -179,6 +180,7 @@ T expm1_imp(T x, const std::integral_constant&, const P& pol) BOOST_MATH_BIG_CONSTANT(T, 64, -0.387922804997682392562e-4), BOOST_MATH_BIG_CONSTANT(T, 64, 0.807473180049193557294e-6) }; + // LCOV_EXCL_STOP T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x); return result; @@ -207,6 +209,7 @@ T expm1_imp(T x, const std::integral_constant&, const P& pol) if(a < tools::epsilon()) return x; + // LCOV_EXCL_START static const float Y = 0.10281276702880859375e1f; static const T n[] = { BOOST_MATH_BIG_CONSTANT(T, 113, -0.28127670288085937499999999999999999854e-1), @@ -233,6 +236,7 @@ T expm1_imp(T x, const std::integral_constant&, const P& pol) BOOST_MATH_BIG_CONSTANT(T, 113, -0.29477341859111589208776402638429026517e-10), BOOST_MATH_BIG_CONSTANT(T, 113, 0.13222065991022301420255904060628100924e-12) }; + // LCOV_EXCL_STOP T result = x * Y + x * tools::evaluate_polynomial(n, x) / tools::evaluate_polynomial(d, x); return result; diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index cce2a8fcf8..42d02caf13 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -164,6 +164,7 @@ test-suite special_fun : [ run ccmath_fma_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : BOOST_MATH_TEST_FLOAT128 "-Bstatic -lquadmath -Bdynamic" ] ] [ run ccmath_signbit_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx17_if_constexpr ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : BOOST_MATH_TEST_FLOAT128 "-Bstatic -lquadmath -Bdynamic" ] ] [ run log1p_expm1_test.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] + [ run test_expm1_extra.cpp ../../test/build//boost_unit_test_framework : : : [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : BOOST_MATH_TEST_FLOAT128 "-Bstatic -lquadmath -Bdynamic" ] ] [ run powm1_sqrtp1m1_test.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] [ run git_issue_705.cpp ../../test/build//boost_unit_test_framework ] [ run git_issue_810.cpp ../../test/build//boost_unit_test_framework ] diff --git a/test/test_expm1_extra.cpp b/test/test_expm1_extra.cpp new file mode 100644 index 0000000000..02e1c59ea3 --- /dev/null +++ b/test/test_expm1_extra.cpp @@ -0,0 +1,67 @@ +// Copyright John Maddock 2024. +// 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_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#define BOOST_TEST_MAIN +#include +#include +#include +#include + +#define SC_(x) static_cast(BOOST_STRINGIZE(x)) + +#include "log1p_expm1_test.hpp" + +// +// DESCRIPTION: +// ~~~~~~~~~~~~ +// +// This file tests the functions log1p and expm1 for +// simple multiprecision types which are able to utilize +// our own rational approximations. This tests those +// approximations even when float/double/long double +// are not using them. +// + +void expected_results() +{ + // + // Define the max and mean errors expected for + // various compilers and platforms. + // + + // + // Catch all cases come last: + // + add_expected_result( + ".*", // compiler + ".*", // stdlib + ".*", // platform + ".*", // test type(s) + ".*", // test data group + ".*", // test function + 4, // Max Peek error + 3); // Max mean error + + // + // Finish off by printing out the compiler/stdlib/platform names, + // we do this to make it easier to mark up expected error rates. + // + std::cout << "Tests run with " << BOOST_COMPILER << ", " + << BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl; +} + + +BOOST_AUTO_TEST_CASE( test_main ) +{ + expected_results(); + + test(boost::multiprecision::cpp_bin_float_double(0), "UDT double"); + test(boost::multiprecision::cpp_bin_float_double_extended(0), "UDT 80-bit long double"); + test(boost::multiprecision::cpp_bin_float_quad(0), "UDT quad-float"); + test(boost::multiprecision::number >(0), "UDT 36 digits"); +} + From 47414e94202992fef0419eeb19061881c8689bd2 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 8 Mar 2024 17:28:28 +0000 Subject: [PATCH 04/27] Tidy up expint coverage. --- .../boost/math/special_functions/expint.hpp | 41 +++++++++++++------ test/test_expint.hpp | 8 ++++ 2 files changed, 36 insertions(+), 13 deletions(-) diff --git a/include/boost/math/special_functions/expint.hpp b/include/boost/math/special_functions/expint.hpp index 9b2b420552..cee9c562e3 100644 --- a/include/boost/math/special_functions/expint.hpp +++ b/include/boost/math/special_functions/expint.hpp @@ -530,7 +530,7 @@ T expint_i_imp(T z, const Policy& pol, const Tag& tag) if(z < 0) return -expint_imp(1, T(-z), pol, tag); if(z == 0) - return -policies::raise_overflow_error(function, nullptr, pol); + return -policies::raise_overflow_error(function, nullptr, pol); // LCOV_EXCL_LINE confirmed covered by real_concept tests return expint_i_as_series(z, pol); } @@ -766,6 +766,7 @@ T expint_i_imp(T z, const Policy& pol, const std::integral_constant& ta // Expected Error Term: 3.883e-21 // Max Error found at long double precision = Poly: 3.344801e-19 Cheb: 4.989937e-19 + // LCOV_EXCL_START static const T P[11] = { BOOST_MATH_BIG_CONSTANT(T, 64, 2.98677224343598593764), BOOST_MATH_BIG_CONSTANT(T, 64, 0.25891613550886736592), @@ -796,6 +797,8 @@ T expint_i_imp(T z, const Policy& pol, const std::integral_constant& ta static const T r1 = c1 / c2; static const T r2 = BOOST_MATH_BIG_CONSTANT(T, 64, 0.131401834143860282009280387409357165515556574352422001206362e-16); static const T r = static_cast(BOOST_MATH_BIG_CONSTANT(T, 64, 0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392)); + // LCOV_EXCL_STOP + T t = (z / 3) - 1; result = tools::evaluate_polynomial(P, t) / tools::evaluate_polynomial(Q, t); @@ -815,7 +818,7 @@ T expint_i_imp(T z, const Policy& pol, const std::integral_constant& ta // Maximum Deviation Found: 2.622e-21 // Expected Error Term: -2.622e-21 // Max Error found at long double precision = Poly: 1.208328e-20 Cheb: 1.073723e-20 - + // LCOV_EXCL_START static const T Y = 1.158985137939453125F; static const T P[9] = { BOOST_MATH_BIG_CONSTANT(T, 64, 0.00139324086199409049399), @@ -840,6 +843,7 @@ T expint_i_imp(T z, const Policy& pol, const std::integral_constant& ta BOOST_MATH_BIG_CONSTANT(T, 64, 0.204339282037446434827e-5), BOOST_MATH_BIG_CONSTANT(T, 64, 0.146951181174930425744e-7) }; + // LCOV_EXCL_STOP T t = z / 2 - 4; result = Y + tools::evaluate_polynomial(P, t) / tools::evaluate_polynomial(Q, t); @@ -852,7 +856,7 @@ T expint_i_imp(T z, const Policy& pol, const std::integral_constant& ta // Expected Error Term: 3.220e-20 // Max Error found at long double precision = Poly: 7.696841e-20 Cheb: 6.205163e-20 - + // LCOV_EXCL_START static const T Y = 1.0869731903076171875F; static const T P[10] = { BOOST_MATH_BIG_CONSTANT(T, 64, -0.00893891094356946995368), @@ -878,6 +882,7 @@ T expint_i_imp(T z, const Policy& pol, const std::integral_constant& ta BOOST_MATH_BIG_CONSTANT(T, 64, 0.000577048986213535829925), BOOST_MATH_BIG_CONSTANT(T, 64, 0.290976943033493216793e-4) }; + // LCOV_EXCL_STOP T t = z / 5 - 3; result = Y + tools::evaluate_polynomial(P, t) / tools::evaluate_polynomial(Q, t); @@ -889,7 +894,7 @@ T expint_i_imp(T z, const Policy& pol, const std::integral_constant& ta // Maximum Deviation Found: 2.940e-21 // Expected Error Term: -2.938e-21 // Max Error found at long double precision = Poly: 3.419893e-19 Cheb: 3.359874e-19 - + // LCOV_EXCL_START static const T Y = 1.03937530517578125F; static const T P[12] = { BOOST_MATH_BIG_CONSTANT(T, 64, -0.00356165148914447278177), @@ -916,6 +921,7 @@ T expint_i_imp(T z, const Policy& pol, const std::integral_constant& ta BOOST_MATH_BIG_CONSTANT(T, 64, 0.0280128013584653182994), BOOST_MATH_BIG_CONSTANT(T, 64, 0.00182034930799902922549) }; + // LCOV_EXCL_STOP T t = z / 10 - 3; result = Y + tools::evaluate_polynomial(P, t) / tools::evaluate_polynomial(Q, t); @@ -929,7 +935,7 @@ T expint_i_imp(T z, const Policy& pol, const std::integral_constant& ta { // Maximum Deviation Found: 3.536e-20 // Max Error found at long double precision = Poly: 1.310671e-19 Cheb: 8.630943e-11 - + // LCOV_EXCL_START static const T exp40 = static_cast(BOOST_MATH_BIG_CONSTANT(T, 64, 2.35385266837019985407899910749034804508871617254555467236651e17)); static const T Y= 1.013065338134765625F; static const T P[9] = { @@ -954,6 +960,7 @@ T expint_i_imp(T z, const Policy& pol, const std::integral_constant& ta BOOST_MATH_BIG_CONSTANT(T, 64, 384647824.678554961174), BOOST_MATH_BIG_CONSTANT(T, 64, -166288297.874583961493) }; + // LCOV_EXCL_STOP T t = 1 / z; result = Y + tools::evaluate_polynomial(P, t) / tools::evaluate_polynomial(Q, t); @@ -993,7 +1000,7 @@ void expint_i_imp_113a(T& result, const T& z, const Policy& pol) // Expected Error Term: -1.230e-36 // Max Error found at long double precision = Poly: 4.355299e-34 Cheb: 7.512581e-34 - + // LCOV_EXCL_START static const T P[15] = { BOOST_MATH_BIG_CONSTANT(T, 113, 2.98677224343598593765287235997328555), BOOST_MATH_BIG_CONSTANT(T, 113, -0.333256034674702967028780537349334037), @@ -1038,6 +1045,7 @@ void expint_i_imp_113a(T& result, const T& z, const Policy& pol) static const T r2 = c3 / c4 / c5; static const T r3 = static_cast(BOOST_MATH_BIG_CONSTANT(T, 113, 0.283806480836357377069325311780969887585024578164571984232357e-31)); static const T r = static_cast(BOOST_MATH_BIG_CONSTANT(T, 113, 0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392)); + // LCOV_EXCL_STOP T t = (z / 3) - 1; result = tools::evaluate_polynomial(P, t) / tools::evaluate_polynomial(Q, t); @@ -1060,7 +1068,7 @@ void expint_i_113b(T& result, const T& z) // Maximum Deviation Found: 7.779e-36 // Expected Error Term: -7.779e-36 // Max Error found at long double precision = Poly: 2.576723e-35 Cheb: 1.236001e-34 - + // LCOV_EXCL_START static const T Y = 1.158985137939453125F; static const T P[15] = { BOOST_MATH_BIG_CONSTANT(T, 113, 0.00139324086199409049282472239613554817), @@ -1096,6 +1104,7 @@ void expint_i_113b(T& result, const T& z) BOOST_MATH_BIG_CONSTANT(T, 113, 0.139007266881450521776529705677086902e-9), BOOST_MATH_BIG_CONSTANT(T, 113, 0.234715286125516430792452741830364672e-11) }; + // LCOV_EXCL_STOP T t = z / 2 - 4; result = Y + tools::evaluate_polynomial(P, t) / tools::evaluate_polynomial(Q, t); @@ -1111,7 +1120,7 @@ void expint_i_113c(T& result, const T& z) // Expected Error Term: 1.080e-34 // Max Error found at long double precision = Poly: 1.958294e-34 Cheb: 2.472261e-34 - + // LCOV_EXCL_START static const T Y = 1.091579437255859375F; static const T P[17] = { BOOST_MATH_BIG_CONSTANT(T, 113, -0.00685089599550151282724924894258520532), @@ -1149,6 +1158,7 @@ void expint_i_113c(T& result, const T& z) BOOST_MATH_BIG_CONSTANT(T, 113, 0.233593219218823384508105943657387644e-7), BOOST_MATH_BIG_CONSTANT(T, 113, 0.554900353169148897444104962034267682e-9) }; + // LCOV_EXCL_STOP T t = z / 4 - 3.5; result = Y + tools::evaluate_polynomial(P, t) / tools::evaluate_polynomial(Q, t); @@ -1163,7 +1173,7 @@ void expint_i_113d(T& result, const T& z) // Maximum Deviation Found: 3.163e-35 // Expected Error Term: 3.163e-35 // Max Error found at long double precision = Poly: 4.158110e-35 Cheb: 5.385532e-35 - + // LCOV_EXCL_START static const T Y = 1.051731109619140625F; static const T P[14] = { BOOST_MATH_BIG_CONSTANT(T, 113, -0.00144552494420652573815404828020593565), @@ -1197,6 +1207,7 @@ void expint_i_113d(T& result, const T& z) BOOST_MATH_BIG_CONSTANT(T, 113, 0.127552010539733113371132321521204458e-7), BOOST_MATH_BIG_CONSTANT(T, 113, 0.25737310826983451144405899970774587e-9) }; + // LCOV_EXCL_STOP T t = z / 4 - 5.5; result = Y + tools::evaluate_polynomial(P, t) / tools::evaluate_polynomial(Q, t); @@ -1214,7 +1225,7 @@ void expint_i_113e(T& result, const T& z) // Maximum Deviation Found: 7.972e-36 // Expected Error Term: 7.962e-36 // Max Error found at long double precision = Poly: 1.711721e-34 Cheb: 3.100018e-34 - + // LCOV_EXCL_START static const T Y = 1.032726287841796875F; static const T P[15] = { BOOST_MATH_BIG_CONSTANT(T, 113, -0.00141056919297307534690895009969373233), @@ -1251,6 +1262,7 @@ void expint_i_113e(T& result, const T& z) BOOST_MATH_BIG_CONSTANT(T, 113, 0.322153582559488797803027773591727565e-7), BOOST_MATH_BIG_CONSTANT(T, 113, -0.161635525318683508633792845159942312e-16) }; + // LCOV_EXCL_STOP T t = z / 8 - 4.25; result = Y + tools::evaluate_polynomial(P, t) / tools::evaluate_polynomial(Q, t); @@ -1268,7 +1280,7 @@ void expint_i_113f(T& result, const T& z) // Maximum Deviation Found: 4.469e-36 // Expected Error Term: 4.468e-36 // Max Error found at long double precision = Poly: 1.288958e-35 Cheb: 2.304586e-35 - + // LCOV_EXCL_START static const T Y = 1.0216197967529296875F; static const T P[12] = { BOOST_MATH_BIG_CONSTANT(T, 113, -0.000322999116096627043476023926572650045), @@ -1298,6 +1310,7 @@ void expint_i_113f(T& result, const T& z) BOOST_MATH_BIG_CONSTANT(T, 113, 0.233740058688179614344680531486267142e-5), BOOST_MATH_BIG_CONSTANT(T, 113, 0.498800627828842754845418576305379469e-7) }; + // LCOV_EXCL_STOP T t = z / 7 - 7; result = Y + tools::evaluate_polynomial(P, t) / tools::evaluate_polynomial(Q, t); @@ -1315,7 +1328,7 @@ void expint_i_113g(T& result, const T& z) // Maximum Deviation Found: 5.588e-35 // Expected Error Term: -5.566e-35 // Max Error found at long double precision = Poly: 9.976345e-35 Cheb: 8.358865e-35 - + // LCOV_EXCL_START static const T Y = 1.015148162841796875F; static const T P[11] = { BOOST_MATH_BIG_CONSTANT(T, 113, -0.000435714784725086961464589957142615216), @@ -1343,6 +1356,7 @@ void expint_i_113g(T& result, const T& z) BOOST_MATH_BIG_CONSTANT(T, 113, 0.000167479843750859222348869769094711093), BOOST_MATH_BIG_CONSTANT(T, 113, 0.475673638665358075556452220192497036e-5) }; + // LCOV_EXCL_STOP T t = z / 14 - 5; result = Y + tools::evaluate_polynomial(P, t) / tools::evaluate_polynomial(Q, t); @@ -1360,7 +1374,7 @@ void expint_i_113h(T& result, const T& z) // Maximum Deviation Found: 4.448e-36 // Expected Error Term: 4.445e-36 // Max Error found at long double precision = Poly: 2.058532e-35 Cheb: 2.165465e-27 - + // LCOV_EXCL_START static const T Y= 1.00849151611328125F; static const T P[9] = { BOOST_MATH_BIG_CONSTANT(T, 113, -0.0084915161132812500000001440233607358), @@ -1385,6 +1399,7 @@ void expint_i_113h(T& result, const T& z) BOOST_MATH_BIG_CONSTANT(T, 113, 8354144.67882768405803322344185185517), BOOST_MATH_BIG_CONSTANT(T, 113, 355076.853106511136734454134915432571) }; + // LCOV_EXCL_STOP T t = 1 / z; result = Y + tools::evaluate_polynomial(P, t) / tools::evaluate_polynomial(Q, t); diff --git a/test/test_expint.hpp b/test/test_expint.hpp index d398fa1620..fb2f138c26 100644 --- a/test/test_expint.hpp +++ b/test/test_expint.hpp @@ -203,7 +203,15 @@ void test_spots(T, const char* t) BOOST_CHECK_EQUAL(boost::math::expint(1, T(0)), std::numeric_limits::infinity()); BOOST_CHECK_EQUAL(boost::math::expint(T(0)), -std::numeric_limits::infinity()); BOOST_CHECK_EQUAL(boost::math::expint(boost::math::tools::log_max_value() * 2), std::numeric_limits::infinity()); + if (boost::math::tools::log_max_value() < boost::math::tools::log_max_value()) + { + BOOST_CHECK_EQUAL(boost::math::expint(boost::math::tools::log_max_value() * 2), std::numeric_limits::infinity()); + } BOOST_CHECK_EQUAL(boost::math::expint(boost::math::tools::log_max_value() + T(38)), std::numeric_limits::infinity()); + if (boost::math::tools::log_max_value() < boost::math::tools::log_max_value()) + { + BOOST_CHECK_EQUAL(boost::math::expint(boost::math::tools::log_max_value() + T(38)), std::numeric_limits::infinity()); + } } else { From e945458aa46c8387156a5a741dcd03655971160b Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 8 Mar 2024 17:51:21 +0000 Subject: [PATCH 05/27] Factorial coverage. --- test/test_factorials.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/test/test_factorials.cpp b/test/test_factorials.cpp index 6b413242c6..9741448dfe 100644 --- a/test/test_factorials.cpp +++ b/test/test_factorials.cpp @@ -309,6 +309,13 @@ void test_spots(T) BOOST_CHECK_CLOSE(::boost::math::falling_factorial(ldexp(static_cast(1), -300), 200), static_cast(-1.93579759151806711025267355739174942986011285920860098569075e282L), 10 * tolerance); } + // for coverage: + BOOST_MATH_IF_CONSTEXPR(std::numeric_limits::has_infinity && (std::numeric_limits::max_exponent <= std::numeric_limits::max_exponent)) + { + BOOST_CHECK_EQUAL(::boost::math::falling_factorial(boost::math::tools::epsilon(), 200), -std::numeric_limits::infinity()); + } + + tolerance = boost::math::tools::epsilon() * 100 * 20; // 20 eps as a percent. unsigned i = boost::math::max_factorial::value; From bb9b10e4dd62e69099fd934e44575d6fa0ca4c6e Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 8 Mar 2024 18:42:06 +0000 Subject: [PATCH 06/27] Correct expint test case. --- test/test_expint.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_expint.hpp b/test/test_expint.hpp index fb2f138c26..a5a98d6efb 100644 --- a/test/test_expint.hpp +++ b/test/test_expint.hpp @@ -205,12 +205,12 @@ void test_spots(T, const char* t) BOOST_CHECK_EQUAL(boost::math::expint(boost::math::tools::log_max_value() * 2), std::numeric_limits::infinity()); if (boost::math::tools::log_max_value() < boost::math::tools::log_max_value()) { - BOOST_CHECK_EQUAL(boost::math::expint(boost::math::tools::log_max_value() * 2), std::numeric_limits::infinity()); + BOOST_CHECK_EQUAL(boost::math::expint(static_cast(boost::math::tools::log_max_value() * 2)), std::numeric_limits::infinity()); } BOOST_CHECK_EQUAL(boost::math::expint(boost::math::tools::log_max_value() + T(38)), std::numeric_limits::infinity()); if (boost::math::tools::log_max_value() < boost::math::tools::log_max_value()) { - BOOST_CHECK_EQUAL(boost::math::expint(boost::math::tools::log_max_value() + T(38)), std::numeric_limits::infinity()); + BOOST_CHECK_EQUAL(boost::math::expint(static_cast(boost::math::tools::log_max_value() + T(38))), std::numeric_limits::infinity()); } } else From 718004ed27489622819255e8a3dde68a00f27ff8 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sun, 10 Mar 2024 19:29:11 +0000 Subject: [PATCH 07/27] fpclassify coverage. --- include/boost/math/special_functions/fpclassify.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/math/special_functions/fpclassify.hpp b/include/boost/math/special_functions/fpclassify.hpp index 2c504d7ac8..6f29f80510 100644 --- a/include/boost/math/special_functions/fpclassify.hpp +++ b/include/boost/math/special_functions/fpclassify.hpp @@ -163,7 +163,7 @@ inline int fpclassify_imp BOOST_NO_MACRO_EXPAND(T t, const generic_tag&) // whenever possible check for Nan's first: #if defined(BOOST_HAS_FPCLASSIFY) && !defined(BOOST_MATH_DISABLE_STD_FPCLASSIFY) if(::boost::math_detail::is_nan_helper(t, typename std::is_floating_point::type())) - return FP_NAN; + return FP_NAN; // LCOV_EXCL_LINE only called in UDT contexts (excluded from coverage checks). #elif defined(isnan) if(boost::math_detail::is_nan_helper(t, typename std::is_floating_point::type())) return FP_NAN; From 0da61dece03070a5e3722fe9d1b7b17edb1c6425 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sun, 10 Mar 2024 19:29:36 +0000 Subject: [PATCH 08/27] Gamma coverage. --- .../boost/math/special_functions/gamma.hpp | 169 +++++------------- test/test_gamma.hpp | 24 +++ 2 files changed, 67 insertions(+), 126 deletions(-) diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index a58ea3e693..535cd6b14e 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -87,7 +87,7 @@ T sinpx(T z) z = -z; } T fl = floor(z); - T dist; + T dist; // LCOV_EXCL_LINE if(is_odd(fl)) { fl += 1; @@ -128,17 +128,27 @@ T gamma_imp(T z, const Policy& pol, const Lanczos& l) { if(floor(z) == z) return policies::raise_pole_error(function, "Evaluation of tgamma at a negative integer %1%.", z, pol); - if(z <= -20) + if (z <= -20) { - result = gamma_imp(T(-z), pol, l) * sinpx(z); + try + { + result = gamma_imp(T(-z), pol, l) * sinpx(z); + } + catch (const std::overflow_error&) + { + return policies::raise_underflow_error(function, "Result of tgamma is too small to represent.", pol); + } BOOST_MATH_INSTRUMENT_VARIABLE(result); - if((fabs(result) < 1) && (tools::max_value() * fabs(result) < boost::math::constants::pi())) - return -boost::math::sign(result) * policies::raise_overflow_error(function, "Result of tgamma is too large to represent.", pol); + // Result can never be small: tgamma[-z] is always larger than sinpx[z] is small: + BOOST_MATH_ASSERT((fabs(result) > 1) || (tools::max_value() * fabs(result) > boost::math::constants::pi())); result = -boost::math::constants::pi() / result; - if(result == 0) + if (result == 0) return policies::raise_underflow_error(function, "Result of tgamma is too small to represent.", pol); + /* + * Result can never be subnormal as we have a value > 1 in the numerator: if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL) return policies::raise_denorm_error(function, "Result of tgamma is denormalized.", result, pol); + */ BOOST_MATH_INSTRUMENT_VARIABLE(result); return result; } @@ -435,14 +445,16 @@ T scaled_tgamma_no_lanczos(const T& z, const Policy& pol, bool islog = false) break; } if (n > number_of_bernoullis_b2n) - return policies::raise_evaluation_error("scaled_tgamma_no_lanczos<%1%>()", "Exceeded maximum series iterations without reaching convergence, best approximation was %1%", T(exp(sum) * half_ln_two_pi_over_z), pol); + // Safety net, we hope to never get here: + return policies::raise_evaluation_error("scaled_tgamma_no_lanczos<%1%>()", "Exceeded maximum series iterations without reaching convergence, best approximation was %1%", T(exp(sum) * half_ln_two_pi_over_z), pol); // LCOV_EXCL_LINE sum += term; // Sanity check for divergence: T fterm = fabs(term); if(fterm > last_term) - return policies::raise_evaluation_error("scaled_tgamma_no_lanczos<%1%>()", "Series became divergent without reaching convergence, best approximation was %1%", T(exp(sum) * half_ln_two_pi_over_z), pol); + // Safety net, we hope to never get here: + return policies::raise_evaluation_error("scaled_tgamma_no_lanczos<%1%>()", "Series became divergent without reaching convergence, best approximation was %1%", T(exp(sum) * half_ln_two_pi_over_z), pol); // LCOV_EXCL_LINE last_term = fterm; } @@ -513,16 +525,16 @@ T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&) if (!n_recur) { if (zz > tools::log_max_value()) - return policies::raise_overflow_error(function, nullptr, pol); + return b_neg ? policies::raise_underflow_error(function, nullptr, pol) : policies::raise_overflow_error(function, nullptr, pol); if (log(zz) * zz / 2 > tools::log_max_value()) - return policies::raise_overflow_error(function, nullptr, pol); + return b_neg ? policies::raise_underflow_error(function, nullptr, pol) : policies::raise_overflow_error(function, nullptr, pol); } T gamma_value = scaled_tgamma_no_lanczos(zz, pol); T power_term = pow(zz, zz / 2); T exp_term = exp(-zz); gamma_value *= (power_term * exp_term); - if(!n_recur && (tools::max_value() / power_term < gamma_value)) - return policies::raise_overflow_error(function, nullptr, pol); + if (!n_recur && (tools::max_value() / power_term < gamma_value)) + return b_neg ? policies::raise_underflow_error(function, nullptr, pol) : policies::raise_overflow_error(function, nullptr, pol); gamma_value *= power_term; // Rescale the result using downward recursion if necessary. @@ -553,24 +565,23 @@ T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&) if(floor_of_z_is_equal_to_z) return policies::raise_pole_error(function, "Evaluation of tgamma at a negative integer %1%.", z, pol); - gamma_value *= sinpx(z); + T s = sinpx(z); + if ((gamma_value > 1) && (tools::max_value() / gamma_value < fabs(s))) + return policies::raise_underflow_error(function, nullptr, pol); + gamma_value *= s; BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value); - - const bool result_is_too_large_to_represent = ( (abs(gamma_value) < 1) - && ((tools::max_value() * abs(gamma_value)) < boost::math::constants::pi())); - - if(result_is_too_large_to_represent) - return policies::raise_overflow_error(function, "Result of tgamma is too large to represent.", pol); + // + // Result can never overflow, since sinpx(z) can never be smaller than machine epsilon and gamma_value > 1. + // + BOOST_MATH_ASSERT( (abs(gamma_value) > 1) || ((tools::max_value() * abs(gamma_value)) > boost::math::constants::pi())); gamma_value = -boost::math::constants::pi() / gamma_value; BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value); - - if(gamma_value == 0) - return policies::raise_underflow_error(function, "Result of tgamma is too small to represent.", pol); - - if((boost::math::fpclassify)(gamma_value) == static_cast(FP_SUBNORMAL)) - return policies::raise_denorm_error(function, "Result of tgamma is denormalized.", gamma_value, pol); + // + // We can never underflow since the numerator > 1 above: + // + BOOST_MATH_ASSERT(gamma_value != 0); } return gamma_value; @@ -586,18 +597,21 @@ inline T log_gamma_near_1(const T& z, Policy const& pol) // BOOST_MATH_STD_USING // ADL of std names - BOOST_MATH_ASSERT(fabs(z) < 1); + // For some reason, several lines aren't triggered for coverage even though + // adjacent lines are... weird! + + BOOST_MATH_ASSERT(fabs(z) < 1); // LCOV_EXCL_LINE T result = -constants::euler() * z; T power_term = z * z / 2; - int n = 2; - T term = 0; + int n = 2; // LCOV_EXCL_LINE + T term = 0; // LCOV_EXCL_LINE do { term = power_term * boost::math::polygamma(n - 1, T(1), pol); - result += term; + result += term; // LCOV_EXCL_LINE ++n; power_term *= z / n; } while (fabs(result) * tools::epsilon() < fabs(term)); @@ -1832,93 +1846,6 @@ inline typename tools::promote_args::type return policies::checked_narrowing_cast(detail::gamma_imp(static_cast(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma<%1%>(%1%)"); } -template -struct igamma_initializer -{ - struct init - { - init() - { - typedef typename policies::precision::type precision_type; - - typedef std::integral_constant tag_type; - - do_init(tag_type()); - } - template - static void do_init(const std::integral_constant&) - { - // If std::numeric_limits::digits is zero, we must not call - // our initialization code here as the precision presumably - // varies at runtime, and will not have been set yet. Plus the - // code requiring initialization isn't called when digits == 0. - if(std::numeric_limits::digits) - { - boost::math::gamma_p(static_cast(400), static_cast(400), Policy()); - } - } - static void do_init(const std::integral_constant&){} - void force_instantiate()const{} - }; - static const init initializer; - static void force_instantiate() - { - initializer.force_instantiate(); - } -}; - -template -const typename igamma_initializer::init igamma_initializer::initializer; - -template -struct lgamma_initializer -{ - struct init - { - init() - { - typedef typename policies::precision::type precision_type; - typedef std::integral_constant tag_type; - - do_init(tag_type()); - } - static void do_init(const std::integral_constant&) - { - boost::math::lgamma(static_cast(2.5), Policy()); - boost::math::lgamma(static_cast(1.25), Policy()); - boost::math::lgamma(static_cast(1.75), Policy()); - } - static void do_init(const std::integral_constant&) - { - boost::math::lgamma(static_cast(2.5), Policy()); - boost::math::lgamma(static_cast(1.25), Policy()); - boost::math::lgamma(static_cast(1.5), Policy()); - boost::math::lgamma(static_cast(1.75), Policy()); - } - static void do_init(const std::integral_constant&) - { - } - void force_instantiate()const{} - }; - static const init initializer; - static void force_instantiate() - { - initializer.force_instantiate(); - } -}; - -template -const typename lgamma_initializer::init lgamma_initializer::initializer; - template inline tools::promote_args_t tgamma(T1 a, T2 z, const Policy&, const std::false_type) @@ -1934,8 +1861,6 @@ inline tools::promote_args_t policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; - igamma_initializer::force_instantiate(); - return policies::checked_narrowing_cast( detail::gamma_incomplete_imp(static_cast(a), static_cast(z), false, true, @@ -1974,8 +1899,6 @@ inline typename tools::promote_args::type policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; - detail::lgamma_initializer::force_instantiate(); - return policies::checked_narrowing_cast(detail::lgamma_imp(static_cast(z), forwarding_policy(), evaluation_type(), sign), "boost::math::lgamma<%1%>(%1%)"); } @@ -2065,8 +1988,6 @@ inline tools::promote_args_t policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; - detail::igamma_initializer::force_instantiate(); - return policies::checked_narrowing_cast( detail::gamma_incomplete_imp(static_cast(a), static_cast(z), false, false, @@ -2096,8 +2017,6 @@ inline tools::promote_args_t policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; - detail::igamma_initializer::force_instantiate(); - return policies::checked_narrowing_cast( detail::gamma_incomplete_imp(static_cast(a), static_cast(z), true, true, @@ -2127,8 +2046,6 @@ inline tools::promote_args_t policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; - detail::igamma_initializer::force_instantiate(); - return policies::checked_narrowing_cast( detail::gamma_incomplete_imp(static_cast(a), static_cast(z), true, false, diff --git a/test/test_gamma.hpp b/test/test_gamma.hpp index c21573dac6..5efa96d2f3 100644 --- a/test/test_gamma.hpp +++ b/test/test_gamma.hpp @@ -346,5 +346,29 @@ void test_spots(T, const char* name) value /= 2; } } + // + // Extra coverage: + // + BOOST_CHECK_THROW(boost::math::tgamma(T(-3)), std::domain_error); + BOOST_CHECK_THROW(boost::math::lgamma(T(-3)), std::domain_error); + BOOST_CHECK_THROW(boost::math::lgamma(T(0)), std::domain_error); + if(boost::math::tools::log_max_value() <= 11356) + { + BOOST_CHECK_GE(boost::math::tgamma(T(11360.0f)), boost::math::tools::max_value()); + BOOST_CHECK_EQUAL(boost::math::tgamma(T(-11360.5f)), T(0)); + BOOST_CHECK_GE(boost::math::tgamma(T(1755.75)), boost::math::tools::max_value()); + BOOST_CHECK_EQUAL(boost::math::tgamma(T(-1755.5)), T(0)); + } + if (boost::math::tools::log_max_value() <= 709) + { + BOOST_CHECK_GE(boost::math::tgamma(T(173)), boost::math::tools::max_value()); + // + // There is an error here, between -173 and -178.5 we should not underflow + // for type double, but do so spuriously because tgamma(-z) overflows. + // There is no easy fix for this, as it's next to impossible to figure out + // ahead of time when you're going to be "in the zone". + // + BOOST_CHECK_EQUAL(boost::math::tgamma(T(-178.5)), T(0)); + } } From c206a8d01983bba94276d8b52f62e120572ea101 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Mon, 11 Mar 2024 11:06:48 +0000 Subject: [PATCH 09/27] Allow for no-exceptions. --- include/boost/math/special_functions/gamma.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 535cd6b14e..1da88fc609 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -130,14 +130,18 @@ T gamma_imp(T z, const Policy& pol, const Lanczos& l) return policies::raise_pole_error(function, "Evaluation of tgamma at a negative integer %1%.", z, pol); if (z <= -20) { +#ifndef BOOST_MATH_NO_EXCEPTIONS try +#endif { result = gamma_imp(T(-z), pol, l) * sinpx(z); } +#ifndef BOOST_MATH_NO_EXCEPTIONS catch (const std::overflow_error&) { return policies::raise_underflow_error(function, "Result of tgamma is too small to represent.", pol); } +#endif BOOST_MATH_INSTRUMENT_VARIABLE(result); // Result can never be small: tgamma[-z] is always larger than sinpx[z] is small: BOOST_MATH_ASSERT((fabs(result) > 1) || (tools::max_value() * fabs(result) > boost::math::constants::pi())); From b90cdedeab35c06e12bddae99f9aebacc3d7a6f0 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 28 Jun 2024 19:39:37 +0100 Subject: [PATCH 10/27] Remove unneeded initializers. Clean up Bessel_ik. --- .../special_functions/detail/bessel_ik.hpp | 7 +-- .../special_functions/detail/bessel_j1.hpp | 27 ---------- .../special_functions/detail/bessel_k0.hpp | 35 ------------- .../special_functions/detail/bessel_k1.hpp | 36 -------------- .../special_functions/detail/polygamma.hpp | 24 +-------- .../detail/unchecked_factorial.hpp | 22 --------- .../boost/math/special_functions/log1p.hpp | 30 ------------ include/boost/math/special_functions/next.hpp | 27 ---------- .../boost/math/special_functions/owens_t.hpp | 30 ------------ .../boost/math/special_functions/trigamma.hpp | 32 ------------ include/boost/math/special_functions/zeta.hpp | 49 ------------------- 11 files changed, 5 insertions(+), 314 deletions(-) diff --git a/include/boost/math/special_functions/detail/bessel_ik.hpp b/include/boost/math/special_functions/detail/bessel_ik.hpp index 0c653b4753..7a99037085 100644 --- a/include/boost/math/special_functions/detail/bessel_ik.hpp +++ b/include/boost/math/special_functions/detail/bessel_ik.hpp @@ -326,6 +326,9 @@ int bessel_ik(T v, T x, T* result_I, T* result_K, int kind, const Policy& pol) T scale = 1; T scale_sign = 1; + n = iround(v, pol); + u = v - n; // -1/2 <= u < 1/2 + if (((kind & need_i) == 0) && (fabs(4 * v * v - 25) / (8 * x) < tools::forth_root_epsilon())) { // A&S 9.7.2 @@ -337,8 +340,6 @@ int bessel_ik(T v, T x, T* result_I, T* result_K, int kind, const Policy& pol) } else { - n = iround(v, pol); - u = v - n; // -1/2 <= u < 1/2 BOOST_MATH_INSTRUMENT_VARIABLE(n); BOOST_MATH_INSTRUMENT_VARIABLE(u); @@ -412,7 +413,7 @@ int bessel_ik(T v, T x, T* result_I, T* result_K, int kind, const Policy& pol) else Iv = std::numeric_limits::quiet_NaN(); // any value will do } - if (reflect) + if (reflect && (kind & need_i)) { T z = (u + n % 2); T fact = (2 / pi()) * (boost::math::sin_pi(z, pol) * Kv); diff --git a/include/boost/math/special_functions/detail/bessel_j1.hpp b/include/boost/math/special_functions/detail/bessel_j1.hpp index 6d354dcce7..90f26bcc84 100644 --- a/include/boost/math/special_functions/detail/bessel_j1.hpp +++ b/include/boost/math/special_functions/detail/bessel_j1.hpp @@ -34,36 +34,9 @@ namespace boost { namespace math{ namespace detail{ template T bessel_j1(T x); -template -struct bessel_j1_initializer -{ - struct init - { - init() - { - do_init(); - } - static void do_init() - { - bessel_j1(T(1)); - } - void force_instantiate()const{} - }; - static const init initializer; - static void force_instantiate() - { - initializer.force_instantiate(); - } -}; - -template -const typename bessel_j1_initializer::init bessel_j1_initializer::initializer; - template T bessel_j1(T x) { - bessel_j1_initializer::force_instantiate(); - static const T P1[] = { static_cast(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4258509801366645672e+11)), static_cast(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6781041261492395835e+09)), diff --git a/include/boost/math/special_functions/detail/bessel_k0.hpp b/include/boost/math/special_functions/detail/bessel_k0.hpp index f29ffa75c4..69a48b779f 100644 --- a/include/boost/math/special_functions/detail/bessel_k0.hpp +++ b/include/boost/math/special_functions/detail/bessel_k0.hpp @@ -46,40 +46,6 @@ namespace boost { namespace math { namespace detail{ template T bessel_k0(const T& x); -template -struct bessel_k0_initializer -{ - struct init - { - init() - { - do_init(tag()); - } - static void do_init(const std::integral_constant&) - { - bessel_k0(T(0.5)); - bessel_k0(T(1.5)); - } - static void do_init(const std::integral_constant&) - { - bessel_k0(T(0.5)); - bessel_k0(T(1.5)); - } - template - static void do_init(const U&){} - void force_instantiate()const{} - }; - static const init initializer; - static void force_instantiate() - { - initializer.force_instantiate(); - } -}; - -template -const typename bessel_k0_initializer::init bessel_k0_initializer::initializer; - - template T bessel_k0_imp(const T&, const std::integral_constant&) { @@ -505,7 +471,6 @@ inline T bessel_k0(const T& x) 113 : -1 > tag_type; - bessel_k0_initializer::force_instantiate(); return bessel_k0_imp(x, tag_type()); } diff --git a/include/boost/math/special_functions/detail/bessel_k1.hpp b/include/boost/math/special_functions/detail/bessel_k1.hpp index bd37f90215..b7259e467c 100644 --- a/include/boost/math/special_functions/detail/bessel_k1.hpp +++ b/include/boost/math/special_functions/detail/bessel_k1.hpp @@ -46,41 +46,6 @@ namespace boost { namespace math { namespace detail{ template T bessel_k1(const T&); - template - struct bessel_k1_initializer - { - struct init - { - init() - { - do_init(tag()); - } - static void do_init(const std::integral_constant&) - { - bessel_k1(T(0.5)); - bessel_k1(T(2)); - bessel_k1(T(6)); - } - static void do_init(const std::integral_constant&) - { - bessel_k1(T(0.5)); - bessel_k1(T(6)); - } - template - static void do_init(const U&) {} - void force_instantiate()const {} - }; - static const init initializer; - static void force_instantiate() - { - initializer.force_instantiate(); - } - }; - - template - const typename bessel_k1_initializer::init bessel_k1_initializer::initializer; - - template inline T bessel_k1_imp(const T&, const std::integral_constant&) { @@ -547,7 +512,6 @@ namespace boost { namespace math { namespace detail{ 113 : -1 > tag_type; - bessel_k1_initializer::force_instantiate(); return bessel_k1_imp(x, tag_type()); } diff --git a/include/boost/math/special_functions/detail/polygamma.hpp b/include/boost/math/special_functions/detail/polygamma.hpp index 8ca3663147..a608f22c7f 100644 --- a/include/boost/math/special_functions/detail/polygamma.hpp +++ b/include/boost/math/special_functions/detail/polygamma.hpp @@ -470,34 +470,12 @@ namespace boost { namespace math { namespace detail{ return exp(power_terms) * ((s < 0) && ((n + 1) & 1) ? -1 : 1) * boost::math::sign(sum); } - template - struct polygamma_initializer - { - struct init - { - init() - { - // Forces initialization of our table of coefficients and mutex: - boost::math::polygamma(30, T(-2.5f), Policy()); - } - void force_instantiate()const{} - }; - static const init initializer; - static void force_instantiate() - { - initializer.force_instantiate(); - } - }; - - template - const typename polygamma_initializer::init polygamma_initializer::initializer; - template inline T polygamma_imp(const int n, T x, const Policy &pol) { BOOST_MATH_STD_USING static const char* function = "boost::math::polygamma<%1%>(int, %1%)"; - polygamma_initializer::initializer.force_instantiate(); + if(n < 0) return policies::raise_domain_error(function, "Order must be >= 0, but got %1%", static_cast(n), pol); if(x < 0) diff --git a/include/boost/math/special_functions/detail/unchecked_factorial.hpp b/include/boost/math/special_functions/detail/unchecked_factorial.hpp index b528a24fe9..8428f7400c 100644 --- a/include/boost/math/special_functions/detail/unchecked_factorial.hpp +++ b/include/boost/math/special_functions/detail/unchecked_factorial.hpp @@ -1271,28 +1271,6 @@ struct max_factorial #endif -template -struct unchecked_factorial_initializer -{ - struct init - { - init() - { - boost::math::unchecked_factorial(3); - } - void force_instantiate()const {} - }; - static const init initializer; - static void force_instantiate() - { - initializer.force_instantiate(); - } -}; - -template -const typename unchecked_factorial_initializer::init unchecked_factorial_initializer::initializer; - - template inline T unchecked_factorial_imp(unsigned i, const std::integral_constant&) { diff --git a/include/boost/math/special_functions/log1p.hpp b/include/boost/math/special_functions/log1p.hpp index 9b8a8e0eb7..1fdd6dc52e 100644 --- a/include/boost/math/special_functions/log1p.hpp +++ b/include/boost/math/special_functions/log1p.hpp @@ -263,34 +263,6 @@ T log1p_imp(T const& x, const Policy& pol, const std::integral_constant return result; } -template -struct log1p_initializer -{ - struct init - { - init() - { - do_init(tag()); - } - template - static void do_init(const std::integral_constant&){} - static void do_init(const std::integral_constant&) - { - boost::math::log1p(static_cast(0.25), Policy()); - } - void force_instantiate()const{} - }; - static const init initializer; - static void force_instantiate() - { - initializer.force_instantiate(); - } -}; - -template -const typename log1p_initializer::init log1p_initializer::initializer; - - } // namespace detail template @@ -312,8 +284,6 @@ inline typename tools::promote_args::type log1p(T x, const Policy&) precision_type::value <= 64 ? 64 : 0 > tag_type; - detail::log1p_initializer::force_instantiate(); - return policies::checked_narrowing_cast( detail::log1p_imp(static_cast(x), forwarding_policy(), tag_type()), "boost::math::log1p<%1%>(%1%)"); } diff --git a/include/boost/math/special_functions/next.hpp b/include/boost/math/special_functions/next.hpp index 02a208e4eb..ee22ce78cd 100644 --- a/include/boost/math/special_functions/next.hpp +++ b/include/boost/math/special_functions/next.hpp @@ -121,31 +121,6 @@ inline bool has_denorm_now() { template T get_min_shift_value(); -template -struct min_shift_initializer -{ - struct init - { - init() - { - do_init(); - } - static void do_init() - { - get_min_shift_value(); - } - void force_instantiate()const{} - }; - static const init initializer; - static void force_instantiate() - { - initializer.force_instantiate(); - } -}; - -template -const typename min_shift_initializer::init min_shift_initializer::initializer; - template inline T calc_min_shifted(const std::true_type&) { @@ -166,8 +141,6 @@ template inline T get_min_shift_value() { static const T val = calc_min_shifted(std::integral_constant::is_specialized || std::numeric_limits::radix == 2>()); - min_shift_initializer::force_instantiate(); - return val; } diff --git a/include/boost/math/special_functions/owens_t.hpp b/include/boost/math/special_functions/owens_t.hpp index c794f0cacb..e09608b397 100644 --- a/include/boost/math/special_functions/owens_t.hpp +++ b/include/boost/math/special_functions/owens_t.hpp @@ -1026,34 +1026,6 @@ namespace boost return val; } // RealType owens_t(RealType h, RealType a) - template - struct owens_t_initializer - { - struct init - { - init() - { - do_init(tag()); - } - template - static void do_init(const std::integral_constant&){} - static void do_init(const std::integral_constant&) - { - boost::math::owens_t(static_cast(7), static_cast(0.96875), Policy()); - boost::math::owens_t(static_cast(2), static_cast(0.5), Policy()); - } - void force_instantiate()const{} - }; - static const init initializer; - static void force_instantiate() - { - initializer.force_instantiate(); - } - }; - - template - const typename owens_t_initializer::init owens_t_initializer::initializer; - } // namespace detail template @@ -1067,8 +1039,6 @@ namespace boost precision_type::value <= 64 ? 64 : 65 > tag_type; - detail::owens_t_initializer::force_instantiate(); - return policies::checked_narrowing_cast(detail::owens_t(static_cast(h), static_cast(a), pol), "boost::math::owens_t<%1%>(%1%,%1%)"); } diff --git a/include/boost/math/special_functions/trigamma.hpp b/include/boost/math/special_functions/trigamma.hpp index f74b43db1f..d40be39cbd 100644 --- a/include/boost/math/special_functions/trigamma.hpp +++ b/include/boost/math/special_functions/trigamma.hpp @@ -394,35 +394,6 @@ T trigamma_imp(T x, const std::integral_constant*, const Policy& pol) { return polygamma_imp(1, x, pol); } -// -// Initializer: ensure all our constants are initialized prior to the first call of main: -// -template -struct trigamma_initializer -{ - struct init - { - init() - { - typedef typename policies::precision::type precision_type; - do_init(std::integral_constant()); - } - void do_init(const std::true_type&) - { - boost::math::trigamma(T(2.5), Policy()); - } - void do_init(const std::false_type&){} - void force_instantiate()const{} - }; - static const init initializer; - static void force_instantiate() - { - initializer.force_instantiate(); - } -}; - -template -const typename trigamma_initializer::init trigamma_initializer::initializer; } // namespace detail @@ -446,9 +417,6 @@ inline typename tools::promote_args::type policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; - // Force initialization of constants: - detail::trigamma_initializer::force_instantiate(); - return policies::checked_narrowing_cast(detail::trigamma_imp( static_cast(x), static_cast(nullptr), forwarding_policy()), "boost::math::trigamma<%1%>(%1%)"); diff --git a/include/boost/math/special_functions/zeta.hpp b/include/boost/math/special_functions/zeta.hpp index 1a3bba7485..de74d3d05c 100644 --- a/include/boost/math/special_functions/zeta.hpp +++ b/include/boost/math/special_functions/zeta.hpp @@ -1010,53 +1010,6 @@ T zeta_imp(T s, T sc, const Policy& pol, const Tag& tag) return result; } -template -struct zeta_initializer -{ - struct init - { - init() - { - do_init(tag()); - } - static void do_init(const std::integral_constant&){ boost::math::zeta(static_cast(5), Policy()); } - static void do_init(const std::integral_constant&){ boost::math::zeta(static_cast(5), Policy()); } - static void do_init(const std::integral_constant&) - { - boost::math::zeta(static_cast(0.5), Policy()); - boost::math::zeta(static_cast(1.5), Policy()); - boost::math::zeta(static_cast(3.5), Policy()); - boost::math::zeta(static_cast(6.5), Policy()); - boost::math::zeta(static_cast(14.5), Policy()); - boost::math::zeta(static_cast(40.5), Policy()); - - boost::math::zeta(static_cast(5), Policy()); - } - static void do_init(const std::integral_constant&) - { - boost::math::zeta(static_cast(0.5), Policy()); - boost::math::zeta(static_cast(1.5), Policy()); - boost::math::zeta(static_cast(3.5), Policy()); - boost::math::zeta(static_cast(5.5), Policy()); - boost::math::zeta(static_cast(9.5), Policy()); - boost::math::zeta(static_cast(16.5), Policy()); - boost::math::zeta(static_cast(25.5), Policy()); - boost::math::zeta(static_cast(70.5), Policy()); - - boost::math::zeta(static_cast(5), Policy()); - } - void force_instantiate()const{} - }; - static const init initializer; - static void force_instantiate() - { - initializer.force_instantiate(); - } -}; - -template -const typename zeta_initializer::init zeta_initializer::initializer; - } // detail template @@ -1078,8 +1031,6 @@ inline typename tools::promote_args::type zeta(T s, const Policy&) precision_type::value <= 113 ? 113 : 0 > tag_type; - detail::zeta_initializer::force_instantiate(); - return policies::checked_narrowing_cast(detail::zeta_imp( static_cast(s), static_cast(1 - static_cast(s)), From 0c6522fc29fb88424d335e7b22428445c1cb0aea Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sat, 29 Jun 2024 10:18:36 +0100 Subject: [PATCH 11/27] Mark up last missing erf.hpp line: it's covered. --- include/boost/math/special_functions/erf.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/math/special_functions/erf.hpp b/include/boost/math/special_functions/erf.hpp index 54964617bf..2b6d5603eb 100644 --- a/include/boost/math/special_functions/erf.hpp +++ b/include/boost/math/special_functions/erf.hpp @@ -716,7 +716,7 @@ T erf_imp(T z, bool invert, const Policy& pol, const std::integral_constant Date: Sat, 29 Jun 2024 11:50:50 +0100 Subject: [PATCH 12/27] tgamma coverage #1. --- test/Jamfile.v2 | 1 + test/test_tgamma_eh.cpp | 16 ++++++++++++++++ 2 files changed, 17 insertions(+) create mode 100644 test/test_tgamma_eh.cpp diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index d52604aa2f..a9a14fca21 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -230,6 +230,7 @@ test-suite special_fun : [ run test_factorials.cpp pch ../../test/build//boost_unit_test_framework ] [ run test_gamma.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] [ run test_gamma_edge.cpp ] + [ run test_tgamma_eh.cpp ] [ run test_gamma_mp.cpp ../../test/build//boost_unit_test_framework : : : release TEST=1 [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] : test_gamma_mp_1 ] [ run test_gamma_mp.cpp ../../test/build//boost_unit_test_framework : : : release TEST=2 [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] : test_gamma_mp_2 ] [ run test_gamma_mp.cpp ../../test/build//boost_unit_test_framework : : : release TEST=3 [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] : test_gamma_mp_3 ] diff --git a/test/test_tgamma_eh.cpp b/test/test_tgamma_eh.cpp new file mode 100644 index 0000000000..b34ec4bf8b --- /dev/null +++ b/test/test_tgamma_eh.cpp @@ -0,0 +1,16 @@ +// (C) Copyright John Maddock 2024. +// 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) + +#include "math_unit_test.hpp" +#include +#include +#include + + +int main() +{ + CHECK_EQUAL(boost::math::tgamma(-200.5), 0.0); // triggers internal exception handling + return boost::math::test::report_errors(); +} From f60b35f217d6a5a0a7039f383a3c0863433b6682 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sat, 29 Jun 2024 19:29:51 +0100 Subject: [PATCH 13/27] More tgamma coverage. --- .../boost/math/special_functions/gamma.hpp | 30 +++++++++---------- test/test_gamma.hpp | 5 ++++ 2 files changed, 19 insertions(+), 16 deletions(-) diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 1da88fc609..919a15c4c7 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -529,16 +529,16 @@ T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&) if (!n_recur) { if (zz > tools::log_max_value()) - return b_neg ? policies::raise_underflow_error(function, nullptr, pol) : policies::raise_overflow_error(function, nullptr, pol); + return b_neg ? policies::raise_underflow_error(function, nullptr, pol) : policies::raise_overflow_error(function, nullptr, pol); // LCOV_EXCL_LINE MP only if (log(zz) * zz / 2 > tools::log_max_value()) - return b_neg ? policies::raise_underflow_error(function, nullptr, pol) : policies::raise_overflow_error(function, nullptr, pol); + return b_neg ? policies::raise_underflow_error(function, nullptr, pol) : policies::raise_overflow_error(function, nullptr, pol); // LCOV_EXCL_LINE MP only } T gamma_value = scaled_tgamma_no_lanczos(zz, pol); T power_term = pow(zz, zz / 2); T exp_term = exp(-zz); gamma_value *= (power_term * exp_term); if (!n_recur && (tools::max_value() / power_term < gamma_value)) - return b_neg ? policies::raise_underflow_error(function, nullptr, pol) : policies::raise_overflow_error(function, nullptr, pol); + return b_neg ? policies::raise_underflow_error(function, nullptr, pol) : policies::raise_overflow_error(function, nullptr, pol); // LCOV_EXCL_LINE MP only gamma_value *= power_term; // Rescale the result using downward recursion if necessary. @@ -571,21 +571,21 @@ T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&) T s = sinpx(z); if ((gamma_value > 1) && (tools::max_value() / gamma_value < fabs(s))) - return policies::raise_underflow_error(function, nullptr, pol); - gamma_value *= s; + return policies::raise_underflow_error(function, nullptr, pol); // LCOV_EXCL_LINE MP only + gamma_value *= s; // LCOV_EXCL_LINE MP only BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value); // // Result can never overflow, since sinpx(z) can never be smaller than machine epsilon and gamma_value > 1. // - BOOST_MATH_ASSERT( (abs(gamma_value) > 1) || ((tools::max_value() * abs(gamma_value)) > boost::math::constants::pi())); + BOOST_MATH_ASSERT( (abs(gamma_value) > 1) || ((tools::max_value() * abs(gamma_value)) > boost::math::constants::pi())); // LCOV_EXCL_LINE MP only gamma_value = -boost::math::constants::pi() / gamma_value; - BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value); + BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value); // LCOV_EXCL_LINE MP only // // We can never underflow since the numerator > 1 above: // - BOOST_MATH_ASSERT(gamma_value != 0); + BOOST_MATH_ASSERT(gamma_value != 0); // LCOV_EXCL_LINE MP only } return gamma_value; @@ -648,7 +648,7 @@ T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sig if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial::value)) { if (sign) - *sign = 1; + *sign = 1; // LCOV_EXCL_LINE return log(boost::math::unchecked_factorial(itrunc(z) - 1)); } @@ -665,7 +665,7 @@ T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sig // inefficient, but simple. The rationale is that the argument here // is relatively small and overflow is not expected to be likely. if (sign) - * sign = 1; + * sign = 1; // // LCOV_EXCL_LINE if(fabs(z - 1) < 0.25) { log_gamma_value = log_gamma_near_1(T(zz - 1), pol); @@ -679,7 +679,7 @@ T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sig // Reflection formula may fail if z is very close to zero, let the series // expansion for tgamma close to zero do the work: if (sign) - *sign = z < 0 ? -1 : 1; + *sign = z < 0 ? -1 : 1; // LCOV_EXCL_LINE return log(abs(gamma_imp(z, pol, lanczos::undefined_lanczos()))); } else @@ -710,7 +710,7 @@ T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sig // Check if the argument of lgamma is exactly equal to a negative integer. if(floor_of_z_is_equal_to_z) - return policies::raise_pole_error(function, "Evaluation of lgamma at a negative integer %1%.", z, pol); + return policies::raise_pole_error(function, "Evaluation of lgamma at a negative integer %1%.", z, pol); // LCOV_EXCL_LINE MP only T t = sinpx(z); @@ -720,12 +720,10 @@ T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sig } else { - sign_of_result = -sign_of_result; + sign_of_result = -sign_of_result; // LCOV_EXCL_LINE MP only } - log_gamma_value = - log_gamma_value - + log(boost::math::constants::pi()) - - log(t); + log_gamma_value = - log_gamma_value + log(boost::math::constants::pi()) - log(t); } if(sign != static_cast(nullptr)) { *sign = sign_of_result; } diff --git a/test/test_gamma.hpp b/test/test_gamma.hpp index 5efa96d2f3..3efc7e9604 100644 --- a/test/test_gamma.hpp +++ b/test/test_gamma.hpp @@ -370,5 +370,10 @@ void test_spots(T, const char* name) // BOOST_CHECK_EQUAL(boost::math::tgamma(T(-178.5)), T(0)); } + // + // Extra coverage for tgamma1pm1: + // + BOOST_CHECK_CLOSE(::boost::math::tgamma1pm1(T(-2.5)), boost::math::tgamma(T(-1.5)) - 1, tolerance); + BOOST_CHECK_CLOSE(::boost::math::tgamma1pm1(T(2.5)), boost::math::tgamma(T(3.5)) - 1, tolerance); } From 87f242f707b61dfe1e4a9e938de2cce4854f012b Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sun, 30 Jun 2024 12:04:52 +0100 Subject: [PATCH 14/27] gamma.hpp coverage. --- .../boost/math/special_functions/gamma.hpp | 60 ++++++++++++------- test/test_igamma.cpp | 2 +- test/test_igamma.hpp | 16 ++++- test/test_tgamma_eh.cpp | 1 + test/test_tgamma_ratio.hpp | 6 ++ 5 files changed, 61 insertions(+), 24 deletions(-) diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 919a15c4c7..523e1bc84d 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -567,7 +567,7 @@ T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&) // Check if the argument of tgamma is exactly equal to a negative integer. if(floor_of_z_is_equal_to_z) - return policies::raise_pole_error(function, "Evaluation of tgamma at a negative integer %1%.", z, pol); + return policies::raise_pole_error(function, "Evaluation of tgamma at a negative integer %1%.", z, pol); // LCOV_EXCL_LINE MP only T s = sinpx(z); if ((gamma_value > 1) && (tools::max_value() / gamma_value < fabs(s))) @@ -689,7 +689,7 @@ T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sig T g = gamma_imp(zz, pol, lanczos::undefined_lanczos()); if (sign) { - *sign = g < 0 ? -1 : 1; + *sign = g < 0 ? -1 : 1; // LCOV_EXCL_LINE MP only } log_gamma_value = log(abs(g)); } @@ -748,7 +748,7 @@ T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& l) precision_type::value <= 113 ? 113 : 0 > tag_type; - T result; + T result{}; if(dz < 0) { if(dz < T(-0.5)) @@ -848,7 +848,7 @@ T full_igamma_prefix(T a, T z, const Policy& pol) } else { - prefix = exp(alz - z); + prefix = exp(alz - z); // LCOV_EXCL_LINE defensive programming, can probably never get here? } } else @@ -857,6 +857,8 @@ T full_igamma_prefix(T a, T z, const Policy& pol) { prefix = pow(z, a) * exp(-z); } + // LCOV_EXCL_START + // Defensive programming, can probably never get here, very hard to prove though! else if(z/a < tools::log_max_value()) { prefix = pow(T(z / exp(z/a)), a); @@ -865,13 +867,15 @@ T full_igamma_prefix(T a, T z, const Policy& pol) { prefix = exp(alz - z); } + // LCOV_EXCL_STOP } // // This error handling isn't very good: it happens after the fact // rather than before it... + // Typically though this method is used when the result is small, we should probably not overflow here... // if((boost::math::fpclassify)(prefix) == (int)FP_INFINITE) - return policies::raise_overflow_error("boost::math::detail::full_igamma_prefix<%1%>(%1%, %1%)", "Result of incomplete gamma function is too large to represent.", pol); + return policies::raise_overflow_error("boost::math::detail::full_igamma_prefix<%1%>(%1%, %1%)", "Result of incomplete gamma function is too large to represent.", pol); // LCOV_EXCL_LINE return prefix; } @@ -886,7 +890,7 @@ T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l) if (z >= tools::max_value()) return 0; T agh = a + static_cast(Lanczos::g()) - T(0.5); - T prefix; + T prefix{}; T d = ((z - a) - static_cast(Lanczos::g()) + T(0.5)) / agh; if(a < 1) @@ -1102,7 +1106,7 @@ T finite_half_gamma_q(T a, T x, T* p_derivative, const Policy& pol) { T term = exp(-x) / sqrt(constants::pi() * x); term *= x; - static const T half = T(1) / 2; + static const T half = T(1) / 2; // LCOV_EXCL_LINE term /= half; T sum = term; for(unsigned n = 2; n < a; ++n) @@ -1190,16 +1194,24 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, { // This is method 4 below, done in logs: result = a * log(x) - x; + /* + * We do not compute derivatives in the non-regularized case. if(p_derivative) *p_derivative = exp(result); + */ + BOOST_MATH_ASSERT(nullptr == p_derivative); result += log(upper_gamma_fraction(a, x, policies::get_epsilon())); } else if(!invert && (a > 4 * x)) { // This is method 2 below, done in logs: result = a * log(x) - x; + /* + * We do not compute derivatives in the non-regularized case. if(p_derivative) *p_derivative = exp(result); + */ + BOOST_MATH_ASSERT(nullptr == p_derivative); T init_value = 0; result += log(detail::lower_gamma_series(a, x, pol, init_value) / a); } @@ -1213,8 +1225,12 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, // Try http://functions.wolfram.com/06.06.06.0039.01 result = 1 + 1 / (12 * a) + 1 / (288 * a * a); result = log(result) - a + (a - 0.5f) * log(a) + log(boost::math::constants::root_two_pi()); + /* + * We do not compute derivatives in the non-regularized case. if(p_derivative) *p_derivative = exp(a * log(x) - x); + */ + BOOST_MATH_ASSERT(nullptr == p_derivative); } else { @@ -1222,8 +1238,12 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, // range of this method, but since the result is almost certainly // infinite, we should probably be OK: result = a * log(x) - x; + /* + * We do not compute derivatives in the non-regularized case. if(p_derivative) *p_derivative = exp(result); + */ + BOOST_MATH_ASSERT(nullptr == p_derivative); T init_value = 0; result += log(detail::lower_gamma_series(a, x, pol, init_value) / a); } @@ -1434,7 +1454,7 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, { // Compute Q: invert = !invert; - T g; + T g{}; result = tgamma_small_upper_part(a, x, pol, &g, invert, p_derivative); invert = false; if(normalised) @@ -1572,21 +1592,19 @@ T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos& } } T zgh = static_cast(z + T(Lanczos::g()) - constants::half()); - T result; + T result{}; if(z + delta == z) { - if (fabs(delta / zgh) < boost::math::tools::epsilon()) - { - // We have: - // result = exp((constants::half() - z) * boost::math::log1p(delta / zgh, pol)); - // 0.5 - z == -z - // log1p(delta / zgh) = delta / zgh = delta / z - // multiplying we get -delta. - result = exp(-delta); - } - else - // from the pow formula below... but this may actually be wrong, we just can't really calculate it :( - result = 1; + // Given delta < z * eps + // and zgh > z + // Then this must follow: + BOOST_MATH_ASSERT(fabs(delta / zgh) < boost::math::tools::epsilon()); + // We have: + // result = exp((constants::half() - z) * boost::math::log1p(delta / zgh, pol)); + // 0.5 - z == -z + // log1p(delta / zgh) = delta / zgh = delta / z + // multiplying we get -delta. + result = exp(-delta); } else { diff --git a/test/test_igamma.cpp b/test/test_igamma.cpp index 8e80c772c4..fc6f133e65 100644 --- a/test/test_igamma.cpp +++ b/test/test_igamma.cpp @@ -281,7 +281,7 @@ void expected_results() // // Large exponent range causes more extreme test cases to be evaluated: // - if(std::numeric_limits::max_exponent > std::numeric_limits::max_exponent) + BOOST_IF_CONSTEXPR(std::numeric_limits::max_exponent > std::numeric_limits::max_exponent) { add_expected_result( "[^|]*", // compiler diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp index b434f727ee..2a24a94b27 100644 --- a/test/test_igamma.hpp +++ b/test/test_igamma.hpp @@ -202,7 +202,7 @@ void test_spots(T) //typedef boost::math::policies::policy > throw_policy; - if(std::numeric_limits::max_exponent <= 1024 && std::numeric_limits::has_infinity) + BOOST_IF_CONSTEXPR(std::numeric_limits::max_exponent <= 1024 && std::numeric_limits::has_infinity) { BOOST_CHECK_EQUAL(::boost::math::tgamma(static_cast(176), static_cast(100)), std::numeric_limits::infinity()); //BOOST_MATH_CHECK_THROW(::boost::math::tgamma(static_cast(176), static_cast(100), throw_policy()), std::overflow_error); @@ -212,7 +212,7 @@ void test_spots(T) BOOST_CHECK_EQUAL(::boost::math::tgamma(static_cast(740.5), static_cast(2500)), std::numeric_limits::infinity()); BOOST_CHECK_EQUAL(::boost::math::tgamma_lower(static_cast(10000.0f), static_cast(10000.0f / 4)), std::numeric_limits::infinity()); } - if(std::numeric_limits::max_exponent >= 1024) + BOOST_IF_CONSTEXPR(std::numeric_limits::max_exponent >= 1024) { BOOST_CHECK_CLOSE(::boost::math::tgamma(static_cast(170), static_cast(165)), static_cast(2.737338337642022829223832094019477918166996032112404370e304L), (std::numeric_limits::digits > 100 ? 10 : 3) * tolerance); BOOST_CHECK_CLOSE(::boost::math::tgamma_lower(static_cast(170), static_cast(165)), static_cast(1.531729671362682445715419794880088619901822603944331733e304L), 3 * tolerance); @@ -254,5 +254,17 @@ void test_spots(T) // BOOST_CHECK_EQUAL(::boost::math::gamma_q(static_cast(1770), static_cast(1e-12)), 1); BOOST_CHECK_EQUAL(::boost::math::gamma_p(static_cast(1770), static_cast(1e-12)), 0); + // + // Coverage: + // + BOOST_CHECK_THROW(boost::math::gamma_p(static_cast(-1), static_cast(2)), std::domain_error); + BOOST_CHECK_THROW(boost::math::gamma_q(static_cast(-1), static_cast(2)), std::domain_error); + BOOST_CHECK_THROW(boost::math::gamma_p(static_cast(1), static_cast(-2)), std::domain_error); + BOOST_CHECK_THROW(boost::math::gamma_q(static_cast(1), static_cast(-2)), std::domain_error); + BOOST_CHECK_THROW(boost::math::gamma_p(static_cast(0), static_cast(2)), std::domain_error); + BOOST_CHECK_THROW(boost::math::gamma_q(static_cast(0), static_cast(2)), std::domain_error); + BOOST_CHECK_THROW(boost::math::gamma_p_derivative(static_cast(-1), static_cast(2)), std::domain_error); + BOOST_CHECK_THROW(boost::math::gamma_p_derivative(static_cast(1), static_cast(-2)), std::domain_error); + BOOST_CHECK_THROW(boost::math::gamma_p_derivative(static_cast(0), static_cast(2)), std::domain_error); } diff --git a/test/test_tgamma_eh.cpp b/test/test_tgamma_eh.cpp index b34ec4bf8b..d27e7e5c1e 100644 --- a/test/test_tgamma_eh.cpp +++ b/test/test_tgamma_eh.cpp @@ -12,5 +12,6 @@ int main() { CHECK_EQUAL(boost::math::tgamma(-200.5), 0.0); // triggers internal exception handling + CHECK_EQUAL(boost::math::gamma_p(500.125, 1e-50), 0.0); // triggers internal exception handling return boost::math::test::report_errors(); } diff --git a/test/test_tgamma_ratio.hpp b/test/test_tgamma_ratio.hpp index bbcf9ea67d..86901371ad 100644 --- a/test/test_tgamma_ratio.hpp +++ b/test/test_tgamma_ratio.hpp @@ -180,4 +180,10 @@ void test_spots(T, const char*) BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_delta_ratio(T(200), T(ldexp(T(1), -1074))), T(1), tol); } } + // + // Coverage: + // + BOOST_CHECK_EQUAL(boost::math::tgamma_delta_ratio(T(2), T(0)), T(1)) + BOOST_CHECK_EQUAL(boost::math::tgamma_delta_ratio(T(200), T(0)), T(1)) + BOOST_CHECK_EQUAL(boost::math::tgamma_delta_ratio(T(2000), T(0)), T(1)) } From d441001a5eaedf2a2a9fbaa9472937075ac3b435 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sun, 30 Jun 2024 12:49:07 +0100 Subject: [PATCH 15/27] Correct test. --- test/test_tgamma_ratio.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_tgamma_ratio.hpp b/test/test_tgamma_ratio.hpp index 86901371ad..32a109290b 100644 --- a/test/test_tgamma_ratio.hpp +++ b/test/test_tgamma_ratio.hpp @@ -183,7 +183,7 @@ void test_spots(T, const char*) // // Coverage: // - BOOST_CHECK_EQUAL(boost::math::tgamma_delta_ratio(T(2), T(0)), T(1)) - BOOST_CHECK_EQUAL(boost::math::tgamma_delta_ratio(T(200), T(0)), T(1)) - BOOST_CHECK_EQUAL(boost::math::tgamma_delta_ratio(T(2000), T(0)), T(1)) + BOOST_CHECK_EQUAL(boost::math::tgamma_delta_ratio(T(2), T(0)), T(1)); + BOOST_CHECK_EQUAL(boost::math::tgamma_delta_ratio(T(200), T(0)), T(1)); + BOOST_CHECK_EQUAL(boost::math::tgamma_delta_ratio(T(2000), T(0)), T(1)); } From d8a4900e59560025fe83b0b540e1d2abd073db2a Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Mon, 1 Jul 2024 12:25:15 +0100 Subject: [PATCH 16/27] Mark up zeta constants for coverage. --- include/boost/math/special_functions/zeta.hpp | 52 ++++++++++++++++--- 1 file changed, 44 insertions(+), 8 deletions(-) diff --git a/include/boost/math/special_functions/zeta.hpp b/include/boost/math/special_functions/zeta.hpp index de74d3d05c..453e59b878 100644 --- a/include/boost/math/special_functions/zeta.hpp +++ b/include/boost/math/special_functions/zeta.hpp @@ -209,6 +209,7 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant(0.24339294433593750202L), static_cast(-0.49092470516353571651L), @@ -225,6 +226,7 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant(0.00024978985622317935355L), static_cast(-0.101855788418564031874e-4L), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc); result -= 1.2433929443359375F; result += (sc); @@ -234,6 +236,7 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant(0.577215664901532860516L), static_cast(0.243210646940107164097L), @@ -250,6 +253,7 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant(0.000255784226140488490982L), static_cast(0.10991819782396112081e-4L), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, T(-sc)) / tools::evaluate_polynomial(Q, T(-sc)); result += 1 / (-sc); } @@ -257,6 +261,7 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant(-0.0537258300023595030676L), @@ -275,6 +280,7 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant(0.106951867532057341359e-4L), static_cast(0.236276623974978646399e-7L), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, T(s - 2)) / tools::evaluate_polynomial(Q, T(s - 2)); result += Y + 1 / (-sc); } @@ -283,7 +289,7 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant(-2.49710190602259410021L), static_cast(-2.60013301809475665334L), @@ -303,6 +309,7 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant(0.718833729365459760664e-8L), static_cast(-0.1129200113474947419e-9L), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, T(s - 4)) / tools::evaluate_polynomial(Q, T(s - 4)); result = 1 + exp(result); } @@ -311,6 +318,7 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant(-4.78558028495135619286L), static_cast(-1.89197364881972536382L), @@ -331,6 +339,7 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant(-0.833378440625385520576e-10L), static_cast(0.699841545204845636531e-12L), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, T(s - 7)) / tools::evaluate_polynomial(Q, T(s - 7)); result = 1 + exp(result); } @@ -338,6 +347,7 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant(-10.3948950573308896825L), static_cast(-2.85827219671106697179L), @@ -358,6 +368,7 @@ inline T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant(0.118507153474022900583e-7L), static_cast(0.222609483627352615142e-14L), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, T(s - 15)) / tools::evaluate_polynomial(Q, T(s - 15)); result = 1 + exp(result); } @@ -383,6 +394,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant& // Maximum Deviation Found: 3.099e-20 // Expected Error Term: 3.099e-20 // Max error found at long double precision: 5.890498e-20 + // LCOV_EXCL_START static const T P[6] = { BOOST_MATH_BIG_CONSTANT(T, 64, 0.243392944335937499969), BOOST_MATH_BIG_CONSTANT(T, 64, -0.496837806864865688082), @@ -400,6 +412,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant& BOOST_MATH_BIG_CONSTANT(T, 64, -0.159600883054550987633e-4), BOOST_MATH_BIG_CONSTANT(T, 64, 0.339770279812410586032e-6), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc); result -= 1.2433929443359375F; result += (sc); @@ -410,7 +423,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant& // Maximum Deviation Found: 1.059e-21 // Expected Error Term: 1.059e-21 // Max error found at long double precision: 1.626303e-19 - + // LCOV_EXCL_START static const T P[6] = { BOOST_MATH_BIG_CONSTANT(T, 64, 0.577215664901532860605), BOOST_MATH_BIG_CONSTANT(T, 64, 0.222537368917162139445), @@ -428,6 +441,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant& BOOST_MATH_BIG_CONSTANT(T, 64, 0.635994377921861930071e-5), BOOST_MATH_BIG_CONSTANT(T, 64, 0.226583954978371199405e-7), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, T(-sc)) / tools::evaluate_polynomial(Q, T(-sc)); result += 1 / (-sc); } @@ -435,6 +449,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant& { // Maximum Deviation Found: 5.946e-22 // Expected Error Term: -5.946e-22 + // LCOV_EXCL_START static const float Y = 0.6986598968505859375; static const T P[7] = { BOOST_MATH_BIG_CONSTANT(T, 64, -0.053725830002359501027), @@ -455,12 +470,14 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant& BOOST_MATH_BIG_CONSTANT(T, 64, 0.278090318191657278204e-6), BOOST_MATH_BIG_CONSTANT(T, 64, -0.19683620233222028478e-8), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, T(s - 2)) / tools::evaluate_polynomial(Q, T(s - 2)); result += Y + 1 / (-sc); } else if(s <= 7) { // Max error found at long double precision: 8.132216e-19 + // LCOV_EXCL_START static const T P[8] = { BOOST_MATH_BIG_CONSTANT(T, 64, -2.49710190602259407065), BOOST_MATH_BIG_CONSTANT(T, 64, -3.36664913245960625334), @@ -482,6 +499,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant& BOOST_MATH_BIG_CONSTANT(T, 64, -0.927884739284359700764e-8), BOOST_MATH_BIG_CONSTANT(T, 64, 0.119810501805618894381e-9), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, T(s - 4)) / tools::evaluate_polynomial(Q, T(s - 4)); result = 1 + exp(result); } @@ -489,6 +507,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant& { // Max error in interpolated form: 1.133e-18 // Max error found at long double precision: 2.183198e-18 + // LCOV_EXCL_START static const T P[9] = { BOOST_MATH_BIG_CONSTANT(T, 64, -4.78558028495135548083), BOOST_MATH_BIG_CONSTANT(T, 64, -3.23873322238609358947), @@ -511,6 +530,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant& BOOST_MATH_BIG_CONSTANT(T, 64, 0.117957556472335968146e-7), BOOST_MATH_BIG_CONSTANT(T, 64, -0.193432300973017671137e-12), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, T(s - 7)) / tools::evaluate_polynomial(Q, T(s - 7)); result = 1 + exp(result); } @@ -518,6 +538,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant& { // Max error in interpolated form: 1.668e-17 // Max error found at long double precision: 1.669714e-17 + // LCOV_EXCL_START static const T P[9] = { BOOST_MATH_BIG_CONSTANT(T, 64, -10.3948950573308861781), BOOST_MATH_BIG_CONSTANT(T, 64, -2.82646012777913950108), @@ -541,6 +562,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant& BOOST_MATH_BIG_CONSTANT(T, 64, -0.939798249922234703384e-16), BOOST_MATH_BIG_CONSTANT(T, 64, 0.264584017421245080294e-18), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, T(s - 15)) / tools::evaluate_polynomial(Q, T(s - 15)); result = 1 + exp(result); } @@ -566,7 +588,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant // Maximum Deviation Found: 9.493e-37 // Expected Error Term: 9.492e-37 // Max error found at long double precision: 7.281332e-31 - + // LCOV_EXCL_START static const T P[10] = { BOOST_MATH_BIG_CONSTANT(T, 113, -1.0), BOOST_MATH_BIG_CONSTANT(T, 113, -0.0353008629988648122808504280990313668), @@ -592,6 +614,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant BOOST_MATH_BIG_CONSTANT(T, 113, -0.217062446168217797598596496310953025e-9), BOOST_MATH_BIG_CONSTANT(T, 113, 0.315823200002384492377987848307151168e-11), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc); result += (sc); result /= (sc); @@ -600,7 +623,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant { // Maximum Deviation Found: 1.616e-37 // Expected Error Term: -1.615e-37 - + // LCOV_EXCL_START static const T P[10] = { BOOST_MATH_BIG_CONSTANT(T, 113, 0.577215664901532860606512090082402431), BOOST_MATH_BIG_CONSTANT(T, 113, 0.255597968739771510415479842335906308), @@ -626,6 +649,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant BOOST_MATH_BIG_CONSTANT(T, 113, -0.377105263588822468076813329270698909e-11), BOOST_MATH_BIG_CONSTANT(T, 113, -0.581926559304525152432462127383600681e-13), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, T(-sc)) / tools::evaluate_polynomial(Q, T(-sc)); result += 1 / (-sc); } @@ -634,7 +658,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant // Maximum Deviation Found: 1.891e-36 // Expected Error Term: -1.891e-36 // Max error found: 2.171527e-35 - + // LCOV_EXCL_START static const float Y = 0.6986598968505859375; static const T P[11] = { BOOST_MATH_BIG_CONSTANT(T, 113, -0.0537258300023595010275848333539748089), @@ -663,6 +687,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant BOOST_MATH_BIG_CONSTANT(T, 113, 0.105677416606909614301995218444080615e-11), BOOST_MATH_BIG_CONSTANT(T, 113, 0.547223964564003701979951154093005354e-15), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, T(s - 2)) / tools::evaluate_polynomial(Q, T(s - 2)); result += Y + 1 / (-sc); } @@ -670,7 +695,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant { // Max error in interpolated form: 1.510e-37 // Max error found at long double precision: 2.769266e-34 - + // LCOV_EXCL_START static const T Y = 3.28348541259765625F; static const T P[13] = { @@ -704,6 +729,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant BOOST_MATH_BIG_CONSTANT(T, 113, -0.275363878344548055574209713637734269e-13), BOOST_MATH_BIG_CONSTANT(T, 113, 0.221564186807357535475441900517843892e-15), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, T(s - 4)) / tools::evaluate_polynomial(Q, T(s - 4)); result -= Y; result = 1 + exp(result); @@ -712,7 +738,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant { // Max error in interpolated form: 1.999e-34 // Max error found at long double precision: 2.156186e-33 - + // LCOV_EXCL_START static const T P[13] = { BOOST_MATH_BIG_CONSTANT(T, 113, -4.0545627381873738086704293881227365), BOOST_MATH_BIG_CONSTANT(T, 113, -4.70088348734699134347906176097717782), @@ -744,6 +770,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant BOOST_MATH_BIG_CONSTANT(T, 113, 0.294670713571839023181857795866134957e-16), BOOST_MATH_BIG_CONSTANT(T, 113, -0.147003914536437243143096875069813451e-18), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, T(s - 6)) / tools::evaluate_polynomial(Q, T(s - 6)); result = 1 + exp(result); } @@ -751,6 +778,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant { // Max error in interpolated form: 1.641e-32 // Max error found at long double precision: 1.696121e-32 + // LCOV_EXCL_START static const T P[13] = { BOOST_MATH_BIG_CONSTANT(T, 113, -6.91319491921722925920883787894829678), BOOST_MATH_BIG_CONSTANT(T, 113, -3.65491257639481960248690596951049048), @@ -782,6 +810,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant BOOST_MATH_BIG_CONSTANT(T, 113, 0.887948682401000153828241615760146728e-19), BOOST_MATH_BIG_CONSTANT(T, 113, -0.34980761098820347103967203948619072e-21), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, T(s - 10)) / tools::evaluate_polynomial(Q, T(s - 10)); result = 1 + exp(result); } @@ -789,7 +818,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant { // Max error in interpolated form: 1.563e-31 // Max error found at long double precision: 1.562725e-31 - + // LCOV_EXCL_START static const T P[13] = { BOOST_MATH_BIG_CONSTANT(T, 113, -11.7824798233959252791987402769438322), BOOST_MATH_BIG_CONSTANT(T, 113, -4.36131215284987731928174218354118102), @@ -821,6 +850,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant BOOST_MATH_BIG_CONSTANT(T, 113, -0.291354445847552426900293580511392459e-22), BOOST_MATH_BIG_CONSTANT(T, 113, 0.73614324724785855925025452085443636e-25), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, T(s - 17)) / tools::evaluate_polynomial(Q, T(s - 17)); result = 1 + exp(result); } @@ -828,6 +858,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant { // Max error in interpolated form: 2.311e-27 // Max error found at long double precision: 2.297544e-27 + // LCOV_EXCL_START static const T P[14] = { BOOST_MATH_BIG_CONSTANT(T, 113, -20.7944102007844314586649688802236072), BOOST_MATH_BIG_CONSTANT(T, 113, -4.95759941987499442499908748130192187), @@ -862,6 +893,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant BOOST_MATH_BIG_CONSTANT(T, 113, -0.557103423021951053707162364713587374e-31), BOOST_MATH_BIG_CONSTANT(T, 113, 0.618708773442584843384712258199645166e-34), }; + // LCOV_EXCL_END result = tools::evaluate_polynomial(P, T(s - 30)) / tools::evaluate_polynomial(Q, T(s - 30)); result = 1 + exp(result); } @@ -879,9 +911,11 @@ T zeta_imp_prec(T s, T sc, const Policy&, const std::integral_constant template T zeta_imp_odd_integer(int s, const T&, const Policy&, const std::true_type&) { + // LCOV_EXCL_START static const T results[] = { BOOST_MATH_BIG_CONSTANT(T, 113, 1.2020569031595942853997381615114500), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0369277551433699263313654864570342), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0083492773819228268397975498497968), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0020083928260822144178527692324121), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0004941886041194645587022825264699), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0001227133475784891467518365263574), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000305882363070204935517285106451), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000076371976378997622736002935630), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000019082127165539389256569577951), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000004769329867878064631167196044), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000001192199259653110730677887189), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000298035035146522801860637051), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000074507117898354294919810042), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000018626597235130490064039099), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000004656629065033784072989233), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000001164155017270051977592974), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000291038504449709968692943), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000072759598350574810145209), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000018189896503070659475848), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000004547473783042154026799), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000001136868407680227849349), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000284217097688930185546), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000071054273952108527129), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000017763568435791203275), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000004440892103143813364), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000001110223025141066134), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000277555756213612417), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000069388939045441537), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000017347234760475766), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000004336808690020650), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000001084202172494241), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000271050543122347), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000067762635780452), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000016940658945098), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000004235164736273), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000001058791184068), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000264697796017), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000066174449004), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000016543612251), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000004135903063), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000001033975766), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000258493941), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000064623485), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000016155871), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000004038968), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000001009742), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000252435), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000063109), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000015777), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000003944), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000986), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000247), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000062), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000015), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000004), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000001), }; + // LCOV_EXCL_END return s > 113 ? 1 : results[(s - 3) / 2]; } @@ -891,9 +925,11 @@ T zeta_imp_odd_integer(int s, const T& sc, const Policy& pol, const std::false_t #ifdef BOOST_MATH_NO_THREAD_LOCAL_WITH_NON_TRIVIAL_TYPES static_assert(std::is_trivially_destructible::value, "Your platform does not support thread_local with non-trivial types, last checked with Mingw-x64-8.1, Jan 2021. Please try a Mingw build with the POSIX threading model, see https://sourceforge.net/p/mingw-w64/bugs/527/"); #endif + // LCOV_EXCL_START static BOOST_MATH_THREAD_LOCAL bool is_init = false; static BOOST_MATH_THREAD_LOCAL T results[50] = {}; static BOOST_MATH_THREAD_LOCAL int digits = tools::digits(); + // LCOV_EXCL_END int current_digits = tools::digits(); if(digits != current_digits) { From 24ede6fcdd73945c0d57962ca62088e99d557f72 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Wed, 3 Jul 2024 19:52:07 +0100 Subject: [PATCH 17/27] gamma.hpp coverage. --- include/boost/math/special_functions/gamma.hpp | 18 ++++++++++++------ test/test_tgamma_eh.cpp | 6 ++++++ test/test_tgamma_ratio.hpp | 16 ++++++++++++++++ 3 files changed, 34 insertions(+), 6 deletions(-) diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 523e1bc84d..3ca4e4c635 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -1436,7 +1436,7 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, optimised_invert = true; } else - init_value = 0; + init_value = 0; // LCOV_EXCL_LINE Unreachable for any "sensible" floating point type. } else init_value = 0; @@ -1546,15 +1546,18 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, } if(p_derivative) { - // - // Need to convert prefix term to derivative: - // + /* + * We can never reach this: if((x < 1) && (tools::max_value() * x < *p_derivative)) { // overflow, just return an arbitrarily large value: *p_derivative = tools::max_value() / 2; } - + */ + BOOST_MATH_ASSERT((x >= 1) || (tools::max_value() * x >= *p_derivative)); + // + // Need to convert prefix term to derivative: + // *p_derivative /= x; } @@ -1800,7 +1803,10 @@ T tgamma_ratio_imp(T x, T y, const Policy& pol) // // Result will almost certainly overflow, try logs just in case: // - return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol)); + prefix = boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol); + if (prefix > boost::math::tools::log_max_value()) + return policies::raise_overflow_error("tgamma_ratio", nullptr, pol); + return exp(prefix); } // // Regular case, x and y both large and similar in magnitude: diff --git a/test/test_tgamma_eh.cpp b/test/test_tgamma_eh.cpp index d27e7e5c1e..1ebb4e4ce8 100644 --- a/test/test_tgamma_eh.cpp +++ b/test/test_tgamma_eh.cpp @@ -3,6 +3,8 @@ // 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_MATH_PROMOTE_DOUBLE_POLICY false + #include "math_unit_test.hpp" #include #include @@ -13,5 +15,9 @@ int main() { CHECK_EQUAL(boost::math::tgamma(-200.5), 0.0); // triggers internal exception handling CHECK_EQUAL(boost::math::gamma_p(500.125, 1e-50), 0.0); // triggers internal exception handling + + // Lines that can only be hit when promotion to 80-bit reals is turned off + CHECK_ULP_CLOSE(boost::math::tgamma(44.0, 0.000001), 6.04152630633738356373551320685139975072645120000000000000000e52, 10); + CHECK_ULP_CLOSE(boost::math::gamma_p(1.0001, boost::math::gamma_p_inv(1.0001, 1e-200)), 1e-200, 10); return boost::math::test::report_errors(); } diff --git a/test/test_tgamma_ratio.hpp b/test/test_tgamma_ratio.hpp index 32a109290b..37910223a1 100644 --- a/test/test_tgamma_ratio.hpp +++ b/test/test_tgamma_ratio.hpp @@ -179,6 +179,16 @@ void test_spots(T, const char*) BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_delta_ratio(T(ldexp(T(1), -1074)), T(200)), T(5.13282785052571536804189023927976812551830809667482691717029e-50L), tol * 10); BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_delta_ratio(T(200), T(ldexp(T(1), -1074))), T(1), tol); } +#if LDBL_MAX_EXP >= 16384 + if(0 != ldexp(T(1), -16445)) + { + // This is denorm_min at 80 bit long double precision: + BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_ratio(T(ldexp(T(1), -16445)), T(200)), T(6.956968705776973348795055087648739671964943862403874276047e4577L), tol * 50); + BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_ratio(T(200), T(ldexp(T(1), -16445))), T(1.43740764446678254374798600325709882852078655978124426386e-4578L), tol * 10); + BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_delta_ratio(T(ldexp(T(1), -16445)), T(200)), T(6.956968705776973348795055087648739671964943862403874276047e4577L), tol * 10); + BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_delta_ratio(T(200), T(ldexp(T(1), -16445))), T(1), tol); + } +#endif } // // Coverage: @@ -186,4 +196,10 @@ void test_spots(T, const char*) BOOST_CHECK_EQUAL(boost::math::tgamma_delta_ratio(T(2), T(0)), T(1)); BOOST_CHECK_EQUAL(boost::math::tgamma_delta_ratio(T(200), T(0)), T(1)); BOOST_CHECK_EQUAL(boost::math::tgamma_delta_ratio(T(2000), T(0)), T(1)); + + BOOST_CHECK_EQUAL(boost::math::tgamma_ratio(T(0.5), T(100000)), T(0)); + BOOST_IF_CONSTEXPR(std::numeric_limits::has_infinity) + { + BOOST_CHECK_EQUAL(boost::math::tgamma_ratio(T(100000), T(0.5)), std::numeric_limits::infinity()); + } } From e322c3143f67e3c5d80d27d4de9245deda6d9667 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Thu, 4 Jul 2024 19:18:01 +0100 Subject: [PATCH 18/27] Hopefully complete gamma coverage. --- .../special_functions/detail/igamma_large.hpp | 9 +++++++ .../special_functions/detail/lgamma_small.hpp | 26 ++++++++++++++++++- .../boost/math/special_functions/gamma.hpp | 10 ++++++- 3 files changed, 43 insertions(+), 2 deletions(-) diff --git a/include/boost/math/special_functions/detail/igamma_large.hpp b/include/boost/math/special_functions/detail/igamma_large.hpp index 5483b53fb6..c05225efce 100644 --- a/include/boost/math/special_functions/detail/igamma_large.hpp +++ b/include/boost/math/special_functions/detail/igamma_large.hpp @@ -88,6 +88,7 @@ T igamma_temme_large(T a, T x, const Policy& pol, std::integral_constant(workspace, 1/a); @@ -293,6 +296,7 @@ T igamma_temme_large(T a, T x, const Policy& pol, std::integral_constant(-0.33333333333333333L), static_cast(0.083333333333333333L), @@ -406,6 +410,7 @@ T igamma_temme_large(T a, T x, const Policy& pol, std::integral_constant(0.00083949872067208728L), static_cast(-0.00043829709854172101L), }; + // LCOV_EXCL_END workspace[8] = tools::evaluate_polynomial(C8, z); workspace[9] = static_cast(-0.00059676129019274625L); @@ -435,6 +440,7 @@ T igamma_temme_large(T a, T x, const Policy& pol, std::integral_constant(-0.333333333L), static_cast(0.0833333333L), @@ -461,6 +467,7 @@ T igamma_temme_large(T a, T x, const Policy& pol, std::integral_constant(0.000771604938L), }; workspace[2] = tools::evaluate_polynomial(C2, z); + // LCOV_EXCL_END T result = tools::evaluate_polynomial(workspace, 1/a); result *= exp(-y) / sqrt(2 * constants::pi() * a); @@ -478,6 +485,7 @@ T igamma_temme_large(T a, T x, const Policy& pol, std::integral_constant T igamma_temme_large(T a, T x, const Policy& pol, std::integral_constant const *) { @@ -768,6 +776,7 @@ T igamma_temme_large(T a, T x, const Policy& pol, std::integral_constant&, co BOOST_MATH_STD_USING // for ADL of std names T result = 0; + + BOOST_ASSERT(z >= tools::root_epsilon()); + /* + * Can not be reached: + * if(z < tools::epsilon()) { result = -log(z); } + */ else if((zm1 == 0) || (zm2 == 0)) { // nothing to do, result is zero.... @@ -87,6 +93,7 @@ T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, co // At long double: Max error found: 1.987e-21 // Maximum Deviation Found (approximation error): 5.900e-24 // + // LCOV_EXCL_START static const T P[] = { static_cast(BOOST_MATH_BIG_CONSTANT(T, 64, -0.180355685678449379109e-1)), static_cast(BOOST_MATH_BIG_CONSTANT(T, 64, 0.25126649619989678683e-1)), @@ -108,6 +115,7 @@ T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, co }; static const float Y = 0.158963680267333984375e0f; + // LCOV_EXCL_END T r = zm2 * (z + 1); T R = tools::evaluate_polynomial(P, zm2); @@ -151,7 +159,7 @@ T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, co // Maximum Deviation Found: 3.139e-021 // Expected Error Term: 3.139e-021 - // + // LCOV_EXCL_START static const float Y = 0.52815341949462890625f; static const T P[] = { @@ -172,6 +180,7 @@ T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, co static_cast(BOOST_MATH_BIG_CONSTANT(T, 64, 0.577039722690451849648e-1)), static_cast(BOOST_MATH_BIG_CONSTANT(T, 64, 0.195768102601107189171e-2)) }; + // LCOV_EXCL_END T r = tools::evaluate_polynomial(P, zm1) / tools::evaluate_polynomial(Q, zm1); T prefix = zm1 * zm2; @@ -197,6 +206,7 @@ T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, co // Maximum Deviation Found: 2.151e-021 // Expected Error Term: 2.150e-021 // + // LCOV_EXCL_START static const float Y = 0.452017307281494140625f; static const T P[] = { @@ -216,6 +226,8 @@ T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, co static_cast(BOOST_MATH_BIG_CONSTANT(T, 64, -0.100666795539143372762e-2)), static_cast(BOOST_MATH_BIG_CONSTANT(T, 64, -0.827193521891290553639e-6)) }; + // LCOV_EXCL_END + T r = zm2 * zm1; T R = tools::evaluate_polynomial(P, T(-zm2)) / tools::evaluate_polynomial(Q, T(-zm2)); @@ -224,6 +236,11 @@ T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, co } return result; } + +// +// 128-bit floats aren't directly tested in our coverage tests (takes too long) +// LCOV_EXCL_START + template T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, const Policy& /* l */, const Lanczos&) { @@ -481,6 +498,8 @@ T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, c BOOST_MATH_INSTRUMENT_CODE(result); return result; } +// LCOV_EXCL_END + template T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, const Policy& pol, const Lanczos& l) { @@ -493,10 +512,15 @@ T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, con // BOOST_MATH_STD_USING // for ADL of std names T result = 0; + + BOOST_MATH_ASSERT(z >= tools::root_epsilon()); + /* + * Not reachable: if(z < tools::epsilon()) { result = -log(z); } + */ else if(z < 0.5) { // taking the log of tgamma reduces the error, no danger of overflow here: diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 3ca4e4c635..02100f32f7 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -1806,7 +1806,11 @@ T tgamma_ratio_imp(T x, T y, const Policy& pol) prefix = boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol); if (prefix > boost::math::tools::log_max_value()) return policies::raise_overflow_error("tgamma_ratio", nullptr, pol); - return exp(prefix); + // + // This is unreachable, unless max_factorial is small compared to the exponent + // range of the type, ie multiprecision types only here... + // + return exp(prefix); // LCOV_EXCL_LINE } // // Regular case, x and y both large and similar in magnitude: @@ -1838,11 +1842,15 @@ T gamma_p_derivative_imp(T a, T x, const Policy& pol) // typedef typename lanczos::lanczos::type lanczos_type; T f1 = detail::regularised_gamma_prefix(a, x, pol, lanczos_type()); + /* + * Derivative goes to zero as x -> 0, this should be unreachable: + * if((x < 1) && (tools::max_value() * x < f1)) { // overflow: return policies::raise_overflow_error("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", nullptr, pol); } + */ if(f1 == 0) { // Underflow in calculation, use logs instead: From 67dc6ec33502ac5ab47f108b048637bb5dc1058d Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Thu, 4 Jul 2024 19:22:54 +0100 Subject: [PATCH 19/27] Correct standalone failure. --- include/boost/math/special_functions/detail/lgamma_small.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/math/special_functions/detail/lgamma_small.hpp b/include/boost/math/special_functions/detail/lgamma_small.hpp index 608327b550..02c6835aba 100644 --- a/include/boost/math/special_functions/detail/lgamma_small.hpp +++ b/include/boost/math/special_functions/detail/lgamma_small.hpp @@ -46,7 +46,7 @@ T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, co BOOST_MATH_STD_USING // for ADL of std names T result = 0; - BOOST_ASSERT(z >= tools::root_epsilon()); + BOOST_MATH_ASSERT(z >= tools::root_epsilon()); /* * Can not be reached: * From 1c6c602c6708a2a4638c5f19091575dae1655db7 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Thu, 4 Jul 2024 19:27:03 +0100 Subject: [PATCH 20/27] Corrections for gamma coverage. --- .../math/special_functions/detail/lgamma_small.hpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/include/boost/math/special_functions/detail/lgamma_small.hpp b/include/boost/math/special_functions/detail/lgamma_small.hpp index 02c6835aba..0c2a1ad505 100644 --- a/include/boost/math/special_functions/detail/lgamma_small.hpp +++ b/include/boost/math/special_functions/detail/lgamma_small.hpp @@ -55,7 +55,7 @@ T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, co result = -log(z); } */ - else if((zm1 == 0) || (zm2 == 0)) + if((zm1 == 0) || (zm2 == 0)) { // nothing to do, result is zero.... } @@ -251,12 +251,16 @@ T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, c // BOOST_MATH_STD_USING // for ADL of std names T result = 0; + BOOST_MATH_ASSERT(z >= tools::root_epsilon()); + /* + * Can not be reached: if(z < tools::epsilon()) { result = -log(z); BOOST_MATH_INSTRUMENT_CODE(result); } - else if((zm1 == 0) || (zm2 == 0)) + */ + if((zm1 == 0) || (zm2 == 0)) { // nothing to do, result is zero.... } @@ -521,7 +525,7 @@ T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, con result = -log(z); } */ - else if(z < 0.5) + if(z < 0.5) { // taking the log of tgamma reduces the error, no danger of overflow here: result = log(gamma_imp(z, pol, Lanczos())); From 0848e3bb770d5b95c3a950343d3897db9da9cd17 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 5 Jul 2024 10:01:41 +0100 Subject: [PATCH 21/27] Exclude owen's T tables from coverage. --- include/boost/math/special_functions/owens_t.hpp | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/include/boost/math/special_functions/owens_t.hpp b/include/boost/math/special_functions/owens_t.hpp index e09608b397..f53effc309 100644 --- a/include/boost/math/special_functions/owens_t.hpp +++ b/include/boost/math/special_functions/owens_t.hpp @@ -70,6 +70,7 @@ namespace boost template inline unsigned short owens_t_compute_code(const RealType h, const RealType a) { + // LCOV_EXCL_START static const RealType hrange[] = { 0.02f, 0.06f, 0.09f, 0.125f, 0.26f, 0.4f, 0.6f, 1.6f, 1.7f, 2.33f, 2.4f, 3.36f, 3.4f, 4.8f }; @@ -97,6 +98,7 @@ namespace boost 1 , 2 , 3 , 3 , 5 , 5 , 7 , 7 , 16 , 16 , 16 , 16 , 16 , 11 , 11, 1 , 2 , 3 , 3 , 5 , 5 , 17 , 17 , 17 , 17 , 16 , 16 , 16 , 11 , 11 }; + // LCOV_EXCL_END unsigned short ihint = 14, iaint = 7; for(unsigned short i = 0; i != 14; i++) @@ -125,7 +127,9 @@ namespace boost template inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const std::integral_constant&) { + // LCOV_EXCL_START static const unsigned short ord[] = {2, 3, 4, 5, 7, 10, 12, 18, 10, 20, 30, 0, 4, 7, 8, 20, 0, 0}; // 18 entries + // LCOV_EXCL_END BOOST_MATH_ASSERT(icode<18); @@ -136,7 +140,9 @@ namespace boost inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const std::integral_constant&) { // method ================>>> {1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 4, 4, 5, 6} - static const unsigned short ord[] = {3, 4, 5, 6, 8, 11, 13, 19, 10, 20, 30, 0, 7, 10, 11, 23, 0, 0}; // 18 entries + // LCOV_EXCL_START + static const unsigned short ord[] = {3, 4, 5, 6, 8, 11, 13, 19, 10, 20, 30, 0, 7, 10, 11, 23, 0, 0}; // 18 entries + // LCOV_EXCL_END BOOST_MATH_ASSERT(icode<18); @@ -233,6 +239,7 @@ namespace boost const unsigned short m = 20; + // LCOV_EXCL_START static const RealType c2[] = { static_cast(0.99999999999999987510), @@ -247,6 +254,7 @@ namespace boost static_cast(-0.82813631607004984866E-01), static_cast(0.24167984735759576523E-01), static_cast(-0.44676566663971825242E-02), static_cast(0.39141169402373836468E-03) }; + // LCOV_EXCL_END const RealType as = a*a; const RealType hs = h*h; @@ -285,6 +293,7 @@ namespace boost const unsigned short m = 30; + // LCOV_EXCL_START static const RealType c2[] = { BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999999999999999729978162447266851932041876728736094298092917625009873), @@ -319,6 +328,7 @@ namespace boost BOOST_MATH_BIG_CONSTANT(RealType, 260, -1.489155613350368934073453260689881330166342484405529981510694514036264969925132e-4), BOOST_MATH_BIG_CONSTANT(RealType, 260, 9.072354320794357587710929507988814669454281514268844884841547607134260303118208e-6) }; + // LCOV_EXCL_END const RealType as = a*a; const RealType hs = h*h; @@ -404,6 +414,7 @@ namespace boost */ const unsigned short m = 13; + // LCOV_EXCL_START static const RealType pts[] = { static_cast(0.35082039676451715489E-02), static_cast(0.31279042338030753740E-01), static_cast(0.85266826283219451090E-01), @@ -423,6 +434,7 @@ namespace boost const RealType as = a*a; const RealType hs = -h*h*boost::math::constants::half(); + // LCOV_EXCL_END RealType val = 0; for(unsigned short i = 0; i < m; ++i) @@ -450,6 +462,7 @@ namespace boost */ const unsigned short m = 19; + // LCOV_EXCL_START static const RealType pts[] = { BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0016634282895983227941), BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.014904509242697054183), @@ -492,6 +505,7 @@ namespace boost BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0018483371329504443947), BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.00079623320100438873578) }; + // LCOV_EXCL_END const RealType as = a*a; const RealType hs = -h*h*boost::math::constants::half(); From fc0b406af4ed142acb51c84c0d937a399f9e77d6 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 5 Jul 2024 10:19:36 +0100 Subject: [PATCH 22/27] Hankel error handling tests. --- test/test_hankel.cpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/test/test_hankel.cpp b/test/test_hankel.cpp index f8bd173da8..1d537e1b0a 100644 --- a/test/test_hankel.cpp +++ b/test/test_hankel.cpp @@ -80,6 +80,18 @@ void test_hankel(T, const char* name) boost::math::sph_hankel_2(data[i][0].real(), data[i][1].real())); } } + + BOOST_IF_CONSTEXPR(std::numeric_limits::has_infinity) + { + std::complex r = boost::math::cyl_hankel_1(T(1), T(0)); + BOOST_CHECK_EQUAL(r, std::complex(std::numeric_limits::infinity(), std::numeric_limits::infinity())); + r = boost::math::cyl_hankel_2(T(1), T(0)); + BOOST_CHECK_EQUAL(r, std::complex(std::numeric_limits::infinity(), -std::numeric_limits::infinity())); + r = boost::math::cyl_hankel_1(T(0), T(0)); + BOOST_CHECK_EQUAL(r, std::complex(1, -std::numeric_limits::infinity())); + r = boost::math::cyl_hankel_2(T(0), T(0)); + BOOST_CHECK_EQUAL(r, std::complex(1, std::numeric_limits::infinity())); + } } // From 6440f731b18326663dba97ea6bfd8a4fe99e4dfb Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 5 Jul 2024 10:26:55 +0100 Subject: [PATCH 23/27] Heumann Lambda coverage. --- test/test_heuman_lambda.hpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/test/test_heuman_lambda.hpp b/test/test_heuman_lambda.hpp index 23720b2d02..cdcbbba51b 100644 --- a/test/test_heuman_lambda.hpp +++ b/test/test_heuman_lambda.hpp @@ -71,5 +71,12 @@ void test_spots(T, const char* type_name) #include "heuman_lambda_data.ipp" do_test_heuman_lambda(heuman_lambda_data, type_name, "Elliptic Integral Heuman Lambda: Random Data"); + + // + // Special cases for coverage: + // + BOOST_CHECK_THROW(boost::math::heuman_lambda(T(1.1), T(0.5)), std::domain_error); + BOOST_CHECK_THROW(boost::math::heuman_lambda(static_cast(1e-100), T(2.5)), std::domain_error); + } From ba6838ce073547f5d6938299ab9d340a8ffc41ba Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sat, 13 Jul 2024 20:04:46 +0100 Subject: [PATCH 24/27] Improve 0F1, 1F0 and 1F1 coverage. --- .../special_functions/hypergeometric_0F1.hpp | 4 +- .../special_functions/hypergeometric_1F0.hpp | 9 +- .../special_functions/hypergeometric_1F1.hpp | 38 ++- test/test_0F1.hpp | 6 + test/test_1F1.hpp | 218 ++++++++++-------- test/test_1F1_log.hpp | 14 ++ test/test_1F1_regularized.cpp | 5 - test/test_1F1_regularized.hpp | 16 ++ test/test_2F0.hpp | 5 + 9 files changed, 188 insertions(+), 127 deletions(-) diff --git a/include/boost/math/special_functions/hypergeometric_0F1.hpp b/include/boost/math/special_functions/hypergeometric_0F1.hpp index 5f0dedd115..a721e7e462 100644 --- a/include/boost/math/special_functions/hypergeometric_0F1.hpp +++ b/include/boost/math/special_functions/hypergeometric_0F1.hpp @@ -64,9 +64,7 @@ namespace boost { namespace math { namespace detail { return T(1); if ((b <= 0) && (b == floor(b))) - return policies::raise_pole_error( - function, - "Evaluation of 0f1 with nonpositive integer b = %1%.", b, pol); + return policies::raise_pole_error(function, "Evaluation of 0f1 with nonpositive integer b = %1%.", b, pol); if (z < -5 && b > -5) { diff --git a/include/boost/math/special_functions/hypergeometric_1F0.hpp b/include/boost/math/special_functions/hypergeometric_1F0.hpp index c78903af18..3fbf81ba8d 100644 --- a/include/boost/math/special_functions/hypergeometric_1F0.hpp +++ b/include/boost/math/special_functions/hypergeometric_1F0.hpp @@ -25,16 +25,11 @@ inline T hypergeometric_1F0_imp(const T& a, const T& z, const Policy& pol) BOOST_MATH_STD_USING // pow if (z == 1) - return policies::raise_pole_error( - function, - "Evaluation of 1F0 with z = %1%.", - z, - pol); + return policies::raise_pole_error(function, "Evaluation of 1F0 with z = %1%.", z, pol); if (1 - z < 0) { if (floor(a) != a) - return policies::raise_domain_error(function, - "Result is complex when a is non-integral and z > 1, but got z = %1%", z, pol); + return policies::raise_domain_error(function, "Result is complex when a is non-integral and z > 1, but got z = %1%", z, pol); } // more naive and convergent method than series return pow(T(1 - z), T(-a)); diff --git a/include/boost/math/special_functions/hypergeometric_1F1.hpp b/include/boost/math/special_functions/hypergeometric_1F1.hpp index fd53cddbab..35d087ed6f 100644 --- a/include/boost/math/special_functions/hypergeometric_1F1.hpp +++ b/include/boost/math/special_functions/hypergeometric_1F1.hpp @@ -80,7 +80,8 @@ namespace boost { namespace math { namespace detail { if (a < 0) { if ((b < a) && (z < -b / 4)) - return hypergeometric_1F1_from_function_ratio_negative_ab(a, b, z, pol, log_scaling); + // Defensive programming: it is *almost* certain that we can never get here, proving that is hard though... + return hypergeometric_1F1_from_function_ratio_negative_ab(a, b, z, pol, log_scaling); // LCOV_EXCL_LINE else { // @@ -123,6 +124,7 @@ namespace boost { namespace math { namespace detail { { if (boost::math::policies::digits() <= 64) return hypergeometric_1F1_from_function_ratio_negative_b_forwards(a, b, z, pol, log_scaling); + // LCOV_EXCL_START, what follows is multiprecision only #ifndef BOOST_MATH_NO_EXCEPTIONS try #endif @@ -138,6 +140,7 @@ namespace boost { namespace math { namespace detail { return hypergeometric_1F1_from_function_ratio_negative_b_forwards(a, b, z, pol, log_scaling); } #endif + // LCOV_EXCL_END } // // We could fall back to Tricomi's approximation if we're in the transition zone @@ -160,6 +163,8 @@ namespace boost { namespace math { namespace detail { return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling); } +#if 0 + // Archived, not used, see comments at call site. template bool is_convergent_negative_z_series(const T& a, const T& b, const T& z, const T& b_minus_a) { @@ -180,7 +185,7 @@ namespace boost { namespace math { namespace detail { // Double check for divergence when we cross the origin on a and b: if (a < 0) { - T n = 300 - floor(a); + T n = 3 - floor(a); if (fabs((a + n) * z / ((b + n) * n)) < 1) { if (b < 0) @@ -225,7 +230,7 @@ namespace boost { namespace math { namespace detail { } return false; } - +#endif template inline T cyl_bessel_i_shrinkage_rate(const T& z) { @@ -313,11 +318,7 @@ namespace boost { namespace math { namespace detail { // undefined result: if (!detail::check_hypergeometric_1F1_parameters(a, b)) - return policies::raise_domain_error( - function, - "Function is indeterminate for negative integer b = %1%.", - b, - pol); + return policies::raise_domain_error(function, "Function is indeterminate for negative integer b = %1%.", b, pol); // other checks: if (a == -1) @@ -444,6 +445,12 @@ namespace boost { namespace math { namespace detail { { if (a == 1) return detail::hypergeometric_1F1_pade(b, z, pol); +#if 0 + // + // Commented out: is_convergent_negative_z_series is fine so far as it goes + // but there appear to be no cases that use it, and in extremis, we will + // fall through to the series evaluation anyway. + // if (is_convergent_negative_z_series(a, b, z, b_minus_a)) { if ((boost::math::sign(b_minus_a) == boost::math::sign(b)) && ((b > 0) || (b < -200))) @@ -462,6 +469,7 @@ namespace boost { namespace math { namespace detail { return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling); } } +#endif if ((b < 0) && (floor(b) == b)) { // Negative integer b, so a must be a negative integer too. @@ -673,16 +681,26 @@ namespace boost { namespace math { namespace detail { while (scale > max_scaling) { + if((fabs(result) > 1) && (fabs(tools::max_value()) / result <= max_scale_factor)) + return policies::raise_overflow_error("hypergeometric_1F1_regularized", nullptr, pol); + // This is *probably* unreachable: + // LCOV_EXCL_START result *= max_scale_factor; scale -= max_scaling; + // LCOV_EXCL_STOP } while (scale < -max_scaling) { result /= max_scale_factor; - scale += max_scaling; + scale += max_scaling; } if (scale != 0) - result *= exp(scale); + { + scale = exp(scale); + if ((scale > 1) && (fabs(result) > 1) && (fabs(tools::max_value() / result) <= scale)) + return policies::raise_overflow_error("hypergeometric_1F1_regularized", nullptr, pol); + result *= scale; + } return result * result_sign; } diff --git a/test/test_0F1.hpp b/test/test_0F1.hpp index 020ae20201..05bd4d8d25 100644 --- a/test/test_0F1.hpp +++ b/test/test_0F1.hpp @@ -215,5 +215,11 @@ void test_spots(T, const char* type_name) do_test_0F1(hypergeometric_0F1_integer_data, type_name, "Integer Arguments"); do_test_0F1(hypergeometric_0F1_real_data, type_name, "Real Arguments"); do_test_0F1(hypergeometric_0F1_large_data, type_name, "Large Arguments"); + + // + // Coverage: + // + BOOST_CHECK_THROW(boost::math::hypergeometric_0F1(static_cast(-1), static_cast(0.5)), std::domain_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_0F1(static_cast(-20), static_cast(0.5)), std::domain_error); } diff --git a/test/test_1F1.hpp b/test/test_1F1.hpp index e4ab8151cb..1d4b9b307e 100644 --- a/test/test_1F1.hpp +++ b/test/test_1F1.hpp @@ -162,94 +162,94 @@ void test_spots5(T, const char* type_name) template void test_spots6(T, const char* type_name) { - static const std::array, 186> hypergeometric_1F1_bugs = { { - { { static_cast(17955.561660766602), static_cast(9.6968994205831605e-09), static_cast(-82.406154185533524), SC_(6.98056008378736714088730927132364938220428678e-11) }}, - { { static_cast(17955.561660766602), static_cast(-9.6968994205831605e-09), static_cast(-82.406154185533524), SC_(-6.98055306629610746072607353939306734740549551e-11) }}, - { { static_cast(-17955.561660766602), static_cast(-9.6968994205831605e-09), static_cast(82.406154185533524), SC_(-42897094853118832762870100.8669248353530950866) }} , - { { static_cast(17955.561660766602), static_cast(17956.061660766602), static_cast(82.406154185533524), SC_(613117565438499794408370861624072730.553215432) }}, - { { static_cast(2.9127331452327709e-07), static_cast(-0.99999970872668542), static_cast(0.15018942760070786), SC_(0.987526018990506843793601092932108059727149508) }}, - { { static_cast(-2.9127331452327709e-07), static_cast(-1.0000002912733146), static_cast(0.15018942760070786), SC_(0.987526120661366412484942089372497015837368389) }}, - { { static_cast(6.7191087900739423e-13), static_cast(-0.99999999999932809), static_cast(0.0011913633891253994), SC_(0.999999289758605006762757201699750974296453229) }}, - { { static_cast(6.7191087900739423e-13), static_cast(-0.99999999999932809), static_cast(-0.0011913633891253994), SC_(0.999999290885918468326416221021126912154021802) }}, - { { static_cast(-6.7191087900739423e-13), static_cast(-1.0000000000006719), static_cast(0.0011913633891253994), SC_(0.999999289758606609651292394510404091049823243) }}, - { { static_cast(-6.7191087900739423e-13), static_cast(-1.0000000000006719), static_cast(-0.0011913633891253994), SC_(0.999999290885916869252591036674587894145399498) }}, - { { static_cast(1.2860067365774887e-17), static_cast(6.2442285664031425e-16), static_cast(-2539.60133934021), SC_(0.979404874070484696999110600576068012417904384) }}, - { { static_cast(1.2860067365774887e-17), static_cast(-6.2442285664031425e-16), static_cast(-2539.60133934021), SC_(1.0205951259295150865252112924093487321207727) }}, - { { static_cast(-1.2860067365774887e-17), static_cast(6.2442285664031425e-16), static_cast(-2539.60133934021), SC_(1.02059512592951530745923325071510441026202975) }}, - { { static_cast(-1.2860067365774887e-17), static_cast(-6.2442285664031425e-16), static_cast(-2539.60133934021), SC_(0.979404874070484909016444856299500644331897735) }}, - { { static_cast(1.2860067365774887e-17), static_cast(1), static_cast(-2539.60133934021), SC_(0.999999999999999891757095137551552220860540801) }}, - { { static_cast(-1.2860067365774887e-17), static_cast(1), static_cast(-2539.60133934021), SC_(1.00000000000000010824290486244845922375479178) }}, - { { static_cast(1.2860067365774887e-17), static_cast(0.5), static_cast(-2539.60133934021), SC_(0.999999999999999873931788919689096760455570214) }}, - { { static_cast(-1.2860067365774887e-17), static_cast(0.5), static_cast(-2539.60133934021), SC_(1.0000000000000001260682110803109183167444166) }}, - { { static_cast(1.2860067365774887e-17), static_cast(-0.5), static_cast(-2539.60133934021), SC_(0.999999999999999899656990458526368219886894767) }}, - { { static_cast(-1.2860067365774887e-17), static_cast(-0.5), static_cast(-2539.60133934021), SC_(1.00000000000000010034300954147364037131355735) }}, - { { static_cast(1.9561377367172441e-13), static_cast(-0.99999999999980438), static_cast(0.53720525559037924), SC_(0.791950585963666119273677451162365759080483409) }}, - { { static_cast(1.9561377367172441e-13), static_cast(-0.99999999999980438), static_cast(-0.53720525559037924), SC_(0.898314630992769591673208399706587643905527327) }}, - { { static_cast(-1.9561377367172441e-13), static_cast(-1.0000000000001956), static_cast(0.53720525559037924), SC_(0.791950585964025761367113514279915403442035074) }}, - { { static_cast(-1.9561377367172441e-13), static_cast(-1.0000000000001956), static_cast(-0.53720525559037924), SC_(0.898314630992646771749564140770704893561753597) }}, - { { static_cast(5.1851756946064858e-12), static_cast(-0.99999999999481481), static_cast(-774.06985878944397), SC_(1.91306610467163858324476828831735612399803649e-06) }}, - { { static_cast(-5.1851756946064858e-12), static_cast(-1.0000000000051852), static_cast(-774.06985878944397), SC_(1.91306610479516297551035931150910859922270467e-06) }}, - - {{ static_cast(4.782769898853794e-15), static_cast(1.0000000000000049), static_cast(43.289540141820908), SC_(715.678254892476818206948251991084031658534788) }}, - { { static_cast(-4.782769898853794e-15), static_cast(0.99999999999999523), static_cast(43.289540141820908), SC_(-713.67825489247727251051792450091274703212426) }}, - { { static_cast(4.782769898853794e-15), static_cast(0.50000000000000477), static_cast(43.289540141820908), SC_(8235.578376364917373771471380274179857713986) }}, - { { static_cast(-4.782769898853794e-15), static_cast(0.49999999999999523), static_cast(43.289540141820908), SC_(-8233.57837636502669085205930058992320862281194) }}, - { { static_cast(4.782769898853794e-15), static_cast(-0.49999999999999523), static_cast(43.289540141820908), SC_(-696269.800378137841948029488304613132151506346) }}, - { { static_cast(-4.782769898853794e-15), static_cast(-0.50000000000000477), static_cast(43.289540141820908), SC_(696271.8003781336298001417038674968935893361) }}, - { { static_cast(8.1104991963343309e-05), static_cast(-0.99991889500803666), static_cast(-289.12455415725708), SC_(7.89625448009377635153307897651433007437615965e-124) }}, - { { static_cast(-8.1104991963343309e-05), static_cast(-1.0000811049919633), static_cast(-289.12455415725708), SC_(7.8949781467741574268884621364833028722017032e-124) }}, - - {{ static_cast(-1.98018241448205767), static_cast(1.98450573845762079), static_cast(54.4977916804564302), SC_(2972026581564772.790187123046255523239732028) }}, - - { { static_cast(-7.8238229420435346e-05), static_cast(-0.50007823822942044), static_cast(-1896.0561106204987), SC_(1.00058778866237037053151236215058095904086972) }}, + static const std::array, 188> hypergeometric_1F1_bugs = { { + { { static_cast(17955.561660766602), static_cast(9.6968994205831605e-09), static_cast(-82.406154185533524), SC_(6.98056008378736714088730927132364938220428678e-11) }}, + { { static_cast(17955.561660766602), static_cast(-9.6968994205831605e-09), static_cast(-82.406154185533524), SC_(-6.98055306629610746072607353939306734740549551e-11) }}, + { { static_cast(-17955.561660766602), static_cast(-9.6968994205831605e-09), static_cast(82.406154185533524), SC_(-42897094853118832762870100.8669248353530950866) }} , + { { static_cast(17955.561660766602), static_cast(17956.061660766602), static_cast(82.406154185533524), SC_(613117565438499794408370861624072730.553215432) }}, + { { static_cast(2.9127331452327709e-07), static_cast(-0.99999970872668542), static_cast(0.15018942760070786), SC_(0.987526018990506843793601092932108059727149508) }}, + { { static_cast(-2.9127331452327709e-07), static_cast(-1.0000002912733146), static_cast(0.15018942760070786), SC_(0.987526120661366412484942089372497015837368389) }}, + { { static_cast(6.7191087900739423e-13), static_cast(-0.99999999999932809), static_cast(0.0011913633891253994), SC_(0.999999289758605006762757201699750974296453229) }}, + { { static_cast(6.7191087900739423e-13), static_cast(-0.99999999999932809), static_cast(-0.0011913633891253994), SC_(0.999999290885918468326416221021126912154021802) }}, + { { static_cast(-6.7191087900739423e-13), static_cast(-1.0000000000006719), static_cast(0.0011913633891253994), SC_(0.999999289758606609651292394510404091049823243) }}, + { { static_cast(-6.7191087900739423e-13), static_cast(-1.0000000000006719), static_cast(-0.0011913633891253994), SC_(0.999999290885916869252591036674587894145399498) }}, + { { static_cast(1.2860067365774887e-17), static_cast(6.2442285664031425e-16), static_cast(-2539.60133934021), SC_(0.979404874070484696999110600576068012417904384) }}, + { { static_cast(1.2860067365774887e-17), static_cast(-6.2442285664031425e-16), static_cast(-2539.60133934021), SC_(1.0205951259295150865252112924093487321207727) }}, + { { static_cast(-1.2860067365774887e-17), static_cast(6.2442285664031425e-16), static_cast(-2539.60133934021), SC_(1.02059512592951530745923325071510441026202975) }}, + { { static_cast(-1.2860067365774887e-17), static_cast(-6.2442285664031425e-16), static_cast(-2539.60133934021), SC_(0.979404874070484909016444856299500644331897735) }}, + { { static_cast(1.2860067365774887e-17), static_cast(1), static_cast(-2539.60133934021), SC_(0.999999999999999891757095137551552220860540801) }}, + { { static_cast(-1.2860067365774887e-17), static_cast(1), static_cast(-2539.60133934021), SC_(1.00000000000000010824290486244845922375479178) }}, + { { static_cast(1.2860067365774887e-17), static_cast(0.5), static_cast(-2539.60133934021), SC_(0.999999999999999873931788919689096760455570214) }}, + { { static_cast(-1.2860067365774887e-17), static_cast(0.5), static_cast(-2539.60133934021), SC_(1.0000000000000001260682110803109183167444166) }}, + { { static_cast(1.2860067365774887e-17), static_cast(-0.5), static_cast(-2539.60133934021), SC_(0.999999999999999899656990458526368219886894767) }}, + { { static_cast(-1.2860067365774887e-17), static_cast(-0.5), static_cast(-2539.60133934021), SC_(1.00000000000000010034300954147364037131355735) }}, + { { static_cast(1.9561377367172441e-13), static_cast(-0.99999999999980438), static_cast(0.53720525559037924), SC_(0.791950585963666119273677451162365759080483409) }}, + { { static_cast(1.9561377367172441e-13), static_cast(-0.99999999999980438), static_cast(-0.53720525559037924), SC_(0.898314630992769591673208399706587643905527327) }}, + { { static_cast(-1.9561377367172441e-13), static_cast(-1.0000000000001956), static_cast(0.53720525559037924), SC_(0.791950585964025761367113514279915403442035074) }}, + { { static_cast(-1.9561377367172441e-13), static_cast(-1.0000000000001956), static_cast(-0.53720525559037924), SC_(0.898314630992646771749564140770704893561753597) }}, + { { static_cast(5.1851756946064858e-12), static_cast(-0.99999999999481481), static_cast(-774.06985878944397), SC_(1.91306610467163858324476828831735612399803649e-06) }}, + { { static_cast(-5.1851756946064858e-12), static_cast(-1.0000000000051852), static_cast(-774.06985878944397), SC_(1.91306610479516297551035931150910859922270467e-06) }}, + + {{ static_cast(4.782769898853794e-15), static_cast(1.0000000000000049), static_cast(43.289540141820908), SC_(715.678254892476818206948251991084031658534788) }}, + { { static_cast(-4.782769898853794e-15), static_cast(0.99999999999999523), static_cast(43.289540141820908), SC_(-713.67825489247727251051792450091274703212426) }}, + { { static_cast(4.782769898853794e-15), static_cast(0.50000000000000477), static_cast(43.289540141820908), SC_(8235.578376364917373771471380274179857713986) }}, + { { static_cast(-4.782769898853794e-15), static_cast(0.49999999999999523), static_cast(43.289540141820908), SC_(-8233.57837636502669085205930058992320862281194) }}, + { { static_cast(4.782769898853794e-15), static_cast(-0.49999999999999523), static_cast(43.289540141820908), SC_(-696269.800378137841948029488304613132151506346) }}, + { { static_cast(-4.782769898853794e-15), static_cast(-0.50000000000000477), static_cast(43.289540141820908), SC_(696271.8003781336298001417038674968935893361) }}, + { { static_cast(8.1104991963343309e-05), static_cast(-0.99991889500803666), static_cast(-289.12455415725708), SC_(7.89625448009377635153307897651433007437615965e-124) }}, + { { static_cast(-8.1104991963343309e-05), static_cast(-1.0000811049919633), static_cast(-289.12455415725708), SC_(7.8949781467741574268884621364833028722017032e-124) }}, + + {{ static_cast(-1.98018241448205767), static_cast(1.98450573845762079), static_cast(54.4977916804564302), SC_(2972026581564772.790187123046255523239732028) }}, + + { { static_cast(-7.8238229420435346e-05), static_cast(-0.50007823822942044), static_cast(-1896.0561106204987), SC_(1.00058778866237037053151236215058095904086972) }}, // Unexpected high error : 2.48268e+91 Found : -9.61305e+268 Expected : -1.74382e+193 - { { static_cast(5.9981750131794866e-15), static_cast(-230.70702263712883), static_cast(240.42092034220695), SC_(-1.74381782591884817018404492963109914357365958e+193) }}, + { { static_cast(5.9981750131794866e-15), static_cast(-230.70702263712883), static_cast(240.42092034220695), SC_(-1.74381782591884817018404492963109914357365958e+193) }}, // Unexpected high error : 1.79769313486231570814527423731704356798070568e+308 Found : -9.61305077326281580724540507004198316499661687e+268 Expected : 1.74381782591870724567837900957146707932623893e+193 - { { static_cast(-5.9981750131794866e-15), static_cast(-230.70702263712883), static_cast(240.42092034220695), SC_(1.74381782591870734495565763481520223752372107e+193) }}, - //Unexpected exception : Error in function boost::math::hypergeometric_pFq : Cancellation is so severe that no bits in the result are correct, last result was - 13497.312248525042 - { { static_cast(-0.00023636552788275367), static_cast(0.49976363447211725), static_cast(-55.448519088327885), SC_(1.00141219419064760011631555641142295011268795) }}, - // Unexpected exception: Error in function boost::math::hypergeometric_pFq: Cancellation is so severe that no bits in the result are correct, last result was -13497.312248525042 - {{ static_cast(-0.00023636552788275367), static_cast(-0.50023636552788275), static_cast(-55.448519088327885), SC_(1.00093463146763986302362749764017215184711625) }}, - // Unexpected exception : Error in function boost::math::hypergeometric_pFq : Cancellation is so severe that no bits in the result are correct, last result was - 1.3871133003351527e+47 - { { static_cast(-1.6548533913638905e-10), static_cast(0.49999999983451465), static_cast(-169.20843148231506), SC_(1.00000000117356793527360151094991866549128017) }}, - // Unexpected exception: Error in function boost::math::hypergeometric_pFq: Cancellation is so severe that no bits in the result are correct, last result was -1.3871133003351527e+47 - {{ static_cast(-1.6548533913638905e-10), static_cast(-0.50000000016548529), static_cast(-169.20843148231506), SC_(1.00000000084161045914716192484600809610013447) }}, + { { static_cast(-5.9981750131794866e-15), static_cast(-230.70702263712883), static_cast(240.42092034220695), SC_(1.74381782591870734495565763481520223752372107e+193) }}, + //Unexpected exception : Error in function boost::math::hypergeometric_pFq : Cancellation is so severe that no bits in the result are correct, last result was - 13497.312248525042 + { { static_cast(-0.00023636552788275367), static_cast(0.49976363447211725), static_cast(-55.448519088327885), SC_(1.00141219419064760011631555641142295011268795) }}, + // Unexpected exception: Error in function boost::math::hypergeometric_pFq: Cancellation is so severe that no bits in the result are correct, last result was -13497.312248525042 + {{ static_cast(-0.00023636552788275367), static_cast(-0.50023636552788275), static_cast(-55.448519088327885), SC_(1.00093463146763986302362749764017215184711625) }}, + // Unexpected exception : Error in function boost::math::hypergeometric_pFq : Cancellation is so severe that no bits in the result are correct, last result was - 1.3871133003351527e+47 + { { static_cast(-1.6548533913638905e-10), static_cast(0.49999999983451465), static_cast(-169.20843148231506), SC_(1.00000000117356793527360151094991866549128017) }}, + // Unexpected exception: Error in function boost::math::hypergeometric_pFq: Cancellation is so severe that no bits in the result are correct, last result was -1.3871133003351527e+47 + {{ static_cast(-1.6548533913638905e-10), static_cast(-0.50000000016548529), static_cast(-169.20843148231506), SC_(1.00000000084161045914716192484600809610013447) }}, // Unexpected high error : 17825.7893791562892147339880466461181640625 Found : -0.000253525216373273569459012577453904668800532818 Expected : -0.000253525216374277052779756536082800266740377992 - { { static_cast(-2.0211181797563725e-14), static_cast(-1.0000000000000202), static_cast(-25.653068032115698), SC_(-0.000253525216374277055047768086884693917115210113) }}, + { { static_cast(-2.0211181797563725e-14), static_cast(-1.0000000000000202), static_cast(-25.653068032115698), SC_(-0.000253525216374277055047768086884693917115210113) }}, // Unexpected high error: 1.79769e+308 Found: -inf Expected: -2.63233e-197 - {{ static_cast(235.44106131792068), static_cast(-2.250966744069919e-13), static_cast(-974.28781914710999), SC_(-2.63233018990922939037251029961929844581862228e-197) }}, + {{ static_cast(235.44106131792068), static_cast(-2.250966744069919e-13), static_cast(-974.28781914710999), SC_(-2.63233018990922939037251029961929844581862228e-197) }}, // Unexpected high error : 1.79769313486231570814527423731704356798070568e+308 Found : -inf Expected : -3.53316570137147325345279499243339692001224196e+226 - { { static_cast(-235.44106131792068), static_cast(-2.250966744069919e-13), static_cast(974.28781914710999), SC_(-3.53316570137147343919975579872097464691424847e+226) }}, + { { static_cast(-235.44106131792068), static_cast(-2.250966744069919e-13), static_cast(974.28781914710999), SC_(-3.53316570137147343919975579872097464691424847e+226) }}, // Unexpected high error : 2.48268e+91 Found : -9.61305e+268 Expected : -1.74382e+193 - { { static_cast(5.9981750131794866e-15), static_cast(-230.70702263712883), static_cast(240.42092034220695), SC_(-1.74381782591884817018404492963109914357365958e+193) }}, + { { static_cast(5.9981750131794866e-15), static_cast(-230.70702263712883), static_cast(240.42092034220695), SC_(-1.74381782591884817018404492963109914357365958e+193) }}, // Unexpected high error : 1.79769313486231570814527423731704356798070568e+308 Found : -9.61305077326281580724540507004198316499661687e+268 Expected : 1.74381782591870724567837900957146707932623893e+193 - { { static_cast(-5.9981750131794866e-15), static_cast(-230.70702263712883), static_cast(240.42092034220695), SC_(1.74381782591870734495565763481520223752372107e+193) }}, - // Unexpected exception : Error in function boost::math::hypergeometric_pFq : Cancellation is so severe that no bits in the result are correct, last result was 3.0871891698197084e+73 - { { static_cast(-5.9981750131794866e-15), static_cast(0.499999999999994), static_cast(-240.42092034220695), SC_(1.00000000000004464930530925572237133417488137) }}, - // Unexpected exception : Error in function boost::math::hypergeometric_pFq : Cancellation is so severe that no bits in the result are correct, last result was 3.0871891698197084e+73 - { { static_cast(-5.9981750131794866e-15), static_cast(-0.500000000000006), static_cast(-240.42092034220695), SC_(1.00000000000003262784934420226963147689063665) }}, + { { static_cast(-5.9981750131794866e-15), static_cast(-230.70702263712883), static_cast(240.42092034220695), SC_(1.74381782591870734495565763481520223752372107e+193) }}, + // Unexpected exception : Error in function boost::math::hypergeometric_pFq : Cancellation is so severe that no bits in the result are correct, last result was 3.0871891698197084e+73 + { { static_cast(-5.9981750131794866e-15), static_cast(0.499999999999994), static_cast(-240.42092034220695), SC_(1.00000000000004464930530925572237133417488137) }}, + // Unexpected exception : Error in function boost::math::hypergeometric_pFq : Cancellation is so severe that no bits in the result are correct, last result was 3.0871891698197084e+73 + { { static_cast(-5.9981750131794866e-15), static_cast(-0.500000000000006), static_cast(-240.42092034220695), SC_(1.00000000000003262784934420226963147689063665) }}, // Unexpected high error : 18466.4373304979599197395145893096923828125 Found : 1.32865406167486480872551302123696359558380209e-08 Expected : 1.3286540616694168317751162703647255236560909e-08 - { { static_cast(6.772927684190258e-10), static_cast(-0.99999999932270722), static_cast(-483.69576895236969), SC_(1.32865406166941679958876322759721528297325713e-08) }}, + { { static_cast(6.772927684190258e-10), static_cast(-0.99999999932270722), static_cast(-483.69576895236969), SC_(1.32865406166941679958876322759721528297325713e-08) }}, // Unexpected high error: 1.79769e+308 Found: -nan(ind) Expected: 5.31173e-38 - {{ static_cast(6763.4877452850342), static_cast(3.6834977949762315e-08), static_cast(-210.20976513624191), SC_(5.31173132667573457976877380237496445775181141e-38) }}, + {{ static_cast(6763.4877452850342), static_cast(3.6834977949762315e-08), static_cast(-210.20976513624191), SC_(5.31173132667573457976877380237496445775181141e-38) }}, // Unexpected high error : 1.79769313486231570814527423731704356798070568e+308 Found : -nan(ind) Expected : 1.04274264437409856500364465136386556989276338e+54 - { { static_cast(-6763.4877452850342), static_cast(3.6834977949762315e-08), static_cast(210.20976513624191), SC_(1.04274264437409861991447530452939035771734596e+54) }}, + { { static_cast(-6763.4877452850342), static_cast(3.6834977949762315e-08), static_cast(210.20976513624191), SC_(1.04274264437409861991447530452939035771734596e+54) }}, // Unexpected high error: 3.17219226436543247206316287281668161679098192e+185 Found: 1.00012411189051491970538746574092795078602734e+201 Expected: 14198882672502154063215954231296 - {{ static_cast(76763.042617797852), static_cast(-21.722407214343548), static_cast(-0.60326536209322512), SC_(14198882672502153010712531896984.8126667697959) }}, + {{ static_cast(76763.042617797852), static_cast(-21.722407214343548), static_cast(-0.60326536209322512), SC_(14198882672502153010712531896984.8126667697959) }}, // Unexpected high error: 1.79769313486231570814527423731704356798070568e+308 Found: -2.39521645877904927856730119998375850409649219e+124 Expected: 2.3952164587795095929135248712964248422934629e+124 - {{ static_cast(-1.8857404964801872e-09), static_cast(-226.52341184020042), static_cast(160.86221924424171), SC_(2.39521645877950946848639784331327651190093595e+124) }}, + {{ static_cast(-1.8857404964801872e-09), static_cast(-226.52341184020042), static_cast(160.86221924424171), SC_(2.39521645877950946848639784331327651190093595e+124) }}, // Unexpected high error : 73027.246763920571538619697093963623046875 Found : 0.000111810625893715248580992382976262433658121154 Expected : 0.000111810625895528292111404111697225971511215903 - { { static_cast(-7.5220323642510769e-13), static_cast(-1.0000000000007523), static_cast(-17.948102783411741), SC_(0.00011181062589552829441403260197223627311023229) }}, + { { static_cast(-7.5220323642510769e-13), static_cast(-1.0000000000007523), static_cast(-17.948102783411741), SC_(0.00011181062589552829441403260197223627311023229) }}, // Unexpected high error: 111726.15160028330865316092967987060546875 Found: 0.00856985181006919560786627698689699172973632813 Expected: 0.00856985180985659310282098743982714950107038021 - {{ static_cast(5.6136137469239618e-15), static_cast(-0.99999999999999434), static_cast(-1989.8742001056671), SC_(0.00856985180985659334965068576732515544478559175) }}, + {{ static_cast(5.6136137469239618e-15), static_cast(-0.99999999999999434), static_cast(-1989.8742001056671), SC_(0.00856985180985659334965068576732515544478559175) }}, // Unexpected high error : 10431.000000023717802832834422588348388671875 Found : 0.99999999999772626324556767940521240234375 Expected : 1.00000000000004241051954068097984418272972107 - { { static_cast(-5.6136137469239618e-15), static_cast(-0.50000000000000566), static_cast(-1989.8742001056671), SC_(1.00000000000004243096226509338784935080089269) }}, + { { static_cast(-5.6136137469239618e-15), static_cast(-0.50000000000000566), static_cast(-1989.8742001056671), SC_(1.00000000000004243096226509338784935080089269) }}, // And more from error rate testing: {{ (T)std::ldexp((double)-17079780487168000, -44), (T)std::ldexp((double)9462273998848000, -46), (T)std::ldexp((double)9928190459904000, -48), SC_(7.7358754011357422722746277257633664799903803239195e-72) }}, {{ (T)std::ldexp((double)-16238757384192000, -44), (T)std::ldexp((double)17248812490752000, -44), (T)std::ldexp((double)12549255331840000, -49), SC_(4.7354970214088286546733909450191631190700414608975e-10) }}, - {{ static_cast(-6.8543290253728628), static_cast(607.72073245607316), static_cast(253.26409819535911), SC_(0.024418741483258497441042709681531519387974841769189) }}, + {{ static_cast(-6.8543290253728628), static_cast(607.72073245607316), static_cast(253.26409819535911), SC_(0.024418741483258497441042709681531519387974841769189) }}, {{ (T)std::ldexp((double)-15569844699136000, -52), (T)std::ldexp((double)12855440629760000, -44), (T)std::ldexp((double)12563412279296000, -45), SC_(0.097879401070280078654536987721507669872679020399179) }}, {{ (T)std::ldexp((double)-13521484578816000, -48), (T)std::ldexp((double)11813014388736000, -46), (T)std::ldexp((double)12736881098752000, -48), SC_(9.1262751214688536871555425535678062558805718157237e-08) }}, {{ (T)std::ldexp((double)-13125670141952000, -44), (T)std::ldexp((double)16524914262016000, -44), (T)std::ldexp((double)12270166867968000, -49), SC_(2.0809215788388623809065210261671764534436583442155e-08) }}, @@ -295,36 +295,36 @@ void test_spots6(T, const char* type_name) // // Special cases for a,b equal and negative integers, see: // - { {-1, -1, 0.9999999999990905052982270717620850, SC_(1.9999999999990905052982270717620850)} }, - { { -2, -2, 0.9999999999990905052982270717620850, SC_(2.4999999999981810105964545571144762) } }, - { { -3, -3, 0.9999999999990905052982270717620850, SC_(2.6666666666643929299122351732524916) } }, - { { -4, -4, 0.9999999999990905052982270717620850, SC_(2.7083333333309080141286065586746589) } }, - { { -5, -5, 0.9999999999990905052982270717620850, SC_(2.7166666666642034518493660889297968) } }, - { { -6, -6, 0.9999999999990905052982270717620850, SC_(2.7180555555530847616157402206496325) } }, - { { -7, -7, 0.9999999999990905052982270717620850, SC_(2.7182539682514961968413528410609673) } }, - { { -8, -8, 0.9999999999990905052982270717620850, SC_(2.7182787698387976036876421724036051) } }, - { { -9, -9, 0.9999999999990905052982270717620850, SC_(2.7182815255707199797197951818652022) } }, - { { -10, -10, 0.9999999999990905052982270717620850, SC_(2.7182818011439122170723781245206092) } }, - { { -11, -11, 0.9999999999990905052982270717620850, SC_(2.7182818261960206022634645412810530) } }, - { { -12, -12, 0.9999999999990905052982270717620850, SC_(2.7182818282836963010274896793573756) } }, - { { -13, -13, 0.9999999999990905052982270717620850, SC_(2.7182818284442867393938070953642430) } }, - { { -14, -14, 0.9999999999990905052982270717620850, SC_(2.7182818284557574849913907639252443) } }, - { { -15, -15, 0.9999999999990905052982270717620850, SC_(2.7182818284565222013645623129904880) } }, - { { -16, -16, 0.9999999999990905052982270717620850, SC_(2.7182818284565699961378854913379726) } }, - { { -17, -17, 0.9999999999990905052982270717620850, SC_(2.7182818284565728075951397933896427)} }, - { { -18, -18, 0.9999999999990905052982270717620850, SC_(2.7182818284565729637872094766949018) } }, - { { -19, -19, 0.9999999999990905052982270717620850, SC_(2.7182818284565729720078447231771757) } }, - { { -20, -20, 0.9999999999990905052982270717620850, SC_(2.7182818284565729724188764855009155) } }, - { { -21, -21, 0.9999999999990905052982270717620850, SC_(2.7182818284565729724384494265639330) } }, - { { -22, -22, 0.9999999999990905052982270717620850, SC_(2.7182818284565729724393391057031602) } }, - { { -23, -23, 0.9999999999990905052982270717620850, SC_(2.7182818284565729724393777874048657) } }, - { { -24, -24, 0.9999999999990905052982270717620850, SC_(2.7182818284565729724393793991424368) } }, - { { -25, -25, 0.9999999999990905052982270717620850, SC_(2.7182818284565729724393794636119396) } }, - { { -26, -26, 0.9999999999990905052982270717620850, SC_(2.7182818284565729724393794660915359) } }, - { { -27, -27, 0.9999999999990905052982270717620850, SC_(2.7182818284565729724393794661833728) } }, - { { -28, -28, 0.9999999999990905052982270717620850, SC_(2.7182818284565729724393794661866527) } }, - { { -29, -29, 0.9999999999990905052982270717620850, SC_(2.7182818284565729724393794661867658) } }, - { { -30, -30, 0.9999999999990905052982270717620850, SC_(2.7182818284565729724393794661867696) } }, + { {-1, -1, SC_(0.9999999999990905052982270717620850), SC_(1.9999999999990905052982270717620850)} }, + { { -2, -2, SC_(0.9999999999990905052982270717620850), SC_(2.4999999999981810105964545571144762) } }, + { { -3, -3, SC_(0.9999999999990905052982270717620850), SC_(2.6666666666643929299122351732524916) } }, + { { -4, -4, SC_(0.9999999999990905052982270717620850), SC_(2.7083333333309080141286065586746589) } }, + { { -5, -5, SC_(0.9999999999990905052982270717620850), SC_(2.7166666666642034518493660889297968) } }, + { { -6, -6, SC_(0.9999999999990905052982270717620850), SC_(2.7180555555530847616157402206496325) } }, + { { -7, -7, SC_(0.9999999999990905052982270717620850), SC_(2.7182539682514961968413528410609673) } }, + { { -8, -8, SC_(0.9999999999990905052982270717620850), SC_(2.7182787698387976036876421724036051) } }, + { { -9, -9, SC_(0.9999999999990905052982270717620850), SC_(2.7182815255707199797197951818652022) } }, + { { -10, -10, SC_(0.9999999999990905052982270717620850), SC_(2.7182818011439122170723781245206092) } }, + { { -11, -11, SC_(0.9999999999990905052982270717620850), SC_(2.7182818261960206022634645412810530) } }, + { { -12, -12, SC_(0.9999999999990905052982270717620850), SC_(2.7182818282836963010274896793573756) } }, + { { -13, -13, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284442867393938070953642430) } }, + { { -14, -14, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284557574849913907639252443) } }, + { { -15, -15, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284565222013645623129904880) } }, + { { -16, -16, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284565699961378854913379726) } }, + { { -17, -17, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284565728075951397933896427)} }, + { { -18, -18, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284565729637872094766949018) } }, + { { -19, -19, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284565729720078447231771757) } }, + { { -20, -20, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284565729724188764855009155) } }, + { { -21, -21, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284565729724384494265639330) } }, + { { -22, -22, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284565729724393391057031602) } }, + { { -23, -23, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284565729724393777874048657) } }, + { { -24, -24, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284565729724393793991424368) } }, + { { -25, -25, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284565729724393794636119396) } }, + { { -26, -26, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284565729724393794660915359) } }, + { { -27, -27, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284565729724393794661833728) } }, + { { -28, -28, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284565729724393794661866527) } }, + { { -29, -29, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284565729724393794661867658) } }, + { { -30, -30, SC_(0.9999999999990905052982270717620850), SC_(2.7182818284565729724393794661867696) } }, { {-1, -1, 23.5, 24.500000000000000000000000000000000} }, { { -2, -2, 23.5, SC_(300.62500000000000000000000000000000) } }, @@ -395,14 +395,18 @@ void test_spots6(T, const char* type_name) {{ 13, 1.5f, 61, SC_(1.35508577094765660270265300640877455638098585524020525369044e39)}}, {{ 13, 1.5f - T(1)/128, 61, SC_(1.40067238333701986992154961431485209677766220602448290643906e39)}}, {{ 13, 1.5f + T(1)/128, 61, SC_(1.31105748771677778012064837998217769289913724450105998963999e39)}}, + + // Uncovered code: + { { SC_(3.125), SC_(-4.25), SC_(0.25), SC_(0.8377654489787288706170314748630804839990046835429477943509538691)}}, + { { SC_(3.125), SC_(-4.25), SC_(6.25), SC_(-1.0732328736644883882525050418994056179896661423277961620744e8)}}, } }; static const std::array, 2> hypergeometric_1F1_big_bugs = { { #if DBL_MAX_EXP == LDBL_MAX_EXP - {{ static_cast(7.8238229420435346e-05), static_cast(-5485.3222503662109), static_cast(1896.0561106204987), BOOST_MATH_HUGE_CONSTANT(T, 1000, 4.33129800901478785957996719992774682013355926e+668) }}, - {{ static_cast(-7.8238229420435346e-05), static_cast(-5485.3222503662109), static_cast(1896.0561106204987), BOOST_MATH_HUGE_CONSTANT(T, 1000, -4.3248750673398590673783317624407455467680038e+668) }}, + {{ static_cast(7.8238229420435346e-05), static_cast(-5485.3222503662109), static_cast(1896.0561106204987), BOOST_MATH_HUGE_CONSTANT(T, 1000, 4.33129800901478785957996719992774682013355926e+668) }}, + {{ static_cast(-7.8238229420435346e-05), static_cast(-5485.3222503662109), static_cast(1896.0561106204987), BOOST_MATH_HUGE_CONSTANT(T, 1000, -4.3248750673398590673783317624407455467680038e+668) }}, #else - { { static_cast(7.8238229420435346e-05), static_cast(-5485.3222503662109), static_cast(1896.0561106204987), SC_(4.33129800901478785957996719992774682013355926e+668) } }, - { { static_cast(-7.8238229420435346e-05), static_cast(-5485.3222503662109), static_cast(1896.0561106204987), SC_(-4.3248750673398590673783317624407455467680038e+668) } }, + { { static_cast(7.8238229420435346e-05), static_cast(-5485.3222503662109), static_cast(1896.0561106204987), SC_(4.33129800901478785957996719992774682013355926e+668) } }, + { { static_cast(-7.8238229420435346e-05), static_cast(-5485.3222503662109), static_cast(1896.0561106204987), SC_(-4.3248750673398590673783317624407455467680038e+668) } }, #endif } }; do_test_1F1(hypergeometric_1F1_bugs, type_name, "Bug cases"); @@ -426,6 +430,15 @@ void test_spots7(T, const char* type_name) do_test_1F1(hypergeometric_1f1_neg_int, type_name, "Both parameters negative integers."); } +template +void test_spots8(T, const char* type_name) +{ + BOOST_CHECK_THROW(boost::math::hypergeometric_1F1(static_cast(2), static_cast(-100), static_cast(0.5)), std::domain_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_1F1(static_cast(-100.1), static_cast(-100), static_cast(0.5)), std::domain_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_1F1(static_cast(-102), static_cast(-100), static_cast(0.5)), std::domain_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_1F1(static_cast(-50.25), static_cast(-100), static_cast(0.5)), std::domain_error); +} + template void test_spots(T z, const char* type_name) { @@ -448,6 +461,7 @@ void test_spots(T z, const char* type_name) if(std::numeric_limits::digits >= std::numeric_limits::digits && std::numeric_limits::digits <= 128) test_spots6(z, type_name); test_spots7(z, type_name); + test_spots8(z, type_name); } diff --git a/test/test_1F1_log.hpp b/test/test_1F1_log.hpp index d3b20fe248..6a840cc8dd 100644 --- a/test/test_1F1_log.hpp +++ b/test/test_1F1_log.hpp @@ -77,6 +77,19 @@ void test_spots2(T, const char* type_name) do_test_1F1(hypergeometric_1f1_log_large_unsolved, type_name, "Large random values - log - unsolved"); } +template +void test_sign(T, const char* type_name) +{ +#include "hypergeometric_1F1_small_random.ipp" + + for (unsigned i = 0; i < hypergeometric_1F1_small_random.size(); ++i) + { + int s; + boost::math::log_hypergeometric_1F1(hypergeometric_1F1_small_random[i][0], hypergeometric_1F1_small_random[i][1], hypergeometric_1F1_small_random[i][2], &s); + BOOST_CHECK_EQUAL(s, boost::math::sign(hypergeometric_1F1_small_random[i][3])); + } +} + template void test_spots_bugs(T, const char* type_name) { @@ -99,6 +112,7 @@ void test_spots(T z, const char* type_name) { test_spots1(z, type_name); test_spots_bugs(z, type_name); + test_sign(z, type_name); #ifdef TEST_UNSOLVED test_spots2(z, type_name); #endif diff --git a/test/test_1F1_regularized.cpp b/test/test_1F1_regularized.cpp index e16e32943f..f62e6f86fa 100644 --- a/test/test_1F1_regularized.cpp +++ b/test/test_1F1_regularized.cpp @@ -162,11 +162,6 @@ BOOST_AUTO_TEST_CASE( test_main ) expected_results(); BOOST_MATH_CONTROL_FP; -#if !defined(TEST) || (TEST == 1) - test_hypergeometric_mellin_transform(); - test_hypergeometric_laplace_transform(); -#endif - #ifndef BOOST_MATH_BUGGY_LARGE_FLOAT_CONSTANTS #if !defined(TEST) || (TEST == 2) test_spots(0.0F, "float"); diff --git a/test/test_1F1_regularized.hpp b/test/test_1F1_regularized.hpp index 90415c6de1..e4d2d4e2f9 100644 --- a/test/test_1F1_regularized.hpp +++ b/test/test_1F1_regularized.hpp @@ -77,5 +77,21 @@ template void test_spots(T z, const char* type_name) { test_spots1(z, type_name); + + // + // Special cases for coverage: + // + BOOST_IF_CONSTEXPR(std::numeric_limits::has_infinity) + { + BOOST_CHECK_EQUAL(boost::math::hypergeometric_1F1_regularized(600.25, 102.75, 11512.0), std::numeric_limits::infinity()); + BOOST_IF_CONSTEXPR(std::numeric_limits::max_exponent == 16384) + { + BOOST_CHECK_EQUAL(boost::math::hypergeometric_1F1_regularized(600.25, 102.75, 9985.0), std::numeric_limits::infinity()); + } + else BOOST_IF_CONSTEXPR(std::numeric_limits::max_exponent <= std::numeric_limits::max_exponent) + { + BOOST_CHECK_EQUAL(boost::math::hypergeometric_1F1_regularized(600.25, 102.75, 514.0), std::numeric_limits::infinity()); + } + } } diff --git a/test/test_2F0.hpp b/test/test_2F0.hpp index d0cacf9e38..4113ebd2c5 100644 --- a/test/test_2F0.hpp +++ b/test/test_2F0.hpp @@ -99,4 +99,9 @@ void test_spots(T z, const char* type_name) a1 = -0.5f; BOOST_CHECK((boost::math::isinf)(boost::math::hypergeometric_2F0(a1, a2, z))); } + // + // Coverage: + // + BOOST_CHECK_EQUAL(boost::math::hypergeometric_2F0(T(0), T(20), T(2)), T(1)); + BOOST_CHECK_EQUAL(boost::math::hypergeometric_2F0(T(20), T(0), T(2)), T(1)); } From 2b4512a7cdf3489e769ce087877ca483d8de4e92 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sun, 14 Jul 2024 11:17:21 +0100 Subject: [PATCH 25/27] 1F1 and 2F0 coverage. --- include/boost/math/special_functions/hypergeometric_1F1.hpp | 2 +- test/test_1F1_regularized.hpp | 2 +- test/test_2F0.hpp | 1 + 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/include/boost/math/special_functions/hypergeometric_1F1.hpp b/include/boost/math/special_functions/hypergeometric_1F1.hpp index 35d087ed6f..21e25670c8 100644 --- a/include/boost/math/special_functions/hypergeometric_1F1.hpp +++ b/include/boost/math/special_functions/hypergeometric_1F1.hpp @@ -140,7 +140,7 @@ namespace boost { namespace math { namespace detail { return hypergeometric_1F1_from_function_ratio_negative_b_forwards(a, b, z, pol, log_scaling); } #endif - // LCOV_EXCL_END + // LCOV_EXCL_STOP } // // We could fall back to Tricomi's approximation if we're in the transition zone diff --git a/test/test_1F1_regularized.hpp b/test/test_1F1_regularized.hpp index e4d2d4e2f9..2221e69d5e 100644 --- a/test/test_1F1_regularized.hpp +++ b/test/test_1F1_regularized.hpp @@ -84,7 +84,7 @@ void test_spots(T z, const char* type_name) BOOST_IF_CONSTEXPR(std::numeric_limits::has_infinity) { BOOST_CHECK_EQUAL(boost::math::hypergeometric_1F1_regularized(600.25, 102.75, 11512.0), std::numeric_limits::infinity()); - BOOST_IF_CONSTEXPR(std::numeric_limits::max_exponent == 16384) + BOOST_IF_CONSTEXPR(std::numeric_limits::max_exponent > 16000) { BOOST_CHECK_EQUAL(boost::math::hypergeometric_1F1_regularized(600.25, 102.75, 9985.0), std::numeric_limits::infinity()); } diff --git a/test/test_2F0.hpp b/test/test_2F0.hpp index 4113ebd2c5..7a10f97c82 100644 --- a/test/test_2F0.hpp +++ b/test/test_2F0.hpp @@ -104,4 +104,5 @@ void test_spots(T z, const char* type_name) // BOOST_CHECK_EQUAL(boost::math::hypergeometric_2F0(T(0), T(20), T(2)), T(1)); BOOST_CHECK_EQUAL(boost::math::hypergeometric_2F0(T(20), T(0), T(2)), T(1)); + BOOST_CHECK_EQUAL(boost::math::hypergeometric_2F0(T(20), T(10), T(0)), T(1)); } From 36d46448a9cb3dd6e7c12457c3e233226fef999e Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sun, 14 Jul 2024 11:45:09 +0100 Subject: [PATCH 26/27] hypergeometric_0F1_bessel.hpp coverage. --- .../detail/hypergeometric_0F1_bessel.hpp | 9 ++++----- test/test_0F1.cpp | 2 +- test/test_0F1.hpp | 5 ++++- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/include/boost/math/special_functions/detail/hypergeometric_0F1_bessel.hpp b/include/boost/math/special_functions/detail/hypergeometric_0F1_bessel.hpp index 31ab073766..be91d406c2 100644 --- a/include/boost/math/special_functions/detail/hypergeometric_0F1_bessel.hpp +++ b/include/boost/math/special_functions/detail/hypergeometric_0F1_bessel.hpp @@ -20,12 +20,11 @@ { BOOST_MATH_STD_USING - const bool is_z_nonpositive = z <= 0; + //const bool is_z_nonpositive = z <= 0; + BOOST_MATH_ASSERT(z < 0); // condition used at call site - const T sqrt_z = is_z_nonpositive ? T(sqrt(-z)) : T(sqrt(z)); - const T bessel_mult = is_z_nonpositive ? - boost::math::cyl_bessel_j(b - 1, 2 * sqrt_z, pol) : - boost::math::cyl_bessel_i(b - 1, 2 * sqrt_z, pol) ; + const T sqrt_z = sqrt(-z); + const T bessel_mult = boost::math::cyl_bessel_j(b - 1, 2 * sqrt_z, pol); if (b > boost::math::max_factorial::value) { diff --git a/test/test_0F1.cpp b/test/test_0F1.cpp index 26e1e304b2..f597128c41 100644 --- a/test/test_0F1.cpp +++ b/test/test_0F1.cpp @@ -49,7 +49,7 @@ void expected_results() ".*", // platform largest_type, // test type(s) "Large.*", // test data group - ".*", 400, 100); // test function + ".*", 1100, 300); // test function // diff --git a/test/test_0F1.hpp b/test/test_0F1.hpp index 05bd4d8d25..160234a6a1 100644 --- a/test/test_0F1.hpp +++ b/test/test_0F1.hpp @@ -178,7 +178,7 @@ void test_spots(T, const char* type_name) {SC_(19.75), SC_(-20.25), SC_(0.34910186209886263988734004369646266800659689838760)}, {SC_(19.75), SC_(-16.25), SC_(0.43171851514846723176830271990344447636246814031258)}, {SC_(19.75), SC_(-12.25), SC_(0.53264760691226214051563276361223610040193575209722)}, {SC_(19.75), SC_(-8.25), SC_(0.65570932018638949205413006445398352201545799769694)}, {SC_(19.75), SC_(-4.25), SC_(0.80547695536244789152425915957330496255335464733020)}, {SC_(19.75), SC_(-0.25), SC_(0.98741773517508964436868231216434995093107999788450)}, {SC_(19.75), SC_(3.75), SC_(1.2080586524527299846909832323308446155063281314444)}, {SC_(19.75), SC_(7.75), SC_(1.4751816095060559707610879513647971547078203326785)}, {SC_(19.75), SC_(11.75), SC_(1.7980527870029919079278807774802021351478928537606)}, {SC_(19.75), SC_(15.75), SC_(2.1876919765295108137252526717389200156554872202073)}, {SC_(19.75), SC_(19.75), SC_(2.6571885305131654484163938203293169727636168525176)} } }; - static const std::array, 49> hypergeometric_0F1_large_data = { { + static const std::array, 52> hypergeometric_0F1_large_data = { { {SC_(-3000.25), SC_(-3000.25), SC_(2.7187352281773126038196142076996225709478475910049)}, {SC_(-3000.25), SC_(-2000.25), SC_(1.9479325178881941977364868969702984561194629278005)}, {SC_(-3000.25), SC_(-1000.25), SC_(1.3957158200680510708579318399559946925035296441969)}, @@ -210,6 +210,9 @@ void test_spots(T, const char* type_name) {SC_(999.75), SC_(-3000.25), SC_(0.049513102314849089089531272645546209356947364503693)}, {SC_(999.75), SC_(-2000.25), SC_(0.13496287555692361766808821455578584612232691847368)}, {SC_(999.75), SC_(-1000.25), SC_(0.36751140843960205420082946295266027163679747815171)}, {SC_(999.75), SC_(-0.25), SC_(0.99974996871616159816731662896550362438335384158579)}, {SC_(999.75), SC_(999.75), SC_(2.7169258487486599795213704968075693651612629011137)}, {SC_(999.75), SC_(1999.75), SC_(7.3761834991423536961023917193398034973414374920724)}, {SC_(999.75), SC_(2999.75), SC_(20.005752597249211333476751024563565243417014841294)}, {SC_(1999.75), SC_(-3000.25), SC_(0.22293486720691878909733143630236248153349047769471)}, {SC_(1999.75), SC_(-2000.25), SC_(0.36769546312912945115503617374640562125026416585911)}, {SC_(1999.75), SC_(-1000.25), SC_(0.60637900493362087991135014744274427763909334895653)}, {SC_(1999.75), SC_(-0.25), SC_(0.99987499218326921518475889542511447173965188103381)}, {SC_(1999.75), SC_(999.75), SC_(1.6485152791856256326123376967789594359576461105786)}, {SC_(1999.75), SC_(1999.75), SC_(2.7176030495659989803962821602572269144434394541710)}, {SC_(1999.75), SC_(2999.75), SC_(4.4794520688544123057192034298995619570726142156707)}, {SC_(2999.75), SC_(-3000.25), SC_(0.36775679765931355729803647936785321939707796855752)}, {SC_(2999.75), SC_(-2000.25), SC_(0.51330776829498225175667208437735196999610552588713)}, {SC_(2999.75), SC_(-1000.25), SC_(0.71643841877594009499597675075273584291266593195763)}, {SC_(2999.75), SC_(-0.25), SC_(0.99991666319319078123594326189204686183292938405581)}, {SC_(2999.75), SC_(999.75), SC_(1.3955090626492889292016427682906884819283986108381)}, {SC_(2999.75), SC_(1999.75), SC_(1.9475357570910758403975860806320078849962664462281)}, {SC_(2999.75), SC_(2999.75), SC_(2.7178291334815105222784044221313459883434486746788)}, + {SC_(190.25), SC_(-10000.75), SC_(1.4750057478240259502978056596880193011950123904017370507273e-31)}, + {SC_(192.25), SC_(-10000.75), SC_(1.3702644191119456726760786567669231339553195106433263692267e-29)}, + {SC_(192.25), SC_(-8000.75), SC_(8.9691784861188665237766513631329658977665974807542912974906e-22)}, } }; do_test_0F1(hypergeometric_0F1_integer_data, type_name, "Integer Arguments"); From 977feaf84acba11c93de94712328441c0f035590 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Tue, 6 Aug 2024 12:09:33 +0100 Subject: [PATCH 27/27] 1F1 coverage. --- .../detail/hypergeometric_1F1_bessel.hpp | 27 ++++++++++++++++--- test/test_1F1.hpp | 3 ++- 2 files changed, 26 insertions(+), 4 deletions(-) diff --git a/include/boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp b/include/boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp index 860217b059..07d9dabf89 100644 --- a/include/boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp +++ b/include/boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp @@ -74,8 +74,12 @@ { // We get very limited precision due to rapid denormalisation/underflow of the Bessel values, raise an exception and try something else: policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Underflow in Bessel functions", bessel_cache[cache_size - 1], pol); + // Exceptions are off if we get here, just fill the cache with NaN's and we'll let this method fail and fallback later: + std::fill(bessel_cache.begin(), bessel_cache.end(), std::numeric_limits::quiet_NaN()); + cache_offset = -cache_size; + return; } - if ((term * bessel_cache[cache_size - 1] < tools::min_value() / (tools::epsilon() * tools::epsilon())) || !(boost::math::isfinite)(term) || (!std::numeric_limits::has_infinity && (fabs(term) > tools::max_value()))) + if ((fabs(term * bessel_cache[cache_size - 1]) < tools::min_value() / (tools::epsilon() * tools::epsilon())) || !(boost::math::isfinite)(term) || (!std::numeric_limits::has_infinity && (fabs(term) > tools::max_value()))) { term = -log(fabs(bessel_arg)) * b_minus_1_plus_n / 2; log_scale = lltrunc(term); @@ -88,15 +92,27 @@ if constexpr (std::numeric_limits::has_infinity) { if (!(boost::math::isfinite)(bessel_cache[cache_size - 1])) + { policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Expected finite Bessel function result but got %1%", bessel_cache[cache_size - 1], pol); + // Exceptions are off if we get here, just fill the cache with NaN's and we'll let this method fail and fallback later: + std::fill(bessel_cache.begin(), bessel_cache.end(), std::numeric_limits::quiet_NaN()); + } } else if ((boost::math::isnan)(bessel_cache[cache_size - 1]) || (fabs(bessel_cache[cache_size - 1]) >= tools::max_value())) + { policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Expected finite Bessel function result but got %1%", bessel_cache[cache_size - 1], pol); + // Exceptions are off if we get here, just fill the cache with NaN's and we'll let this method fail and fallback later: + std::fill(bessel_cache.begin(), bessel_cache.end(), std::numeric_limits::quiet_NaN()); + } #else if ((std::numeric_limits::has_infinity && !(boost::math::isfinite)(bessel_cache[cache_size - 1])) || (!std::numeric_limits::has_infinity && ((boost::math::isnan)(bessel_cache[cache_size - 1]) || (fabs(bessel_cache[cache_size - 1]) >= tools::max_value())))) + { policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Expected finite Bessel function result but got %1%", bessel_cache[cache_size - 1], pol); + // Exceptions are off if we get here, just fill the cache with NaN's and we'll let this method fail and fallback later: + std::fill(bessel_cache.begin(), bessel_cache.end(), std::numeric_limits::quiet_NaN()); + } #endif cache_offset = -cache_size; refill_cache(); @@ -108,8 +124,13 @@ // very small (or zero) when b == 2a: // BOOST_MATH_STD_USING + // + // Except in the multiprecision case, we have probably illiminated anything + // would need more than the default 64 Bessel Functions. Anything more + // than that risks becoming a divergent series anyway... + // if(n - 2 - cache_offset >= cache_size) - refill_cache(); + refill_cache(); // LCOV_EXCL_LINE T result = A_minus_2 * term * bessel_cache[n - 2 - cache_offset]; term *= mult; ++n; @@ -122,7 +143,7 @@ if (A_minus_2 != 0) { if (n - 2 - cache_offset >= cache_size) - refill_cache(); + refill_cache(); // LCOV_EXCL_LINE result += A_minus_2 * term * bessel_cache[n - 2 - cache_offset]; } term *= mult; diff --git a/test/test_1F1.hpp b/test/test_1F1.hpp index 1d4b9b307e..b4c67faf1a 100644 --- a/test/test_1F1.hpp +++ b/test/test_1F1.hpp @@ -162,7 +162,7 @@ void test_spots5(T, const char* type_name) template void test_spots6(T, const char* type_name) { - static const std::array, 188> hypergeometric_1F1_bugs = { { + static const std::array, 189> hypergeometric_1F1_bugs = { { { { static_cast(17955.561660766602), static_cast(9.6968994205831605e-09), static_cast(-82.406154185533524), SC_(6.98056008378736714088730927132364938220428678e-11) }}, { { static_cast(17955.561660766602), static_cast(-9.6968994205831605e-09), static_cast(-82.406154185533524), SC_(-6.98055306629610746072607353939306734740549551e-11) }}, { { static_cast(-17955.561660766602), static_cast(-9.6968994205831605e-09), static_cast(82.406154185533524), SC_(-42897094853118832762870100.8669248353530950866) }} , @@ -399,6 +399,7 @@ void test_spots6(T, const char* type_name) // Uncovered code: { { SC_(3.125), SC_(-4.25), SC_(0.25), SC_(0.8377654489787288706170314748630804839990046835429477943509538691)}}, { { SC_(3.125), SC_(-4.25), SC_(6.25), SC_(-1.0732328736644883882525050418994056179896661423277961620744e8)}}, + { { SC_(-1792.0615542826642), SC_(1252.5126953125000), SC_(15.873352050781250), SC_(9.5429538115592332178745238256843491134483360764267909811922e-11)}}, } }; static const std::array, 2> hypergeometric_1F1_big_bugs = { { #if DBL_MAX_EXP == LDBL_MAX_EXP