From aad4f85955459fbae1daa0e44bf491fc88138e0d Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Thu, 24 Nov 2022 18:43:32 +0000 Subject: [PATCH 1/5] Improve powm1 error handling. Makes 0^-n an overflow error (which matches std::pow which returns +INF rather than a NaN). Fixes https://github.com/boostorg/math/issues/781. --- doc/sf/powers.qbk | 2 ++ include/boost/math/special_functions/powm1.hpp | 10 ++++++++-- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/doc/sf/powers.qbk b/doc/sf/powers.qbk index 4f1db7a59d..363c9d90ee 100644 --- a/doc/sf/powers.qbk +++ b/doc/sf/powers.qbk @@ -281,6 +281,8 @@ when T1 and T2 are different types. There are two domains where this is useful: when /y/ is very small, or when /x/ is close to 1. +Note that for invalid input this function may raise a __domain_error or __overflow_error as appropriate. + Implemented in terms of `expm1`. The following graph illustrates the behaviour of powm1: diff --git a/include/boost/math/special_functions/powm1.hpp b/include/boost/math/special_functions/powm1.hpp index de23714370..c81d84ea09 100644 --- a/include/boost/math/special_functions/powm1.hpp +++ b/include/boost/math/special_functions/powm1.hpp @@ -16,6 +16,7 @@ #include #include #include +#include #include namespace boost{ namespace math{ namespace detail{ @@ -40,7 +41,7 @@ inline T powm1_imp(const T x, const T y, const Policy& pol) // fall through.... } } - else if (x < 0) + else if (boost::math::signbit(x)) // Need to error check -0 here as well { // y had better be an integer: if (boost::math::trunc(y) != y) @@ -48,7 +49,12 @@ inline T powm1_imp(const T x, const T y, const Policy& pol) if (boost::math::trunc(y / 2) == y / 2) return powm1_imp(T(-x), y, pol); } - return pow(x, y) - 1; + T result = pow(x, y) - 1; + if((boost::math::isinf)(result)) + return result < 0 ? -boost::math::policies::raise_overflow_error(function, nullptr, pol) : boost::math::policies::raise_overflow_error(function, nullptr, pol); + if((boost::math::isnan)(result)) + return boost::math::policies::raise_domain_error(function, "Result of pow is complex or undefined", x, pol); + return result; } } // detail From 8add2b43739322188594b1cb77821e2c61f6c99f Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 25 Nov 2022 17:08:58 +0000 Subject: [PATCH 2/5] Correct handling of NaN's and infinities in ibeta. Update tests. Fixes https://github.com/boostorg/math/issues/878 --- include/boost/math/special_functions/beta.hpp | 14 +++++++++++ test/test_ibeta.hpp | 18 ++++++++++++++ test/test_ibeta_derivative.hpp | 18 ++++++++++++++ test/test_ibeta_inv.hpp | 17 +++++++++++++ test/test_ibeta_inv_ab.hpp | 24 +++++++++++++++++++ 5 files changed, 91 insertions(+) diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index c950670383..4f0a9b61bc 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -1017,6 +1017,13 @@ T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_de BOOST_MATH_ASSERT((p_derivative == 0) || normalised); + if(!(boost::math::isfinite)(a)) + return policies::raise_domain_error(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol); + if(!(boost::math::isfinite)(b)) + return policies::raise_domain_error(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol); + if(!(boost::math::isfinite)(x)) + return policies::raise_domain_error(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol); + if(p_derivative) *p_derivative = -1; // value not set. @@ -1411,6 +1418,13 @@ T ibeta_derivative_imp(T a, T b, T x, const Policy& pol) // // start with the usual error checks: // + if (!(boost::math::isfinite)(a)) + return policies::raise_domain_error(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol); + if (!(boost::math::isfinite)(b)) + return policies::raise_domain_error(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol); + if (!(boost::math::isfinite)(x)) + return policies::raise_domain_error(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol); + if(a <= 0) return policies::raise_domain_error(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol); if(b <= 0) diff --git a/test/test_ibeta.hpp b/test/test_ibeta.hpp index 9d63ff0543..005cd1e52b 100644 --- a/test/test_ibeta.hpp +++ b/test/test_ibeta.hpp @@ -301,6 +301,24 @@ void test_spots(T) BOOST_MATH_CHECK_THROW(::boost::math::ibetac(static_cast(2), static_cast(2), static_cast(-0.5)), std::domain_error); BOOST_MATH_CHECK_THROW(::boost::math::ibetac(static_cast(2), static_cast(2), static_cast(1.5)), std::domain_error); + if (std::numeric_limits::has_quiet_NaN) + { + T n = std::numeric_limits::quiet_NaN(); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta(n, static_cast(2.125), static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta(static_cast(2.125), n, static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta(static_cast(2.125), static_cast(1.125), n), std::domain_error); + } + if (std::numeric_limits::has_infinity) + { + T n = std::numeric_limits::infinity(); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta(n, static_cast(2.125), static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta(static_cast(2.125), n, static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta(static_cast(2.125), static_cast(1.125), n), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta(-n, static_cast(2.125), static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta(static_cast(2.125), -n, static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta(static_cast(2.125), static_cast(1.125), -n), std::domain_error); + } + // // a = b = 0.5 is a special case: // diff --git a/test/test_ibeta_derivative.hpp b/test/test_ibeta_derivative.hpp index 5bc7be584c..bf820ade3b 100644 --- a/test/test_ibeta_derivative.hpp +++ b/test/test_ibeta_derivative.hpp @@ -130,5 +130,23 @@ void test_spots(T) static_cast(4.5), ldexp(static_cast(1), -557)), static_cast(5.24647512910420109893867082626308082567071751558842352760e-167L), tolerance * 4); + + if (std::numeric_limits::has_quiet_NaN) + { + T n = std::numeric_limits::quiet_NaN(); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_derivative(n, static_cast(2.125), static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_derivative(static_cast(2.125), n, static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_derivative(static_cast(2.125), static_cast(1.125), n), std::domain_error); + } + if (std::numeric_limits::has_infinity) + { + T n = std::numeric_limits::infinity(); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_derivative(n, static_cast(2.125), static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_derivative(static_cast(2.125), n, static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_derivative(static_cast(2.125), static_cast(1.125), n), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_derivative(-n, static_cast(2.125), static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_derivative(static_cast(2.125), -n, static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_derivative(static_cast(2.125), static_cast(1.125), -n), std::domain_error); + } } diff --git a/test/test_ibeta_inv.hpp b/test/test_ibeta_inv.hpp index 5eea61dee3..1e1c841f2f 100644 --- a/test/test_ibeta_inv.hpp +++ b/test/test_ibeta_inv.hpp @@ -288,5 +288,22 @@ void test_spots(T) static_cast(1) / static_cast(10)), static_cast(1), tolerance); } + if (std::numeric_limits::has_quiet_NaN) + { + T n = std::numeric_limits::quiet_NaN(); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_inv(n, static_cast(2.125), static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_inv(static_cast(2.125), n, static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_inv(static_cast(2.125), static_cast(1.125), n), std::domain_error); + } + if (std::numeric_limits::has_infinity) + { + T n = std::numeric_limits::infinity(); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_inv(n, static_cast(2.125), static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_inv(static_cast(2.125), n, static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_inv(static_cast(2.125), static_cast(1.125), n), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_inv(-n, static_cast(2.125), static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_inv(static_cast(2.125), -n, static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_inv(static_cast(2.125), static_cast(1.125), -n), std::domain_error); + } } diff --git a/test/test_ibeta_inv_ab.hpp b/test/test_ibeta_inv_ab.hpp index e270010e16..7982772838 100644 --- a/test/test_ibeta_inv_ab.hpp +++ b/test/test_ibeta_inv_ab.hpp @@ -212,5 +212,29 @@ void test_beta(T, const char* name) test_inverses2(ibeta_inva_data, name, "Inverse incomplete beta"); } #endif + // + // Special spot tests and bug reports: + // + if (std::numeric_limits::has_quiet_NaN) + { + T n = std::numeric_limits::quiet_NaN(); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_inva(n, static_cast(2.125), static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_inva(static_cast(2.125), n, static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_inva(static_cast(2.125), static_cast(1.125), n), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_invb(n, static_cast(2.125), static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_invb(static_cast(2.125), n, static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_invb(static_cast(2.125), static_cast(1.125), n), std::domain_error); + } + if (std::numeric_limits::has_infinity) + { + T n = std::numeric_limits::infinity(); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_invb(n, static_cast(2.125), static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_invb(static_cast(2.125), n, static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_invb(static_cast(2.125), static_cast(1.125), n), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_invb(-n, static_cast(2.125), static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_invb(static_cast(2.125), -n, static_cast(0.125)), std::domain_error); + BOOST_MATH_CHECK_THROW(::boost::math::ibeta_invb(static_cast(2.125), static_cast(1.125), -n), std::domain_error); + } + } From 18fd65aaa1c516c6b3d807d8bd6e5aeeaef4c2e5 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 25 Nov 2022 17:45:01 +0000 Subject: [PATCH 3/5] Make inverses overflow safe. --- .../math/special_functions/detail/ibeta_inverse.hpp | 13 +++++++++++-- include/boost/math/tools/roots.hpp | 10 +++++++++- 2 files changed, 20 insertions(+), 3 deletions(-) diff --git a/include/boost/math/special_functions/detail/ibeta_inverse.hpp b/include/boost/math/special_functions/detail/ibeta_inverse.hpp index b9f0995a6d..90a94e644d 100644 --- a/include/boost/math/special_functions/detail/ibeta_inverse.hpp +++ b/include/boost/math/special_functions/detail/ibeta_inverse.hpp @@ -631,8 +631,17 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py) // Try and compute the easy way first: // T bet = 0; - if(b < 2) - bet = boost::math::beta(a, b, pol); + if (b < 2) + { + try + { + bet = boost::math::beta(a, b, pol); + } + catch (const std::overflow_error&) + { + bet = tools::max_value(); + } + } if(bet != 0) { y = pow(b * q * bet, 1/b); diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index cb2c891239..ccf908f217 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -543,7 +543,15 @@ namespace detail { last_f0 = f0; delta2 = delta1; delta1 = delta; - detail::unpack_tuple(f(result), f0, f1, f2); + try + { + detail::unpack_tuple(f(result), f0, f1, f2); + } + catch (const std::overflow_error&) + { + f0 = max > 0 ? tools::max_value() : -tools::min_value(); + f1 = f2 = 0; + } --count; BOOST_MATH_INSTRUMENT_VARIABLE(f0); From c6fccc286e8c5d2f72dca2ad983bfa2aa7f7d93d Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 25 Nov 2022 19:51:00 +0000 Subject: [PATCH 4/5] Correct concept failures. --- include/boost/math/special_functions/powm1.hpp | 3 +-- test/git_issue_705.cpp | 2 ++ test/std_real_concept_check.cpp | 3 +++ 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/include/boost/math/special_functions/powm1.hpp b/include/boost/math/special_functions/powm1.hpp index c81d84ea09..05b5a5d3fa 100644 --- a/include/boost/math/special_functions/powm1.hpp +++ b/include/boost/math/special_functions/powm1.hpp @@ -26,7 +26,6 @@ inline T powm1_imp(const T x, const T y, const Policy& pol) { BOOST_MATH_STD_USING static const char* function = "boost::math::powm1<%1%>(%1%, %1%)"; - if (x > 0) { if ((fabs(y * (x - 1)) < 0.5) || (fabs(y) < 0.2)) @@ -41,7 +40,7 @@ inline T powm1_imp(const T x, const T y, const Policy& pol) // fall through.... } } - else if (boost::math::signbit(x)) // Need to error check -0 here as well + else if ((boost::math::signbit)(x)) // Need to error check -0 here as well { // y had better be an integer: if (boost::math::trunc(y) != y) diff --git a/test/git_issue_705.cpp b/test/git_issue_705.cpp index f9476aa53f..4b205421a7 100644 --- a/test/git_issue_705.cpp +++ b/test/git_issue_705.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_OVERFLOW_ERROR_POLICY ignore_error + #include "math_unit_test.hpp" #include #include diff --git a/test/std_real_concept_check.cpp b/test/std_real_concept_check.cpp index 4752f9608e..3edccb7381 100644 --- a/test/std_real_concept_check.cpp +++ b/test/std_real_concept_check.cpp @@ -71,6 +71,9 @@ struct numeric_limits static const bool traps = false; static const bool tinyness_before = false; static const float_round_style round_style = round_toward_zero; +#ifndef BOOST_NO_CXX11_NUMERIC_LIMITS + static const int max_digits10 = digits10 + 2; +#endif }; } #endif From ea70672db7f894d256dada5370539fa616bbb7c0 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sat, 26 Nov 2022 09:20:37 +0000 Subject: [PATCH 5/5] Fix gcc -fno-exception build. --- include/boost/math/special_functions/detail/ibeta_inverse.hpp | 4 ++++ include/boost/math/tools/roots.hpp | 4 ++++ 2 files changed, 8 insertions(+) diff --git a/include/boost/math/special_functions/detail/ibeta_inverse.hpp b/include/boost/math/special_functions/detail/ibeta_inverse.hpp index 90a94e644d..9b7bb50b80 100644 --- a/include/boost/math/special_functions/detail/ibeta_inverse.hpp +++ b/include/boost/math/special_functions/detail/ibeta_inverse.hpp @@ -633,14 +633,18 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py) T bet = 0; if (b < 2) { +#ifndef BOOST_NO_EXCEPTIONS try +#endif { bet = boost::math::beta(a, b, pol); } +#ifndef BOOST_NO_EXCEPTIONS catch (const std::overflow_error&) { bet = tools::max_value(); } +#endif } if(bet != 0) { diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index ccf908f217..6493e81228 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -543,15 +543,19 @@ namespace detail { last_f0 = f0; delta2 = delta1; delta1 = delta; +#ifndef BOOST_NO_EXCEPTIONS try +#endif { detail::unpack_tuple(f(result), f0, f1, f2); } +#ifndef BOOST_NO_EXCEPTIONS catch (const std::overflow_error&) { f0 = max > 0 ? tools::max_value() : -tools::min_value(); f1 = f2 = 0; } +#endif --count; BOOST_MATH_INSTRUMENT_VARIABLE(f0);