From 9ca40b9f016bf8f4a7211a8446492ed2aeb06e6c Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Wed, 6 Mar 2024 18:47:48 +0000 Subject: [PATCH] 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 57ff60529..54964617b 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 1475a9a88..9b2b42055 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 eec635603..f53e42ab2 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 56b6cd12f..886e26d25 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 3f44a8091..3db40d70e 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 491db2fcd..d398fa162 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()); + } }