From 89f8dd36b098c49831fdec39fcf9b0f27528a0d0 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 28 Dec 2025 10:33:46 -0800 Subject: [PATCH 1/9] Initial implementation of [skip ci] --- .../math/distributions/non_central_f.hpp | 47 +++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index dedd437144..39fdcf31f9 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -404,7 +404,54 @@ namespace boost Policy()); return (x / (1 - x)) * (c.dist.degrees_of_freedom2() / c.dist.degrees_of_freedom1()); } // quantile complement. + + /* Need to rewrite this according to `find_non_centrality` in + non_central_chi_squared.hpp */ + template + struct non_centrality_finder + { + non_centrality_finder(const RealType x_, const RealType dfn_, const RealType dfd_, const RealType p_, bool c) + : x(x_), dfn(dfn_), dfd(dfd_), p(p_), comp(c) {} + RealType operator()(RealType nc) const + { + non_central_f_distribution d(dfn, dfd, nc); + return comp ? + RealType(p - cdf(complement(d, x))) + : RealType(cdf(d, x) - p); + } + private: + RealType x, dfn, dfd, p; + bool comp; + }; + + template + RealType find_non_centrality(const RealType dfn, const RealType dfd, const RealType p, const RealType x, const Policy& pol) + { + constexpr auto function = "non_central_f<%1%>::find_non_centrality"; + + if ( p =< 0 || p >= 1) { + return policies::raise_domain_error(function, "Can't find non centrality parameter when the probability is <=0 or >=1, only possible answer is %1%", // LCOV_EXCL_LINE + RealType(boost::math::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE + } + + RealType guess = RealType(10); // Starting guess. + RealType factor = 8; // How big steps to take when searching. + boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); + tools::eps_tolerance tol(policies::digits()); + + std::pair result_bracket = tools::bracket_and_solve_root( + non_centrality_finder(x, dfn, dfd, p), guess, factor, + false, tol, max_iter, pol); + + if (max_iter >= policies::get_max_root_iterations()) { + return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE + " or there is no answer to problem. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE + } + + RealType result = result_bracket.first + (result_bracket.second - result_bracket.first)/2; + return result; + } } // namespace math } // namespace boost From ed3289df1a4fcba96cb49f675b117545ce2514c3 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Tue, 30 Dec 2025 20:13:40 -0800 Subject: [PATCH 2/9] Added tests and changed namescapes --- .../math/distributions/non_central_f.hpp | 118 +++++++++++------- test/test_nc_f.cpp | 2 + 2 files changed, 72 insertions(+), 48 deletions(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index 39fdcf31f9..5013c24082 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -22,6 +22,56 @@ namespace boost { namespace math { + namespace detail + { + /* Need to rewrite this according to `find_non_centrality` in + non_central_chi_squared.hpp */ + template + struct non_centrality_finder_f + { + non_centrality_finder_f(const RealType x_, const RealType dfn_, const RealType dfd_, const RealType p_, bool c) + : x(x_), dfn(dfn_), dfd(dfd_), p(p_), comp(c) {} + + RealType operator()(RealType nc) const + { + non_central_f_distribution d(dfn, dfd, nc); + return comp ? + RealType(p - cdf(complement(d, x))) + : RealType(cdf(d, x) - p); + } + private: + RealType x, dfn, dfd, p; + bool comp; + }; + + template + inline RealType find_non_centrality_f(const RealType dfn, const RealType dfd, const RealType p, const RealType x, const Policy& pol) + { + constexpr auto function = "non_central_f<%1%>::find_non_centrality"; + + if ( p <= 0 || p >= 1) { + return policies::raise_domain_error(function, "Can't find non centrality parameter when the probability is <=0 or >=1, only possible answer is %1%", // LCOV_EXCL_LINE + RealType(boost::math::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE + } + + RealType guess = RealType(10); // Starting guess. + RealType factor = 1; // How big steps to take when searching. + boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); + tools::eps_tolerance tol(policies::digits()); + + std::pair result_bracket = tools::bracket_and_solve_root( + detail::non_centrality_finder_f(x, dfn, dfd, p, false), guess, factor, + false, tol, max_iter, pol); + + RealType result = result_bracket.first + (result_bracket.second - result_bracket.first)/2; + if (max_iter >= policies::get_max_root_iterations()) { + return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE + " or there is no answer to problem. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE + } + return result; + } + } // namespace detail + template > class non_central_f_distribution { @@ -58,6 +108,26 @@ namespace boost { // Private data getter function. return ncp; } + static RealType find_non_centrality(const RealType dfn, const RealType dfd, const RealType p, const RealType x) + { + constexpr auto function = "non_central_f_distribution<%1%>::find_non_centrality"; + typedef typename policies::evaluation::type eval_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + eval_type result = detail::find_non_centrality_f( + static_cast(dfn), + static_cast(dfd), + static_cast(p), + static_cast(x), + forwarding_policy()); + return policies::checked_narrowing_cast( + result, + function); + } private: // Data member, initialized by constructor. RealType v1; // alpha. @@ -404,54 +474,6 @@ namespace boost Policy()); return (x / (1 - x)) * (c.dist.degrees_of_freedom2() / c.dist.degrees_of_freedom1()); } // quantile complement. - - /* Need to rewrite this according to `find_non_centrality` in - non_central_chi_squared.hpp */ - template - struct non_centrality_finder - { - non_centrality_finder(const RealType x_, const RealType dfn_, const RealType dfd_, const RealType p_, bool c) - : x(x_), dfn(dfn_), dfd(dfd_), p(p_), comp(c) {} - - RealType operator()(RealType nc) const - { - non_central_f_distribution d(dfn, dfd, nc); - return comp ? - RealType(p - cdf(complement(d, x))) - : RealType(cdf(d, x) - p); - } - private: - RealType x, dfn, dfd, p; - bool comp; - }; - - template - RealType find_non_centrality(const RealType dfn, const RealType dfd, const RealType p, const RealType x, const Policy& pol) - { - constexpr auto function = "non_central_f<%1%>::find_non_centrality"; - - if ( p =< 0 || p >= 1) { - return policies::raise_domain_error(function, "Can't find non centrality parameter when the probability is <=0 or >=1, only possible answer is %1%", // LCOV_EXCL_LINE - RealType(boost::math::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE - } - - RealType guess = RealType(10); // Starting guess. - RealType factor = 8; // How big steps to take when searching. - boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); - tools::eps_tolerance tol(policies::digits()); - - std::pair result_bracket = tools::bracket_and_solve_root( - non_centrality_finder(x, dfn, dfd, p), guess, factor, - false, tol, max_iter, pol); - - if (max_iter >= policies::get_max_root_iterations()) { - return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE - " or there is no answer to problem. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE - } - - RealType result = result_bracket.first + (result_bracket.second - result_bracket.first)/2; - return result; - } } // namespace math } // namespace boost diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index 3a582a29a1..22cc667d5d 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -141,6 +141,8 @@ void test_spot( quantile(dist, P), x, tol * 10); BOOST_CHECK_CLOSE( quantile(complement(dist, Q)), x, tol * 10); + BOOST_CHECK_CLOSE( + dist.find_non_centrality(a, b, P, x), ncp, tol * 10); } if(boost::math::tools::digits() > 50) { From 55e2b17b020e236f6d4b3bac541c08ae5d715aee Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Fri, 2 Jan 2026 10:29:52 -0800 Subject: [PATCH 3/9] Added boost tuple support --- .../math/distributions/non_central_f.hpp | 24 +++++++++++++++++-- test/test_nc_f.cpp | 2 ++ 2 files changed, 24 insertions(+), 2 deletions(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index 5013c24082..63736a2b5e 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -17,6 +17,7 @@ #include #include #include +#include // complements namespace boost { @@ -24,8 +25,6 @@ namespace boost { namespace detail { - /* Need to rewrite this according to `find_non_centrality` in - non_central_chi_squared.hpp */ template struct non_centrality_finder_f { @@ -128,6 +127,27 @@ namespace boost result, function); } + template + BOOST_MATH_GPU_ENABLED static RealType find_non_centrality(const complemented4_type& c) + { + constexpr auto function = "non_central_f_distribution<%1%>::find_non_centrality"; + typedef typename policies::evaluation::type eval_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + eval_type result = detail::find_non_centrality_f( + static_cast(c.dist), + static_cast(c.param1), + static_cast(c.param2), + static_cast(c.param3), + forwarding_policy()); + return policies::checked_narrowing_cast( + result, + function); + } private: // Data member, initialized by constructor. RealType v1; // alpha. diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index 22cc667d5d..064e9eb889 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -143,6 +143,8 @@ void test_spot( quantile(complement(dist, Q)), x, tol * 10); BOOST_CHECK_CLOSE( dist.find_non_centrality(a, b, P, x), ncp, tol * 10); + BOOST_CHECK_CLOSE( + dist.find_non_centrality(boost::math::complement(a, b, P, x)), ncp, tol * 10); } if(boost::math::tools::digits() > 50) { From b330dd7e0cfa2f661373346cc599c3d25200df12 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Fri, 2 Jan 2026 11:09:54 -0800 Subject: [PATCH 4/9] Updated code to match non_central_chi_squared syntax --- .../boost/math/distributions/non_central_f.hpp | 16 +++++++++------- test/test_nc_f.cpp | 4 ++-- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index 63736a2b5e..8bad55c455 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -44,23 +44,23 @@ namespace boost }; template - inline RealType find_non_centrality_f(const RealType dfn, const RealType dfd, const RealType p, const RealType x, const Policy& pol) + inline RealType find_non_centrality_f(RealType x, const RealType dfn, const RealType dfd, const RealType p, const RealType q, const Policy& pol) { constexpr auto function = "non_central_f<%1%>::find_non_centrality"; - if ( p <= 0 || p >= 1) { + if ( p == 0 || q == 0) { return policies::raise_domain_error(function, "Can't find non centrality parameter when the probability is <=0 or >=1, only possible answer is %1%", // LCOV_EXCL_LINE RealType(boost::math::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } - + + non_centrality_finder_f f(x, dfn, dfd, p < q ? p : q, p < q ? false : true); RealType guess = RealType(10); // Starting guess. RealType factor = 1; // How big steps to take when searching. boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); tools::eps_tolerance tol(policies::digits()); std::pair result_bracket = tools::bracket_and_solve_root( - detail::non_centrality_finder_f(x, dfn, dfd, p, false), guess, factor, - false, tol, max_iter, pol); + f, guess, factor, false, tol, max_iter, pol); RealType result = result_bracket.first + (result_bracket.second - result_bracket.first)/2; if (max_iter >= policies::get_max_root_iterations()) { @@ -107,7 +107,7 @@ namespace boost { // Private data getter function. return ncp; } - static RealType find_non_centrality(const RealType dfn, const RealType dfd, const RealType p, const RealType x) + static RealType find_non_centrality(const RealType x, const RealType dfn, const RealType dfd, const RealType p) { constexpr auto function = "non_central_f_distribution<%1%>::find_non_centrality"; typedef typename policies::evaluation::type eval_type; @@ -118,10 +118,11 @@ namespace boost policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; eval_type result = detail::find_non_centrality_f( + static_cast(x), static_cast(dfn), static_cast(dfd), static_cast(p), - static_cast(x), + static_cast(1-p), forwarding_policy()); return policies::checked_narrowing_cast( result, @@ -143,6 +144,7 @@ namespace boost static_cast(c.param1), static_cast(c.param2), static_cast(c.param3), + static_cast(1-c.param3), forwarding_policy()); return policies::checked_narrowing_cast( result, diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index 064e9eb889..6000a051ce 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -142,9 +142,9 @@ void test_spot( BOOST_CHECK_CLOSE( quantile(complement(dist, Q)), x, tol * 10); BOOST_CHECK_CLOSE( - dist.find_non_centrality(a, b, P, x), ncp, tol * 10); + dist.find_non_centrality(x, a, b, P), ncp, tol * 10); BOOST_CHECK_CLOSE( - dist.find_non_centrality(boost::math::complement(a, b, P, x)), ncp, tol * 10); + dist.find_non_centrality(boost::math::complement(x, a, b, P)), ncp, tol * 10); } if(boost::math::tools::digits() > 50) { From 2464bcc8e0c0b0bbdefd8e7ae2e96e1a828af90e Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Fri, 2 Jan 2026 12:04:39 -0800 Subject: [PATCH 5/9] Added documentation; changed parameter naming --- doc/distributions/nc_f.qbk | 19 ++++++++++++++++ .../math/distributions/non_central_f.hpp | 22 +++++++++---------- test/test_nc_f.cpp | 2 +- 3 files changed, 31 insertions(+), 12 deletions(-) diff --git a/doc/distributions/nc_f.qbk b/doc/distributions/nc_f.qbk index d31c8116bc..8611f5c674 100644 --- a/doc/distributions/nc_f.qbk +++ b/doc/distributions/nc_f.qbk @@ -26,6 +26,11 @@ // Accessor to non-centrality parameter lambda: BOOST_MATH_GPU_ENABLED RealType non_centrality()const; + + // Parameter finders: + static RealType find_non_centrality(const RealType x, const RealType v1, const RealType v2, const RealType p); + template + static RealType find_non_centrality(const complemented3_type& c); }; }} // namespaces @@ -53,6 +58,20 @@ for different values of [lambda]: [graph nc_f_pdf] + BOOST_MATH_GPU_ENABLED static RealType find_non_centrality(const RealType x, const RealType v1, const RealType v2, const RealType p); + +This function returns the non centrality parameter /lambda/ such that: + +`cdf(non_central_chi_squared(v1, v2, lambda), x) == p` + + template + BOOST_MATH_GPU_ENABLED static RealType find_non_centrality(const complemented4_type& c); + +When called with argument `boost::math::complement(x, v1, v2, q)` +this function returns the non centrality parameter /lambda/ such that: + +`cdf(complement(non_central_chi_squared(v1, v2, lambda), x)) == q`. + [h4 Member Functions] BOOST_MATH_GPU_ENABLED non_central_f_distribution(RealType v1, RealType v2, RealType lambda); diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index 8bad55c455..dfd9baf0d6 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -28,23 +28,23 @@ namespace boost template struct non_centrality_finder_f { - non_centrality_finder_f(const RealType x_, const RealType dfn_, const RealType dfd_, const RealType p_, bool c) - : x(x_), dfn(dfn_), dfd(dfd_), p(p_), comp(c) {} + non_centrality_finder_f(const RealType x_, const RealType v1_, const RealType v2_, const RealType p_, bool c) + : x(x_), v1(v1_), v2(v2_), p(p_), comp(c) {} RealType operator()(RealType nc) const { - non_central_f_distribution d(dfn, dfd, nc); + non_central_f_distribution d(v1, v2, nc); return comp ? RealType(p - cdf(complement(d, x))) : RealType(cdf(d, x) - p); } private: - RealType x, dfn, dfd, p; + RealType x, v1, v2, p; bool comp; }; template - inline RealType find_non_centrality_f(RealType x, const RealType dfn, const RealType dfd, const RealType p, const RealType q, const Policy& pol) + inline 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"; @@ -53,7 +53,7 @@ namespace boost RealType(boost::math::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } - non_centrality_finder_f f(x, dfn, dfd, p < q ? p : q, p < q ? false : true); + non_centrality_finder_f f(x, v1, v2, p < q ? p : q, p < q ? false : true); RealType guess = RealType(10); // Starting guess. RealType factor = 1; // How big steps to take when searching. boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); @@ -107,7 +107,7 @@ namespace boost { // Private data getter function. return ncp; } - static RealType find_non_centrality(const RealType x, const RealType dfn, const RealType dfd, const RealType p) + static RealType find_non_centrality(const RealType x, const RealType v1, const RealType v2, const RealType p) { constexpr auto function = "non_central_f_distribution<%1%>::find_non_centrality"; typedef typename policies::evaluation::type eval_type; @@ -119,8 +119,8 @@ namespace boost policies::assert_undefined<> >::type forwarding_policy; eval_type result = detail::find_non_centrality_f( static_cast(x), - static_cast(dfn), - static_cast(dfd), + static_cast(v1), + static_cast(v2), static_cast(p), static_cast(1-p), forwarding_policy()); @@ -129,7 +129,7 @@ namespace boost function); } template - BOOST_MATH_GPU_ENABLED static RealType find_non_centrality(const complemented4_type& c) + static RealType find_non_centrality(const complemented4_type& c) { constexpr auto function = "non_central_f_distribution<%1%>::find_non_centrality"; typedef typename policies::evaluation::type eval_type; @@ -143,8 +143,8 @@ namespace boost static_cast(c.dist), static_cast(c.param1), static_cast(c.param2), - static_cast(c.param3), static_cast(1-c.param3), + static_cast(c.param3), forwarding_policy()); return policies::checked_narrowing_cast( result, diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index 6000a051ce..fdd03ab287 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -144,7 +144,7 @@ void test_spot( BOOST_CHECK_CLOSE( dist.find_non_centrality(x, a, b, P), ncp, tol * 10); BOOST_CHECK_CLOSE( - dist.find_non_centrality(boost::math::complement(x, a, b, P)), ncp, tol * 10); + dist.find_non_centrality(boost::math::complement(x, a, b, Q)), ncp, tol * 10); } if(boost::math::tools::digits() > 50) { From b794ab7de91e2222d9cff390b961ffd676a52dce Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Wed, 21 Jan 2026 09:41:24 -0500 Subject: [PATCH 6/9] Add GPU markers and change std::pair to boost::math::pair --- include/boost/math/distributions/non_central_f.hpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index dfd9baf0d6..2acaab2555 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -28,10 +28,10 @@ namespace boost template struct non_centrality_finder_f { - non_centrality_finder_f(const RealType x_, const RealType v1_, const RealType v2_, const RealType p_, bool c) + BOOST_MATH_GPU_ENABLED non_centrality_finder_f(const RealType x_, const RealType v1_, const RealType v2_, const RealType p_, const bool c) : x(x_), v1(v1_), v2(v2_), p(p_), comp(c) {} - RealType operator()(RealType nc) const + BOOST_MATH_GPU_ENABLED RealType operator()(RealType nc) const { non_central_f_distribution d(v1, v2, nc); return comp ? @@ -44,7 +44,7 @@ namespace boost }; template - inline 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 Policy& pol) { constexpr auto function = "non_central_f<%1%>::find_non_centrality"; @@ -59,7 +59,7 @@ namespace boost boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); tools::eps_tolerance tol(policies::digits()); - std::pair result_bracket = tools::bracket_and_solve_root( + boost::math::pair result_bracket = tools::bracket_and_solve_root( f, guess, factor, false, tol, max_iter, pol); RealType result = result_bracket.first + (result_bracket.second - result_bracket.first)/2; @@ -107,7 +107,7 @@ namespace boost { // Private data getter function. return ncp; } - static RealType find_non_centrality(const RealType x, const RealType v1, const RealType v2, const RealType p) + BOOST_MATH_GPU_ENABLED static RealType find_non_centrality(const RealType x, const RealType v1, const RealType v2, const RealType p) { constexpr auto function = "non_central_f_distribution<%1%>::find_non_centrality"; typedef typename policies::evaluation::type eval_type; @@ -129,7 +129,7 @@ namespace boost function); } template - static RealType find_non_centrality(const complemented4_type& c) + BOOST_MATH_GPU_ENABLED static RealType find_non_centrality(const complemented4_type& c) { constexpr auto function = "non_central_f_distribution<%1%>::find_non_centrality"; typedef typename policies::evaluation::type eval_type; From b223789a2f9e9ad124ee8a3ad6225e9b965aa541 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Wed, 21 Jan 2026 09:41:38 -0500 Subject: [PATCH 7/9] Add CUDA tests --- test/cuda_jamfile | 2 + test/test_nc_f_finder_double.cu | 114 ++++++++++++++++++++++++++++++++ test/test_nc_f_finder_float.cu | 114 ++++++++++++++++++++++++++++++++ 3 files changed, 230 insertions(+) create mode 100644 test/test_nc_f_finder_double.cu create mode 100644 test/test_nc_f_finder_float.cu diff --git a/test/cuda_jamfile b/test/cuda_jamfile index 4082057fa5..47ced5f037 100644 --- a/test/cuda_jamfile +++ b/test/cuda_jamfile @@ -165,6 +165,8 @@ run test_nc_f_pdf_double.cu ; run test_nc_f_pdf_float.cu ; run test_nc_f_quan_double.cu ; run test_nc_f_quan_float.cu ; +run test_nc_f_finder_double.cu ; +run test_nc_f_finder_float.cu ; run test_nc_chi_squared_cdf_double.cu ; run test_nc_chi_squared_cdf_float.cu ; diff --git a/test/test_nc_f_finder_double.cu b/test/test_nc_f_finder_double.cu new file mode 100644 index 0000000000..779860d9d1 --- /dev/null +++ b/test/test_nc_f_finder_double.cu @@ -0,0 +1,114 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// 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) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#include +#include +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef double float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + boost::math::non_central_f_distribution dist(5, 5, 5); + out[i] = dist.find_non_centrality(5, 5, 5, in1[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + try{ + + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + boost::random::mt19937 gen; + boost::random::uniform_real_distribution dist(0.1, 0.9); + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = dist(gen); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 256; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + cuda_test<<>>(input_vector1.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + std::cout << "CUDA kernal done in " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + if (err != cudaSuccess) + { + std::cerr << "Failed to launch non_central_f distribution kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + + boost::math::non_central_f_distribution non_central_dist(5, 5, 5); + for(int i = 0; i < numElements; ++i) + { + results.push_back(non_central_dist.find_non_centrality(5, 5, 5, input_vector1[i])); + } + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 100.0) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + std::cerr << "Error rate was: " << boost::math::epsilon_difference(output_vector[i], results[i]) << "eps" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED with calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + } + return 0; +} diff --git a/test/test_nc_f_finder_float.cu b/test/test_nc_f_finder_float.cu new file mode 100644 index 0000000000..73a97c872c --- /dev/null +++ b/test/test_nc_f_finder_float.cu @@ -0,0 +1,114 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// 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) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#include +#include +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef float float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + boost::math::non_central_f_distribution dist(5, 5, 5); + out[i] = dist.find_non_centrality(5, 5, 5, in1[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + try{ + + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + boost::random::mt19937 gen; + boost::random::uniform_real_distribution dist(0.1, 0.9); + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = dist(gen); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 256; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + cuda_test<<>>(input_vector1.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + std::cout << "CUDA kernal done in " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + if (err != cudaSuccess) + { + std::cerr << "Failed to launch non_central_f distribution kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + + boost::math::non_central_f_distribution non_central_dist(5, 5, 5); + for(int i = 0; i < numElements; ++i) + { + results.push_back(non_central_dist.find_non_centrality(5, 5, 5, input_vector1[i])); + } + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 100.0) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + std::cerr << "Error rate was: " << boost::math::epsilon_difference(output_vector[i], results[i]) << "eps" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED with calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + } + return 0; +} From e51d71a474c12d8e55989ed3a4b297b2136a5514 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Wed, 21 Jan 2026 09:43:24 -0500 Subject: [PATCH 8/9] Consistency of GPU markers in the docs --- doc/distributions/nc_f.qbk | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/distributions/nc_f.qbk b/doc/distributions/nc_f.qbk index 8611f5c674..4aa1fa95df 100644 --- a/doc/distributions/nc_f.qbk +++ b/doc/distributions/nc_f.qbk @@ -28,9 +28,9 @@ BOOST_MATH_GPU_ENABLED RealType non_centrality()const; // Parameter finders: - static RealType find_non_centrality(const RealType x, const RealType v1, const RealType v2, const RealType p); + BOOST_MATH_GPU_ENABLED static RealType find_non_centrality(const RealType x, const RealType v1, const RealType v2, const RealType p); template - static RealType find_non_centrality(const complemented3_type& c); + BOOST_MATH_GPU_ENABLED static RealType find_non_centrality(const complemented3_type& c); }; }} // namespaces From 59a4e7a46d41bb9e0bdc4159272e50388815ff5b Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Wed, 21 Jan 2026 15:13:10 -0800 Subject: [PATCH 9/9] Added catch for edge case nc=0 and testing for this --- .../math/distributions/non_central_f.hpp | 13 ++++++++++++- test/test_nc_f.cpp | 19 +++++++++++++++++++ 2 files changed, 31 insertions(+), 1 deletion(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index 2acaab2555..02a75a28d6 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -18,6 +18,8 @@ #include #include #include // complements +#include +#include namespace boost { @@ -53,9 +55,18 @@ namespace boost RealType(boost::math::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } + // Check if nc = 0 (which is just the F-distribution) + // I first tried comparing to boost::math::tools::epsilon(), + // but this was never satisfied for floats or doubles. Only + // long doubles would correctly return 0. + fisher_f_distribution dist(v1, v2); + if (boost::math::relative_difference(pdf(dist, x), p) < 1e-7){ + return 0; + } + non_centrality_finder_f f(x, v1, v2, p < q ? p : q, p < q ? false : true); RealType guess = RealType(10); // Starting guess. - RealType factor = 1; // How big steps to take when searching. + RealType factor = RealType(2); // How big steps to take when searching. boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); tools::eps_tolerance tol(policies::digits()); diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index fdd03ab287..3ccb9e3b43 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -327,6 +327,25 @@ void test_spots(RealType) { BOOST_CHECK_CLOSE(cdf(boost::math::non_central_f_distribution(static_cast(1e-100L), 3.f, 1.5f), static_cast(1e100L)), static_cast(0.6118152873453990639132215575213809716459L), tolerance); } + + // Check find_non_centrality_f edge case handling + // Case when nc=0 + RealType a = 5; + RealType b = 2; + RealType nc = 0; + RealType x_vals[] = { 0.25, 1.25, 10, 100}; + boost::math::non_central_f_distribution dist_no_centrality(a, b, nc); + for (RealType x : x_vals) + { + RealType P = pdf(dist_no_centrality, x); + BOOST_CHECK_CLOSE(dist.find_non_centrality(x, a, b, P), static_cast(0), tolerance); + } + // Case when P=1 or P=0 + BOOST_MATH_CHECK_THROW(dist.find_non_centrality(x, a, b, 1), std::domain_error); + BOOST_MATH_CHECK_THROW(dist.find_non_centrality(x, a, b, 0), std::domain_error); + // Case when Q=1 or Q=0 + BOOST_MATH_CHECK_THROW(dist.find_non_centrality(boost::math::complement(x, a, b, 1)), std::domain_error); + BOOST_MATH_CHECK_THROW(dist.find_non_centrality(boost::math::complement(x, a, b, 0)), std::domain_error); } // template void test_spots(RealType) BOOST_AUTO_TEST_CASE( test_main )