From bf3bc2e6c27fd86bc729e4fe1d2353d96afd3b3d Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sun, 5 Mar 2023 13:18:27 +0000 Subject: [PATCH] Fix ibeta_inv for very small p. (#962) * Fix ibeta_inv for very small p. Change assert's in temme_method_1_ibeta_inverse to corrections when guess goes out of range. Change handling of non-convergence in second_order_root_finder to use bracketing when the end points are many orders of magnitude apart. Fixes: https://github.com/boostorg/math/issues/961. * Add missing copyright. [CI SKIP] --- .../detail/ibeta_inverse.hpp | 19 +++-- include/boost/math/tools/roots.hpp | 24 ++++-- test/Jamfile.v2 | 1 + test/git_issue_961.cpp | 76 +++++++++++++++++++ 4 files changed, 110 insertions(+), 10 deletions(-) create mode 100644 test/git_issue_961.cpp diff --git a/include/boost/math/special_functions/detail/ibeta_inverse.hpp b/include/boost/math/special_functions/detail/ibeta_inverse.hpp index 9f257793d1..79eb95c45f 100644 --- a/include/boost/math/special_functions/detail/ibeta_inverse.hpp +++ b/include/boost/math/special_functions/detail/ibeta_inverse.hpp @@ -118,9 +118,17 @@ T temme_method_1_ibeta_inverse(T a, T b, T z, const Policy& pol) x = 0.5; else x = (1 + eta * sqrt((1 + c) / eta_2)) / 2; - - BOOST_MATH_ASSERT(x >= 0); - BOOST_MATH_ASSERT(x <= 1); + // + // These are post-conditions of the method, but the addition above + // may result in us being out by 1ulp either side of the boundary, + // so just check that we're in bounds and adjust as needed. + // See https://github.com/boostorg/math/issues/961 + // + if (x < 0) + x = 0; + else if (x > 1) + x = 1; + BOOST_MATH_ASSERT(eta * (x - 0.5) >= 0); #ifdef BOOST_INSTRUMENT std::cout << "Estimating x with Temme method 1: " << x << std::endl; @@ -901,6 +909,8 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py) if(x < lower) x = lower; } + std::uintmax_t max_iter = policies::get_max_root_iterations(); + std::uintmax_t max_iter_used = 0; // // Figure out how many digits to iterate towards: // @@ -923,10 +933,9 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py) // Now iterate, we can use either p or q as the target here // depending on which is smaller: // - std::uintmax_t max_iter = policies::get_max_root_iterations(); x = boost::math::tools::halley_iterate( boost::math::detail::ibeta_roots(a, b, (p < q ? p : q), (p < q ? false : true)), x, lower, upper, digits, max_iter); - policies::check_root_iterations("boost::math::ibeta<%1%>(%1%, %1%, %1%)", max_iter, pol); + policies::check_root_iterations("boost::math::ibeta<%1%>(%1%, %1%, %1%)", max_iter + max_iter_used, pol); // // We don't really want these asserts here, but they are useful for sanity // checking that we have the limits right, uncomment if you suspect bugs *only*. diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index 6493e81228..9fd7bee488 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -586,8 +586,8 @@ namespace detail { // we can jump way out of bounds if we're not careful. // See https://svn.boost.org/trac/boost/ticket/8314. delta = f0 / f1; - if (fabs(delta) > 2 * fabs(guess)) - delta = (delta < 0 ? -1 : 1) * 2 * fabs(guess); + if (fabs(delta) > 2 * fabs(result)) + delta = (delta < 0 ? -1 : 1) * 2 * fabs(result); } } else @@ -600,9 +600,19 @@ namespace detail { if ((convergence > 0.8) && (convergence < 2)) { // last two steps haven't converged. - delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2; - if ((result != 0) && (fabs(delta) > result)) - delta = sign(delta) * fabs(result) * 0.9f; // protect against huge jumps! + if (fabs(min) < 1 ? fabs(1000 * min) < fabs(max) : fabs(max / min) > 1000) + { + if(delta > 0) + delta = bracket_root_towards_min(f, result, f0, min, max, count); + else + delta = bracket_root_towards_max(f, result, f0, min, max, count); + } + else + { + delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2; + if ((result != 0) && (fabs(delta) > result)) + delta = sign(delta) * fabs(result) * 0.9f; // protect against huge jumps! + } // reset delta2 so that this branch will *not* be taken on the // next iteration: delta2 = delta * 3; @@ -641,6 +651,8 @@ namespace detail { result = guess - delta; if (result <= min) result = float_next(min); + if (result >= max) + result = float_prior(max); guess = min; continue; } @@ -669,6 +681,8 @@ namespace detail { result = guess - delta; if (result >= max) result = float_prior(max); + if (result <= min) + result = float_next(min); guess = min; continue; } diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 8aa2bdcff0..daff6520e5 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -168,6 +168,7 @@ test-suite special_fun : [ run git_issue_705.cpp ../../test/build//boost_unit_test_framework ] [ run git_issue_810.cpp ../../test/build//boost_unit_test_framework ] [ run git_issue_826.cpp ../../test/build//boost_unit_test_framework ] + [ run git_issue_961.cpp ] [ run special_functions_test.cpp ../../test/build//boost_unit_test_framework ] [ run test_airy.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] [ run test_bessel_j.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] diff --git a/test/git_issue_961.cpp b/test/git_issue_961.cpp new file mode 100644 index 0000000000..7594dcd6fc --- /dev/null +++ b/test/git_issue_961.cpp @@ -0,0 +1,76 @@ +// (C) Copyright John Maddock 2023. +// 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 +#include +#include +#include + +using namespace std; +using namespace boost::math; + +void show_fp_exception_flags() +{ + if (std::fetestexcept(FE_DIVBYZERO)) { + cout << " FE_DIVBYZERO"; + } + // FE_INEXACT is common and not interesting. + // if (std::fetestexcept(FE_INEXACT)) { + // cout << " FE_INEXACT"; + // } + if (std::fetestexcept(FE_INVALID)) { + cout << " FE_INVALID"; + } + if (std::fetestexcept(FE_OVERFLOW)) { + cout << " FE_OVERFLOW"; + } + if (std::fetestexcept(FE_UNDERFLOW)) { + cout << " FE_UNDERFLOW"; + } + cout << endl; +} + +template +int test() +{ + double a = 14.208308325339239; + double b = a; + + double p = 6.4898872103239473e-300; // Throws exception: Assertion `x >= 0' failed. + // double p = 7.8e-307; // No flags set, returns 8.57354094063444939e-23 + // double p = 7.7e-307; // FE_UNDERFLOW set, returns 0.0 + + while (p > (std::numeric_limits::min)()) + { + std::feclearexcept(FE_ALL_EXCEPT); + + try { + + double x = ibeta_inv(a, b, p, Policy()); + + show_fp_exception_flags(); + + std::cout << std::scientific << std::setw(24) + << std::setprecision(17) << x << std::endl; + } + catch (const std::exception& e) + { + std::cout << e.what() << std::endl; + return 1; + } + + p /= 1.25; + } + + return 0; +} + +int main(int argc, char* argv[]) +{ + using namespace boost::math::policies; + if (test>()) + return 1; + return test>>(); +}