From 0e48617447ee5c03023b3351f1dd2bf9d3973f61 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Tue, 10 Feb 2026 11:02:47 -0800 Subject: [PATCH] Reverted to checking the relative error for small nc --- .../math/distributions/non_central_f.hpp | 21 ++++++++----------- 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index a6df68b8f..60d5926a9 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -46,7 +46,8 @@ namespace boost }; template - BOOST_MATH_GPU_ENABLED RealType find_non_centrality_f(const RealType x, const RealType v1, const RealType v2, const RealType p, const RealType q, const RealType p_q_precision, const Policy& pol) { + BOOST_MATH_GPU_ENABLED RealType find_non_centrality_f(const RealType x, const RealType v1, const RealType v2, const RealType p, const RealType q, const Policy& pol) + { constexpr auto function = "non_central_f<%1%>::find_non_centrality"; if ( p == 0 || q == 0) { @@ -54,16 +55,14 @@ namespace boost RealType(boost::math::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } - // Check if nc = 0 (which is just the F-distribution) non_centrality_finder_f f(x, v1, v2, p < q ? p : q, p < q ? false : true); - // This occurs when the root finder would need to find a result smaller than - // tools::min_value (which it cannot do). Note that we have to add in a small - // amount of "tolerance" since the subtraction in our termination condition - // implies a small amount of wobble in the result which should be of the - // order p * eps (note not q * eps, since q is calculated as 1-p). - // Also note that p_q_precision is passed down from our caller as the - // epsilon of the original called values, and not after possible promotion. - if (f(tools::min_value()) <= 3 * p_q_precision * p){ + + // Check if nc = 0 (which is just the F-distribution) + // When the relative difference between the cdf with nc=0 + // and p drops below 1e-7, the function f becomes + // numerically unstable. Thus, we can't find nc below this point. + // See PR 1345 for more details. + if (abs(f(tools::min_value()) / p) <= 1e-7){ return 0; } @@ -136,7 +135,6 @@ namespace boost static_cast(v2), static_cast(p), static_cast(1-p), - static_cast(tools::epsilon()), forwarding_policy()); return policies::checked_narrowing_cast( result, @@ -159,7 +157,6 @@ namespace boost static_cast(c.param2), static_cast(1-c.param3), static_cast(c.param3), - static_cast(tools::epsilon()), forwarding_policy()); return policies::checked_narrowing_cast( result,