From 2a9394a655ee94c4c6c680561cd1b64b59387166 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Thu, 12 Feb 2026 11:18:34 -0800 Subject: [PATCH] Reverted nc=0 condition and changed relaxed test conditions --- .../math/distributions/non_central_f.hpp | 21 +++++++++++-------- test/test_nc_f.cpp | 2 +- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index c0acc8225..0fd33e6da 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -46,8 +46,7 @@ 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 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 RealType p_q_precision, const Policy& pol) { constexpr auto function = "non_central_f<%1%>::find_non_centrality"; if ( p == 0 || q == 0) { @@ -55,14 +54,16 @@ namespace boost RealType(boost::math::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } - non_centrality_finder_f f(x, v1, v2, p < q ? p : q, p < q ? false : true); - // 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-6, 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-6){ + 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()) <= 20 * p_q_precision * p){ return 0; } @@ -135,6 +136,7 @@ namespace boost static_cast(v2), static_cast(p), static_cast(1-p), + static_cast(tools::epsilon()), forwarding_policy()); return policies::checked_narrowing_cast( result, @@ -157,6 +159,7 @@ 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, diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index 5a1cd4da6..6b7c2347a 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -336,7 +336,7 @@ void test_spots(RealType, const char* name = nullptr) for (RealType x : x_vals) { RealType P = cdf(dist_no_centrality, x); - BOOST_CHECK_LE(dist.find_non_centrality(x, a, b, P), boost::math::tools::min_value() * 10); + BOOST_CHECK_LE(dist.find_non_centrality(x, a, b, P), tolerance); } // Case when P=1 or P=0 BOOST_MATH_CHECK_THROW(dist.find_non_centrality(x, a, b, 1), std::domain_error);