From 0909d7cc9c3d06bcc6dadebe598bbfecbe4d16ab Mon Sep 17 00:00:00 2001 From: Ryan Date: Tue, 25 Jul 2023 17:15:56 -0400 Subject: [PATCH] only solve problems if have enough precision --- .../detail/ibeta_inverse.hpp | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/include/boost/math/special_functions/detail/ibeta_inverse.hpp b/include/boost/math/special_functions/detail/ibeta_inverse.hpp index be4bd95399..bb6040e073 100644 --- a/include/boost/math/special_functions/detail/ibeta_inverse.hpp +++ b/include/boost/math/special_functions/detail/ibeta_inverse.hpp @@ -25,23 +25,16 @@ namespace boost{ namespace math{ namespace detail{ template struct temme_root_finder { - temme_root_finder(const T t_, const T a_) : t(t_), a(a_) {} + temme_root_finder(const T t_, const T a_) : t(t_), a(a_) { + const T x_extrema = 1 / (1 + a); + BOOST_MATH_ASSERT(0 < x_extrema && x_extrema < 1); + } boost::math::tuple operator()(T x) { BOOST_MATH_STD_USING // ADL of std names T y = 1 - x; - if(y == 0) - { - T big = tools::max_value() / 4; - return boost::math::make_tuple(static_cast(-big), static_cast(-big)); - } - if(x == 0) - { - T big = tools::max_value() / 4; - return boost::math::make_tuple(static_cast(-big), big); - } T f = log(x) + a * log(y) + t; T f1 = (1 / x) - (a / (y)); return boost::math::make_tuple(f, f1); @@ -410,6 +403,10 @@ T temme_method_3_ibeta_inverse(T a, T b, T p, T q, const Policy& pol) T lower = eta < mu ? cross : 0; T upper = eta < mu ? 1 : cross; T x = (lower + upper) / 2; + + // Early exit for cases with numerical precision issues. + if (cross == 0 || cross == 1) { return cross; } + x = tools::newton_raphson_iterate( temme_root_finder(u, mu), x, lower, upper, policies::digits() / 2); #ifdef BOOST_INSTRUMENT