From 638bfc49dc5283c0290eb97f8dae4caa954a04af Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Tue, 17 Feb 2026 10:44:21 -0800 Subject: [PATCH 01/96] Added nu expansion --- include/boost/math/special_functions/beta.hpp | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index f6327fa220..ca7d6b6bed 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -1142,7 +1142,20 @@ BOOST_MATH_GPU_ENABLED T ibeta_large_ab(T a, T b, T x, T y, bool invert, bool no T x0 = a / (a + b); T y0 = b / (a + b); - T nu = x0 * log(x / x0) + y0 * log(y / y0); + + // Expand nu about x0 + T nu = 0; + for (int i=2; i<5; i++) + { + nu += pow(x-x0, i) / i * (pow(x0, -(i-1)) - pow(x0-1, -(i-1))) * pow(-1, i+1); + } + // Calculate the next term in the series + T remainder = pow(x-x0, 5) / 5 * (pow(x0, -4) - pow(x0-1, 4)); + + // If the remainder is large, then fall back to using the log formula + if (remainder >= tools::forth_root_epsilon()){ + nu = x0 * log(x / x0) + y0 * log(y / y0); + } // // Above compution is unstable, force nu to zero if // something went wrong: @@ -1181,7 +1194,7 @@ BOOST_MATH_GPU_ENABLED T ibeta_large_ab(T a, T b, T x, T y, bool invert, bool no T mul = 1; if (!normalised) mul = boost::math::beta(a, b, pol); - + T log_erf_remainder = -0.5 * log(2 * constants::pi() * (a+b)) + a * log(x / x0) + b * log((1-x) / (1-x0)) + log(abs(1/nu - sqrt(x0 * (1-x0)) / (x-x0))); return mul * ((invert ? (1 + boost::math::erf(-nu * sqrt((a + b) / 2), pol)) / 2 : boost::math::erfc(-nu * sqrt((a + b) / 2), pol) / 2)); } From 3c176f1ee3e7713dea36818c4f3f633084edf840 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Tue, 17 Feb 2026 11:42:06 -0800 Subject: [PATCH 02/96] Added try/catch statement --- include/boost/math/special_functions/beta.hpp | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index ca7d6b6bed..a277c65e8b 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -1194,7 +1194,7 @@ BOOST_MATH_GPU_ENABLED T ibeta_large_ab(T a, T b, T x, T y, bool invert, bool no T mul = 1; if (!normalised) mul = boost::math::beta(a, b, pol); - T log_erf_remainder = -0.5 * log(2 * constants::pi() * (a+b)) + a * log(x / x0) + b * log((1-x) / (1-x0)) + log(abs(1/nu - sqrt(x0 * (1-x0)) / (x-x0))); + // T log_erf_remainder = -0.5 * log(2 * constants::pi() * (a+b)) + a * log(x / x0) + b * log((1-x) / (1-x0)) + log(abs(1/nu - sqrt(x0 * (1-x0)) / (x-x0))); return mul * ((invert ? (1 + boost::math::erf(-nu * sqrt((a + b) / 2), pol)) / 2 : boost::math::erfc(-nu * sqrt((a + b) / 2), pol) / 2)); } @@ -1612,8 +1612,17 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b invert = false; } else - fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative); - + { + try + { + fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative); + } + // If series converges slowly, fall back to erf approximation + catch (boost::math::evaluation_error& e) + { + fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol); + } + } BOOST_MATH_INSTRUMENT_VARIABLE(fract); } } From 8b5b1fdfabe38bf9d41507d5178b417aa4236ac6 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Wed, 18 Feb 2026 15:41:12 -0800 Subject: [PATCH 03/96] Accepted review comments to copy code from ibeta_fraction2 --- include/boost/math/special_functions/beta.hpp | 58 +++++++++---------- 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index a277c65e8b..6940ad929c 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -1584,46 +1584,46 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b else { // a and b both large: - bool use_asym = false; T ma = BOOST_MATH_GPU_SAFE_MAX(a, b); - T xa = ma == a ? x : y; T saddle = ma / (a + b); - T powers = 0; - if ((ma > 1e-5f / tools::epsilon()) && (ma / BOOST_MATH_GPU_SAFE_MIN(a, b) < (xa < saddle ? 2 : 15))) - { - if (a == b) - use_asym = true; - else - { - powers = exp(log(x / (a / (a + b))) * a + log(y / (b / (a + b))) * b); - if (powers < tools::epsilon()) - use_asym = true; - } - } - if(use_asym) + // If ma large and x is close to saddle, use `erf` approximation in `ibeta_large_ab`. + // In this case we don't need to invert + if ((ma > 1e-5f / tools::epsilon()) && (fabs(saddle - x) < 1e-12) && (1 - saddle > 0.125)) { fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol); - if (fract * tools::epsilon() < powers) - { - // Erf approximation failed, correction term is too large, fall back: - fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative); - } - else - invert = false; + invert = false; } else { - try + // This is the same logic that is in `ibeta_fraction2`. We have repeated + // the implementation here because when x is close to saddle, the continued + // fraction will not converge. If this occurs, we fall back to the `erf` + // approximation. Simply using `ibeta_fraction2` will cause an evaluation + // error to be thrown and we won't be able to use the `erf` approximation. + typedef typename lanczos::lanczos::type lanczos_type; + T local_result = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol); + if (p_derivative) { - fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative); + *p_derivative = local_result; + BOOST_MATH_ASSERT(*p_derivative >= 0); } - // If series converges slowly, fall back to erf approximation - catch (boost::math::evaluation_error& e) + if (local_result != 0) { - fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol); + ibeta_fraction2_t f(a, b, x, y); + boost::math::uintmax_t max_terms = boost::math::policies::get_max_series_iterations(); + T local_fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon(), max_terms); + if (max_terms >= boost::math::policies::get_max_series_iterations()) + { + // Continued fraction failed, fall back to asymptotic expansion: + fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol); + invert = false; + } + else + fract = local_result / local_fract; } - } - BOOST_MATH_INSTRUMENT_VARIABLE(fract); + else + fract = 0; + } } } if(p_derivative) From 12047b68d9d7db53e8c7605c45ab76c87efc04a3 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Thu, 19 Feb 2026 12:10:25 -0800 Subject: [PATCH 04/96] Added some sparse testing --- test/test_ibeta.hpp | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/test/test_ibeta.hpp b/test/test_ibeta.hpp index 2f0ea47107..00f82d8be7 100644 --- a/test/test_ibeta.hpp +++ b/test/test_ibeta.hpp @@ -487,5 +487,27 @@ void test_spots(T) } BOOST_CHECK_EQUAL(boost::math::ibeta(static_cast(2), static_cast(1), static_cast(0)), 0); BOOST_CHECK_EQUAL(boost::math::ibeta(static_cast(1), static_cast(2), static_cast(0)), 0); + + // Bug testing for large a,b and x close to a / (a+b). See PR 1363. + // The values for a,b are just too large for floats to handle. + if (!std::is_same::value) + { + T a_values[2] = {10000000272564224, 3.1622776601699636e+16}; + T b_values[2] = {9965820922822656, 3.130654883566682e+18}; + T delta[8] = {0, 1e-15, 1e-14, 1e-13, 1e-12, 1e-11, 1e-10, 1e-9}; + T a, b, x; + for (unsigned int i=0; i<2; i++){ + a = a_values[i]; + b = b_values[i]; + x = a / (a+b); // roughly the median of ibeta + + + for (unsigned int j=0; j < 7; j++) + { + BOOST_CHECK(boost::math::ibeta(a, b, x + delta[j+1]) > boost::math::ibeta(a, b, x + delta[j])); + BOOST_CHECK(boost::math::ibeta(a, b, x - delta[j]) > boost::math::ibeta(a, b, x - delta[j+1])); + } + } + } } From 1611bae3e86811c20c1c01e3cf69fd667a1a5b32 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Fri, 20 Feb 2026 10:21:41 -0800 Subject: [PATCH 05/96] Relaxed tests --- test/test_ibeta.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/test/test_ibeta.hpp b/test/test_ibeta.hpp index 00f82d8be7..7713c74c39 100644 --- a/test/test_ibeta.hpp +++ b/test/test_ibeta.hpp @@ -501,11 +501,10 @@ void test_spots(T) b = b_values[i]; x = a / (a+b); // roughly the median of ibeta - for (unsigned int j=0; j < 7; j++) { - BOOST_CHECK(boost::math::ibeta(a, b, x + delta[j+1]) > boost::math::ibeta(a, b, x + delta[j])); - BOOST_CHECK(boost::math::ibeta(a, b, x - delta[j]) > boost::math::ibeta(a, b, x - delta[j+1])); + BOOST_CHECK(boost::math::ibeta(a, b, x + delta[j+1]) >= boost::math::ibeta(a, b, x + delta[j])); + BOOST_CHECK(boost::math::ibeta(a, b, x - delta[j]) >= boost::math::ibeta(a, b, x - delta[j+1])); } } } From 44cbb533cc0a5cc49ccb236bdf630f16be29ccd2 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Fri, 20 Feb 2026 11:02:57 -0800 Subject: [PATCH 06/96] Initial implementation/test of find_v1 --- .../math/distributions/non_central_f.hpp | 106 +++++++++++++++++- test/test_nc_f.cpp | 12 +- 2 files changed, 113 insertions(+), 5 deletions(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index 783c321980..8efeb6216c 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -44,7 +44,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 RealType p_q_precision, const Policy& pol) + { constexpr auto function = "non_central_f<%1%>::find_non_centrality"; if ( p == 0 || q == 0) { @@ -80,6 +81,66 @@ namespace boost } return result; } + + // Find first degrees of freedom + template + struct f_v1_degrees_of_freedom_finder + { + f_v1_degrees_of_freedom_finder( + RealType x_, RealType v2_, RealType nc_, RealType p_, bool c) + : x(x_), v2(v2_), nc(nc_), p(p_), comp(c) {} + + RealType operator()(const RealType& v1) + { + non_central_f_distribution d(v1, v2, nc); + return comp ? + p - cdf(complement(d, x)) + : cdf(d, x) - p; + } + private: + RealType x; + RealType v2; + RealType nc; + RealType p; + bool comp; + }; + + template + inline RealType find_v1_degrees_of_freedom_f( + const RealType x, const RealType v2, const RealType nc, const RealType p, const RealType q, const Policy& pol) + { + using std::fabs; + const char* function = "non_central_f<%1%>::find_v1_degrees_of_freedom_f"; + if((p == 0) || (q == 0)) + { + // + // Can't find a thing if one of p and q is zero: + // + return policies::raise_evaluation_error(function, "Can't find v1 degrees of freedom when the probability is 0 or 1, only possible answer is %1%", + RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE + } + if (fabs(x) < tools::epsilon()) + { + return policies::raise_evaluation_error(function, "Can't find v1 degrees of freedom when the abscissa value is very close to zero as all degrees of freedom generate the same CDF at x=0: try again further out in the tails!!", + RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE + } + f_v1_degrees_of_freedom_finder f(x, v2, nc, p < q ? p : q, p < q ? false : true); + tools::eps_tolerance tol(policies::digits()); + std::uintmax_t max_iter = policies::get_max_root_iterations(); + // + // Pick an initial guess: + // + RealType guess = 200; + std::pair ir = tools::bracket_and_solve_root( + f, guess, RealType(2), x < 0 ? false : true, tol, max_iter, pol); + RealType result = ir.first + (ir.second - ir.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 > @@ -163,6 +224,49 @@ namespace boost result, function); } + BOOST_MATH_GPU_ENABLED static RealType find_v1(const RealType x, const RealType v2, const RealType nc, const RealType p) + { + constexpr auto function = "non_central_f_distribution<%1%>::find_v1"; + 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_v1_degrees_of_freedom_f( + static_cast(x), + static_cast(v2), + static_cast(nc), + static_cast(p), + static_cast(1-p), + forwarding_policy()); + return policies::checked_narrowing_cast( + result, + function); + } + template + BOOST_MATH_GPU_ENABLED static RealType find_v1(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_v1_degrees_of_freedom_f( + static_cast(c.dist), + static_cast(c.param1), + static_cast(c.param2), + static_cast(1-c.param3), + 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 6b7c2347ab..96664e4d9a 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -138,13 +138,17 @@ void test_spot( BOOST_CHECK_CLOSE( cdf(complement(dist, x)), Q, tol); BOOST_CHECK_CLOSE( - quantile(dist, P), x, tol * 10); + quantile(dist, P), x, tol * 10); BOOST_CHECK_CLOSE( - quantile(complement(dist, Q)), x, tol * 10); + quantile(complement(dist, Q)), x, tol * 10); BOOST_CHECK_CLOSE( - dist.find_non_centrality(x, a, b, P), ncp, tol * 10); + dist.find_non_centrality(x, a, b, P), ncp, tol * 10); BOOST_CHECK_CLOSE( - dist.find_non_centrality(boost::math::complement(x, a, b, Q)), ncp, tol * 10); + dist.find_non_centrality(boost::math::complement(x, a, b, Q)), ncp, tol * 10); + BOOST_CHECK_CLOSE( + dist.find_v1(x, b, ncp, P), a, tol * 10); + BOOST_CHECK_CLOSE( + dist.find_v1(boost::math::complement(x, b, ncp, Q)), a, tol * 10); } if(boost::math::tools::digits() > 50) { From 3eb7561296c5a15cc2c891a68e5c7f7a624d31e4 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Fri, 20 Feb 2026 16:53:26 -0800 Subject: [PATCH 07/96] Added v2 finder --- .../math/distributions/non_central_f.hpp | 83 +++++++++++++++---- test/test_nc_f.cpp | 4 + 2 files changed, 70 insertions(+), 17 deletions(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index 8efeb6216c..3112aeedb4 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -82,16 +82,17 @@ namespace boost return result; } - // Find first degrees of freedom template - struct f_v1_degrees_of_freedom_finder + struct f_degrees_of_freedom_finder { - f_v1_degrees_of_freedom_finder( - RealType x_, RealType v2_, RealType nc_, RealType p_, bool c) - : x(x_), v2(v2_), nc(nc_), p(p_), comp(c) {} + f_degrees_of_freedom_finder( + RealType x_, RealType v_, RealType nc_, bool find_v1_, RealType p_, bool c) + : x(x_), v(v_), nc(nc_), find_v1(find_v1_), p(p_), comp(c) {} - RealType operator()(const RealType& v1) + RealType operator()(const RealType& input_v) { + RealType v1 = find_v1 ? input_v : v; + RealType v2 = find_v1 ? v : input_v; non_central_f_distribution d(v1, v2, nc); return comp ? p - cdf(complement(d, x)) @@ -99,40 +100,41 @@ namespace boost } private: RealType x; - RealType v2; + RealType v; RealType nc; + bool find_v1; RealType p; bool comp; }; template - inline RealType find_v1_degrees_of_freedom_f( - const RealType x, const RealType v2, const RealType nc, const RealType p, const RealType q, const Policy& pol) + inline RealType find_degrees_of_freedom_f( + const RealType x, const RealType v, const RealType nc, const bool find_v1, const RealType p, const RealType q, const Policy& pol) { using std::fabs; - const char* function = "non_central_f<%1%>::find_v1_degrees_of_freedom_f"; + const char* function = "non_central_f<%1%>::find_degrees_of_freedom_f"; if((p == 0) || (q == 0)) { // // Can't find a thing if one of p and q is zero: // - return policies::raise_evaluation_error(function, "Can't find v1 degrees of freedom when the probability is 0 or 1, only possible answer is %1%", + return policies::raise_evaluation_error(function, "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%", RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } if (fabs(x) < tools::epsilon()) { - return policies::raise_evaluation_error(function, "Can't find v1 degrees of freedom when the abscissa value is very close to zero as all degrees of freedom generate the same CDF at x=0: try again further out in the tails!!", + return policies::raise_evaluation_error(function, "Can't find degrees of freedom when the abscissa value is very close to zero as all degrees of freedom generate the same CDF at x=0: try again further out in the tails!!", RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } - f_v1_degrees_of_freedom_finder f(x, v2, nc, p < q ? p : q, p < q ? false : true); + f_degrees_of_freedom_finder f(x, v, nc, find_v1, p < q ? p : q, p < q ? false : true); tools::eps_tolerance tol(policies::digits()); std::uintmax_t max_iter = policies::get_max_root_iterations(); // // Pick an initial guess: // - RealType guess = 200; + RealType guess = 1; std::pair ir = tools::bracket_and_solve_root( - f, guess, RealType(2), x < 0 ? false : true, tol, max_iter, pol); + f, guess, RealType(2), f(guess) < 0 ? true : false, tol, max_iter, pol); RealType result = ir.first + (ir.second - ir.first) / 2; if(max_iter >= policies::get_max_root_iterations()) { @@ -234,10 +236,11 @@ namespace boost policies::promote_double, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; - eval_type result = detail::find_v1_degrees_of_freedom_f( + eval_type result = detail::find_degrees_of_freedom_f( static_cast(x), static_cast(v2), static_cast(nc), + true, static_cast(p), static_cast(1-p), forwarding_policy()); @@ -256,10 +259,56 @@ namespace boost policies::promote_double, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; - eval_type result = detail::find_v1_degrees_of_freedom_f( + eval_type result = detail::find_degrees_of_freedom_f( static_cast(c.dist), static_cast(c.param1), static_cast(c.param2), + true, + static_cast(1-c.param3), + static_cast(c.param3), + forwarding_policy()); + return policies::checked_narrowing_cast( + result, + function); + } + BOOST_MATH_GPU_ENABLED static RealType find_v2(const RealType x, const RealType v2, const RealType nc, const RealType p) + { + constexpr auto function = "non_central_f_distribution<%1%>::find_v1"; + 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_degrees_of_freedom_f( + static_cast(x), + static_cast(v2), + static_cast(nc), + false, + static_cast(p), + static_cast(1-p), + forwarding_policy()); + return policies::checked_narrowing_cast( + result, + function); + } + template + BOOST_MATH_GPU_ENABLED static RealType find_v2(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_degrees_of_freedom_f( + static_cast(c.dist), + static_cast(c.param1), + static_cast(c.param2), + false, static_cast(1-c.param3), static_cast(c.param3), forwarding_policy()); diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index 96664e4d9a..5a61d31947 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -149,6 +149,10 @@ void test_spot( dist.find_v1(x, b, ncp, P), a, tol * 10); BOOST_CHECK_CLOSE( dist.find_v1(boost::math::complement(x, b, ncp, Q)), a, tol * 10); + BOOST_CHECK_CLOSE( + dist.find_v2(x, a, ncp, P), b, tol * 10); + BOOST_CHECK_CLOSE( + dist.find_v2(boost::math::complement(x, a, ncp, Q)), b, tol * 10); } if(boost::math::tools::digits() > 50) { From 48016543b9a0a9545266d7242c3f35c586bf1952 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sat, 21 Feb 2026 10:49:08 -0800 Subject: [PATCH 08/96] Included edge cases for code coverage [ci skip] --- include/boost/math/distributions/non_central_f.hpp | 4 ++-- test/test_nc_f.cpp | 12 ++++++++++++ 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index 3112aeedb4..5831f6622c 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -118,10 +118,10 @@ namespace boost // // Can't find a thing if one of p and q is zero: // - return policies::raise_evaluation_error(function, "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%", + return policies::raise_domain_error(function, "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%", RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } - if (fabs(x) < tools::epsilon()) + if (x < tools::epsilon()) { return policies::raise_evaluation_error(function, "Can't find degrees of freedom when the abscissa value is very close to zero as all degrees of freedom generate the same CDF at x=0: try again further out in the tails!!", RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index 5a61d31947..d96675320d 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -377,6 +377,18 @@ void test_spots(RealType, const char* name = nullptr) } } } + // Check find_v1/v2 edge cases + // Case when P=1 or P=0 + nc = 2; + BOOST_MATH_CHECK_THROW(dist.find_v1(x, b, nc, 1), std::domain_error); + BOOST_MATH_CHECK_THROW(dist.find_v1(x, b, nc, 0), std::domain_error); + // Case when Q=1 or Q=0 + BOOST_MATH_CHECK_THROW(dist.find_v1(boost::math::complement(x, b, nc, 1)), std::domain_error); + BOOST_MATH_CHECK_THROW(dist.find_v1(boost::math::complement(x, b, nc, 0)), std::domain_error); + // Check very small values of x an evaluation error is thrown + x = boost::math::tools::epsilon() / 10; + BOOST_MATH_CHECK_THROW(dist.find_v1(boost::math::complement(x, b, nc, 0.5)), boost::math::evaluation_error); + BOOST_MATH_CHECK_THROW(dist.find_v1(x, b, nc, 0.5), boost::math::evaluation_error); } // template void test_spots(RealType) BOOST_AUTO_TEST_CASE( test_main ) From 92dfc343c3f2f632e301f5b850b2f52618878747 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 22 Feb 2026 14:20:44 -0800 Subject: [PATCH 09/96] Added error message output [apple] --- test/test_ibeta.cpp | 8 ++++---- test/test_ibeta.hpp | 9 ++++++--- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/test/test_ibeta.cpp b/test/test_ibeta.cpp index 563844a06c..0059b46619 100644 --- a/test/test_ibeta.cpp +++ b/test/test_ibeta.cpp @@ -380,18 +380,18 @@ BOOST_AUTO_TEST_CASE( test_main ) gsl_set_error_handler_off(); #endif #ifdef TEST_FLOAT - test_spots(0.0F); + test_spots(0.0F, "float"); #endif #ifdef TEST_DOUBLE - test_spots(0.0); + test_spots(0.0, "double"); #endif #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS #ifdef TEST_LDOUBLE - test_spots(0.0L); + test_spots(0.0L, "long double"); #endif #if !BOOST_WORKAROUND(BOOST_BORLANDC, BOOST_TESTED_AT(0x582)) #if defined(TEST_REAL_CONCEPT) && !defined(BOOST_MATH_NO_REAL_CONCEPT_TESTS) - test_spots(boost::math::concepts::real_concept(0.1)); + test_spots(boost::math::concepts::real_concept(0.1), "real_concept"); #endif #endif #endif diff --git a/test/test_ibeta.hpp b/test/test_ibeta.hpp index 7713c74c39..b6e372ec01 100644 --- a/test/test_ibeta.hpp +++ b/test/test_ibeta.hpp @@ -152,12 +152,13 @@ void test_beta(T, const char* name) } template -void test_spots(T) +void test_spots(T, const char* name) { // // basic sanity checks, tolerance is 30 epsilon expressed as a percentage: // Spot values are from http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=BetaRegularized // using precision of 50 decimal digits. + std::cout << "Testing spot values with type " << name << std::endl; T tolerance = boost::math::tools::epsilon() * 3000; if (boost::math::tools::digits() > 100) tolerance *= 2; @@ -503,8 +504,10 @@ void test_spots(T) for (unsigned int j=0; j < 7; j++) { - BOOST_CHECK(boost::math::ibeta(a, b, x + delta[j+1]) >= boost::math::ibeta(a, b, x + delta[j])); - BOOST_CHECK(boost::math::ibeta(a, b, x - delta[j]) >= boost::math::ibeta(a, b, x - delta[j+1])); + BOOST_CHECK_MESSAGE(boost::math::ibeta(a, b, x + delta[j+1]) > boost::math::ibeta(a, b, x + delta[j]), + "ibeta not monotonically increasing above a/(a+b) for ibeta(" << a << ", " << b << ", " << x << "): [" << boost::math::ibeta(a, b, x - delta[j]) << " >= " << boost::math::ibeta(a, b, x - delta[j+1]) << "]"); + BOOST_CHECK_MESSAGE(boost::math::ibeta(a, b, x - delta[j]) > boost::math::ibeta(a, b, x - delta[j+1]), + "ibeta not monotonically increasing below a/(a+b) for ibeta(" << a << ", " << b << ", " << x << "): [" << boost::math::ibeta(a, b, x - delta[j]) << " >= " << boost::math::ibeta(a, b, x - delta[j+1]) << "]"); } } } From cbf53f88cce9516b9528cfa5a10b67108d62a606 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 22 Feb 2026 15:55:02 -0800 Subject: [PATCH 10/96] Fixed error message in test_ibeta [macos] --- test/test_ibeta.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_ibeta.hpp b/test/test_ibeta.hpp index b6e372ec01..7bef2d9715 100644 --- a/test/test_ibeta.hpp +++ b/test/test_ibeta.hpp @@ -505,9 +505,9 @@ void test_spots(T, const char* name) for (unsigned int j=0; j < 7; j++) { BOOST_CHECK_MESSAGE(boost::math::ibeta(a, b, x + delta[j+1]) > boost::math::ibeta(a, b, x + delta[j]), - "ibeta not monotonically increasing above a/(a+b) for ibeta(" << a << ", " << b << ", " << x << "): [" << boost::math::ibeta(a, b, x - delta[j]) << " >= " << boost::math::ibeta(a, b, x - delta[j+1]) << "]"); + "ibeta not monotonically increasing above a/(a+b) for ibeta(" << a << ", " << b << ", " << x << ") and delta=" << delta[j] << ": [" << boost::math::ibeta(a, b, x + delta[j+1]) << " >= " << boost::math::ibeta(a, b, x + delta[j]) << "]"); BOOST_CHECK_MESSAGE(boost::math::ibeta(a, b, x - delta[j]) > boost::math::ibeta(a, b, x - delta[j+1]), - "ibeta not monotonically increasing below a/(a+b) for ibeta(" << a << ", " << b << ", " << x << "): [" << boost::math::ibeta(a, b, x - delta[j]) << " >= " << boost::math::ibeta(a, b, x - delta[j+1]) << "]"); + "ibeta not monotonically increasing below a/(a+b) for ibeta(" << a << ", " << b << ", " << x << ") and delta=" << delta[j] << ": [" << boost::math::ibeta(a, b, x - delta[j]) << " >= " << boost::math::ibeta(a, b, x - delta[j+1]) << "]"); } } } From bffd9233241e3cf8b871e053621b273a7f44282f Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 1 Mar 2026 14:16:46 -0800 Subject: [PATCH 11/96] Added debugging print statements [apple] --- include/boost/math/special_functions/beta.hpp | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index 6940ad929c..78a1e96c06 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -1210,6 +1210,8 @@ BOOST_MATH_GPU_ENABLED T ibeta_large_ab(T a, T b, T x, T y, bool invert, bool no template BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_derivative) { + std::cout << std::setprecision(64) << "a= " << a << std::endl; + std::cout << "b= " << b << std::endl; constexpr auto function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)"; typedef typename lanczos::lanczos::type lanczos_type; BOOST_MATH_STD_USING // for ADL of std math functions. @@ -1592,6 +1594,7 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b { fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol); invert = false; + std::cout << "Using Erf approximation: " << fract << std::endl; } else { @@ -1617,9 +1620,12 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b // Continued fraction failed, fall back to asymptotic expansion: fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol); invert = false; + std::cout << "Failed continued fractions so falling back to erf" << std::endl; } - else + else{ fract = local_result / local_fract; + std::cout << "Using continued fractions with number of terms: " << max_terms << std::endl; + } } else fract = 0; @@ -1647,6 +1653,12 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b } } } + std::cout << "Fract is: " << fract << std::endl; + std::cout << "invert: " << invert << std::endl; + std::cout << "normalized: " << normalised << std::endl; + std::cout << "Beta(a,b): " << boost::math::beta(a, b, pol) << std::endl; + std::cout << "ibeta(a,b)=" << (invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract) << std::endl; + std::cout << "===========================" << std::endl; return invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract; } // template T ibeta_imp(T a, T b, T x, const Lanczos& l, bool inv, bool normalised) From 486fc145285afa22fc6030a3021cad0cc0ffa9d3 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 1 Mar 2026 14:55:31 -0800 Subject: [PATCH 12/96] Included iostream for testing [apple] --- include/boost/math/special_functions/beta.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index 78a1e96c06..08de2fc55b 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -33,6 +33,7 @@ #include #include #include +#include namespace boost{ namespace math{ From 552284c319badee30bb9a166b3b799a40b5e50fd Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 1 Mar 2026 17:16:21 -0800 Subject: [PATCH 13/96] Added more debugging info --- include/boost/math/special_functions/beta.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index 08de2fc55b..03a6528466 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -1213,6 +1213,7 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b { std::cout << std::setprecision(64) << "a= " << a << std::endl; std::cout << "b= " << b << std::endl; + std::cout << "x= " << x << std::endl; constexpr auto function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)"; typedef typename lanczos::lanczos::type lanczos_type; BOOST_MATH_STD_USING // for ADL of std math functions. @@ -1626,6 +1627,7 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b else{ fract = local_result / local_fract; std::cout << "Using continued fractions with number of terms: " << max_terms << std::endl; + std::cout << "Series converges with epsilon: " << boost::math::policies::get_epsilon() << std::endl; } } else From 4042e6633b7f3bba0dd02380f9a2ac5dc1f96719 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 1 Mar 2026 18:41:48 -0800 Subject: [PATCH 14/96] Added more debugging info for continued fraction --- include/boost/math/special_functions/beta.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index 03a6528466..4162ade158 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -1615,6 +1615,9 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b if (local_result != 0) { ibeta_fraction2_t f(a, b, x, y); + ibeta_fraction2_t g(a, b, x, y); + std::cout << "First (a,b) pair: (" << g().first << ", " << g().second << ")" << std::endl; + std::cout << "Second (a,b) pair: (" << g().first << ", " << g().second << ")" << std::endl; boost::math::uintmax_t max_terms = boost::math::policies::get_max_series_iterations(); T local_fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon(), max_terms); if (max_terms >= boost::math::policies::get_max_series_iterations()) From 89a3455c2b96e0bca1b9cc36840c5ed0fde83e0f Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 1 Mar 2026 19:01:59 -0800 Subject: [PATCH 15/96] Got rid of int upcasting in ibeta_fraction2_t --- include/boost/math/special_functions/beta.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index 4162ade158..49f415a279 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -819,14 +819,14 @@ struct ibeta_fraction2_t { typedef boost::math::pair result_type; - BOOST_MATH_GPU_ENABLED ibeta_fraction2_t(T a_, T b_, T x_, T y_) : a(a_), b(b_), x(x_), y(y_), m(0) {} + BOOST_MATH_GPU_ENABLED ibeta_fraction2_t(T a_, T b_, T x_, T y_) : a(a_), b(b_), x(x_), y(y_), m(static_cast(0)) {} BOOST_MATH_GPU_ENABLED result_type operator()() { T denom = (a + 2 * m - 1); T aN = (m * (a + m - 1) / denom) * ((a + b + m - 1) / denom) * (b - m) * x * x; - T bN = static_cast(m); + T bN = m; bN += (m * (b - m) * x) / (a + 2*m - 1); bN += ((a + m) * (a * y - b * x + 1 + m *(2 - x))) / (a + 2*m + 1); @@ -837,7 +837,7 @@ struct ibeta_fraction2_t private: T a, b, x, y; - int m; + T m; }; // // Evaluate the incomplete beta via the continued fraction representation: From 050a24a8a444dc0f6ef98698b7735a98a9fe655b Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 1 Mar 2026 19:37:05 -0800 Subject: [PATCH 16/96] More debug output for continued fraction --- include/boost/math/special_functions/beta.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index 49f415a279..50deeac59d 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -1631,6 +1631,8 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b fract = local_result / local_fract; std::cout << "Using continued fractions with number of terms: " << max_terms << std::endl; std::cout << "Series converges with epsilon: " << boost::math::policies::get_epsilon() << std::endl; + std::cout << "Local result: " << local_result << std::endl; + std::cout << "Local fract: " << local_fract << std::endl; } } else From 3c6f6ae44395bb5e37c9a4fff21dab768848f491 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 1 Mar 2026 20:27:16 -0800 Subject: [PATCH 17/96] Added debug for lanczos approximation --- include/boost/math/special_functions/beta.hpp | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index 50deeac59d..b98a21ab85 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -476,7 +476,7 @@ BOOST_MATH_GPU_ENABLED T ibeta_power_terms(T a, const char* = "boost::math::ibeta<%1%>(%1%, %1%, %1%)") { BOOST_MATH_STD_USING - + std::cout << "Not using lanczos approximation" << std::endl; if(!normalised) { return prefix * pow(x, a) * pow(y, b); @@ -1616,8 +1616,6 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b { ibeta_fraction2_t f(a, b, x, y); ibeta_fraction2_t g(a, b, x, y); - std::cout << "First (a,b) pair: (" << g().first << ", " << g().second << ")" << std::endl; - std::cout << "Second (a,b) pair: (" << g().first << ", " << g().second << ")" << std::endl; boost::math::uintmax_t max_terms = boost::math::policies::get_max_series_iterations(); T local_fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon(), max_terms); if (max_terms >= boost::math::policies::get_max_series_iterations()) @@ -1629,8 +1627,6 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b } else{ fract = local_result / local_fract; - std::cout << "Using continued fractions with number of terms: " << max_terms << std::endl; - std::cout << "Series converges with epsilon: " << boost::math::policies::get_epsilon() << std::endl; std::cout << "Local result: " << local_result << std::endl; std::cout << "Local fract: " << local_fract << std::endl; } @@ -1662,9 +1658,6 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b } } std::cout << "Fract is: " << fract << std::endl; - std::cout << "invert: " << invert << std::endl; - std::cout << "normalized: " << normalised << std::endl; - std::cout << "Beta(a,b): " << boost::math::beta(a, b, pol) << std::endl; std::cout << "ibeta(a,b)=" << (invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract) << std::endl; std::cout << "===========================" << std::endl; return invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract; From ca4bda357a15ab0a76f0ee19407016c9b76a33fd Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Thu, 5 Mar 2026 19:50:17 -0800 Subject: [PATCH 18/96] Manual type promotion for arcsine distribution --- include/boost/math/distributions/arcsine.hpp | 277 +++++++++++-------- 1 file changed, 165 insertions(+), 112 deletions(-) diff --git a/include/boost/math/distributions/arcsine.hpp b/include/boost/math/distributions/arcsine.hpp index 29feac48fe..b9a09a28d9 100644 --- a/include/boost/math/distributions/arcsine.hpp +++ b/include/boost/math/distributions/arcsine.hpp @@ -164,7 +164,107 @@ namespace boost return check_dist(function, x_min, x_max, result, pol) && check_prob(function, p, result, pol); } // bool check_dist_and_prob + + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_mean(const RealType x_min, const RealType x_max) + { + return (x_min + x_max) / 2; + } + + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_variance(const RealType x_min, const RealType x_max) + { + return (x_max - x_min) * (x_max - x_min) / 8; + } + + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_pdf(const RealType x, const RealType x_min, const RealType x_max) + { + BOOST_MATH_STD_USING + using boost::math::constants::pi; + return 1 / (pi() * sqrt((x - x_min) * (x_max - x))); + } + + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_cdf(const RealType x, const RealType x_min, const RealType x_max) + { + BOOST_MATH_STD_USING + // Special cases: + if (x == x_min) + { + return 0; + } + else if (x == x_max) + { + return 1; + } + using boost::math::constants::pi; + return 2 * asin(sqrt((x - x_min) / (x_max - x_min))) / pi(); + } + + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_ccdf(const RealType x, const RealType x_min, const RealType x_max) + { + BOOST_MATH_STD_USING + // Special cases + if (x == x_min) + { + return 1; + } + else if (x == x_max) + { + return 0; + } + using boost::math::constants::pi; + // Naive version x = 1 - x; + // result = static_cast(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi(); + // is less accurate, so use acos instead of asin for complement. + return 2 * acos(sqrt((x - x_min) / (x_max - x_min))) / pi(); + } + + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_quantile(const RealType x, const RealType p, const RealType x_min, const RealType x_max) + { + BOOST_MATH_STD_USING + // Special cases: + if (p == 0) + { + return 0; + } + if (p == 1) + { + return 1; + } + using boost::math::constants::half_pi; + RealType sin2hpip = sin(half_pi() * p); + RealType sin2hpip2 = sin2hpip * sin2hpip; + result = -x_min * sin2hpip2 + x_min + x_max * sin2hpip2; + + return result; + } + template + BOOST_MATH_GPU_ENABLED inline RealType arcsine_cquantile(const RealType x, const RealType q, const RealType x_min, const RealType x_max) + { + BOOST_MATH_STD_USING + // Special cases: + if (q == 1) + { + return 0; + } + if (q == 0) + { + return 1; + } + // Naive RealType p = 1 - q; result = sin(half_pi() * p); loses accuracy, so use a cos alternative instead. + //result = cos(half_pi() * q); // for arcsine(0,1) + //result = result * result; + // For generalized arcsine: + using boost::math::constants::half_pi; + RealType cos2hpip = cos(half_pi() * q); + RealType cos2hpip2 = cos2hpip * cos2hpip; + result = -x_min * cos2hpip2 + x_min + x_max * cos2hpip2; + return result; } // namespace arcsine_detail template > @@ -230,7 +330,7 @@ namespace boost RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::mean(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -239,7 +339,9 @@ namespace boost { return result; } - return (x_min + x_max) / 2; + + typedef policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_mean(static_cast(x_min), static_cast(x_max))); } // mean template @@ -248,7 +350,7 @@ namespace boost RealType result; RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::variance(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -257,7 +359,9 @@ namespace boost { return result; } - return (x_max - x_min) * (x_max - x_min) / 8; + + typedef policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_variance(static_cast(x_min), static_cast(x_max))); } // variance template @@ -277,7 +381,7 @@ namespace boost RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); RealType result; - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::median(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -286,7 +390,9 @@ namespace boost { return result; } - return (x_min + x_max) / 2; + // The median is the same as the mean + typedef policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_mean(static_cast(x_min), static_cast(x_max))); } template @@ -296,7 +402,7 @@ namespace boost RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::skewness(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -315,7 +421,7 @@ namespace boost RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::kurtosis_excess(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -324,8 +430,8 @@ namespace boost { return result; } - result = -3; - return result / 2; + typedef policies::evaluation_t policy_promoted_type; + return static_cast(static_cast(-3) / static_cast(2)); } // kurtosis_excess template @@ -335,7 +441,7 @@ namespace boost RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist( + if (!arcsine_detail::check_dist( "boost::math::kurtosis(arcsine_distribution<%1%> const&, %1% )", x_min, x_max, @@ -344,75 +450,62 @@ namespace boost { return result; } - - return 3 + kurtosis_excess(dist); + typedef policies::evaluation_t policy_promoted_type; + return static_cast(static_cast(3) / static_cast(2)); } // kurtosis template BOOST_MATH_GPU_ENABLED inline RealType pdf(const arcsine_distribution& dist, const RealType& xx) { // Probability Density/Mass Function arcsine. BOOST_FPU_EXCEPTION_GUARD - BOOST_MATH_STD_USING // For ADL of std functions. - - constexpr auto function = "boost::math::pdf(arcsine_distribution<%1%> const&, %1%)"; - RealType lo = dist.x_min(); - RealType hi = dist.x_max(); + RealType x_min = dist.x_min(); + RealType x_max = dist.x_max(); RealType x = xx; // Argument checks: RealType result = 0; - if (false == arcsine_detail::check_dist_and_x( - function, - lo, hi, x, + if (!arcsine_detail::check_dist_and_x( + "boost::math::pdf(arcsine_distribution<%1%> const&, %1%)", + x_min, + x_max, + x, &result, Policy())) { return result; } - using boost::math::constants::pi; - result = static_cast(1) / (pi() * sqrt((x - lo) * (hi - x))); - return result; + typedef policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_pdf(static_cast(x), + static_cast(x_min), + static_cast(x_max))); } // pdf template BOOST_MATH_GPU_ENABLED inline RealType cdf(const arcsine_distribution& dist, const RealType& x) { // Cumulative Distribution Function arcsine. - BOOST_MATH_STD_USING // For ADL of std functions. - - constexpr auto function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)"; - RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); // Argument checks: RealType result = 0; - if (false == arcsine_detail::check_dist_and_x( - function, - x_min, x_max, x, + if (!arcsine_detail::check_dist_and_x( + "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)", + x_min, + x_max, + x, &result, Policy())) { return result; } - // Special cases: - if (x == x_min) - { - return 0; - } - else if (x == x_max) - { - return 1; - } - using boost::math::constants::pi; - result = static_cast(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi(); - return result; + typedef policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_cdf(static_cast(x), + static_cast(x_min), + static_cast(x_max))); } // arcsine cdf template BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type, RealType>& c) { // Complemented Cumulative Distribution Function arcsine. - BOOST_MATH_STD_USING // For ADL of std functions. - constexpr auto function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)"; - RealType x = c.param; arcsine_distribution const& dist = c.dist; RealType x_min = dist.x_min(); @@ -420,27 +513,19 @@ namespace boost // Argument checks: RealType result = 0; - if (false == arcsine_detail::check_dist_and_x( - function, - x_min, x_max, x, + if (!arcsine_detail::check_dist_and_x( + "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)", + x_min, + x_max, + x, &result, Policy())) { return result; } - if (x == x_min) - { - return 0; - } - else if (x == x_max) - { - return 1; - } - using boost::math::constants::pi; - // Naive version x = 1 - x; - // result = static_cast(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi(); - // is less accurate, so use acos instead of asin for complement. - result = static_cast(2) * acos(sqrt((x - x_min) / (x_max - x_min))) / pi(); - return result; + typedef policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_ccdf(static_cast(x), + static_cast(x_min), + static_cast(x_max))); } // arcsine ccdf template @@ -454,37 +539,23 @@ namespace boost // and return a value such that the probability that a random variable x // will be less than or equal to that value // is whatever probability you supplied as an argument. - BOOST_MATH_STD_USING // For ADL of std functions. - - using boost::math::constants::half_pi; - - constexpr auto function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)"; - RealType result = 0; // of argument checks: RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist_and_prob( - function, - x_min, x_max, p, + if (!arcsine_detail::check_dist_and_prob( + "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)", + x_min, + x_max, + p, &result, Policy())) { return result; } - // Special cases: - if (p == 0) - { - return 0; - } - if (p == 1) - { - return 1; - } - - RealType sin2hpip = sin(half_pi() * p); - RealType sin2hpip2 = sin2hpip * sin2hpip; - result = -x_min * sin2hpip2 + x_min + x_max * sin2hpip2; - - return result; + typedef policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_quantile(static_cast(x), + static_cast(p), + static_cast(x_min), + static_cast(x_max))); } // quantile template @@ -493,19 +564,14 @@ namespace boost // Complement Quantile or Percent Point arcsine function. // Return the number of expected x for a given // complement of the probability q. - BOOST_MATH_STD_USING // For ADL of std functions. - - using boost::math::constants::half_pi; - constexpr auto function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)"; - // Error checks: RealType q = c.param; const arcsine_distribution& dist = c.dist; RealType result = 0; RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); - if (false == arcsine_detail::check_dist_and_prob( - function, + if (!arcsine_detail::check_dist_and_prob( + "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)", x_min, x_max, q, @@ -513,24 +579,11 @@ namespace boost { return result; } - // Special cases: - if (q == 1) - { - return 0; - } - if (q == 0) - { - return 1; - } - // Naive RealType p = 1 - q; result = sin(half_pi() * p); loses accuracy, so use a cos alternative instead. - //result = cos(half_pi() * q); // for arcsine(0,1) - //result = result * result; - // For generalized arcsine: - RealType cos2hpip = cos(half_pi() * q); - RealType cos2hpip2 = cos2hpip * cos2hpip; - result = -x_min * cos2hpip2 + x_min + x_max * cos2hpip2; - - return result; + typedef policies::evaluation_t policy_promoted_type; + return static_cast(arcsine_detail::arcsine_quantile(static_cast(x), + static_cast(q), + static_cast(x_min), + static_cast(x_max))); } // Quantile Complement } // namespace math From d0004eabbd1b6382417a20b81f0f02f4a4b2775f Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Thu, 5 Mar 2026 20:15:02 -0800 Subject: [PATCH 19/96] Added typename to definitions and fixed undeclared result --- include/boost/math/distributions/arcsine.hpp | 29 ++++++++------------ 1 file changed, 12 insertions(+), 17 deletions(-) diff --git a/include/boost/math/distributions/arcsine.hpp b/include/boost/math/distributions/arcsine.hpp index b9a09a28d9..9db778fcfb 100644 --- a/include/boost/math/distributions/arcsine.hpp +++ b/include/boost/math/distributions/arcsine.hpp @@ -238,9 +238,7 @@ namespace boost using boost::math::constants::half_pi; RealType sin2hpip = sin(half_pi() * p); RealType sin2hpip2 = sin2hpip * sin2hpip; - result = -x_min * sin2hpip2 + x_min + x_max * sin2hpip2; - - return result; + return -x_min * sin2hpip2 + x_min + x_max * sin2hpip2; } template BOOST_MATH_GPU_ENABLED inline RealType arcsine_cquantile(const RealType x, const RealType q, const RealType x_min, const RealType x_max) @@ -262,9 +260,7 @@ namespace boost using boost::math::constants::half_pi; RealType cos2hpip = cos(half_pi() * q); RealType cos2hpip2 = cos2hpip * cos2hpip; - result = -x_min * cos2hpip2 + x_min + x_max * cos2hpip2; - - return result; + return -x_min * cos2hpip2 + x_min + x_max * cos2hpip2; } // namespace arcsine_detail template > @@ -340,7 +336,7 @@ namespace boost return result; } - typedef policies::evaluation_t policy_promoted_type; + typedef typename policies::evaluation_t policy_promoted_type; return static_cast(arcsine_detail::arcsine_mean(static_cast(x_min), static_cast(x_max))); } // mean @@ -360,7 +356,7 @@ namespace boost return result; } - typedef policies::evaluation_t policy_promoted_type; + typedef typename policies::evaluation_t policy_promoted_type; return static_cast(arcsine_detail::arcsine_variance(static_cast(x_min), static_cast(x_max))); } // variance @@ -391,7 +387,7 @@ namespace boost return result; } // The median is the same as the mean - typedef policies::evaluation_t policy_promoted_type; + typedef typename policies::evaluation_t policy_promoted_type; return static_cast(arcsine_detail::arcsine_mean(static_cast(x_min), static_cast(x_max))); } @@ -430,7 +426,7 @@ namespace boost { return result; } - typedef policies::evaluation_t policy_promoted_type; + typedef typename policies::evaluation_t policy_promoted_type; return static_cast(static_cast(-3) / static_cast(2)); } // kurtosis_excess @@ -450,7 +446,7 @@ namespace boost { return result; } - typedef policies::evaluation_t policy_promoted_type; + typedef typename policies::evaluation_t policy_promoted_type; return static_cast(static_cast(3) / static_cast(2)); } // kurtosis @@ -474,7 +470,7 @@ namespace boost { return result; } - typedef policies::evaluation_t policy_promoted_type; + typedef typename policies::evaluation_t policy_promoted_type; return static_cast(arcsine_detail::arcsine_pdf(static_cast(x), static_cast(x_min), static_cast(x_max))); @@ -497,7 +493,7 @@ namespace boost { return result; } - typedef policies::evaluation_t policy_promoted_type; + typedef typename policies::evaluation_t policy_promoted_type; return static_cast(arcsine_detail::arcsine_cdf(static_cast(x), static_cast(x_min), static_cast(x_max))); @@ -522,7 +518,7 @@ namespace boost { return result; } - typedef policies::evaluation_t policy_promoted_type; + typedef typename policies::evaluation_t policy_promoted_type; return static_cast(arcsine_detail::arcsine_ccdf(static_cast(x), static_cast(x_min), static_cast(x_max))); @@ -551,7 +547,7 @@ namespace boost { return result; } - typedef policies::evaluation_t policy_promoted_type; + typedef typename policies::evaluation_t policy_promoted_type; return static_cast(arcsine_detail::arcsine_quantile(static_cast(x), static_cast(p), static_cast(x_min), @@ -579,13 +575,12 @@ namespace boost { return result; } - typedef policies::evaluation_t policy_promoted_type; + typedef typename policies::evaluation_t policy_promoted_type; return static_cast(arcsine_detail::arcsine_quantile(static_cast(x), static_cast(q), static_cast(x_min), static_cast(x_max))); } // Quantile Complement - } // namespace math } // namespace boost From 99b6d453a00b8c061c8027f1847b08462b5e4eeb Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Thu, 5 Mar 2026 20:19:08 -0800 Subject: [PATCH 20/96] Added evaluation_t shorthand --- include/boost/math/policies/policy.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/include/boost/math/policies/policy.hpp b/include/boost/math/policies/policy.hpp index ccc33bb9f8..6d80b2b292 100644 --- a/include/boost/math/policies/policy.hpp +++ b/include/boost/math/policies/policy.hpp @@ -766,6 +766,9 @@ struct evaluation using type = typename boost::math::conditional::type; }; +template +using evaluation_t = typename evaluation::type; + template struct precision { From 423874d341b8264bc7b3c1c560e87fea5c3eca5d Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Fri, 6 Mar 2026 10:49:02 -0800 Subject: [PATCH 21/96] A couple of bug fixes --- include/boost/math/distributions/arcsine.hpp | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/include/boost/math/distributions/arcsine.hpp b/include/boost/math/distributions/arcsine.hpp index 9db778fcfb..02bb857ce2 100644 --- a/include/boost/math/distributions/arcsine.hpp +++ b/include/boost/math/distributions/arcsine.hpp @@ -223,7 +223,7 @@ namespace boost } template - BOOST_MATH_GPU_ENABLED inline RealType arcsine_quantile(const RealType x, const RealType p, const RealType x_min, const RealType x_max) + BOOST_MATH_GPU_ENABLED inline RealType arcsine_quantile( const RealType p, const RealType x_min, const RealType x_max) { BOOST_MATH_STD_USING // Special cases: @@ -241,7 +241,7 @@ namespace boost return -x_min * sin2hpip2 + x_min + x_max * sin2hpip2; } template - BOOST_MATH_GPU_ENABLED inline RealType arcsine_cquantile(const RealType x, const RealType q, const RealType x_min, const RealType x_max) + BOOST_MATH_GPU_ENABLED inline RealType arcsine_cquantile(const RealType q, const RealType x_min, const RealType x_max) { BOOST_MATH_STD_USING // Special cases: @@ -460,7 +460,7 @@ namespace boost RealType x = xx; // Argument checks: - RealType result = 0; + RealType result; if (!arcsine_detail::check_dist_and_x( "boost::math::pdf(arcsine_distribution<%1%> const&, %1%)", x_min, @@ -471,9 +471,9 @@ namespace boost return result; } typedef typename policies::evaluation_t policy_promoted_type; - return static_cast(arcsine_detail::arcsine_pdf(static_cast(x), - static_cast(x_min), - static_cast(x_max))); + return static_cast(arcsine_detail::arcsine_pdf(static_cast(x), + static_cast(x_min), + static_cast(x_max))); } // pdf template @@ -548,10 +548,9 @@ namespace boost return result; } typedef typename policies::evaluation_t policy_promoted_type; - return static_cast(arcsine_detail::arcsine_quantile(static_cast(x), - static_cast(p), + return static_cast(arcsine_detail::arcsine_quantile(static_cast(p), static_cast(x_min), - static_cast(x_max))); + static_cast(x_max))); } // quantile template @@ -576,8 +575,7 @@ namespace boost return result; } typedef typename policies::evaluation_t policy_promoted_type; - return static_cast(arcsine_detail::arcsine_quantile(static_cast(x), - static_cast(q), + return static_cast(arcsine_detail::arcsine_quantile(static_cast(q), static_cast(x_min), static_cast(x_max))); } // Quantile Complement From 89b393f22976cb17d79d38a61c0db631dfca48aa Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Fri, 6 Mar 2026 13:27:33 -0800 Subject: [PATCH 22/96] [ci skip] Added missing bracket --- include/boost/math/distributions/arcsine.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/include/boost/math/distributions/arcsine.hpp b/include/boost/math/distributions/arcsine.hpp index 02bb857ce2..3778a4d89a 100644 --- a/include/boost/math/distributions/arcsine.hpp +++ b/include/boost/math/distributions/arcsine.hpp @@ -261,6 +261,7 @@ namespace boost RealType cos2hpip = cos(half_pi() * q); RealType cos2hpip2 = cos2hpip * cos2hpip; return -x_min * cos2hpip2 + x_min + x_max * cos2hpip2; + } } // namespace arcsine_detail template > From 5e7591773a88dbc941a06461d0d758247ec2bd2e Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Fri, 6 Mar 2026 19:07:18 -0800 Subject: [PATCH 23/96] Fixed typo in complement quantile --- include/boost/math/distributions/arcsine.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/math/distributions/arcsine.hpp b/include/boost/math/distributions/arcsine.hpp index 3778a4d89a..2241320634 100644 --- a/include/boost/math/distributions/arcsine.hpp +++ b/include/boost/math/distributions/arcsine.hpp @@ -576,7 +576,7 @@ namespace boost return result; } typedef typename policies::evaluation_t policy_promoted_type; - return static_cast(arcsine_detail::arcsine_quantile(static_cast(q), + return static_cast(arcsine_detail::arcsine_cquantile(static_cast(q), static_cast(x_min), static_cast(x_max))); } // Quantile Complement From 8e919f2380c892f5ea2764d7833113cf9d852b3a Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sat, 7 Mar 2026 10:47:32 -0800 Subject: [PATCH 24/96] Accepted review comments; testing code coverage --- include/boost/math/distributions/arcsine.hpp | 6 ++---- test/test_arcsine.cpp | 8 ++++++++ 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/include/boost/math/distributions/arcsine.hpp b/include/boost/math/distributions/arcsine.hpp index 2241320634..4fe3c7d751 100644 --- a/include/boost/math/distributions/arcsine.hpp +++ b/include/boost/math/distributions/arcsine.hpp @@ -427,8 +427,7 @@ namespace boost { return result; } - typedef typename policies::evaluation_t policy_promoted_type; - return static_cast(static_cast(-3) / static_cast(2)); + return static_cast(-1.5f); } // kurtosis_excess template @@ -447,8 +446,7 @@ namespace boost { return result; } - typedef typename policies::evaluation_t policy_promoted_type; - return static_cast(static_cast(3) / static_cast(2)); + return static_cast(1.5f); } // kurtosis template diff --git a/test/test_arcsine.cpp b/test/test_arcsine.cpp index eed8dc9481..cb5646e627 100644 --- a/test/test_arcsine.cpp +++ b/test/test_arcsine.cpp @@ -437,11 +437,19 @@ void test_spots(RealType) BOOST_CHECK_CLOSE_FRACTION(cdf(as_m2m1, -1.05), static_cast(0.85643370687129372924905811522494428117838480010259L), tolerance); BOOST_CHECK_CLOSE_FRACTION(cdf(as_m2m1, -1.5), static_cast(0.5L), tolerance); BOOST_CHECK_CLOSE_FRACTION(cdf(as_m2m1, -1.95), static_cast(0.14356629312870627075094188477505571882161519989741L), 8 * tolerance); // Not much less accurate. + BOOST_CHECK_EQUAL(cdf(as_m2m1, as_m2m1.x_min()), 0); + BOOST_CHECK_EQUAL(cdf(as_m2m1, as_m2m1.x_max()), 1); + BOOST_CHECK_EQUAL(cdf(complement(as_m2m1, as_m2m1.x_min())), 1); + BOOST_CHECK_EQUAL(cdf(complement(as_m2m1, as_m2m1.x_max())), 0); // Quantile BOOST_CHECK_CLOSE_FRACTION(quantile(as_m2m1, static_cast(0.85643370687129372924905811522494428117838480010259L)), -static_cast(1.05L), 2 * tolerance); // BOOST_CHECK_CLOSE_FRACTION(quantile(as_m2m1, static_cast(0.5L)), -static_cast(1.5L), 2 * tolerance); // BOOST_CHECK_CLOSE_FRACTION(quantile(as_m2m1, static_cast(0.14356629312870627075094188477505571882161519989741L)), -static_cast(1.95L), 4 * tolerance); // + BOOST_CHECK_EQUAL(quantile(as_m2m1, 1), 1); + BOOST_CHECK_EQUAL(quantile(as_m2m1, 0), 0); + BOOST_CHECK_EQUAL(quantile(complement(as_m2m1, 1)), 0); + BOOST_CHECK_EQUAL(quantile(complement(as_m2m1, 0)), 1); BOOST_CHECK_CLOSE_FRACTION(quantile(complement(as_m2m1, static_cast(0.14356629312870627075094188477505571882161519989741L))), -static_cast(1.05L), 2 * tolerance); // BOOST_CHECK_CLOSE_FRACTION(quantile(as_m2m1, static_cast(0.5L)), -static_cast(1.5L), 2 * tolerance); // From d1eff51ed8bd20ae0e0becd8032d69d1c9036157 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 8 Mar 2026 14:33:38 -0700 Subject: [PATCH 25/96] Removed debugging information --- include/boost/math/special_functions/beta.hpp | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index b98a21ab85..3f31bea408 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -33,7 +33,6 @@ #include #include #include -#include namespace boost{ namespace math{ @@ -476,7 +475,6 @@ BOOST_MATH_GPU_ENABLED T ibeta_power_terms(T a, const char* = "boost::math::ibeta<%1%>(%1%, %1%, %1%)") { BOOST_MATH_STD_USING - std::cout << "Not using lanczos approximation" << std::endl; if(!normalised) { return prefix * pow(x, a) * pow(y, b); @@ -1211,9 +1209,6 @@ BOOST_MATH_GPU_ENABLED T ibeta_large_ab(T a, T b, T x, T y, bool invert, bool no template BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_derivative) { - std::cout << std::setprecision(64) << "a= " << a << std::endl; - std::cout << "b= " << b << std::endl; - std::cout << "x= " << x << std::endl; constexpr auto function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)"; typedef typename lanczos::lanczos::type lanczos_type; BOOST_MATH_STD_USING // for ADL of std math functions. @@ -1596,7 +1591,6 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b { fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol); invert = false; - std::cout << "Using Erf approximation: " << fract << std::endl; } else { @@ -1623,12 +1617,9 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b // Continued fraction failed, fall back to asymptotic expansion: fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol); invert = false; - std::cout << "Failed continued fractions so falling back to erf" << std::endl; } else{ fract = local_result / local_fract; - std::cout << "Local result: " << local_result << std::endl; - std::cout << "Local fract: " << local_fract << std::endl; } } else @@ -1657,9 +1648,6 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b } } } - std::cout << "Fract is: " << fract << std::endl; - std::cout << "ibeta(a,b)=" << (invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract) << std::endl; - std::cout << "===========================" << std::endl; return invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract; } // template T ibeta_imp(T a, T b, T x, const Lanczos& l, bool inv, bool normalised) From e1bdcad01783fe8d285c363c948b47173846e734 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 8 Mar 2026 15:32:47 -0700 Subject: [PATCH 26/96] Removed failing tests and added code comments --- test/test_ibeta.hpp | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/test/test_ibeta.hpp b/test/test_ibeta.hpp index 7bef2d9715..13ff157f9d 100644 --- a/test/test_ibeta.hpp +++ b/test/test_ibeta.hpp @@ -490,19 +490,26 @@ void test_spots(T, const char* name) BOOST_CHECK_EQUAL(boost::math::ibeta(static_cast(1), static_cast(2), static_cast(0)), 0); // Bug testing for large a,b and x close to a / (a+b). See PR 1363. - // The values for a,b are just too large for floats to handle. + // The values for a,b are just too large for floats to handle. The tests here only + // check if ibeta is monotonically increasing for a,b fixed and x increasing. Both + // Mathematica and mpmath (in Python) are not able to evaluate ibeta for such large + // values of a,b so spot testing wasn't possible. The accuracy of the fix for + // these large values is very sensitive to the computer architecture. + // Macos with arm64 does the worst but linux and windows pass a larger range of tests. if (!std::is_same::value) { - T a_values[2] = {10000000272564224, 3.1622776601699636e+16}; - T b_values[2] = {9965820922822656, 3.130654883566682e+18}; - T delta[8] = {0, 1e-15, 1e-14, 1e-13, 1e-12, 1e-11, 1e-10, 1e-9}; + // Larger values of a,b become more numerically unstable. Larger/smaller values + // of delta also become unstable. + T a_values[1] = {10000000272564224}; + T b_values[1] = {9965820922822656}; + T delta[7] = {0, 1e-15, 1e-14, 1e-13, 1e-12, 1e-11, 1e-10}; T a, b, x; - for (unsigned int i=0; i<2; i++){ + for (unsigned int i=0; i<1; i++){ a = a_values[i]; b = b_values[i]; x = a / (a+b); // roughly the median of ibeta - for (unsigned int j=0; j < 7; j++) + for (unsigned int j=0; j < 6; j++) { BOOST_CHECK_MESSAGE(boost::math::ibeta(a, b, x + delta[j+1]) > boost::math::ibeta(a, b, x + delta[j]), "ibeta not monotonically increasing above a/(a+b) for ibeta(" << a << ", " << b << ", " << x << ") and delta=" << delta[j] << ": [" << boost::math::ibeta(a, b, x + delta[j+1]) << " >= " << boost::math::ibeta(a, b, x + delta[j]) << "]"); From 8fec9d93cd2a0ab4d8237bd25724a8ed66f43344 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 9 Mar 2026 08:36:02 -0400 Subject: [PATCH 27/96] Fix scalar div for octonions --- include/boost/math/octonion.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/include/boost/math/octonion.hpp b/include/boost/math/octonion.hpp index 5d147df584..4239140564 100644 --- a/include/boost/math/octonion.hpp +++ b/include/boost/math/octonion.hpp @@ -942,6 +942,10 @@ namespace boost b /= rhs; \ c /= rhs; \ d /= rhs; \ + e /= rhs; \ + f /= rhs; \ + g /= rhs; \ + h /= rhs; \ \ return(*this); \ } From 0074609cda2ebfbe26bfc4f76a0317e96a63ffc8 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 9 Mar 2026 08:36:15 -0400 Subject: [PATCH 28/96] Add division test suite mirroring multiplication --- test/octonion_test.cpp | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/test/octonion_test.cpp b/test/octonion_test.cpp index 96f43103c5..00ebf535ff 100644 --- a/test/octonion_test.cpp +++ b/test/octonion_test.cpp @@ -699,6 +699,34 @@ BOOST_TEST_CASE_TEMPLATE_FUNCTION(multiplication_test, T) } } +BOOST_TEST_CASE_TEMPLATE_FUNCTION(division_test, T) +{ + #if BOOST_WORKAROUND(__GNUC__, < 3) + #else /* BOOST_WORKAROUND(__GNUC__, < 3) */ + using ::std::numeric_limits; + + using ::boost::math::abs; + #endif /* BOOST_WORKAROUND(__GNUC__, < 3) */ + + + BOOST_TEST_MESSAGE("Testing division for " + << string_type_name::_() << "."); + + BOOST_REQUIRE_PREDICATE(::std::less_equal(), + (abs(::boost::math::octonion(1,0,0,0,0,0,0,0)/ + ::boost::math::octonion(1,0,0,0,0,0,0,0)- + static_cast(1))) + (numeric_limits::epsilon())); + + for (int idx = 1; idx < 8; ++idx) + { + ::boost::math::octonion toto = index_i_element(idx); + + BOOST_REQUIRE_PREDICATE(::std::less_equal(), + (abs(toto/toto-static_cast(1))) + (numeric_limits::epsilon())); + } +} BOOST_TEST_CASE_TEMPLATE_FUNCTION(exp_test, T) { @@ -752,6 +780,7 @@ boost::unit_test::test_suite * init_unit_test_suite(int, char *[]) #define BOOST_OCTONION_TEST \ BOOST_OCTONION_COMMON_GENERATOR(multiplication) \ + BOOST_OCTONION_COMMON_GENERATOR(division) \ BOOST_OCTONION_COMMON_GENERATOR_NEAR_EPS(exp) From a51619294964b2ef395be6c3c94ecf51faa99ed7 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 9 Mar 2026 09:10:19 -0400 Subject: [PATCH 29/96] Make general port of patch file from vcpkg --- CMakeLists.txt | 110 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 110 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index b789eb22aa..ba46c65bdf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -46,6 +46,116 @@ endif() target_compile_features(boost_math INTERFACE cxx_std_14) +# Legacy C99/TR1 compiled libraries +option(BOOST_MATH_BUILD_WITH_LEGACY_FUNCTIONS "Build the C99 and TR1 compiled libraries" OFF) + +if(BOOST_MATH_BUILD_WITH_LEGACY_FUNCTIONS) + +include(CheckCXXSourceCompiles) +set(CMAKE_REQUIRED_LIBRARIES Boost::config) +set(CMAKE_REQUIRED_INCLUDES "${CMAKE_CURRENT_SOURCE_DIR}/include") +check_cxx_source_compiles("#include <${CMAKE_CURRENT_SOURCE_DIR}/config/has_long_double_support.cpp> \n int main() { return 0;}" BOOST_MATH_HAS_LONG_DOUBLE) +unset(CMAKE_REQUIRED_LIBRARIES) +unset(CMAKE_REQUIRED_INCLUDES) + +set(C99_SOURCES + acosh + asinh + atanh + cbrt + copysign + erfc + erf + expm1 + fmax + fmin + fpclassify + hypot + lgamma + llround + log1p + lround + nextafter + nexttoward + round + tgamma + trunc +) + +set(TR1_SOURCES + assoc_laguerre + assoc_legendre + beta + comp_ellint_1 + comp_ellint_2 + comp_ellint_3 + cyl_bessel_i + cyl_bessel_j + cyl_bessel_k + cyl_neumann + ellint_1 + ellint_2 + ellint_3 + expint + hermite + laguerre + legendre + riemann_zeta + sph_bessel + sph_legendre + sph_neumann +) + +list(TRANSFORM C99_SOURCES PREPEND "src/tr1/") +list(TRANSFORM TR1_SOURCES PREPEND "src/tr1/") + +list(TRANSFORM C99_SOURCES APPEND "f.cpp" OUTPUT_VARIABLE C99_SOURCESf) +list(TRANSFORM TR1_SOURCES APPEND "f.cpp" OUTPUT_VARIABLE TR1_SOURCESf) + +set(types "" f) + +if(BOOST_MATH_HAS_LONG_DOUBLE) + list(TRANSFORM C99_SOURCES APPEND "l.cpp" OUTPUT_VARIABLE C99_SOURCESl) + list(TRANSFORM TR1_SOURCES APPEND "l.cpp" OUTPUT_VARIABLE TR1_SOURCESl) + list(APPEND types l) +endif() + +list(TRANSFORM C99_SOURCES APPEND ".cpp") +list(TRANSFORM TR1_SOURCES APPEND ".cpp") + +foreach(type IN LISTS types) + add_library(boost_math_tr1${type} ${TR1_SOURCES${type}}) + add_library(Boost::math_tr1${type} ALIAS boost_math_tr1${type}) + target_link_libraries(boost_math_tr1${type} PUBLIC Boost::config) + target_include_directories(boost_math_tr1${type} PRIVATE src/tr1) + target_include_directories(boost_math_tr1${type} PRIVATE include) + + add_library(boost_math_c99${type} ${C99_SOURCES${type}}) + add_library(Boost::math_c99${type} ALIAS boost_math_c99${type}) + target_link_libraries(boost_math_c99${type} PUBLIC Boost::config) + target_include_directories(boost_math_c99${type} PRIVATE src/tr1) + target_include_directories(boost_math_c99${type} PRIVATE include) + + if(BUILD_SHARED_LIBS) + target_compile_definitions(boost_math_tr1${type} PUBLIC BOOST_MATH_TR1_DYN_LINK=1) + target_compile_definitions(boost_math_c99${type} PUBLIC BOOST_MATH_TR1_DYN_LINK=1) + if(MSVC) + target_compile_definitions(boost_math_tr1${type} PRIVATE "BOOST_SYMBOL_EXPORT=__declspec(dllexport)" BOOST_ALL_NO_LIB) + target_compile_definitions(boost_math_c99${type} PRIVATE "BOOST_SYMBOL_EXPORT=__declspec(dllexport)" BOOST_ALL_NO_LIB) + endif() + endif() + + target_compile_features(boost_math_tr1${type} PUBLIC cxx_std_14) + target_compile_features(boost_math_c99${type} PUBLIC cxx_std_14) + + if(APPLE) + target_compile_definitions(boost_math_tr1${type} PRIVATE _DARWIN_C_SOURCE) + target_compile_definitions(boost_math_c99${type} PRIVATE _DARWIN_C_SOURCE) + endif() +endforeach() + +endif() + if(BUILD_TESTING AND EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/test/CMakeLists.txt") add_subdirectory(test) From 99048089124fb1dded8f5251272ec5e638f1b999 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 9 Mar 2026 09:25:07 -0400 Subject: [PATCH 30/96] Fix linking of config and installation --- CMakeLists.txt | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index ba46c65bdf..e2d6c055d3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -52,10 +52,14 @@ option(BOOST_MATH_BUILD_WITH_LEGACY_FUNCTIONS "Build the C99 and TR1 compiled li if(BOOST_MATH_BUILD_WITH_LEGACY_FUNCTIONS) include(CheckCXXSourceCompiles) -set(CMAKE_REQUIRED_LIBRARIES Boost::config) -set(CMAKE_REQUIRED_INCLUDES "${CMAKE_CURRENT_SOURCE_DIR}/include") +get_target_property(_config_type Boost::config TYPE) +if(_config_type STREQUAL "INTERFACE_LIBRARY") + get_target_property(_config_inc Boost::config INTERFACE_INCLUDE_DIRECTORIES) +else() + set(_config_inc "") +endif() +set(CMAKE_REQUIRED_INCLUDES "${CMAKE_CURRENT_SOURCE_DIR}/include" ${_config_inc}) check_cxx_source_compiles("#include <${CMAKE_CURRENT_SOURCE_DIR}/config/has_long_double_support.cpp> \n int main() { return 0;}" BOOST_MATH_HAS_LONG_DOUBLE) -unset(CMAKE_REQUIRED_LIBRARIES) unset(CMAKE_REQUIRED_INCLUDES) set(C99_SOURCES @@ -154,6 +158,14 @@ foreach(type IN LISTS types) endif() endforeach() +if(BOOST_SUPERPROJECT_VERSION AND NOT CMAKE_VERSION VERSION_LESS 3.13) + set(boost_math_legacy_targets "") + foreach(type IN LISTS types) + list(APPEND boost_math_legacy_targets boost_math_c99${type} boost_math_tr1${type}) + endforeach() + boost_install(TARGETS ${boost_math_legacy_targets} VERSION "${BOOST_SUPERPROJECT_VERSION}") +endif() + endif() if(BUILD_TESTING AND EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/test/CMakeLists.txt") From 7043cf3552a0a15caf87403110ba7001e7ea0f3b Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 9 Mar 2026 09:26:10 -0400 Subject: [PATCH 31/96] Add quick test for cmake install and subdir --- test/quick.cpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) create mode 100644 test/quick.cpp diff --git a/test/quick.cpp b/test/quick.cpp new file mode 100644 index 0000000000..30ac1b896b --- /dev/null +++ b/test/quick.cpp @@ -0,0 +1,17 @@ +// Copyright 2026 Matt Borland +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#include +#include + +int main() +{ + auto result = boost::math::tgamma(2.0); + // tgamma(2) == 1! == 1.0 + if (result < 0.99 || result > 1.01) + { + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +} From c3a3040a947242ac7a00c9e2d5f89c1b84b9d74f Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 9 Mar 2026 09:26:19 -0400 Subject: [PATCH 32/96] Add cmake install and subdir tests --- test/cmake_install_test/CMakeLists.txt | 17 +++++++++++++++++ test/cmake_subdir_test/CMakeLists.txt | 18 ++++++++++++++++++ 2 files changed, 35 insertions(+) create mode 100644 test/cmake_install_test/CMakeLists.txt create mode 100644 test/cmake_subdir_test/CMakeLists.txt diff --git a/test/cmake_install_test/CMakeLists.txt b/test/cmake_install_test/CMakeLists.txt new file mode 100644 index 0000000000..010d5e9fe3 --- /dev/null +++ b/test/cmake_install_test/CMakeLists.txt @@ -0,0 +1,17 @@ +# Copyright 2025 Matt Borland +# Distributed under the Boost Software License, Version 1.0. +# https://www.boost.org/LICENSE_1_0.txt + +cmake_minimum_required(VERSION 3.8...3.16) + +project(cmake_install_test LANGUAGES CXX) + +find_package(boost_math REQUIRED) + +add_executable(quick ../quick.cpp) +target_link_libraries(quick Boost::math) + +enable_testing() +add_test(quick quick) + +add_custom_target(check COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure -C $) diff --git a/test/cmake_subdir_test/CMakeLists.txt b/test/cmake_subdir_test/CMakeLists.txt new file mode 100644 index 0000000000..09d57446d7 --- /dev/null +++ b/test/cmake_subdir_test/CMakeLists.txt @@ -0,0 +1,18 @@ +# Copyright 2025 Matt Borland +# Distributed under the Boost Software License, Version 1.0. +# https://www.boost.org/LICENSE_1_0.txt + +cmake_minimum_required(VERSION 3.8...3.20) + +project(cmake_subdir_test LANGUAGES CXX) + +# Math is header-only and standalone by default (no deps needed) +add_subdirectory(../.. boostorg/math) + +add_executable(quick ../quick.cpp) +target_link_libraries(quick Boost::math) + +enable_testing() +add_test(quick quick) + +add_custom_target(check COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure -C $) From cca3167fecac7cf33c76cc040b12f9af22143c4c Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 9 Mar 2026 09:26:32 -0400 Subject: [PATCH 33/96] Add legacy c file intall and subdir tests --- test/cmake_install_test_legacy/CMakeLists.txt | 19 ++++++++++++++ test/cmake_subdir_test_legacy/CMakeLists.txt | 23 ++++++++++++++++ test/quick_tr1.cpp | 26 +++++++++++++++++++ 3 files changed, 68 insertions(+) create mode 100644 test/cmake_install_test_legacy/CMakeLists.txt create mode 100644 test/cmake_subdir_test_legacy/CMakeLists.txt create mode 100644 test/quick_tr1.cpp diff --git a/test/cmake_install_test_legacy/CMakeLists.txt b/test/cmake_install_test_legacy/CMakeLists.txt new file mode 100644 index 0000000000..b5db787104 --- /dev/null +++ b/test/cmake_install_test_legacy/CMakeLists.txt @@ -0,0 +1,19 @@ +# Copyright 2025 Matt Borland +# Distributed under the Boost Software License, Version 1.0. +# https://www.boost.org/LICENSE_1_0.txt + +cmake_minimum_required(VERSION 3.8...3.16) + +project(cmake_install_test_legacy LANGUAGES CXX) + +find_package(boost_math REQUIRED) +find_package(boost_math_c99 REQUIRED) +find_package(boost_math_tr1 REQUIRED) + +add_executable(quick_tr1 ../quick_tr1.cpp) +target_link_libraries(quick_tr1 Boost::math Boost::math_c99 Boost::math_tr1) + +enable_testing() +add_test(quick_tr1 quick_tr1) + +add_custom_target(check COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure -C $) diff --git a/test/cmake_subdir_test_legacy/CMakeLists.txt b/test/cmake_subdir_test_legacy/CMakeLists.txt new file mode 100644 index 0000000000..c606705fa0 --- /dev/null +++ b/test/cmake_subdir_test_legacy/CMakeLists.txt @@ -0,0 +1,23 @@ +# Copyright 2025 Matt Borland +# Distributed under the Boost Software License, Version 1.0. +# https://www.boost.org/LICENSE_1_0.txt + +cmake_minimum_required(VERSION 3.8...3.20) + +project(cmake_subdir_test_legacy LANGUAGES CXX) + +# config must be added before math (compiled targets depend on Boost::config) +add_subdirectory(../../../config boostorg/config) + +# Enable legacy compiled targets before adding math +set(BOOST_MATH_BUILD_WITH_LEGACY_FUNCTIONS ON CACHE BOOL "" FORCE) + +add_subdirectory(../.. boostorg/math) + +add_executable(quick_tr1 ../quick_tr1.cpp) +target_link_libraries(quick_tr1 Boost::math Boost::math_c99 Boost::math_tr1) + +enable_testing() +add_test(quick_tr1 quick_tr1) + +add_custom_target(check COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure -C $) diff --git a/test/quick_tr1.cpp b/test/quick_tr1.cpp new file mode 100644 index 0000000000..33e21a1495 --- /dev/null +++ b/test/quick_tr1.cpp @@ -0,0 +1,26 @@ +// Copyright 2026 Matt Borland +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#include +#include + +int main() +{ + // boost_cbrt from boost_math_c99 + auto cbrt_result = boost::math::tr1::boost_cbrt(8.0); + if (cbrt_result < 1.99 || cbrt_result > 2.01) + { + return EXIT_FAILURE; + } + + // boost_riemann_zeta from boost_math_tr1 + // zeta(2) == pi^2/6 ~= 1.6449 + auto zeta_result = boost::math::tr1::boost_riemann_zeta(2.0); + if (zeta_result < 1.64 || zeta_result > 1.65) + { + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +} From 43ba867230b4d8dcb9fd34659b2b17efb81eb8bc Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 9 Mar 2026 09:26:59 -0400 Subject: [PATCH 34/96] Add cmake subdir and install tests to CI --- .github/workflows/ci.yml | 169 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 169 insertions(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index fee18d22c7..0750f95e63 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -486,6 +486,175 @@ jobs: cd ../boost-root/__build__ cmake --build . --target tests + posix-cmake-subdir: + strategy: + fail-fast: false + matrix: + include: + - os: ubuntu-22.04 + - os: ubuntu-24.04 + - os: macos-15 + + runs-on: ${{matrix.os}} + + steps: + - uses: actions/checkout@v4 + + - name: Use library with add_subdirectory (header-only) + run: | + cd test/cmake_subdir_test + mkdir __build__ && cd __build__ + cmake .. + cmake --build . + ctest --output-on-failure --no-tests=error + + posix-cmake-subdir-legacy: + strategy: + fail-fast: false + matrix: + include: + - os: ubuntu-22.04 + - os: ubuntu-24.04 + - os: macos-15 + + runs-on: ${{matrix.os}} + + steps: + - uses: actions/checkout@v4 + + - name: Setup Boost + run: | + echo GITHUB_REPOSITORY: $GITHUB_REPOSITORY + LIBRARY=${GITHUB_REPOSITORY#*/} + echo LIBRARY: $LIBRARY + echo "LIBRARY=$LIBRARY" >> $GITHUB_ENV + echo GITHUB_BASE_REF: $GITHUB_BASE_REF + echo GITHUB_REF: $GITHUB_REF + REF=${GITHUB_BASE_REF:-$GITHUB_REF} + REF=${REF#refs/heads/} + echo REF: $REF + BOOST_BRANCH=develop && [ "$REF" == "master" ] && BOOST_BRANCH=master || true + echo BOOST_BRANCH: $BOOST_BRANCH + cd .. + git clone -b $BOOST_BRANCH --depth 1 https://github.com/boostorg/boost.git boost-root + cd boost-root + mkdir -p libs/$LIBRARY + cp -r $GITHUB_WORKSPACE/* libs/$LIBRARY + git submodule update --init tools/boostdep + python tools/boostdep/depinst/depinst.py --git_args "--jobs 3" $LIBRARY + + - name: Use library with add_subdirectory (legacy compiled) + run: | + cd ../boost-root/libs/$LIBRARY/test/cmake_subdir_test_legacy + mkdir __build__ && cd __build__ + cmake .. + cmake --build . + ctest --output-on-failure --no-tests=error + + posix-cmake-install: + strategy: + fail-fast: false + matrix: + include: + - os: macos-15 + + runs-on: ${{matrix.os}} + + steps: + - uses: actions/checkout@v4 + + - name: Setup Boost + run: | + echo GITHUB_REPOSITORY: $GITHUB_REPOSITORY + LIBRARY=${GITHUB_REPOSITORY#*/} + echo LIBRARY: $LIBRARY + echo "LIBRARY=$LIBRARY" >> $GITHUB_ENV + echo GITHUB_BASE_REF: $GITHUB_BASE_REF + echo GITHUB_REF: $GITHUB_REF + REF=${GITHUB_BASE_REF:-$GITHUB_REF} + REF=${REF#refs/heads/} + echo REF: $REF + BOOST_BRANCH=develop && [ "$REF" == "master" ] && BOOST_BRANCH=master || true + echo BOOST_BRANCH: $BOOST_BRANCH + cd .. + git clone -b $BOOST_BRANCH --depth 1 https://github.com/boostorg/boost.git boost-root + cd boost-root + mkdir -p libs/$LIBRARY + cp -r $GITHUB_WORKSPACE/* libs/$LIBRARY + git submodule update --init tools/boostdep + python tools/boostdep/depinst/depinst.py --git_args "--jobs 3" $LIBRARY + + - name: Configure + run: | + cd ../boost-root + mkdir __build__ && cd __build__ + cmake -DBOOST_INCLUDE_LIBRARIES=$LIBRARY -DCMAKE_INSTALL_PREFIX=~/.local .. + + - name: Install + run: | + cd ../boost-root/__build__ + cmake --build . --target install + + - name: Use the installed library (header-only) + run: | + cd ../boost-root/libs/$LIBRARY/test/cmake_install_test + mkdir __build__ && cd __build__ + cmake -DCMAKE_INSTALL_PREFIX=~/.local .. + cmake --build . + ctest --output-on-failure --no-tests=error + + posix-cmake-install-legacy: + strategy: + fail-fast: false + matrix: + include: + - os: macos-15 + + runs-on: ${{matrix.os}} + + steps: + - uses: actions/checkout@v4 + + - name: Setup Boost + run: | + echo GITHUB_REPOSITORY: $GITHUB_REPOSITORY + LIBRARY=${GITHUB_REPOSITORY#*/} + echo LIBRARY: $LIBRARY + echo "LIBRARY=$LIBRARY" >> $GITHUB_ENV + echo GITHUB_BASE_REF: $GITHUB_BASE_REF + echo GITHUB_REF: $GITHUB_REF + REF=${GITHUB_BASE_REF:-$GITHUB_REF} + REF=${REF#refs/heads/} + echo REF: $REF + BOOST_BRANCH=develop && [ "$REF" == "master" ] && BOOST_BRANCH=master || true + echo BOOST_BRANCH: $BOOST_BRANCH + cd .. + git clone -b $BOOST_BRANCH --depth 1 https://github.com/boostorg/boost.git boost-root + cd boost-root + mkdir -p libs/$LIBRARY + cp -r $GITHUB_WORKSPACE/* libs/$LIBRARY + git submodule update --init tools/boostdep + python tools/boostdep/depinst/depinst.py --git_args "--jobs 3" $LIBRARY + + - name: Configure + run: | + cd ../boost-root + mkdir __build__ && cd __build__ + cmake -DBOOST_INCLUDE_LIBRARIES=$LIBRARY -DBOOST_MATH_BUILD_WITH_LEGACY_FUNCTIONS=ON -DCMAKE_INSTALL_PREFIX=~/.local .. + + - name: Install + run: | + cd ../boost-root/__build__ + cmake --build . --target install + + - name: Use the installed library (legacy compiled) + run: | + cd ../boost-root/libs/$LIBRARY/test/cmake_install_test_legacy + mkdir __build__ && cd __build__ + cmake -DCMAKE_INSTALL_PREFIX=~/.local .. + cmake --build . + ctest --output-on-failure --no-tests=error + sycl-cmake-test: strategy: fail-fast: false From 6f59f3852acd161d56eaa04540bcc0b59b5bdd4f Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Wed, 11 Mar 2026 16:39:39 +0000 Subject: [PATCH 35/96] Correct ellint_2 for N * PI / 2 phi argument. Fixes https://github.com/boostorg/math/issues/1377. --- include/boost/math/special_functions/ellint_2.hpp | 12 +++++++++--- test/test_ellint_2.hpp | 6 +++++- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/include/boost/math/special_functions/ellint_2.hpp b/include/boost/math/special_functions/ellint_2.hpp index 0cc1fa0944..76aff168de 100644 --- a/include/boost/math/special_functions/ellint_2.hpp +++ b/include/boost/math/special_functions/ellint_2.hpp @@ -136,9 +136,15 @@ BOOST_MATH_GPU_ENABLED T ellint_e_imp(T phi, T k, const Policy& pol) return policies::raise_domain_error("boost::math::ellint_2<%1%>(%1%, %1%)", "The parameter k is out of range, got k = %1%", k, pol); } T cosp = cos(rphi); - T c = 1 / (sinp * sinp); - T cm1 = cosp * cosp / (sinp * sinp); // c - 1 - result = s * ((1 - k2) * ellint_rf_imp(cm1, T(c - k2), c, pol) + k2 * (1 - k2) * ellint_rd(cm1, c, T(c - k2), pol) / 3 + k2 * sqrt(cm1 / (c * (c - k2)))); + T c = sinp * sinp; + if (c > tools::min_value()) + { + c = 1 / c; + T cm1 = cosp * cosp / (sinp * sinp); // c - 1 + result = s * ((1 - k2) * ellint_rf_imp(cm1, T(c - k2), c, pol) + k2 * (1 - k2) * ellint_rd(cm1, c, T(c - k2), pol) / 3 + k2 * sqrt(cm1 / (c * (c - k2)))); + } + else + result = 0; } if (m != 0) { diff --git a/test/test_ellint_2.hpp b/test/test_ellint_2.hpp index 29a73c9961..3b565a493d 100644 --- a/test/test_ellint_2.hpp +++ b/test/test_ellint_2.hpp @@ -97,7 +97,7 @@ void test_spots(T, const char* type_name) // called a second time. // #ifndef TEST_UDT - static const std::array, 28> data1 = {{ + static const std::array, 31> data1 = {{ {{ SC_(0.0), SC_(0.0), SC_(0.0) }}, {{ SC_(-10.0), SC_(0.0), SC_(-10.0) }}, {{ SC_(-1.0), SC_(-1.0), SC_(-0.84147098480789650665250232163029899962256306079837) }}, @@ -127,6 +127,10 @@ void test_spots(T, const char* type_name) {{ -2 * boost::math::constants::pi(), SC_(1.0), SC_(-4.0) }}, {{ -(20 * boost::math::constants::pi()) / 21, SC_(1.0), SC_(-1.85095773382382555307064528472278244309033056100177750424044959888147151556780639619867349409562910808166870200808) }}, {{ -(20 * boost::math::constants::pi()) / 19, SC_(1.0), SC_(-2.16459459028073389414365205908793841951217248335965412335147127404907690285247170629778731438009399303864259295133) }}, + // Test at n * pi / 2, see https://github.com/boostorg/math/issues/1377 + {{ boost::math::constants::half_pi(), SC_(0.25), boost::math::ellint_2(static_cast(0.25)) }}, + {{ 3 * boost::math::constants::half_pi(), SC_(0.25), 3 * boost::math::ellint_2(static_cast(0.25)) }}, + {{ 6 * boost::math::constants::half_pi(), SC_(0.25), 6 * boost::math::ellint_2(static_cast(0.25)) }}, }}; do_test_ellint_e2(data1, type_name, "Elliptic Integral E: Mathworld Data"); From c793e1cd1b4ad84dfdb5bcd551857a951ca4556d Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Thu, 12 Mar 2026 10:22:45 +0000 Subject: [PATCH 36/96] Ellint2: Simplify fix. --- include/boost/math/special_functions/ellint_2.hpp | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/include/boost/math/special_functions/ellint_2.hpp b/include/boost/math/special_functions/ellint_2.hpp index 76aff168de..501d36b66b 100644 --- a/include/boost/math/special_functions/ellint_2.hpp +++ b/include/boost/math/special_functions/ellint_2.hpp @@ -122,7 +122,7 @@ BOOST_MATH_GPU_ENABLED T ellint_e_imp(T phi, T k, const Policy& pol) rphi = constants::half_pi() - rphi; } T k2 = k * k; - if(boost::math::pow<3>(rphi) * k2 / 6 < tools::epsilon() * fabs(rphi)) + if(boost::math::pow<3>(rphi) * k2 / 6 <= tools::epsilon() * fabs(rphi)) { // See http://functions.wolfram.com/EllipticIntegrals/EllipticE2/06/01/03/0001/ result = s * rphi; @@ -136,15 +136,9 @@ BOOST_MATH_GPU_ENABLED T ellint_e_imp(T phi, T k, const Policy& pol) return policies::raise_domain_error("boost::math::ellint_2<%1%>(%1%, %1%)", "The parameter k is out of range, got k = %1%", k, pol); } T cosp = cos(rphi); - T c = sinp * sinp; - if (c > tools::min_value()) - { - c = 1 / c; - T cm1 = cosp * cosp / (sinp * sinp); // c - 1 - result = s * ((1 - k2) * ellint_rf_imp(cm1, T(c - k2), c, pol) + k2 * (1 - k2) * ellint_rd(cm1, c, T(c - k2), pol) / 3 + k2 * sqrt(cm1 / (c * (c - k2)))); - } - else - result = 0; + T c = 1 / (sinp * sinp); + T cm1 = cosp * cosp / (sinp * sinp); // c - 1 + result = s * ((1 - k2) * ellint_rf_imp(cm1, T(c - k2), c, pol) + k2 * (1 - k2) * ellint_rd(cm1, c, T(c - k2), pol) / 3 + k2 * sqrt(cm1 / (c * (c - k2)))); } if (m != 0) { From 372c10b46d87a9862a6b6f723ccbf6d1eb092b47 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sun, 15 Mar 2026 13:47:42 +0000 Subject: [PATCH 37/96] Remove dead code and fix conceptual error. --- include/boost/math/special_functions/beta.hpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index 3f31bea408..c201736ddd 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -19,7 +19,6 @@ #include #include #include -#include #include #include #include @@ -828,7 +827,7 @@ struct ibeta_fraction2_t bN += (m * (b - m) * x) / (a + 2*m - 1); bN += ((a + m) * (a * y - b * x + 1 + m *(2 - x))) / (a + 2*m + 1); - ++m; + m += 1; return boost::math::make_pair(aN, bN); } @@ -1609,7 +1608,6 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b if (local_result != 0) { ibeta_fraction2_t f(a, b, x, y); - ibeta_fraction2_t g(a, b, x, y); boost::math::uintmax_t max_terms = boost::math::policies::get_max_series_iterations(); T local_fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon(), max_terms); if (max_terms >= boost::math::policies::get_max_series_iterations()) From 918bb51d81e74e9db1a5fcc12b7972732e4649ac Mon Sep 17 00:00:00 2001 From: Jacob Hass <72838561+JacobHass8@users.noreply.github.com> Date: Sun, 15 Mar 2026 13:25:13 -0700 Subject: [PATCH 38/96] Reverted changes to ibeta_fraction2_t to use int for m --- include/boost/math/special_functions/beta.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index c201736ddd..5c3a762ca0 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -816,25 +816,25 @@ struct ibeta_fraction2_t { typedef boost::math::pair result_type; - BOOST_MATH_GPU_ENABLED ibeta_fraction2_t(T a_, T b_, T x_, T y_) : a(a_), b(b_), x(x_), y(y_), m(static_cast(0)) {} + BOOST_MATH_GPU_ENABLED ibeta_fraction2_t(T a_, T b_, T x_, T y_) : a(a_), b(b_), x(x_), y(y_), m(0) {} BOOST_MATH_GPU_ENABLED result_type operator()() { T denom = (a + 2 * m - 1); T aN = (m * (a + m - 1) / denom) * ((a + b + m - 1) / denom) * (b - m) * x * x; - T bN = m; + T bN = static_cast(m); bN += (m * (b - m) * x) / (a + 2*m - 1); bN += ((a + m) * (a * y - b * x + 1 + m *(2 - x))) / (a + 2*m + 1); - m += 1; + ++m; return boost::math::make_pair(aN, bN); } private: T a, b, x, y; - T m; + int m; }; // // Evaluate the incomplete beta via the continued fraction representation: From fd70189643a3480928043ae65a1f7ca266498b23 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Tue, 17 Mar 2026 12:06:47 -0500 Subject: [PATCH 39/96] Add user facing domain access --- include/boost/math/interpolators/pchip.hpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/include/boost/math/interpolators/pchip.hpp b/include/boost/math/interpolators/pchip.hpp index 6697166df0..a4ad5a3579 100644 --- a/include/boost/math/interpolators/pchip.hpp +++ b/include/boost/math/interpolators/pchip.hpp @@ -121,6 +121,11 @@ class pchip { } } + std::pair domain() const + { + return impl_->domain(); + } + private: std::shared_ptr> impl_; }; From 9f8cf5caa8adbd4d38ec2f956bc90d490e49bc48 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Tue, 17 Mar 2026 12:06:55 -0500 Subject: [PATCH 40/96] Test domain function --- test/pchip_test.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/pchip_test.cpp b/test/pchip_test.cpp index 6e8ea235bf..5e8fb3c2ac 100644 --- a/test/pchip_test.cpp +++ b/test/pchip_test.cpp @@ -38,6 +38,10 @@ void test_constant() auto pchip_spline = pchip(std::move(x_copy), std::move(y_copy)); //std::cout << "Constant value pchip spline = " << pchip_spline << "\n"; + auto [lo, hi] = pchip_spline.domain(); + CHECK_ULP_CLOSE(lo, 0, 1); + CHECK_ULP_CLOSE(hi, 81, 1); + for (Real t = x[0]; t <= x.back(); t += Real(0.25)) { CHECK_ULP_CLOSE(Real(7), pchip_spline(t), 2); CHECK_ULP_CLOSE(Real(0), pchip_spline.prime(t), 2); From 883b7a1c8efeecf271426a8b81f0d6c7315e29fa Mon Sep 17 00:00:00 2001 From: Philipp Remy Date: Wed, 18 Mar 2026 18:04:50 +0100 Subject: [PATCH 41/96] [CI] Test legacy libraries with static and shared linkage Signed-off-by: Philipp Remy --- .github/workflows/ci.yml | 54 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 53 insertions(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 0750f95e63..578fb5e527 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -603,7 +603,7 @@ jobs: cmake --build . ctest --output-on-failure --no-tests=error - posix-cmake-install-legacy: + posix-cmake-install-legacy-static: strategy: fail-fast: false matrix: @@ -655,6 +655,58 @@ jobs: cmake --build . ctest --output-on-failure --no-tests=error + posix-cmake-install-legacy-shared: + strategy: + fail-fast: false + matrix: + include: + - os: macos-15 + + runs-on: ${{matrix.os}} + + steps: + - uses: actions/checkout@v4 + + - name: Setup Boost + run: | + echo GITHUB_REPOSITORY: $GITHUB_REPOSITORY + LIBRARY=${GITHUB_REPOSITORY#*/} + echo LIBRARY: $LIBRARY + echo "LIBRARY=$LIBRARY" >> $GITHUB_ENV + echo GITHUB_BASE_REF: $GITHUB_BASE_REF + echo GITHUB_REF: $GITHUB_REF + REF=${GITHUB_BASE_REF:-$GITHUB_REF} + REF=${REF#refs/heads/} + echo REF: $REF + BOOST_BRANCH=develop && [ "$REF" == "master" ] && BOOST_BRANCH=master || true + echo BOOST_BRANCH: $BOOST_BRANCH + cd .. + git clone -b $BOOST_BRANCH --depth 1 https://github.com/boostorg/boost.git boost-root + cd boost-root + mkdir -p libs/$LIBRARY + cp -r $GITHUB_WORKSPACE/* libs/$LIBRARY + git submodule update --init tools/boostdep + python tools/boostdep/depinst/depinst.py --git_args "--jobs 3" $LIBRARY + + - name: Configure + run: | + cd ../boost-root + mkdir __build__ && cd __build__ + cmake -DBOOST_INCLUDE_LIBRARIES=$LIBRARY -DBOOST_MATH_BUILD_WITH_LEGACY_FUNCTIONS=ON -DCMAKE_INSTALL_PREFIX=~/.local -DBUILD_SHARED_LIBS=ON .. + + - name: Install + run: | + cd ../boost-root/__build__ + cmake --build . --target install + + - name: Use the installed library (legacy compiled) + run: | + cd ../boost-root/libs/$LIBRARY/test/cmake_install_test_legacy + mkdir __build__ && cd __build__ + cmake -DCMAKE_INSTALL_PREFIX=~/.local .. + cmake --build . + ctest --output-on-failure --no-tests=error + sycl-cmake-test: strategy: fail-fast: false From 4aaeeacdbbfeee26ad0c3bb31d80c46ee1c382f0 Mon Sep 17 00:00:00 2001 From: Philipp Remy Date: Wed, 18 Mar 2026 18:07:56 +0100 Subject: [PATCH 42/96] [legacy] Allow building the legacy libraries in standalone mode Signed-off-by: Philipp Remy --- CMakeLists.txt | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e2d6c055d3..2b9f1b949b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -52,11 +52,12 @@ option(BOOST_MATH_BUILD_WITH_LEGACY_FUNCTIONS "Build the C99 and TR1 compiled li if(BOOST_MATH_BUILD_WITH_LEGACY_FUNCTIONS) include(CheckCXXSourceCompiles) -get_target_property(_config_type Boost::config TYPE) -if(_config_type STREQUAL "INTERFACE_LIBRARY") - get_target_property(_config_inc Boost::config INTERFACE_INCLUDE_DIRECTORIES) -else() - set(_config_inc "") +set(_config_inc "") +if(NOT BOOST_MATH_STANDALONE) + get_target_property(_config_type Boost::config TYPE) + if(_config_type STREQUAL "INTERFACE_LIBRARY") + get_target_property(_config_inc Boost::config INTERFACE_INCLUDE_DIRECTORIES) + endif() endif() set(CMAKE_REQUIRED_INCLUDES "${CMAKE_CURRENT_SOURCE_DIR}/include" ${_config_inc}) check_cxx_source_compiles("#include <${CMAKE_CURRENT_SOURCE_DIR}/config/has_long_double_support.cpp> \n int main() { return 0;}" BOOST_MATH_HAS_LONG_DOUBLE) @@ -130,13 +131,17 @@ list(TRANSFORM TR1_SOURCES APPEND ".cpp") foreach(type IN LISTS types) add_library(boost_math_tr1${type} ${TR1_SOURCES${type}}) add_library(Boost::math_tr1${type} ALIAS boost_math_tr1${type}) - target_link_libraries(boost_math_tr1${type} PUBLIC Boost::config) + if(NOT BOOST_MATH_STANDALONE) + target_link_libraries(boost_math_tr1${type} PUBLIC Boost::config) + endif() target_include_directories(boost_math_tr1${type} PRIVATE src/tr1) target_include_directories(boost_math_tr1${type} PRIVATE include) add_library(boost_math_c99${type} ${C99_SOURCES${type}}) add_library(Boost::math_c99${type} ALIAS boost_math_c99${type}) - target_link_libraries(boost_math_c99${type} PUBLIC Boost::config) + if(NOT BOOST_MATH_STANDALONE) + target_link_libraries(boost_math_c99${type} PUBLIC Boost::config) + endif() target_include_directories(boost_math_c99${type} PRIVATE src/tr1) target_include_directories(boost_math_c99${type} PRIVATE include) From b1523d96b132f11a3dfd72683fe0fd9ac0e29c01 Mon Sep 17 00:00:00 2001 From: Philipp Remy Date: Wed, 18 Mar 2026 18:09:00 +0100 Subject: [PATCH 43/96] [legacy] Handle MinGW and fix build with shared linkage (1) Change MSVC to WIN32: The symbol export rules should also apply if a MinGW toolchain is used. (2) On targets other than MSVC, a build with BUILD_SHARED_LIBS=ON failed because the macro BOOST_SYMBOL_EXPORT was not properly defined. Define it to default visibility for proper symbol export in shared libraries. Signed-off-by: Philipp Remy --- CMakeLists.txt | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2b9f1b949b..7d5f6d270b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -148,9 +148,12 @@ foreach(type IN LISTS types) if(BUILD_SHARED_LIBS) target_compile_definitions(boost_math_tr1${type} PUBLIC BOOST_MATH_TR1_DYN_LINK=1) target_compile_definitions(boost_math_c99${type} PUBLIC BOOST_MATH_TR1_DYN_LINK=1) - if(MSVC) + if(WIN32) target_compile_definitions(boost_math_tr1${type} PRIVATE "BOOST_SYMBOL_EXPORT=__declspec(dllexport)" BOOST_ALL_NO_LIB) target_compile_definitions(boost_math_c99${type} PRIVATE "BOOST_SYMBOL_EXPORT=__declspec(dllexport)" BOOST_ALL_NO_LIB) + else() + target_compile_definitions(boost_math_tr1${type} PRIVATE "BOOST_SYMBOL_EXPORT=__attribute__((visibility(\"default\")))" BOOST_ALL_NO_LIB) + target_compile_definitions(boost_math_c99${type} PRIVATE "BOOST_SYMBOL_EXPORT=__attribute__((visibility(\"default\")))" BOOST_ALL_NO_LIB) endif() endif() From c936a099fc90fc82b8452f92c44044ee939a8e32 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 30 Mar 2026 10:02:40 -0400 Subject: [PATCH 44/96] User must explicitly enable CUDA support via CMake --- include/boost/math/tools/config.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/boost/math/tools/config.hpp b/include/boost/math/tools/config.hpp index 4829ebe36e..26a3a7f73b 100644 --- a/include/boost/math/tools/config.hpp +++ b/include/boost/math/tools/config.hpp @@ -11,7 +11,7 @@ #pragma once #endif -#ifndef __CUDACC_RTC__ +#if !(defined(__CUDACC_RTC__) && defined(BOOST_MATH_ENABLE_NVRTC)) #include @@ -678,7 +678,7 @@ namespace boost{ namespace math{ // CUDA support: // -#ifdef __CUDACC__ +#if defined(__CUDACC__) && defined(BOOST_MATH_ENABLE_CUDA) // We have to get our include order correct otherwise you get compilation failures #include From f6ca42f320e2828e5f1e6afd6a6225a8b8606d1b Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 30 Mar 2026 10:02:52 -0400 Subject: [PATCH 45/96] Fix naked NVRTC macro check --- include/boost/math/special_functions/sign.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/boost/math/special_functions/sign.hpp b/include/boost/math/special_functions/sign.hpp index 4f76522654..5d2dfc2e23 100644 --- a/include/boost/math/special_functions/sign.hpp +++ b/include/boost/math/special_functions/sign.hpp @@ -14,7 +14,7 @@ #pragma once #endif -#ifndef __CUDACC_RTC__ +#ifndef BOOST_MATH_HAS_NVRTC #include #include @@ -234,7 +234,7 @@ BOOST_MATH_GPU_ENABLED T sign(T z) } // namespace math } // namespace boost -#endif // __CUDACC_RTC__ +#endif // BOOST_MATH_HAS_NVRTC #endif // BOOST_MATH_TOOLS_SIGN_HPP From d6368ee45ed1c15eb8cdae059aab2a15fb7edc33 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 30 Mar 2026 10:05:19 -0400 Subject: [PATCH 46/96] Replace naked NVCC checks --- include/boost/math/special_functions/lanczos.hpp | 2 +- include/boost/math/special_functions/next.hpp | 2 +- include/boost/math/tools/config.hpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/include/boost/math/special_functions/lanczos.hpp b/include/boost/math/special_functions/lanczos.hpp index 0ec24bddbf..409e14c5e1 100644 --- a/include/boost/math/special_functions/lanczos.hpp +++ b/include/boost/math/special_functions/lanczos.hpp @@ -2751,7 +2751,7 @@ struct lanczos } // namespace math } // namespace boost -#if !defined(_CRAYC) && !defined(__CUDACC__) && (!defined(__GNUC__) || (__GNUC__ > 3) || ((__GNUC__ == 3) && (__GNUC_MINOR__ > 3))) +#if !defined(_CRAYC) && !defined(BOOST_MATH_ENABLE_CUDA) && (!defined(__GNUC__) || (__GNUC__ > 3) || ((__GNUC__ == 3) && (__GNUC_MINOR__ > 3))) #if ((defined(_M_IX86_FP) && (_M_IX86_FP >= 2)) || defined(__SSE2__) || defined(_M_AMD64) || defined(_M_X64)) && !defined(_MANAGED) && !defined(BOOST_MATH_HAS_GPU_SUPPORT) #include #endif diff --git a/include/boost/math/special_functions/next.hpp b/include/boost/math/special_functions/next.hpp index 74c34f06ad..f73cb2f6a5 100644 --- a/include/boost/math/special_functions/next.hpp +++ b/include/boost/math/special_functions/next.hpp @@ -25,7 +25,7 @@ #include -#if !defined(_CRAYC) && !defined(__CUDACC__) && (!defined(__GNUC__) || (__GNUC__ > 3) || ((__GNUC__ == 3) && (__GNUC_MINOR__ > 3))) +#if !defined(_CRAYC) && !defined(BOOST_MATH_ENABLE_CUDA) && (!defined(__GNUC__) || (__GNUC__ > 3) || ((__GNUC__ == 3) && (__GNUC_MINOR__ > 3))) #if (defined(_M_IX86_FP) && (_M_IX86_FP >= 2)) || defined(__SSE2__) #include "xmmintrin.h" #define BOOST_MATH_CHECK_SSE2 diff --git a/include/boost/math/tools/config.hpp b/include/boost/math/tools/config.hpp index 26a3a7f73b..c97dbdede2 100644 --- a/include/boost/math/tools/config.hpp +++ b/include/boost/math/tools/config.hpp @@ -168,7 +168,7 @@ # define BOOST_MATH_NOINLINE __declspec(noinline) # elif defined(__GNUC__) && __GNUC__ > 3 // Clang also defines __GNUC__ (as 4) -# if defined(__CUDACC__) +# if defined(__CUDACC__) && defined(BOOST_MATH_ENABLE_CUDA) // nvcc doesn't always parse __noinline__, // see: https://svn.boost.org/trac/boost/ticket/9392 # define BOOST_MATH_NOINLINE __attribute__ ((noinline)) From b0cef13e8e63d83363d1cd4c92740ae76512227d Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 30 Mar 2026 10:09:55 -0400 Subject: [PATCH 47/96] Add test set --- test/cuda_jamfile | 3 ++ test/github_issue_1383.cu | 111 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 114 insertions(+) create mode 100644 test/github_issue_1383.cu diff --git a/test/cuda_jamfile b/test/cuda_jamfile index e4d7505da8..b9ebf8d049 100644 --- a/test/cuda_jamfile +++ b/test/cuda_jamfile @@ -9,6 +9,9 @@ project : requirements [ requires cxx14_decltype_auto cxx14_generic_lambdas cxx14_return_type_deduction cxx14_variable_templates cxx14_constexpr ] ; +# Github Issues +run github_issue_1383.cu ; + # Quad run test_exp_sinh_quad_float.cu ; run test_exp_sinh_quad_double.cu ; diff --git a/test/github_issue_1383.cu b/test/github_issue_1383.cu new file mode 100644 index 0000000000..4ec82e785f --- /dev/null +++ b/test/github_issue_1383.cu @@ -0,0 +1,111 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024 - 2026. +// 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) + +#ifdef BOOST_MATH_ENABLE_CUDA +# undef BOOST_MATH_ENABLE_CUDA +#endif // BOOST_MATH_ENABLE_CUDA + +// Purposefully pull in headers that caused errors in the linked issue +#include +#include +#include + +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +#ifdef BOOST_MATH_ENABLE_CUDA +# error "We should not be enabling this ourselves" +#endif // BOOST_MATH_ENABLE_CUDA + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef double float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in, float_type *out, int numElements) +{ + using std::cos; + + if (i < numElements) + { + out[i] = cos(in[i]); + } +} + +/** + * Host main routine + */ +int main() +{ + // 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_vector(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector[i] = rand()/(float_type)RAND_MAX; + } + + // 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_vector.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 vectorAdd 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(); + for(int i = 0; i < numElements; ++i) + results.push_back(std::cos(input_vector[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]) > 10) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED, normal calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + + return 0; +} From 61e689d98435a8999c09eeabc0fb4358bd2b9530 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 30 Mar 2026 10:15:38 -0400 Subject: [PATCH 48/96] Fix deprecated find package in CML --- test/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4715cf379b..9d673adc7d 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -10,13 +10,13 @@ if(HAVE_BOOST_TEST) message(STATUS "Building boost.math with CUDA") - find_package(CUDA REQUIRED) enable_language(CUDA) + find_package(CUDAToolkit REQUIRED) set(CMAKE_CUDA_EXTENSIONS OFF) enable_testing() - boost_test_jamfile(FILE cuda_jamfile LINK_LIBRARIES Boost::math Boost::assert Boost::concept_check Boost::config Boost::core Boost::integer Boost::lexical_cast Boost::multiprecision Boost::predef Boost::random Boost::throw_exception Boost::unit_test_framework ${CUDA_LIBRARIES} INCLUDE_DIRECTORIES ${CUDA_INCLUDE_DIRS} ) + boost_test_jamfile(FILE cuda_jamfile LINK_LIBRARIES Boost::math Boost::assert Boost::concept_check Boost::config Boost::core Boost::integer Boost::lexical_cast Boost::multiprecision Boost::predef Boost::random Boost::throw_exception Boost::unit_test_framework CUDA::cudart COMPILE_DEFINITIONS BOOST_MATH_ENABLE_CUDA=1 ) elseif (BOOST_MATH_ENABLE_NVRTC) From e45257161113be0f0e1fc8a239c5999303702cfa Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 30 Mar 2026 10:17:06 -0400 Subject: [PATCH 49/96] Update CI script --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 578fb5e527..5e91faa944 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -820,7 +820,7 @@ jobs: run: | cd ../boost-root mkdir __build__ && cd __build__ - cmake -DBOOST_INCLUDE_LIBRARIES=$LIBRARY -DBUILD_TESTING=ON -DCMAKE_CUDA_COMPILER=/usr/local/cuda/bin/nvcc -DBOOST_MATH_ENABLE_CUDA=1 -DCMAKE_CUDA_ARCHITECTURES=86 -DCUDA_TOOLKIT_ROOT_DIR=/usr/local/cuda-12.8 .. + cmake -DBOOST_INCLUDE_LIBRARIES=$LIBRARY -DBUILD_TESTING=ON -DCMAKE_CUDA_COMPILER=/usr/local/cuda/bin/nvcc -DBOOST_MATH_ENABLE_CUDA=1 -DCMAKE_CUDA_ARCHITECTURES="75;86" -DCMAKE_CUDA_STANDARD=17 .. - name: Build tests run: | cd ../boost-root/__build__ From 55d82b06d479e8f30a50c69ab2a746c2cd56b46c Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 30 Mar 2026 10:20:17 -0400 Subject: [PATCH 50/96] Fix copy-paste error --- test/github_issue_1383.cu | 1 + 1 file changed, 1 insertion(+) diff --git a/test/github_issue_1383.cu b/test/github_issue_1383.cu index 4ec82e785f..b3544077de 100644 --- a/test/github_issue_1383.cu +++ b/test/github_issue_1383.cu @@ -36,6 +36,7 @@ typedef double float_type; __global__ void cuda_test(const float_type *in, float_type *out, int numElements) { using std::cos; + const int i = blockDim.x * blockIdx.x + threadIdx.x; if (i < numElements) { From ecf9fc1afb8c4c1c834e1e63878776e8b96427de Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 30 Mar 2026 10:39:46 -0400 Subject: [PATCH 51/96] Use more internal functions --- test/github_issue_1383.cu | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/test/github_issue_1383.cu b/test/github_issue_1383.cu index b3544077de..119f5d1d5e 100644 --- a/test/github_issue_1383.cu +++ b/test/github_issue_1383.cu @@ -12,11 +12,13 @@ #include #include #include +#include #include #include #include #include +#include #include "cuda_managed_ptr.hpp" #include "stopwatch.hpp" @@ -63,9 +65,15 @@ int main() cuda_managed_ptr output_vector(numElements); // Initialize the input vectors + // Check some of our numeric_limits for viability + std::mt19937_64 rng {42}; + std::uniform_real_distribution dist(0, boost::math::constants::pi()); + static_assert(boost::math::numeric_limits::is_specialized, "Should be since it's a double"); + static_assert(boost::math::numeric_limits::is_signed, "Should be since it's a double"); + for (int i = 0; i < numElements; ++i) { - input_vector[i] = rand()/(float_type)RAND_MAX; + input_vector[i] = dist(rng); } // Launch the Vector Add CUDA Kernel From 936df78a477e2e72969bac5686420b5811f72646 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 30 Mar 2026 10:42:46 -0400 Subject: [PATCH 52/96] Add additional test this time with math's CUDA support --- test/cuda_jamfile | 1 + test/github_issue_1383_pt_2.cu | 116 +++++++++++++++++++++++++++++++++ 2 files changed, 117 insertions(+) create mode 100644 test/github_issue_1383_pt_2.cu diff --git a/test/cuda_jamfile b/test/cuda_jamfile index b9ebf8d049..f2f2c9aa82 100644 --- a/test/cuda_jamfile +++ b/test/cuda_jamfile @@ -11,6 +11,7 @@ project : requirements # Github Issues run github_issue_1383.cu ; +run github_issue_1383_pt_2.cu ; # Quad run test_exp_sinh_quad_float.cu ; diff --git a/test/github_issue_1383_pt_2.cu b/test/github_issue_1383_pt_2.cu new file mode 100644 index 0000000000..8e746abfce --- /dev/null +++ b/test/github_issue_1383_pt_2.cu @@ -0,0 +1,116 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024 - 2026. +// 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) + +// Purposefully pull in headers that caused errors in the linked issue +#include +#include +#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 *in, float_type *out, int numElements) +{ + using std::cos; + const int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = cos(in[i]); + if (out[i] > boost::math::numeric_limits::max() || !boost::math::numeric_limits::is_signed) + { + __trap(); + } + } +} + +/** + * Host main routine + */ +int main() +{ + // 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_vector(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + // Initialize the input vectors + // Check some of our numeric_limits for viability + std::mt19937_64 rng {42}; + std::uniform_real_distribution dist(0, boost::math::constants::pi()); + static_assert(boost::math::numeric_limits::is_specialized, "Should be since it's a double"); + static_assert(boost::math::numeric_limits::is_signed, "Should be since it's a double"); + + for (int i = 0; i < numElements; ++i) + { + input_vector[i] = dist(rng); + } + + // 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_vector.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 vectorAdd 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(); + for(int i = 0; i < numElements; ++i) + results.push_back(std::cos(input_vector[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]) > 10) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED, normal calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + + return 0; +} From 1d9e9279da978bc3cd68f5f4f25d3d5cb9f69945 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Mon, 30 Mar 2026 10:45:00 -0400 Subject: [PATCH 53/96] Change macro definition --- include/boost/math/tools/config.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/math/tools/config.hpp b/include/boost/math/tools/config.hpp index c97dbdede2..40346c9b0f 100644 --- a/include/boost/math/tools/config.hpp +++ b/include/boost/math/tools/config.hpp @@ -774,7 +774,7 @@ BOOST_MATH_GPU_ENABLED constexpr T gpu_safe_max(const T& a, const T& b) { return # define BOOST_MATH_STATIC_LOCAL_VARIABLE # else # define BOOST_MATH_INLINE_CONSTEXPR constexpr -# define BOOST_MATH_STATIC constexpr +# define BOOST_MATH_STATIC static # define BOOST_MATH_STATIC_LOCAL_VARIABLE static # endif #endif From 727a07708400c995488d72d06817cab6483d2fe0 Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Sun, 19 Apr 2026 19:37:55 +0200 Subject: [PATCH 54/96] Update gh actions checkout version --- .github/workflows/ci.yml | 34 +++++++++++++++++----------------- .github/workflows/codecov.yml | 4 ++-- 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 5e91faa944..8a6bbcf460 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -28,7 +28,7 @@ jobs: standard: [ c++20 ] suite: [ github_ci_block_1, github_ci_block_2 ] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Set TOOLSET @@ -86,7 +86,7 @@ jobs: compiler: [ g++-13, clang++-19 ] standard: [ c++14, c++17, c++20, c++23 ] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Set TOOLSET @@ -145,7 +145,7 @@ jobs: standard: [ 14, 17, 20 ] suite: [ github_ci_block_1, github_ci_block_2 ] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Checkout main boost @@ -188,7 +188,7 @@ jobs: standard: [ 14 ] suite: [ github_ci_block_1, github_ci_block_2 ] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Checkout main boost @@ -231,7 +231,7 @@ jobs: standard: [ 14, 17, latest ] suite: [ github_ci_block_1, github_ci_block_2 ] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Checkout main boost @@ -271,7 +271,7 @@ jobs: env: TOOLSET: gcc steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Install Cygwin @@ -305,7 +305,7 @@ jobs: matrix: compiler: [ g++-13 ] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Add repository @@ -345,7 +345,7 @@ jobs: matrix: compiler: [ clang++-19 ] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Add repository @@ -387,7 +387,7 @@ jobs: standard: [ c++14, c++17, c++20, c++23 ] suite: [ github_ci_block_1, github_ci_block_2 ] steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: fetch-depth: '0' - name: Set TOOLSET @@ -448,7 +448,7 @@ jobs: runs-on: ${{matrix.os}} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 - name: Install packages if: matrix.install @@ -498,7 +498,7 @@ jobs: runs-on: ${{matrix.os}} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 - name: Use library with add_subdirectory (header-only) run: | @@ -520,7 +520,7 @@ jobs: runs-on: ${{matrix.os}} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 - name: Setup Boost run: | @@ -561,7 +561,7 @@ jobs: runs-on: ${{matrix.os}} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 - name: Setup Boost run: | @@ -613,7 +613,7 @@ jobs: runs-on: ${{matrix.os}} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 - name: Setup Boost run: | @@ -665,7 +665,7 @@ jobs: runs-on: ${{matrix.os}} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 - name: Setup Boost run: | @@ -733,7 +733,7 @@ jobs: printenv >> $GITHUB_ENV - name: checkout project code - uses: actions/checkout@v4 + uses: actions/checkout@v6 - name: Install Packages run: | @@ -791,7 +791,7 @@ jobs: echo "Installed cuda version is: ${{steps.cuda-toolkit.outputs.cuda}}"+ echo "Cuda install location: ${{steps.cuda-toolkit.outputs.CUDA_PATH}}" nvcc -V - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 - name: Install Packages run: | diff --git a/.github/workflows/codecov.yml b/.github/workflows/codecov.yml index 6c76d1f238..f867c38311 100644 --- a/.github/workflows/codecov.yml +++ b/.github/workflows/codecov.yml @@ -92,7 +92,7 @@ jobs: fi git config --global pack.threads 0 - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 with: # For coverage builds fetch the whole history, else only 1 commit using a 'fake ternary' fetch-depth: ${{ matrix.coverage && '0' || '1' }} @@ -106,7 +106,7 @@ jobs: restore-keys: ${{matrix.os}}-${{matrix.container}}-${{matrix.compiler}}- - name: Fetch Boost.CI - uses: actions/checkout@v4 + uses: actions/checkout@v6 with: repository: boostorg/boost-ci ref: master From e1108f6af935fa9b97a04f4be0b5269aad451ba6 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Thu, 23 Apr 2026 11:10:54 -0500 Subject: [PATCH 55/96] Update boost math version --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 7d5f6d270b..616c21e27f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -5,7 +5,7 @@ cmake_minimum_required(VERSION 3.8...3.16) -project(boost_math VERSION 1.91.0 LANGUAGES CXX) +project(boost_math VERSION 1.92.0 LANGUAGES CXX) add_library(boost_math INTERFACE) From 2c4580c119cc40ca5341d634bb579a77f6d88ce6 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Tue, 28 Apr 2026 11:27:23 -0400 Subject: [PATCH 56/96] Allow fallthrough of Android to 128-bit float --- include/boost/math/special_functions/detail/fp_traits.hpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/include/boost/math/special_functions/detail/fp_traits.hpp b/include/boost/math/special_functions/detail/fp_traits.hpp index 7cdb12b188..368df69f3b 100644 --- a/include/boost/math/special_functions/detail/fp_traits.hpp +++ b/include/boost/math/special_functions/detail/fp_traits.hpp @@ -305,9 +305,12 @@ template<> struct fp_traits_non_native // long double (>64 bits), x86 and x64 ----------------------------------------- -#elif defined(__i386) || defined(__i386__) || defined(_M_IX86) \ +// Andorid is a known exception because on x64 it reports LDBL_MANT_DIG == 113 +// See: https://github.com/boostorg/math/issues/1389 + +#elif (defined(__i386) || defined(__i386__) || defined(_M_IX86) \ || defined(__amd64) || defined(__amd64__) || defined(_M_AMD64) \ - || defined(__x86_64) || defined(__x86_64__) || defined(_M_X64) + || defined(__x86_64) || defined(__x86_64__) || defined(_M_X64)) && !defined(__ANDROID__) // Intel extended double precision format (80 bits) From 8d8f4e4b22a02943023de47283d4948e10e20d9e Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Tue, 28 Apr 2026 13:28:33 -0400 Subject: [PATCH 57/96] Match drone.sh to boost.sh --- .drone/drone.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.drone/drone.sh b/.drone/drone.sh index 86dbf129f5..7fc9d61b11 100644 --- a/.drone/drone.sh +++ b/.drone/drone.sh @@ -51,8 +51,8 @@ echo '==================================> BEFORE_SCRIPT' echo '==================================> SCRIPT' echo "using $TOOLSET : : $COMPILER : -std=$CXXSTD $OPTIONS ;" > ~/user-config.jam -(cd libs/config/test && ../../../b2 config_info_travis_install toolset=$TOOLSET && ./config_info_travis) -(cd libs/math/test && ../../../b2 toolset=$TOOLSET $TEST_SUITE define=$DEFINE) +(cd libs/config/test && ../../../b2 print_config_info print_math_info toolset=$TOOLSET) +(cd libs/math/test && ../../../b2 -j3 toolset=$TOOLSET $TEST_SUITE) echo '==================================> AFTER_SUCCESS' From 74406e5074fb8868d8496f6a478b28915cec8c19 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Tue, 28 Apr 2026 16:02:01 -0400 Subject: [PATCH 58/96] Simplify logic that we are asserting anyway --- .../special_functions/detail/fp_traits.hpp | 19 ++----------------- 1 file changed, 2 insertions(+), 17 deletions(-) diff --git a/include/boost/math/special_functions/detail/fp_traits.hpp b/include/boost/math/special_functions/detail/fp_traits.hpp index 368df69f3b..7362d0fada 100644 --- a/include/boost/math/special_functions/detail/fp_traits.hpp +++ b/include/boost/math/special_functions/detail/fp_traits.hpp @@ -270,13 +270,7 @@ template<> struct fp_traits_non_native // long double (64 bits) ------------------------------------------------------- -#if defined(BOOST_NO_INT64_T) || defined(BOOST_NO_INCLASS_MEMBER_INITIALIZATION)\ - || defined(BOOST_BORLANDC) || defined(__CODEGEAR__) || (defined(__APPLE__) && defined(__aarch64__)) || defined(_MSC_VER)\ - || (defined(__GNUC__) && defined(__aarch64__) && defined(_WIN32))\ - || (defined(__GNUC__) && (defined(__arm__) || defined(__thumb__)))\ - || defined(__SYCL_DEVICE_ONLY__) || (LDBL_MANT_DIG == 53) - -static_assert(LDBL_MANT_DIG == 53, "Oops, assumption that long double is a 64-bit quantity is incorrect!!"); +#if LDBL_MANT_DIG == 53 || defined(__SYCL_DEVICE_ONLY__) template<> struct fp_traits_non_native { @@ -305,17 +299,10 @@ template<> struct fp_traits_non_native // long double (>64 bits), x86 and x64 ----------------------------------------- -// Andorid is a known exception because on x64 it reports LDBL_MANT_DIG == 113 -// See: https://github.com/boostorg/math/issues/1389 - -#elif (defined(__i386) || defined(__i386__) || defined(_M_IX86) \ - || defined(__amd64) || defined(__amd64__) || defined(_M_AMD64) \ - || defined(__x86_64) || defined(__x86_64__) || defined(_M_X64)) && !defined(__ANDROID__) +#elif LDBL_MANT_DIG == 64 // Intel extended double precision format (80 bits) -static_assert(LDBL_MANT_DIG == 64, "Oops, assumption that long double is an 80-bit quantity is incorrect!!"); - template<> struct fp_traits_non_native { @@ -443,8 +430,6 @@ struct fp_traits_non_native // IEEE extended double precision format with 15 exponent bits (128 bits) -static_assert(LDBL_MANT_DIG == 113, "Oops, assumption that long double is a 128-bit quantity is incorrect!!"); - template<> struct fp_traits_non_native { From 03ba15aa427ab0b7a445d1fbc4fafff8f292760c Mon Sep 17 00:00:00 2001 From: Charles Bicari Date: Thu, 30 Apr 2026 15:49:42 -0400 Subject: [PATCH 59/96] constants: Protect F which is defined as a macro by some frameworks For instance, ESP32 platform globally defines: esp32/WString.h: #define F(string_literal) (FPSTR(PSTR(string_literal))) --- include/boost/math/constants/constants.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/math/constants/constants.hpp b/include/boost/math/constants/constants.hpp index bdc0a74e27..d162af7475 100644 --- a/include/boost/math/constants/constants.hpp +++ b/include/boost/math/constants/constants.hpp @@ -129,7 +129,7 @@ namespace boost{ namespace math { initializer() { - F(); + (F)(); } void force_instantiate()const{} }; From f5b16da1628310df9e20f65427584cf71c05aff8 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Wed, 6 May 2026 08:18:01 -0400 Subject: [PATCH 60/96] Update drone runners OSes to not require use of toolchain PPA --- .drone.star | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/.drone.star b/.drone.star index 691246e337..1878d18723 100644 --- a/.drone.star +++ b/.drone.star @@ -47,22 +47,22 @@ def main(ctx): result.append(linux_cxx("Ubuntu g++-8 " + cxx + " " + suite, "g++-8", packages="g++-8", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-8', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) result.append(linux_cxx("Ubuntu g++-9 " + cxx + " " + suite, "g++-9", packages="g++-9", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-9', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) for cxx in clang_6_stds: - result.append(linux_cxx("Ubuntu clang++-6 " + cxx + " " + suite, "clang++-6.0", packages="clang-6.0", llvm_os="xenial", llvm_ver="6.0", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-6.0', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) - result.append(linux_cxx("Ubuntu clang++-7 " + cxx + " " + suite, "clang++-7", packages="clang-7", llvm_os="xenial", llvm_ver="7", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-7', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) - result.append(linux_cxx("Ubuntu clang++-8 " + cxx + " " + suite, "clang++-8", packages="clang-8", llvm_os="xenial", llvm_ver="8", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-8', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) - result.append(linux_cxx("Ubuntu clang++-9 " + cxx + " " + suite, "clang++-9", packages="clang-9", llvm_os="xenial", llvm_ver="9", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-9', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu clang++-6 " + cxx + " " + suite, "clang++-6.0", packages="clang-6.0", llvm_os="focal", llvm_ver="6.0", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-6.0', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu clang++-7 " + cxx + " " + suite, "clang++-7", packages="clang-7", llvm_os="focal", llvm_ver="7", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-7', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu clang++-8 " + cxx + " " + suite, "clang++-8", packages="clang-8", llvm_os="focal", llvm_ver="8", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-8', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu clang++-9 " + cxx + " " + suite, "clang++-9", packages="clang-9", llvm_os="focal", llvm_ver="9", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-9', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) for cxx in gnu_9_stds: - result.append(linux_cxx("Ubuntu g++-10 " + cxx + " " + suite, "g++-10", packages="g++-10", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-10', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) - result.append(linux_cxx("Ubuntu g++-11 " + cxx + " " + suite, "g++-11", packages="g++-11", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-11', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu g++-10 " + cxx + " " + suite, "g++-10", packages="g++-10", buildtype="boost", image="cppalliance/droneubuntu2204:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-10', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu g++-11 " + cxx + " " + suite, "g++-11", packages="g++-11", buildtype="boost", image="cppalliance/droneubuntu2204:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-11', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) result.append(linux_cxx("Ubuntu g++-12 " + cxx + " " + suite, "g++-12", packages="g++-12", buildtype="boost", image="cppalliance/droneubuntu2204:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-12', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) result.append(linux_cxx("Ubuntu clang++-14 " + cxx + " " + suite, "clang++-14", packages="clang-14", llvm_os="jammy", llvm_ver="14", buildtype="boost", image="cppalliance/droneubuntu2204:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-14', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) result.append(linux_cxx("Ubuntu clang++-15 " + cxx + " " + suite, "clang++-15", packages="clang-15", llvm_os="jammy", llvm_ver="15", buildtype="boost", image="cppalliance/droneubuntu2204:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-15', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) for cxx in clang_10_stds: result.append(linux_cxx("Ubuntu clang++-10 " + cxx + " " + suite, "clang++-10", packages="clang-10", llvm_os="xenial", llvm_ver="10", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-10', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) for cxx in gnu_non_native: - result.append(linux_cxx("Ubuntu g++ s390s " + cxx + " " + suite, "g++", packages="g++", buildtype="boost", image="cppalliance/droneubuntu2204:multiarch", arch="s390x", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu g++ s390s " + cxx + " " + suite, "g++", packages="g++", buildtype="boost", image="cppalliance/droneubuntu2404:multiarch", arch="s390x", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) for cxx in gnu_non_native: - result.append(linux_cxx("Ubuntu g++ ARM64" + cxx + " " + suite, "g++", packages="g++", buildtype="boost", image="cppalliance/droneubuntu2204:multiarch", arch="arm64", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu g++ ARM64" + cxx + " " + suite, "g++", packages="g++", buildtype="boost", image="cppalliance/droneubuntu2404:multiarch", arch="arm64", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) for cxx in gnu_non_native: result.append(osx_cxx("M1 Clang " + cxx + " " + suite, "clang++", buildscript="drone", buildtype="boost", xcode_version="14.1", environment={'TOOLSET': 'clang', 'CXXSTD': cxx, 'TEST_SUITE': suite, 'DEFINE': 'BOOST_MATH_NO_REAL_CONCEPT_TESTS,BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS,BOOST_MATH_MULTI_ARCH_CI_RUN', }, globalenv=globalenv)) for suite in gcc13_things_to_test: From a23376aba314af2dc65b906b00029a15c3515fa2 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Wed, 6 May 2026 11:41:42 -0400 Subject: [PATCH 61/96] Try bionic for Clangs 6,7,8 --- .drone.star | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.drone.star b/.drone.star index 1878d18723..17b2006d7a 100644 --- a/.drone.star +++ b/.drone.star @@ -47,9 +47,9 @@ def main(ctx): result.append(linux_cxx("Ubuntu g++-8 " + cxx + " " + suite, "g++-8", packages="g++-8", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-8', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) result.append(linux_cxx("Ubuntu g++-9 " + cxx + " " + suite, "g++-9", packages="g++-9", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-9', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) for cxx in clang_6_stds: - result.append(linux_cxx("Ubuntu clang++-6 " + cxx + " " + suite, "clang++-6.0", packages="clang-6.0", llvm_os="focal", llvm_ver="6.0", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-6.0', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) - result.append(linux_cxx("Ubuntu clang++-7 " + cxx + " " + suite, "clang++-7", packages="clang-7", llvm_os="focal", llvm_ver="7", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-7', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) - result.append(linux_cxx("Ubuntu clang++-8 " + cxx + " " + suite, "clang++-8", packages="clang-8", llvm_os="focal", llvm_ver="8", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-8', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu clang++-6 " + cxx + " " + suite, "clang++-6.0", packages="clang-6.0", llvm_os="bionic", llvm_ver="6.0", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-6.0', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu clang++-7 " + cxx + " " + suite, "clang++-7", packages="clang-7", llvm_os="bionic", llvm_ver="7", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-7', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) + result.append(linux_cxx("Ubuntu clang++-8 " + cxx + " " + suite, "clang++-8", packages="clang-8", llvm_os="bionic", llvm_ver="8", buildtype="boost", image="cppalliance/droneubuntu1804:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-8', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) result.append(linux_cxx("Ubuntu clang++-9 " + cxx + " " + suite, "clang++-9", packages="clang-9", llvm_os="focal", llvm_ver="9", buildtype="boost", image="cppalliance/droneubuntu2004:1", environment={'TOOLSET': 'clang', 'COMPILER': 'clang++-9', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) for cxx in gnu_9_stds: result.append(linux_cxx("Ubuntu g++-10 " + cxx + " " + suite, "g++-10", packages="g++-10", buildtype="boost", image="cppalliance/droneubuntu2204:1", environment={'TOOLSET': 'gcc', 'COMPILER': 'g++-10', 'CXXSTD': cxx, 'TEST_SUITE': suite, }, globalenv=globalenv)) From f49ab98b02d4cbc8145f6395172074780ae263cd Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sat, 9 May 2026 17:06:29 +0100 Subject: [PATCH 62/96] InverseGaussian: Improve quantile behaviour. (#1395) * InverseGaussian: Improve quantile behaviour. Provide a better initial guess for the root, and bracket the root to avoid shooting off to infinity. Fixes: https://github.com/boostorg/math/issues/1392 Refs: https://github.com/scipy/scipy/issues/25096 * InverseGaussian: Make bracket_root_towards_min/max non-recursive on CUDA. --- .../math/distributions/inverse_gaussian.hpp | 63 +++++++++++-------- include/boost/math/tools/roots.hpp | 10 +++ test/test_inverse_gaussian.cpp | 7 +++ 3 files changed, 54 insertions(+), 26 deletions(-) diff --git a/include/boost/math/distributions/inverse_gaussian.hpp b/include/boost/math/distributions/inverse_gaussian.hpp index 20d3b6bdd5..795ad23d4a 100644 --- a/include/boost/math/distributions/inverse_gaussian.hpp +++ b/include/boost/math/distributions/inverse_gaussian.hpp @@ -52,6 +52,7 @@ #include #include #include // for erf/erfc. +#include // for erf/erfc. #include #include #include @@ -298,7 +299,7 @@ struct inverse_gaussian_quantile_complement_functor namespace detail { template - BOOST_MATH_GPU_ENABLED inline RealType guess_ig(RealType p, RealType mu = 1, RealType lambda = 1) + BOOST_MATH_GPU_ENABLED inline RealType guess_ig(RealType p, RealType q, RealType mu, RealType lambda) { // guess at random variate value x for inverse gaussian quantile. BOOST_MATH_STD_USING using boost::math::policies::policy; @@ -323,27 +324,16 @@ namespace detail x = mu * exp(quantile(n01, p) / sqrt(phi) - 1/(2 * phi)); } else - { // phi < 2 so much less symmetrical with long tail, - // so use gamma distribution as an approximation. - using boost::math::gamma_distribution; - - // Define the distribution, using gamma_nooverflow: - using gamma_nooverflow = gamma_distribution; - - gamma_nooverflow g(static_cast(0.5), static_cast(1.)); - - // R qgamma(0.2, 0.5, 1) = 0.0320923 - RealType qg = quantile(complement(g, p)); - x = lambda / (qg * 2); - // - if (x > mu/2) // x > mu /2? - { // x too large for the gamma approximation to work well. - //x = qgamma(p, 0.5, 1.0); // qgamma(0.270614, 0.5, 1) = 0.05983807 - RealType q = quantile(g, p); - // x = mu * exp(q * static_cast(0.1)); // Said to improve at high p - // x = mu * x; // Improves at high p? - x = mu * exp(q / sqrt(phi) - 1/(2 * phi)); - } + { + // Construct a gamma distribution with the same mean and standard deviation, and hope it + // get's us in more or less the right ballpark: + inverse_gaussian_distribution ig(mu, lambda); + RealType m = mean(ig); + RealType v = standard_deviation(ig); + RealType k = m * m / v; + RealType theta = v / m; + gamma_distribution n01(k, theta); + x = p < q ? quantile(n01, p) : quantile(complement(n01, q)); } return x; } // guess_ig @@ -379,7 +369,7 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const inverse_gaussian_distribut return result; // infinity; } - RealType guess = detail::guess_ig(p, dist.mean(), dist.scale()); + RealType guess = detail::guess_ig(p, RealType(1 - p), dist.mean(), dist.scale()); using boost::math::tools::max_value; RealType min = static_cast(0); // Minimum possible value is bottom of range of distribution. @@ -390,8 +380,19 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const inverse_gaussian_distribut int get_digits = policies::digits();// get digits from policy, boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); // and max iterations. using boost::math::tools::newton_raphson_iterate; + inverse_gaussian_quantile_functor func(dist, p); + // + // Our guess is often not very good, so lets bracket the root just in case: + // + RealType f0 = boost::math::get<0>(func(guess)); + if (f0 < 0) + tools::detail::bracket_root_towards_max(func, guess, f0, min, max, max_iter); + else + tools::detail::bracket_root_towards_min(func, guess, f0, min, max, max_iter); + if((guess < min) || (guess > max)) + guess = (min + max) / 2; result = - newton_raphson_iterate(inverse_gaussian_quantile_functor(dist, p), guess, min, max, get_digits, max_iter); + newton_raphson_iterate(func, guess, min, max, get_digits, max_iter); 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 @@ -455,7 +456,7 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type(); boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); using boost::math::tools::newton_raphson_iterate; - result = newton_raphson_iterate(inverse_gaussian_quantile_complement_functor(c.dist, q), guess, min, max, get_digits, max_iter); + inverse_gaussian_quantile_complement_functor func(c.dist, q); + RealType f0 = boost::math::get<0>(func(guess)); + if (f0 > 0) + tools::detail::bracket_root_towards_max(func, guess, f0, min, max, max_iter); + else + tools::detail::bracket_root_towards_min(func, guess, f0, min, max, max_iter); + if ((guess < min) || (guess > max)) + guess = (min + max) / 2; + result = + newton_raphson_iterate(func, guess, min, max, get_digits, max_iter); + result = newton_raphson_iterate(func, guess, min, max, get_digits, max_iter); 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 diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index b0b0fc246c..af9b76d998 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -456,8 +456,13 @@ namespace detail { if (count) { max = guess; + // + // We can't have this extra recursive tidy up step on CUDA: + // +#if !defined(BOOST_MATH_ENABLE_CUDA) && !defined(BOOST_MATH_ENABLE_SYCL) if (multiplier > 16) return (guess0 - guess) + bracket_root_towards_min(f, guess, f_current, min, max, count); +#endif } return guess0 - (max + min) / 2; } @@ -520,8 +525,13 @@ namespace detail { if (count) { min = guess; + // + // We can't have this extra recursive tidy up step on CUDA: + // +#if !defined(BOOST_MATH_ENABLE_CUDA) && !defined(BOOST_MATH_ENABLE_SYCL) if (multiplier > 16) return (guess0 - guess) + bracket_root_towards_max(f, guess, f_current, min, max, count); +#endif } return guess0 - (max + min) / 2; } diff --git a/test/test_inverse_gaussian.cpp b/test/test_inverse_gaussian.cpp index 3825da9397..0fe64223dc 100644 --- a/test/test_inverse_gaussian.cpp +++ b/test/test_inverse_gaussian.cpp @@ -252,6 +252,13 @@ BOOST_AUTO_TEST_CASE( test_main ) // but better than R that gives up completely! // R Error in SuppDists::qinverse_gaussian(4.87914430108515e-219, 1, 1) : Infinite value in NewtonRoot() + inverse_gaussian w_big(66.99652081); + BOOST_CHECK_CLOSE_FRACTION( + quantile(w_big, 0.97969), 591.567880739988823, 15 * tolfeweps); + BOOST_CHECK_CLOSE_FRACTION( + quantile(complement(w_big, 1 - 0.97969)), 591.567880739988823, 10 * tolfeweps); + + BOOST_CHECK_CLOSE_FRACTION( pdf(w11, 0.5), static_cast(0.87878257893544476), tolfeweps); // pdf BOOST_CHECK_CLOSE_FRACTION( From 51efb17098ddfd209788ae985f3716057fc07e24 Mon Sep 17 00:00:00 2001 From: Charles Bicari Date: Mon, 11 May 2026 15:50:00 -0400 Subject: [PATCH 63/96] constants: Protect F which is defined as a macro by some frameworks - part 2 This is the sibling of : https://github.com/boostorg/math/pull/1391 The same issue occurs with unprotected macro. Had forgotten this one in previous correction. --- include/boost/math/constants/constants.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/math/constants/constants.hpp b/include/boost/math/constants/constants.hpp index d162af7475..3daf2d266e 100644 --- a/include/boost/math/constants/constants.hpp +++ b/include/boost/math/constants/constants.hpp @@ -151,7 +151,7 @@ namespace boost{ namespace math { initializer() { - F(); + (F)(); } void force_instantiate()const{} }; From f724aec8286f68ed870d506f632d5e76046b6ec1 Mon Sep 17 00:00:00 2001 From: Daniel Schmitz <40656107+dschmitz89@users.noreply.github.com> Date: Fri, 22 May 2026 18:49:30 +0200 Subject: [PATCH 64/96] ENH: add parameter finder for degrees or freedom for students_t distribution (#1385) * ENH: add parameter finder for degrees or freedom * Add early exit if approximation is precise to machine precision * Loosen tolerance and add tests in the tails * Modularize and simplify * Use overflow error instead of inf * Relax all test tolerances * Should use a linter .. * Relax test tolerance for float precision and remove ill conditioned test with df=1e6 * Fix indentation * Loosen tolerance for very high degrees of freedom * Test if approximation is a good hint, otherwise fall back to 0.01 * Add test coverage of all cases for approximation * Do not use symmetry for better numerical accuracy in the tails * Add edge case tests * Exclude float from testing for extreme values * Update test comment * Replace non ASCII character.. * Simplify using boost's own relative error function * Add early exits for df=1 and df=inf Co-authored-by: Copilot * Try to fix CI carnage Co-authored-by: Copilot * Simplify analytical cases and raise overflow error for normal case * Rename to find_degrees_of_freedom --------- Co-authored-by: Copilot --- doc/distributions/students_t.qbk | 13 ++ .../boost/math/distributions/students_t.hpp | 168 +++++++++++++++++- test/test_students_t.cpp | 144 +++++++++++++++ 3 files changed, 316 insertions(+), 9 deletions(-) diff --git a/doc/distributions/students_t.qbk b/doc/distributions/students_t.qbk index 3396048c5f..2ff2812cbb 100644 --- a/doc/distributions/students_t.qbk +++ b/doc/distributions/students_t.qbk @@ -29,6 +29,11 @@ RealType beta, RealType sd, RealType hint = 100); + + // degrees of freedom inversion from a quantile and probability: + BOOST_MATH_GPU_ENABLED static RealType find_degrees_of_freedom( + RealType t, + RealType p); }; }} // namespaces @@ -106,6 +111,13 @@ For more information on this function see the [@http://www.itl.nist.gov/div898/handbook/prc/section2/prc222.htm NIST Engineering Statistics Handbook]. + BOOST_MATH_GPU_ENABLED static RealType find_degrees_of_freedom( + RealType t, + RealType p); + +Returns the degrees of freedom [nu] such that CDF(/x/; [nu]) = /p/. +Requires 0 < /p/ < 1 and /x/ [ne] 0, otherwise calls __domain_error. + [h4 Non-member Accessors] All the [link math_toolkit.dist_ref.nmp usual non-member accessor functions] that are generic to all @@ -164,6 +176,7 @@ without the subtraction implied above.]] [[skewness][if (v > 3) 0 else NaN ]] [[kurtosis][if (v > 4) 3 * (v - 2) \/ (v - 4) else NaN]] [[kurtosis excess][if (v > 4) 6 \/ (df - 4) else NaN]] +[[find_degrees_of_freedom(t, p)][Uses a 2nd-order Edgeworth expansion as initial guess for a numerical root finder]] ] If the moment index /k/ is less than /v/, then the moment is undefined. diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index 39f20d6e41..bc1c7139eb 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -19,9 +19,11 @@ #include #include // for ibeta(a, b, x). #include +#include #include #include -#include +#include +#include #include #ifdef _MSC_VER @@ -58,6 +60,10 @@ class students_t_distribution RealType sd, RealType hint = 100); + BOOST_MATH_GPU_ENABLED static RealType find_degrees_of_freedom( + RealType t, + RealType p); + private: // Data member: RealType df_; // degrees of freedom is a real number > 0 or +infinity. @@ -281,6 +287,18 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type +BOOST_MATH_GPU_ENABLED inline bool analytical_df_if_cdf_matches(const Distribution& dist, RealType x, RealType p) +{ + RealType cdf_val = cdf(dist, x); + return epsilon_difference(cdf_val, p) < 4; +} + +// Minimum degrees-of-freedom used as the warm-start fallback when the +// Edgeworth approximation yields no valid positive root or is inaccurate +constexpr double df_hint_fallback = 0.01; + // // Functors for finding degrees of freedom: // @@ -308,6 +326,93 @@ struct sample_size_func RealType alpha, beta, ratio; }; +template +struct cdf_to_df_func +{ + BOOST_MATH_GPU_ENABLED cdf_to_df_func(RealType x_val, RealType p_val, bool c) + : x(x_val), p(p_val), comp(c) {} + + BOOST_MATH_GPU_ENABLED RealType operator()(const RealType& df) + { + students_t_distribution t(df); + return comp ? + RealType(p - cdf(complement(t, x))) + : RealType(cdf(t, x) - p); + } + RealType x, p; + bool comp; +}; + +// +// Shared root-finding helper used by find_degrees_of_freedom. +// +template +BOOST_MATH_GPU_ENABLED RealType solve_for_degrees_of_freedom( + Func f, + RealType hint, + bool rising, + const char* function, + Policy const&) +{ + tools::eps_tolerance tol(policies::digits()); + boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); + boost::math::pair r = tools::bracket_and_solve_root( + f, hint, RealType(2), rising, tol, max_iter, Policy()); + RealType result = r.first + (r.second - r.first) / 2; + if (max_iter >= policies::get_max_root_iterations()) + { + return policies::raise_evaluation_error(function, // LCOV_EXCL_LINE + "Unable to locate solution in a reasonable time: either there is no answer to " // LCOV_EXCL_LINE + "the degrees of freedom or the answer is infinite. Current best guess is %1%", // LCOV_EXCL_LINE + result, Policy()); // LCOV_EXCL_LINE + } + return result; +} + +// +// Edgeworth warm-start for find_degrees_of_freedom(t, p). +// On success writes a df estimate into 'result' and returns true. +// Returns false when no positive root is found; 'result' is left unchanged. +// +template +BOOST_MATH_GPU_ENABLED bool approximate_df_with_edgeworth_expansion(RealType x, RealType p, RealType& result) +{ + BOOST_MATH_STD_USING + // F(x; nu) ~ cdf_normal(x) - pdf_normal(x)*(x + x^3)/(4*nu) + pdf_normal(x)*(3x + 5x^3 + 7x^5 - 3x^7)/(96*nu^2) + // Substituting u = 1/nu and setting F(x; nu) = p gives b*u^2 - a*u + c = 0, where: + // a = pdf_normal(x) * (x + x^3) / 4 + // b = pdf_normal(x) * (3x + 5x^3 + 7x^5 - 3x^7) / 96 + // c = cdf_normal(x) - p + normal_distribution std_normal(0, 1); + RealType pdf_val = pdf(std_normal, x); + RealType cdf_val = cdf(std_normal, x); + + RealType x2 = x * x; + RealType x3 = x2 * x; + RealType x5 = x3 * x2; + RealType x7 = x5 * x2; + + RealType a = pdf_val * (x + x3) / 4; + RealType b = pdf_val * (3 * x + 5 * x3 + 7 * x5 - 3 * x7) / 96; + RealType c = cdf_val - p; + + RealType discriminant = a * a - 4 * b * c; + if (discriminant >= 0 && b != 0) + { + RealType sqrt_disc = sqrt(discriminant); + // The two roots of b*u^2 - a*u + c = 0. Pick the smallest positive u (= largest df = 1/u). + RealType u = (a - sqrt_disc) / (2 * b); + if (u <= 0) + u = (a + sqrt_disc) / (2 * b); + if (u > 0) + { + result = 1 / u; + return true; + } + } + return false; +} + } // namespace detail template @@ -332,16 +437,61 @@ BOOST_MATH_GPU_ENABLED RealType students_t_distribution::find_ hint = 1; detail::sample_size_func f(alpha, beta, sd, difference_from_mean); - tools::eps_tolerance tol(policies::digits()); - boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); - boost::math::pair r = tools::bracket_and_solve_root(f, hint, RealType(2), false, tol, max_iter, Policy()); - RealType result = r.first + (r.second - r.first) / 2; - if(max_iter >= policies::get_max_root_iterations()) + return detail::solve_for_degrees_of_freedom(f, hint, false, function, Policy()); +} + +template +BOOST_MATH_GPU_ENABLED RealType students_t_distribution::find_degrees_of_freedom( + RealType t, + RealType p) +{ + BOOST_MATH_STD_USING // for ADL of std functions + constexpr auto function = "boost::math::students_t_distribution<%1%>::find_degrees_of_freedom"; + // + // Check for domain errors: + // + RealType error_result; + if (false == detail::check_probability(function, p, &error_result, Policy())) + return error_result; + + if (t == 0) { - return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time: either there is no answer to how many degrees of freedom are required" // LCOV_EXCL_LINE - " or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE + // CDF(0; df) = 0.5 for all df; only solvable when p == 0.5. + if (p == static_cast(0.5)) + return policies::raise_overflow_error(function, nullptr, Policy()); + return policies::raise_domain_error( + function, + "No degrees of freedom can satisfy CDF(0; df) == %1% (must be 0.5).", + p, Policy()); } - return result; + + // Analytical cases: df = infinity (normal), df = 1 (cauchy) + boost::math::normal_distribution norm(0, 1); + if (detail::analytical_df_if_cdf_matches(norm, t, p)) + return policies::raise_overflow_error(function, nullptr, Policy()); + boost::math::cauchy_distribution cauchy(0, 1); + if (detail::analytical_df_if_cdf_matches(cauchy, t, p)) + return RealType(1); + + // Edgeworth warm start: compute a df estimate; fall back to df_hint_fallback if it fails + // or is too inaccurate + RealType hint = static_cast(detail::df_hint_fallback); + if (detail::approximate_df_with_edgeworth_expansion(t, p, hint)) + { + // Check that approximation is at least somewhat close; + // for small degrees of freedom it does not fail but is very inaccurate. + RealType relative_error_threshold = static_cast(0.1); + students_t_distribution t_approx(hint); + RealType p_approx = cdf(t_approx, t); + if (relative_difference(p_approx, p) > relative_error_threshold) + hint = static_cast(detail::df_hint_fallback); + } + // Root-find on f(df) = CDF(t; df) - p. + // CDF(t; df) is strictly increasing in df for t > 0, decreasing for t < 0. + // Use the smaller of p and q=1-p to avoid cancellation near 1. + RealType q = 1 - p; + detail::cdf_to_df_func f(t, p < q ? p : q, p < q ? false : true); + return detail::solve_for_degrees_of_freedom(f, hint, t > 0, function, Policy()); } template diff --git a/test/test_students_t.cpp b/test/test_students_t.cpp index a08cf75a01..3c63d3d639 100644 --- a/test/test_students_t.cpp +++ b/test/test_students_t.cpp @@ -548,6 +548,150 @@ void test_spots(RealType) static_cast(1.0))), 9); + // Tests for find_degrees_of_freedom(t, p) overload. + // Each case is derived from the CDF spot tests above: the exact df is known, + // so we verify that inverting CDF(x; df) = p recovers df to tight tolerance. + // float has ~7 significant digits; large-df inversion is ill-conditioned at that precision, + // so use a looser tolerance for single precision. + RealType tol_inv_df = std::is_same::value + ? static_cast(0.1) // float: 0.1% + : static_cast(0.01); // double+: 0.01% + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(-6.96455673428326), + static_cast(0.01)), + static_cast(2), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(-3.36492999890721), + static_cast(0.01)), + static_cast(5), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(-0.559429644), + static_cast(0.3)), + static_cast(5), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(1.475884049), + static_cast(0.9)), + static_cast(5), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(-1.475884049), + static_cast(0.1)), + static_cast(5), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(-5.2410429995425), + static_cast(0.00001)), + static_cast(25), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(6.96455673428326), + static_cast(0.99)), + static_cast(2), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(2), + static_cast(0.610822886098362)), + static_cast(0.1), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(2), + static_cast(0.777242554908434)), + static_cast(0.5), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(2), + static_cast(0.822925875908677)), + static_cast(0.75), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(2), + static_cast(0.977114826753374)), + static_cast(1000), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(1), + static_cast(0.8413326478347855)), + static_cast(1e4), + tol_inv_df * static_cast(5)); // Looser tolerance for ill-conditioned problem with very large df + // Small df edge cases where approximation becomes problematic: would need different values for float + // For now just test double and long double to have test overage of these code paths + if (!std::is_same::value) + { + // Small df case 1: Edgeworth expansion breaks down for small degrees of freedom, use fallback + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(1e20), + static_cast(0.5244796002843015)), + static_cast(1e-3), + tol_inv_df); + // Small df case 2: Edgeworth expansion succeeds but is inaccurate, use fallback + BOOST_CHECK_CLOSE( + students_t_distribution::find_degrees_of_freedom( + static_cast(2.0), + static_cast(0.500000000644961)), + static_cast(1e-10), + tol_inv_df); + } + // Analytical test: df=1 (Cauchy) case + { + boost::math::cauchy_distribution cauchy(0, 1); + RealType x = static_cast(1.0); + RealType p = cdf(cauchy, x); + RealType df_result = students_t_distribution::find_degrees_of_freedom(x, p); + BOOST_CHECK_EQUAL(df_result, static_cast(1.0)); + } + + + // Domain error: p outside (0,1) +#ifndef BOOST_NO_EXCEPTIONS + BOOST_MATH_CHECK_THROW( + students_t_distribution::find_degrees_of_freedom( + static_cast(2), + static_cast(-0.1)), + std::domain_error); + BOOST_MATH_CHECK_THROW( + students_t_distribution::find_degrees_of_freedom( + static_cast(2), + static_cast(1.1)), + std::domain_error); + // x == 0, p != 0.5: no df can satisfy CDF(0; df) == p -> domain error + BOOST_MATH_CHECK_THROW( + students_t_distribution::find_degrees_of_freedom( + static_cast(0), + static_cast(0.3)), + std::domain_error); + // x == 0, p == 0.5: CDF(0; df) == 0.5 for all df -> overflow (infinite solutions) + BOOST_MATH_CHECK_THROW( + students_t_distribution::find_degrees_of_freedom( + static_cast(0), + static_cast(0.5)), + std::overflow_error); + { + // Analytical test: df=infinity (Normal) case returns overflow error + boost::math::normal_distribution norm(0, 1); + RealType x = static_cast(1.0); + RealType p = cdf(norm, x); + BOOST_MATH_CHECK_THROW( + students_t_distribution::find_degrees_of_freedom(x, p), + std::overflow_error); + } +#endif + // Test for large degrees of freedom when should be same as normal. RealType inf = std::numeric_limits::infinity(); RealType nan = std::numeric_limits::quiet_NaN(); From 31197cace95de7fc2fd02b8b9dd0ba3a690db1d5 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 22 May 2026 18:06:26 +0100 Subject: [PATCH 65/96] Fix edge case short circuit logic in binomial distro inverse. Fixes: https://github.com/boostorg/math/issues/1397 --- include/boost/math/distributions/binomial.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/math/distributions/binomial.hpp b/include/boost/math/distributions/binomial.hpp index b17893e422..93b148eac8 100644 --- a/include/boost/math/distributions/binomial.hpp +++ b/include/boost/math/distributions/binomial.hpp @@ -226,7 +226,7 @@ namespace boost // but zero is the best we can do: return 0; } - if(p == 1 || success_fraction == 1) + if(q == 0 || success_fraction == 1) { // Probability of n or fewer successes is always one, // so n is the most sensible answer here: return trials; From 04eace43b14ae206f96fa8e5d2e94c16f21b4f3a Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Sun, 24 May 2026 11:11:05 -0700 Subject: [PATCH 66/96] Added checks for multiple degrees of freedom --- .../math/distributions/non_central_f.hpp | 15 +++++++++ test/test_nc_f.cpp | 33 ++++++++++++++----- 2 files changed, 39 insertions(+), 9 deletions(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index 5831f6622c..ccfa18f65c 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -111,6 +111,7 @@ namespace boost inline RealType find_degrees_of_freedom_f( const RealType x, const RealType v, const RealType nc, const bool find_v1, const RealType p, const RealType q, const Policy& pol) { + BOOST_MATH_STD_USING using std::fabs; const char* function = "non_central_f<%1%>::find_degrees_of_freedom_f"; if((p == 0) || (q == 0)) @@ -127,6 +128,20 @@ namespace boost RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } f_degrees_of_freedom_finder f(x, v, nc, find_v1, p < q ? p : q, p < q ? false : true); + + // There are times when f has two roots - thus, two degrees of freedom can + // be found. We find this case by checking if the sign of f for large and + // small values of v have the same sign. If the sign is the same, then there + // are an even number of roots. If the signs differ, there is only one root + // and we can safely find the root. + RealType vLarge = sqrt(boost::math::tools::max_value()); + RealType vSmall = 1 / vLarge; + + if ((f(vLarge) < 0) == (f(vSmall) < 0)){ + return policies::raise_evaluation_error(function, "Can't find degrees of freedom because two degrees of freedom can be found using the given parameters", + RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE + } + tools::eps_tolerance tol(policies::digits()); std::uintmax_t max_iter = policies::get_max_root_iterations(); // diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index d96675320d..78d4eea701 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -145,14 +145,6 @@ void test_spot( dist.find_non_centrality(x, a, b, P), ncp, tol * 10); BOOST_CHECK_CLOSE( dist.find_non_centrality(boost::math::complement(x, a, b, Q)), ncp, tol * 10); - BOOST_CHECK_CLOSE( - dist.find_v1(x, b, ncp, P), a, tol * 10); - BOOST_CHECK_CLOSE( - dist.find_v1(boost::math::complement(x, b, ncp, Q)), a, tol * 10); - BOOST_CHECK_CLOSE( - dist.find_v2(x, a, ncp, P), b, tol * 10); - BOOST_CHECK_CLOSE( - dist.find_v2(boost::math::complement(x, a, ncp, Q)), b, tol * 10); } if(boost::math::tools::digits() > 50) { @@ -377,8 +369,22 @@ void test_spots(RealType, const char* name = nullptr) } } } + + // Quick spot check for finding degrees of freedom + RealType v1 = 5; + RealType v2 = 2; + nc = 1; + x = 6; + boost::math::non_central_f_distribution ref(v1, v2, nc); + RealType P = cdf(ref, x); + BOOST_CHECK_CLOSE(ref.find_v1(x, v2, nc, P), v1, tolerance); + + // Check case where two degrees of freedom solve the inversion problem + BOOST_MATH_CHECK_THROW(dist.find_v1(RealType(1.5), RealType(2.0), RealType(1.0), RealType(0.49845842011686358665786775091245664L)), boost::math::evaluation_error); + BOOST_MATH_CHECK_THROW(dist.find_v1(RealType(3.51), RealType(5), RealType(0), RealType(0.85802971653663762108266155337333L)), boost::math::evaluation_error); + // Check find_v1/v2 edge cases - // Case when P=1 or P=0 + // Case when P=1 or P=0 nc = 2; BOOST_MATH_CHECK_THROW(dist.find_v1(x, b, nc, 1), std::domain_error); BOOST_MATH_CHECK_THROW(dist.find_v1(x, b, nc, 0), std::domain_error); @@ -389,6 +395,15 @@ void test_spots(RealType, const char* name = nullptr) x = boost::math::tools::epsilon() / 10; BOOST_MATH_CHECK_THROW(dist.find_v1(boost::math::complement(x, b, nc, 0.5)), boost::math::evaluation_error); BOOST_MATH_CHECK_THROW(dist.find_v1(x, b, nc, 0.5), boost::math::evaluation_error); + + BOOST_MATH_CHECK_THROW(dist.find_v2(x, b, nc, 1), std::domain_error); + BOOST_MATH_CHECK_THROW(dist.find_v2(x, b, nc, 0), std::domain_error); + // Case when Q=1 or Q=0 + BOOST_MATH_CHECK_THROW(dist.find_v2(boost::math::complement(x, b, nc, 1)), std::domain_error); + BOOST_MATH_CHECK_THROW(dist.find_v2(boost::math::complement(x, b, nc, 0)), std::domain_error); + // Check very small values of x an evaluation error is thrown + BOOST_MATH_CHECK_THROW(dist.find_v2(boost::math::complement(x, b, nc, 0.5)), boost::math::evaluation_error); + BOOST_MATH_CHECK_THROW(dist.find_v2(x, b, nc, 0.5), boost::math::evaluation_error); } // template void test_spots(RealType) BOOST_AUTO_TEST_CASE( test_main ) From d43d08c92d8cce7916d69cb9a3e3b2d9065ec52d Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Mon, 25 May 2026 12:18:04 -0700 Subject: [PATCH 67/96] Improved test coverage --- test/test_nc_f.cpp | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index 78d4eea701..c836ddc5fb 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -370,14 +370,21 @@ void test_spots(RealType, const char* name = nullptr) } } - // Quick spot check for finding degrees of freedom - RealType v1 = 5; - RealType v2 = 2; - nc = 1; - x = 6; - boost::math::non_central_f_distribution ref(v1, v2, nc); - RealType P = cdf(ref, x); - BOOST_CHECK_CLOSE(ref.find_v1(x, v2, nc, P), v1, tolerance); + // Quick spot check for finding degrees of freedom. When checking for two degrees + // of freedom for real_concept types, the cdf at large/small v2 can be greater than 1 + // or less than 0. + if (!std::is_same::value){ + RealType v1 = 10; + RealType v2 = 5; + nc = 1; + x = 6; + boost::math::non_central_f_distribution ref(v1, v2, nc); + RealType P = cdf(ref, x); + BOOST_CHECK_CLOSE(ref.find_v2(x, v1, nc, P), v2, tolerance); + BOOST_CHECK_CLOSE(ref.find_v2(boost::math::complement(x, v1, nc, 1-P)), v2, tolerance); + BOOST_CHECK_CLOSE(ref.find_v1(x, v2, nc, P), v1, tolerance); + BOOST_CHECK_CLOSE(ref.find_v1(boost::math::complement(x, v2, nc, 1-P)), v1, tolerance); + } // Check case where two degrees of freedom solve the inversion problem BOOST_MATH_CHECK_THROW(dist.find_v1(RealType(1.5), RealType(2.0), RealType(1.0), RealType(0.49845842011686358665786775091245664L)), boost::math::evaluation_error); From fbffcf0700392f85c0fe7c2feb392f90a181e630 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Mon, 25 May 2026 16:24:42 -0700 Subject: [PATCH 68/96] Completed test coverage for arcsine distribution --- test/test_arcsine.cpp | 36 +++++++++++++++++++++++++++++++----- 1 file changed, 31 insertions(+), 5 deletions(-) diff --git a/test/test_arcsine.cpp b/test/test_arcsine.cpp index cb5646e627..dc2fba919e 100644 --- a/test/test_arcsine.cpp +++ b/test/test_arcsine.cpp @@ -105,6 +105,7 @@ void test_ignore_policy(RealType) //std::cout << "pdf(ignore_error_arcsine(-1, +1), std::numeric_limits::infinity()) = " << pdf(ignore_error_arcsine(-1, +1), std::numeric_limits::infinity()) << std::endl; // Outputs: pdf(ignore_error_arcsine(-1, +1), std::numeric_limits::infinity()) = 1.#QNAN } + // PDF BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_arcsine(0, 1), std::numeric_limits::infinity()))); // x == infinity BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_arcsine(-1, 1), std::numeric_limits::infinity()))); // x == infinity BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_arcsine(0, 1), static_cast (-2)))); // x < xmin @@ -120,9 +121,31 @@ void test_ignore_policy(RealType) BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(0, 1), static_cast (+2)))); // x > x_max BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_arcsine(-1, 1), static_cast (+2)))); // x > x_max + // CDF + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_arcsine(0, 1), std::numeric_limits::infinity()))); // x == infinity + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_arcsine(-1, 1), std::numeric_limits::infinity()))); // x == infinity + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_arcsine(0, 1), static_cast (-2)))); // x < xmin + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_arcsine(-1, 1), static_cast (-2)))); // x < xmin + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_arcsine(0, 1), static_cast (+2)))); // x > x_max + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_arcsine(-1, 1), static_cast (+2)))); // x > x_max + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_arcsine(0, 1), std::numeric_limits::infinity())))); // x == infinity + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_arcsine(0, 1), static_cast (-2))))); // x < xmin + + // Quantile + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_arcsine(0, 1), std::numeric_limits::infinity()))); // p == infinity + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_arcsine(0, 1), static_cast(-1)))); // p < 0 + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_arcsine(0, 1), static_cast(2)))); // p > 1 + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_arcsine(0, 1), std::numeric_limits::infinity())))); // q == infinity + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_arcsine(0, 1), static_cast(-1))))); // q < 0 + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_arcsine(0, 1), static_cast(2))))); // q > 1 + // Mean BOOST_CHECK((boost::math::isnan)(mean(ignore_error_arcsine(-nan, 0)))); BOOST_CHECK((boost::math::isnan)(mean(ignore_error_arcsine(+nan, 0)))); + + // Median + BOOST_CHECK((boost::math::isnan)(median(ignore_error_arcsine(-nan, 0)))); + BOOST_CHECK((boost::math::isnan)(median(ignore_error_arcsine(+nan, 0)))); if (std::numeric_limits::has_infinity) { @@ -276,7 +299,8 @@ void test_spots(RealType) BOOST_CHECK_EQUAL(variance(arcsine_01), 0.125); // 1/8 = 0.125 BOOST_CHECK_CLOSE_FRACTION(standard_deviation(arcsine_01), one_div_root_two() / 2, tolerance); // 1/ sqrt(s) = 0.35355339059327379 BOOST_CHECK_EQUAL(skewness(arcsine_01), 0); // - BOOST_CHECK_EQUAL(kurtosis_excess(arcsine_01), -1.5); // 3/2 + BOOST_CHECK_EQUAL(kurtosis_excess(arcsine_01), -1.5); // -3/2 + BOOST_CHECK_EQUAL(kurtosis(arcsine_01), 1.5); // 3/2 BOOST_CHECK_EQUAL(support(arcsine_01).first, 0); // BOOST_CHECK_EQUAL(range(arcsine_01).first, 0); // BOOST_CHECK_THROW(mode(arcsine_01), std::domain_error); // Two modes at x_min and x_max, so throw instead. @@ -374,8 +398,8 @@ void test_spots(RealType) BOOST_CHECK_EQUAL(variance(as_m11), 0.5); // 1 - (-1) = 2 ^ 2 = 4 /8 = 0.5 BOOST_CHECK_EQUAL(skewness(as_m11), 0); // - BOOST_CHECK_EQUAL(kurtosis_excess(as_m11), -1.5); // 3/2 - + BOOST_CHECK_EQUAL(kurtosis_excess(as_m11), -1.5); // -3/2 + BOOST_CHECK_EQUAL(kurtosis(as_m11), 1.5); // 3/2 BOOST_CHECK_CLOSE_FRACTION(pdf(as_m11, 0.05), static_cast(0.31870852113797122803869876869296281629727218095644L), tolerance); BOOST_CHECK_CLOSE_FRACTION(pdf(as_m11, 0.5), static_cast(0.36755259694786136634088433220864629426492432024443L), tolerance); @@ -428,7 +452,8 @@ void test_spots(RealType) BOOST_CHECK_EQUAL(median(as_m2m1), -1.5); // 1 / (1 + 1) = 1/2 exactly. BOOST_CHECK_EQUAL(variance(as_m2m1), 0.125); BOOST_CHECK_EQUAL(skewness(as_m2m1), 0); // - BOOST_CHECK_EQUAL(kurtosis_excess(as_m2m1), -1.5); // 3/2 + BOOST_CHECK_EQUAL(kurtosis(as_m2m1), 1.5); // 3/2 + BOOST_CHECK_EQUAL(kurtosis_excess(as_m2m1), -1.5); // -3/2 BOOST_CHECK_CLOSE_FRACTION(pdf(as_m2m1, -1.95), static_cast(1.4605059227421865250256574657088244053723856445614L), 4 * tolerance); BOOST_CHECK_CLOSE_FRACTION(pdf(as_m2m1, -1.5), static_cast(0.63661977236758134307553505349005744813783858296183L), tolerance); @@ -615,7 +640,8 @@ void test_spots(RealType) BOOST_CHECK_EQUAL(variance(as), 0.125); //0.125 BOOST_CHECK_CLOSE_FRACTION(standard_deviation(as), one_div_root_two() / 2, std::numeric_limits::epsilon()); // 0.353553 BOOST_CHECK_EQUAL(skewness(as), 0); // - BOOST_CHECK_EQUAL(kurtosis_excess(as), -1.5); // 3/2 + BOOST_CHECK_EQUAL(kurtosis(as), 1.5); // 3/2 + BOOST_CHECK_EQUAL(kurtosis_excess(as), -1.5); // -3/2 BOOST_CHECK_EQUAL(support(as).first, 0); // BOOST_CHECK_EQUAL(range(as).first, 0); // BOOST_CHECK_THROW(mode(as), std::domain_error); // Two modes at x_min and x_max, so throw instead. From 0d7c2459b9724c6efa129692aef9a6951e0c8459 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Mon, 25 May 2026 17:53:49 -0700 Subject: [PATCH 69/96] Added bernoulli distribution tests --- test/test_bernoulli.cpp | 45 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/test/test_bernoulli.cpp b/test/test_bernoulli.cpp index 5f78b0afc5..256a897783 100644 --- a/test/test_bernoulli.cpp +++ b/test/test_bernoulli.cpp @@ -114,6 +114,11 @@ void test_spots(RealType) bernoulli_distribution(static_cast(0.4L)), static_cast(0)), static_cast(0.6L), tolerance); // Expect 1- p. + BOOST_CHECK_CLOSE_FRACTION( + pdf( // OK k (or n) + bernoulli_distribution(static_cast(0.4L)), static_cast(1)), + static_cast(0.4L), tolerance); // Expect p. + BOOST_CHECK_EQUAL( mean(bernoulli_distribution(static_cast(0.5L))), static_cast(0.5L)); @@ -280,6 +285,46 @@ void test_spots(RealType) BOOST_MATH_CHECK_THROW(quantile(complement(w, +inf)), std::domain_error); // p = + inf } // has_infinity #endif + + using boost::math::policies::policy; + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + typedef bernoulli_distribution ignore_error_bernoulli; + + if (std::numeric_limits::has_quiet_NaN) + { + // PDF + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_bernoulli(0.5), std::numeric_limits::infinity()))); // k == infinity + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_bernoulli(0.5), 2))); // k > 1 + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_bernoulli(0.5), static_cast(0.5)))); // k != 0, 1 + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_bernoulli(0.5), static_cast(-0.5)))); // k != 0, 1 + + // CDF + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_bernoulli(0.5), std::numeric_limits::infinity()))); // k == infinity + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_bernoulli(0.5), 2))); // k > 1 + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_bernoulli(0.5), static_cast(0.5)))); // k != 0, 1 + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_bernoulli(0.5), static_cast(-0.5)))); // k != 0, 1 + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_bernoulli(0.5), std::numeric_limits::infinity())))); // k == infinity + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_bernoulli(0.5), 2)))); // k > 1 + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_bernoulli(0.5), static_cast(0.5))))); // k != 0, 1 + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_bernoulli(0.5), static_cast(-0.5))))); // k != 0, 1 + + // // Quantile + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_bernoulli(0.5), std::numeric_limits::infinity()))); // p == infinity + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_bernoulli(0.5), static_cast(-1)))); // p < 0 + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_bernoulli(0.5), static_cast(2)))); // p > 1 + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_bernoulli(0.5), std::numeric_limits::infinity())))); // q == infinity + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_bernoulli(0.5), static_cast(-1))))); // q < 0 + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_bernoulli(0.5), static_cast(2))))); // q > 1 + + } // has_quiet_NaN } // template void test_spots(RealType) BOOST_AUTO_TEST_CASE( test_main ) From 264487e3543d86ee3d877d7691225cb8bcd7e98f Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Mon, 25 May 2026 20:11:23 -0700 Subject: [PATCH 70/96] Added beta distribution --- test/test_beta_dist.cpp | 68 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 66 insertions(+), 2 deletions(-) diff --git a/test/test_beta_dist.cpp b/test/test_beta_dist.cpp index d989197b6b..3983f526ff 100644 --- a/test/test_beta_dist.cpp +++ b/test/test_beta_dist.cpp @@ -65,6 +65,18 @@ using std::numeric_limits; # include #endif +template +void test_find_alpha_beta(RealType mean, RealType variance){ + BOOST_MATH_CHECK_THROW(beta_distribution::find_alpha(mean, variance), std::domain_error); + BOOST_MATH_CHECK_THROW(beta_distribution::find_beta(mean, variance), std::domain_error); +} + +template +void test_find_alpha_beta(RealType param, RealType x, RealType p){ + BOOST_MATH_CHECK_THROW(beta_distribution::find_alpha(param, x, p), std::domain_error); + BOOST_MATH_CHECK_THROW(beta_distribution::find_beta(param, x, p), std::domain_error); +} + template void test_spot( RealType a, // alpha a @@ -175,7 +187,9 @@ void test_spots(RealType) using ::boost::math::pdf; // Tests that should throw: - BOOST_MATH_CHECK_THROW(mode(beta_distribution(static_cast(1), static_cast(1))), std::domain_error); + BOOST_MATH_CHECK_THROW(mode(beta_distribution(static_cast(1), static_cast(1))), std::domain_error); // alpha = beta = 1 + BOOST_MATH_CHECK_THROW(mode(beta_distribution(static_cast(2), static_cast(1))), std::domain_error); // beta = 1 + BOOST_MATH_CHECK_THROW(mode(beta_distribution(static_cast(1), static_cast(2))), std::domain_error); // alpha = 1 // mode is undefined, and throws domain_error! // BOOST_MATH_CHECK_THROW(median(beta_distribution(static_cast(1), static_cast(1))), std::domain_error); @@ -497,7 +511,19 @@ void test_spots(RealType) BOOST_MATH_CHECK_THROW(quantile(dist, -1), std::domain_error); BOOST_MATH_CHECK_THROW(quantile(complement(dist, -1)), std::domain_error); - // No longer allow any parameter to be NaN or inf, so all these tests should throw. + BOOST_MATH_CHECK_THROW(pdf(beta_distribution(static_cast(0.5), static_cast(1.5)), static_cast(0)), std::overflow_error); // alpha < 1, x = 0 + BOOST_MATH_CHECK_THROW(pdf(beta_distribution(static_cast(1.5), static_cast(0.5)), static_cast(1)), std::overflow_error); // beta < 1, x = 1 + + // There aren't any upper bounds on the mean/variance. Should we add these? + test_find_alpha_beta(static_cast(-1), static_cast(1)); // mean < 0 + test_find_alpha_beta(static_cast(-1), static_cast(-1)); // var < 0 + + test_find_alpha_beta(static_cast(1), static_cast(0.5), static_cast(-0.5)); // p < 0 + test_find_alpha_beta(static_cast(1), static_cast(0.5), static_cast(1.5)); // p > 0 + test_find_alpha_beta(static_cast(1), static_cast(1.5), static_cast(0.5)); // x > 0 + test_find_alpha_beta(static_cast(1), static_cast(-0.5), static_cast(0.5)); // x < 0 + + // No longer allow any parameter to be NaN or inf, so all these tests should throw. if (std::numeric_limits::has_quiet_NaN) { // Attempt to construct from non-finite should throw. @@ -545,6 +571,37 @@ void test_spots(RealType) BOOST_MATH_CHECK_THROW(cdf(complement(w, +inf)), std::domain_error); // x = + inf BOOST_MATH_CHECK_THROW(quantile(w, +inf), std::domain_error); // p = + inf BOOST_MATH_CHECK_THROW(quantile(complement(w, +inf)), std::domain_error); // p = + inf + + test_find_alpha_beta(static_cast(1), static_cast(0.5), inf); // p = inf + test_find_alpha_beta(inf, static_cast(0.5)); // mean = inf + test_find_alpha_beta(static_cast(0.5), inf); // var = inf + test_find_alpha_beta(static_cast(1), inf, static_cast(0.5)); // x = inf + + // Check isnan policies + using boost::math::policies::policy; + + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + typedef beta_distribution ignore_error_beta; + + BOOST_CHECK((boost::math::isnan)(ignore_error_beta::find_alpha(inf, static_cast(0.5)))); // mean = inf + BOOST_CHECK((boost::math::isnan)(ignore_error_beta::find_beta(inf, static_cast(0.5)))); // mean = inf + BOOST_CHECK((boost::math::isnan)(ignore_error_beta::find_alpha(inf, static_cast(0.5), static_cast(0.5)))); // beta = inf + BOOST_CHECK((boost::math::isnan)(ignore_error_beta::find_beta(inf, static_cast(0.5), static_cast(0.5)))); // alpha = inf + BOOST_CHECK((boost::math::isnan)(mode(ignore_error_beta(static_cast(2), static_cast(1))))); // beta = 1 + BOOST_CHECK((boost::math::isnan)(mode(ignore_error_beta(static_cast(1), static_cast(2))))); // alpha = 1 + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_beta(static_cast(1), static_cast(1)), 2))); // x > 1 + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_beta(static_cast(1), static_cast(1)), 2))); // x > 1 + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_beta(static_cast(1), static_cast(1)), 2)))); // x > 1 + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_beta(static_cast(1), static_cast(1)), 2))); // p > 1 + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_beta(static_cast(1), static_cast(1)), 2)))); // p > 1 } // has_infinity // Error handling checks: @@ -600,6 +657,8 @@ BOOST_AUTO_TEST_CASE( test_main ) BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta11, 0.5), 0.5, 2 * std::numeric_limits::epsilon()); BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta11, 0.9), 0.9, 2 * std::numeric_limits::epsilon()); BOOST_CHECK_EQUAL(cdf(mybeta11, 1), 1.); // Exact unity expected. + BOOST_CHECK_EQUAL(cdf(complement(mybeta11, 1)), 0.0); // Exact 0 expected. + BOOST_CHECK_EQUAL(cdf(complement(mybeta11, 0)), 1.0); // Exact unity expected. double tol = std::numeric_limits::epsilon() * 10; BOOST_CHECK_EQUAL(pdf(mybeta22, 1), 0); // is dome shape. @@ -627,6 +686,11 @@ BOOST_AUTO_TEST_CASE( test_main ) // quantile. BOOST_CHECK_CLOSE_FRACTION(quantile(mybeta22, 0.028), 0.1, tol); BOOST_CHECK_CLOSE_FRACTION(quantile(complement(mybeta22, 1 - 0.028)), 0.1, tol); + BOOST_CHECK_EQUAL(quantile(mybeta22, 0), 0.); + BOOST_CHECK_EQUAL(quantile(mybeta22, 1), 1.); + BOOST_CHECK_EQUAL(quantile(complement(mybeta22, 0)), 1.); + BOOST_CHECK_EQUAL(quantile(complement(mybeta22, 1)), 0.); + BOOST_CHECK_EQUAL(kurtosis(mybeta11), 3+ kurtosis_excess(mybeta11)); // Check kurtosis_excess = kurtosis - 3; BOOST_CHECK_CLOSE_FRACTION(variance(mybeta22), 0.05, tol); BOOST_CHECK_CLOSE_FRACTION(mean(mybeta22), 0.5, tol); From 96f7df8e25d93b0e74dde886e3007427e1dd30da Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Tue, 26 May 2026 11:43:58 -0700 Subject: [PATCH 71/96] Added som basic binomial tests --- include/boost/math/distributions/binomial.hpp | 20 +++--- test/test_binomial.cpp | 67 +++++++++++++++++++ 2 files changed, 77 insertions(+), 10 deletions(-) diff --git a/include/boost/math/distributions/binomial.hpp b/include/boost/math/distributions/binomial.hpp index 93b148eac8..f2a800aeac 100644 --- a/include/boost/math/distributions/binomial.hpp +++ b/include/boost/math/distributions/binomial.hpp @@ -88,7 +88,7 @@ #include // error checks #include // isnan. #include // for root finding. - +#include #include namespace boost @@ -157,7 +157,7 @@ namespace boost template BOOST_MATH_CUDA_ENABLED inline bool check_dist_and_prob(const char* function, const RealType& N, RealType p, RealType prob, RealType* result, const Policy& pol) { - if((check_dist(function, N, p, result, pol) && detail::check_probability(function, prob, result, pol)) == false) + if(!(check_dist(function, N, p, result, pol) && detail::check_probability(function, prob, result, pol))) return false; return true; } @@ -321,11 +321,11 @@ namespace boost BOOST_MATH_STATIC const char* function = "boost::math::binomial_distribution<%1%>::find_lower_bound_on_p"; // Error checks: RealType result = 0; - if(false == binomial_detail::check_dist_and_k( + if(!(binomial_detail::check_dist_and_k( function, trials, RealType(0), successes, &result, Policy()) && binomial_detail::check_dist_and_prob( - function, trials, RealType(0), probability, &result, Policy())) + function, trials, RealType(0), probability, &result, Policy()))) { return result; } if(successes == 0) @@ -346,11 +346,11 @@ namespace boost BOOST_MATH_STATIC const char* function = "boost::math::binomial_distribution<%1%>::find_upper_bound_on_p"; // Error checks: RealType result = 0; - if(false == binomial_detail::check_dist_and_k( + if(!(binomial_detail::check_dist_and_k( function, trials, RealType(0), successes, &result, Policy()) && binomial_detail::check_dist_and_prob( - function, trials, RealType(0), probability, &result, Policy())) + function, trials, RealType(0), probability, &result, Policy()))) { return result; } if(trials == successes) @@ -373,11 +373,11 @@ namespace boost BOOST_MATH_STATIC const char* function = "boost::math::binomial_distribution<%1%>::find_minimum_number_of_trials"; // Error checks: RealType result = 0; - if(false == binomial_detail::check_dist_and_k( + if(!(binomial_detail::check_dist_and_k( function, k, p, k, &result, Policy()) && binomial_detail::check_dist_and_prob( - function, k, p, alpha, &result, Policy())) + function, k, p, alpha, &result, Policy()))) { return result; } result = ibetac_invb(k + 1, p, alpha, Policy()); // returns n - k @@ -392,11 +392,11 @@ namespace boost BOOST_MATH_STATIC const char* function = "boost::math::binomial_distribution<%1%>::find_maximum_number_of_trials"; // Error checks: RealType result = 0; - if(false == binomial_detail::check_dist_and_k( + if(!(binomial_detail::check_dist_and_k( function, k, p, k, &result, Policy()) && binomial_detail::check_dist_and_prob( - function, k, p, alpha, &result, Policy())) + function, k, p, alpha, &result, Policy()))) { return result; } result = ibeta_invb(k + 1, p, alpha, Policy()); // returns n - k diff --git a/test/test_binomial.cpp b/test/test_binomial.cpp index ef7f171723..45027aa823 100644 --- a/test/test_binomial.cpp +++ b/test/test_binomial.cpp @@ -599,7 +599,74 @@ void test_spots(RealType T) binomial_distribution(static_cast(8), static_cast(1)), static_cast(7)), static_cast(0) ); + BOOST_CHECK_EQUAL( + cdf( + complement(binomial_distribution(static_cast(8), static_cast(0)), + static_cast(7))), static_cast(0) + ); + BOOST_CHECK_EQUAL( + cdf( + complement(binomial_distribution(static_cast(8), static_cast(1)), + static_cast(7))), static_cast(1) + ); + // Check Error handling and edge cases + BOOST_CHECK_EQUAL(binomial_distribution::find_lower_bound_on_p(static_cast(5), + static_cast(0), + static_cast(0.5)), + static_cast(0)); // success = 0 + BOOST_CHECK_THROW(binomial_distribution::find_lower_bound_on_p(static_cast(5), + static_cast(2), + static_cast(-0.5)), + std::domain_error); // probability < 0 + BOOST_CHECK_THROW(binomial_distribution::find_lower_bound_on_p(static_cast(5), + static_cast(2), + static_cast(1.5)), + std::domain_error); // probability > 1 + BOOST_CHECK_EQUAL(binomial_distribution::find_upper_bound_on_p(static_cast(5), + static_cast(5), + static_cast(0.5)), + static_cast(1)); // trials = successes + BOOST_CHECK_THROW(binomial_distribution::find_upper_bound_on_p(static_cast(5), + static_cast(2), + static_cast(-0.5)), + std::domain_error); // probability < 0 + BOOST_CHECK_THROW(binomial_distribution::find_upper_bound_on_p(static_cast(5), + static_cast(2), + static_cast(1.5)), + std::domain_error); // probability > 1 + BOOST_CHECK_THROW(binomial_distribution::find_minimum_number_of_trials(static_cast(10), + static_cast(1.5), + static_cast(0.05)), + std::domain_error); // probability > 1 + BOOST_CHECK_THROW(binomial_distribution::find_minimum_number_of_trials(static_cast(10), + static_cast(-1), + static_cast(0.05)), + std::domain_error); // probability < 0 + BOOST_CHECK_THROW(binomial_distribution::find_minimum_number_of_trials(static_cast(10), + static_cast(0.5), + static_cast(-0.05)), + std::domain_error); // alpha < 0 + BOOST_CHECK_THROW(binomial_distribution::find_minimum_number_of_trials(static_cast(10), + static_cast(0.5), + static_cast(1.314)), + std::domain_error); // alpha > 1 + BOOST_CHECK_THROW(binomial_distribution::find_maximum_number_of_trials(static_cast(10), + static_cast(1.5), + static_cast(0.05)), + std::domain_error); // probability > 1 + BOOST_CHECK_THROW(binomial_distribution::find_maximum_number_of_trials(static_cast(10), + static_cast(-1), + static_cast(0.05)), + std::domain_error); // probability < 0 + BOOST_CHECK_THROW(binomial_distribution::find_maximum_number_of_trials(static_cast(10), + static_cast(0.5), + static_cast(-0.05)), + std::domain_error); // alpha < 0 + BOOST_CHECK_THROW(binomial_distribution::find_maximum_number_of_trials(static_cast(10), + static_cast(0.5), + static_cast(1.5)), + std::domain_error); // alpha > 1 #endif { From f2b43c29df778f10f9e3a9633b58df9d98f1e565 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Tue, 26 May 2026 12:48:04 -0700 Subject: [PATCH 72/96] Added nan handling for binomial distribution --- include/boost/math/distributions/binomial.hpp | 2 +- test/test_binomial.cpp | 37 +++++++++++++++++++ 2 files changed, 38 insertions(+), 1 deletion(-) diff --git a/include/boost/math/distributions/binomial.hpp b/include/boost/math/distributions/binomial.hpp index f2a800aeac..262823170f 100644 --- a/include/boost/math/distributions/binomial.hpp +++ b/include/boost/math/distributions/binomial.hpp @@ -88,7 +88,7 @@ #include // error checks #include // isnan. #include // for root finding. -#include + #include namespace boost diff --git a/test/test_binomial.cpp b/test/test_binomial.cpp index 45027aa823..b76c8e03e4 100644 --- a/test/test_binomial.cpp +++ b/test/test_binomial.cpp @@ -667,6 +667,43 @@ void test_spots(RealType T) static_cast(0.5), static_cast(1.5)), std::domain_error); // alpha > 1 + + if (std::numeric_limits::has_infinity) + { + // Check isnan policies + using boost::math::policies::policy; + + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + typedef binomial_distribution ignore_error_binomial; + + BOOST_CHECK((boost::math::isnan)(ignore_error_binomial::find_lower_bound_on_p(static_cast(5), + static_cast(2), + static_cast(-0.5)))); // probability < 0 + BOOST_CHECK((boost::math::isnan)(ignore_error_binomial::find_upper_bound_on_p(static_cast(5), + static_cast(2), + static_cast(-0.5)))); // probability < 0 + BOOST_CHECK((boost::math::isnan)(ignore_error_binomial::find_minimum_number_of_trials(static_cast(10), + static_cast(-1), + static_cast(0.05)))); // probability < 0 + BOOST_CHECK((boost::math::isnan)(ignore_error_binomial::find_maximum_number_of_trials(static_cast(10), + static_cast(-1), + static_cast(0.05)))); // probability < 0 + + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_binomial(static_cast(-8), static_cast(0.25)), static_cast(8)))); + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_binomial(static_cast(8), static_cast(-0.25)), static_cast(8)))); + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_binomial(static_cast(8), static_cast(0.25)), static_cast(9)))); + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_binomial(static_cast(8), static_cast(0.25)), static_cast(9)))); + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_binomial(static_cast(8), static_cast(0.25)), static_cast(9))))); + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_binomial(static_cast(8), static_cast(0.25)), static_cast(-1)))); + } #endif { From 2b07b7a3ecbb44fb78da985fd287e464b0b0b46e Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Tue, 26 May 2026 18:06:34 -0700 Subject: [PATCH 73/96] Added cauchy tests --- include/boost/math/distributions/cauchy.hpp | 18 +++---- test/test_cauchy.cpp | 56 ++++++++++++++++++++- 2 files changed, 64 insertions(+), 10 deletions(-) diff --git a/include/boost/math/distributions/cauchy.hpp b/include/boost/math/distributions/cauchy.hpp index 00a9ac9ac9..5a19f499a3 100644 --- a/include/boost/math/distributions/cauchy.hpp +++ b/include/boost/math/distributions/cauchy.hpp @@ -57,11 +57,11 @@ BOOST_MATH_GPU_ENABLED RealType cdf_imp(const cauchy_distribution((complement) ? 1 : 0); } #endif - if(false == detail::check_x(function, x, &result, Policy())) + if(!detail::check_x(function, x, &result, Policy())) { // Catches x == NaN return result; } @@ -111,15 +111,15 @@ BOOST_MATH_GPU_ENABLED RealType quantile_imp( RealType result = 0; RealType location = dist.location(); RealType scale = dist.scale(); - if(false == detail::check_location(function, location, &result, Policy())) + if(!detail::check_location(function, location, &result, Policy())) { return result; } - if(false == detail::check_scale(function, scale, &result, Policy())) + if(!detail::check_scale(function, scale, &result, Policy())) { return result; } - if(false == detail::check_probability(function, p, &result, Policy())) + if(!detail::check_probability(function, p, &result, Policy())) { return result; } @@ -224,11 +224,11 @@ BOOST_MATH_GPU_ENABLED inline RealType pdf(const cauchy_distribution >(0, 1); // (All) valid constructor parameter values. - + if (std::numeric_limits::has_infinity) + { + BOOST_CHECK_EQUAL( + boost::math::pdf( + cauchy_distribution(-2, 0.25), + std::numeric_limits::infinity()), // x + static_cast(0) // probability + ); + BOOST_CHECK_EQUAL( + boost::math::cdf( + cauchy_distribution(-2, 0.25), + std::numeric_limits::infinity()), // x + static_cast(1) // probability + ); + BOOST_CHECK_EQUAL( + boost::math::cdf( + cauchy_distribution(-2, 0.25), + -std::numeric_limits::infinity()), // x + static_cast(0) // probability + ); + BOOST_CHECK_EQUAL( + boost::math::quantile( + cauchy_distribution(-2, 0.25), + static_cast(0.5)), // p + static_cast(-2) // location + ); + + using boost::math::policies::policy; + + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + typedef boost::math::cauchy_distribution ignore_error_cauchy; + + // PDF + BOOST_CHECK((boost::math::isnan)(boost::math::pdf(ignore_error_cauchy(static_cast(0), static_cast(-1)), static_cast(0)))); + BOOST_CHECK((boost::math::isnan)(boost::math::pdf(ignore_error_cauchy(std::numeric_limits::infinity(), static_cast(1)), static_cast(0)))); + BOOST_CHECK((boost::math::isnan)(boost::math::pdf(ignore_error_cauchy(static_cast(0), static_cast(1)), std::numeric_limits::quiet_NaN()))); + + // CDF + BOOST_CHECK((boost::math::isnan)(boost::math::cdf(ignore_error_cauchy(static_cast(0), static_cast(-1)), static_cast(0)))); + BOOST_CHECK((boost::math::isnan)(boost::math::cdf(ignore_error_cauchy(std::numeric_limits::infinity(), static_cast(1)), static_cast(0)))); + BOOST_CHECK((boost::math::isnan)(boost::math::cdf(ignore_error_cauchy(static_cast(0), static_cast(1)), std::numeric_limits::quiet_NaN()))); + + // Quantile + BOOST_CHECK((boost::math::isnan)(boost::math::quantile(ignore_error_cauchy(static_cast(0), static_cast(-1)), static_cast(0.25)))); + BOOST_CHECK((boost::math::isnan)(boost::math::quantile(ignore_error_cauchy(std::numeric_limits::infinity(), static_cast(1)), static_cast(0.25)))); + BOOST_CHECK((boost::math::isnan)(boost::math::quantile(ignore_error_cauchy(static_cast(0), static_cast(1)), static_cast(-0.25)))); + } } // template void test_spots(RealType) From 5abd8a2a37885fb01400e101a3391e8b6da14007 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Tue, 26 May 2026 18:46:05 -0700 Subject: [PATCH 74/96] Added chi_squared distribution --- .../boost/math/distributions/chi_squared.hpp | 14 ++++----- test/test_chi_squared.cpp | 29 +++++++++++++++++++ 2 files changed, 35 insertions(+), 8 deletions(-) diff --git a/include/boost/math/distributions/chi_squared.hpp b/include/boost/math/distributions/chi_squared.hpp index 3944569e89..9823a2d555 100644 --- a/include/boost/math/distributions/chi_squared.hpp +++ b/include/boost/math/distributions/chi_squared.hpp @@ -103,7 +103,7 @@ BOOST_MATH_GPU_ENABLED RealType pdf(const chi_squared_distribution::find constexpr auto function = "boost::math::chi_squared_distribution<%1%>::find_degrees_of_freedom(%1%,%1%,%1%,%1%,%1%)"; // Check for domain errors: RealType error_result; - if(false == - detail::check_probability(function, alpha, &error_result, Policy()) - && detail::check_probability(function, beta, &error_result, Policy())) + if(!(detail::check_probability(function, alpha, &error_result, Policy()) + && detail::check_probability(function, beta, &error_result, Policy()))) { // Either probability is outside 0 to 1. return error_result; } diff --git a/test/test_chi_squared.cpp b/test/test_chi_squared.cpp index a64926ca58..360f6947de 100644 --- a/test/test_chi_squared.cpp +++ b/test/test_chi_squared.cpp @@ -546,6 +546,9 @@ void test_spots(RealType T) chi_squared_distribution(static_cast(8)), static_cast(1.1))), std::domain_error ); + BOOST_MATH_CHECK_THROW( + mode(chi_squared_distribution(static_cast(1))), std::domain_error + ); // This first test value is taken from an example here: // http://www.itl.nist.gov/div898/handbook/prc/section2/prc232.htm @@ -566,6 +569,32 @@ void test_spots(RealType T) check_out_of_range >(1); // (All) valid constructor parameter values. + // NaN handling + if (std::numeric_limits::has_infinity) + { + using boost::math::policies::policy; + + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + typedef boost::math::chi_squared_distribution ignore_error_chi_squared; + + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_chi_squared(static_cast(-1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_chi_squared(static_cast(-1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_chi_squared(static_cast(-1)), static_cast(1))))); + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_chi_squared(static_cast(-1)), static_cast(0.5)))); + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_chi_squared(static_cast(-1)), static_cast(0.5))))); + BOOST_CHECK((boost::math::isnan)(mode(ignore_error_chi_squared(static_cast(1))))); + BOOST_CHECK((boost::math::isnan)(ignore_error_chi_squared::find_degrees_of_freedom(10, -1, 0.01f, 100))); + BOOST_CHECK((boost::math::isnan)(ignore_error_chi_squared::find_degrees_of_freedom(10, 0.05f, -1, 100))); + + } } // template void test_spots(RealType) BOOST_AUTO_TEST_CASE( test_main ) From bfcc51ce5988def3a7c46dd187bd45f6e4513dcf Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Wed, 27 May 2026 14:23:31 -0700 Subject: [PATCH 75/96] Added exponential distribution tests --- test/test_exponential_dist.cpp | 42 ++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/test/test_exponential_dist.cpp b/test/test_exponential_dist.cpp index 124e2929ad..ecf5106d01 100644 --- a/test/test_exponential_dist.cpp +++ b/test/test_exponential_dist.cpp @@ -362,6 +362,48 @@ void test_spots(RealType T) BOOST_CHECK_EQUAL(pdf(exponential_distribution(2), inf), 0); BOOST_CHECK_EQUAL(cdf(exponential_distribution(2), inf), 1); BOOST_CHECK_EQUAL(cdf(complement(exponential_distribution(2), inf)), 0); + + // Test NaN policies + using boost::math::policies::policy; + + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + typedef boost::math::exponential_distribution ignore_error_exponential; + + // PDF + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_exponential(static_cast(-1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_exponential(static_cast(1)), static_cast(-1)))); + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_exponential(static_cast(-1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_exponential(static_cast(1)), static_cast(-1)))); + + // CDF + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_exponential(static_cast(-1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_exponential(static_cast(1)), static_cast(-1)))); + BOOST_CHECK((boost::math::isnan)(logcdf(ignore_error_exponential(static_cast(-1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(logcdf(ignore_error_exponential(static_cast(1)), static_cast(-1)))); + + // Complement CDF + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_exponential(static_cast(-1)), static_cast(1))))); + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_exponential(static_cast(1)), static_cast(-1))))); + BOOST_CHECK((boost::math::isnan)(logcdf(complement(ignore_error_exponential(static_cast(-1)), static_cast(1))))); + BOOST_CHECK((boost::math::isnan)(logcdf(complement(ignore_error_exponential(static_cast(1)), static_cast(-1))))); + + // Quantile + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_exponential(static_cast(-1)), static_cast(0.5)))); + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_exponential(static_cast(1)), static_cast(-0.5)))); + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_exponential(static_cast(-1)), static_cast(0.5))))); + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_exponential(static_cast(1)), static_cast(-0.5))))); + + // Mean and Standard Deviation + BOOST_CHECK((boost::math::isnan)(mean(ignore_error_exponential(static_cast(-1))))); + BOOST_CHECK((boost::math::isnan)(standard_deviation(ignore_error_exponential(static_cast(-1))))); } } // template void test_spots(RealType) From 3aa42a3f82f989ed746ed6050021f4d86465f810 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Wed, 27 May 2026 20:53:01 -0700 Subject: [PATCH 76/96] Added extreme value distribution --- .../math/distributions/extreme_value.hpp | 4 -- test/test_extreme_value.cpp | 59 +++++++++++++++++++ 2 files changed, 59 insertions(+), 4 deletions(-) diff --git a/include/boost/math/distributions/extreme_value.hpp b/include/boost/math/distributions/extreme_value.hpp index 73454d29d4..7049ac6a15 100644 --- a/include/boost/math/distributions/extreme_value.hpp +++ b/include/boost/math/distributions/extreme_value.hpp @@ -173,8 +173,6 @@ BOOST_MATH_GPU_ENABLED inline RealType cdf(const extreme_value_distribution&, %1%)", x, &result, Policy())) return result; @@ -199,8 +197,6 @@ BOOST_MATH_GPU_ENABLED inline RealType logcdf(const extreme_value_distribution&, %1%)", x, &result, Policy())) return result; diff --git a/test/test_extreme_value.cpp b/test/test_extreme_value.cpp index 255777914f..dd0c32ae47 100644 --- a/test/test_extreme_value.cpp +++ b/test/test_extreme_value.cpp @@ -253,6 +253,65 @@ void test_spots(RealType) BOOST_CHECK_EQUAL(cdf(complement(extreme_value_distribution(), inf)), 0); BOOST_CHECK_EQUAL(logcdf(extreme_value_distribution(), -inf), 0); BOOST_CHECK_EQUAL(logcdf(extreme_value_distribution(), inf), 1); + + // Test NaN policy + using boost::math::policies::policy; + + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + typedef boost::math::extreme_value_distribution ignore_error_extreme; + + // PDF + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_extreme(static_cast(0), static_cast(-1)), static_cast(0)))); + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_extreme(inf, static_cast(1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(pdf(ignore_error_extreme(static_cast(0), static_cast(1)), std::numeric_limits::quiet_NaN()))); + + // log(PDF) + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_extreme(static_cast(0), static_cast(-1)), static_cast(0)))); + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_extreme(inf, static_cast(1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(logpdf(ignore_error_extreme(static_cast(0), static_cast(1)), std::numeric_limits::quiet_NaN()))); + + // CDF + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_extreme(static_cast(0), static_cast(-1)), static_cast(0)))); + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_extreme(inf, static_cast(1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(cdf(ignore_error_extreme(static_cast(0), static_cast(1)), std::numeric_limits::quiet_NaN()))); + + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_extreme(static_cast(0), static_cast(-1)), static_cast(0))))); + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_extreme(inf, static_cast(1)), static_cast(1))))); + BOOST_CHECK((boost::math::isnan)(cdf(complement(ignore_error_extreme(static_cast(0), static_cast(1)), std::numeric_limits::quiet_NaN())))); + + // log(CDF) + BOOST_CHECK((boost::math::isnan)(logcdf(ignore_error_extreme(static_cast(0), static_cast(-1)), static_cast(0)))); + BOOST_CHECK((boost::math::isnan)(logcdf(ignore_error_extreme(inf, static_cast(1)), static_cast(1)))); + BOOST_CHECK((boost::math::isnan)(logcdf(ignore_error_extreme(static_cast(0), static_cast(1)), std::numeric_limits::quiet_NaN()))); + + BOOST_CHECK((boost::math::isnan)(logcdf(complement(ignore_error_extreme(static_cast(0), static_cast(-1)), static_cast(0))))); + BOOST_CHECK((boost::math::isnan)(logcdf(complement(ignore_error_extreme(inf, static_cast(1)), static_cast(1))))); + BOOST_CHECK((boost::math::isnan)(logcdf(complement(ignore_error_extreme(static_cast(0), static_cast(1)), std::numeric_limits::quiet_NaN())))); + + // Quantile + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_extreme(static_cast(0), static_cast(-1)), static_cast(0.5)))); + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_extreme(inf, static_cast(1)), static_cast(0.5)))); + BOOST_CHECK((boost::math::isnan)(quantile(ignore_error_extreme(static_cast(0), static_cast(1)), static_cast(-1)))); + + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_extreme(static_cast(0), static_cast(-1)), static_cast(0.5))))); + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_extreme(inf, static_cast(1)), static_cast(0.5))))); + BOOST_CHECK((boost::math::isnan)(quantile(complement(ignore_error_extreme(static_cast(0), static_cast(1)), static_cast(-1))))); + + // Mean + BOOST_CHECK((boost::math::isnan)(mean(ignore_error_extreme(static_cast(0), static_cast(-1))))); + BOOST_CHECK((boost::math::isnan)(mean(ignore_error_extreme(inf, static_cast(1))))); + + // Standard deviation + BOOST_CHECK((boost::math::isnan)(standard_deviation(ignore_error_extreme(static_cast(0), static_cast(-1))))); + BOOST_CHECK((boost::math::isnan)(standard_deviation(ignore_error_extreme(inf, static_cast(1))))); } // // Bug reports: From ba0f27463e0772231794bcbda634a54d89a71690 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Thu, 28 May 2026 09:53:28 -0700 Subject: [PATCH 77/96] Added find_location tests --- test/test_find_location.cpp | 36 ++++++++++++++++++++++++++++++++---- 1 file changed, 32 insertions(+), 4 deletions(-) diff --git a/test/test_find_location.cpp b/test/test_find_location.cpp index 7e0fb543a5..e16aaf5ae3 100644 --- a/test/test_find_location.cpp +++ b/test/test_find_location.cpp @@ -74,18 +74,18 @@ void test_spots(RealType) BOOST_MATH_CHECK_THROW(find_location(0., 2., 0.), std::domain_error); // p above 0 to 1. BOOST_MATH_CHECK_THROW(find_location(numeric_limits::infinity(), 0.5, 0.), std::domain_error); // z not finite. - BOOST_MATH_CHECK_THROW(find_location(numeric_limits::quiet_NaN(), -1., 0.), + BOOST_MATH_CHECK_THROW(find_location(numeric_limits::quiet_NaN(), 0.5, 0.), std::domain_error); // z not finite - BOOST_MATH_CHECK_THROW(find_location(0., -1., numeric_limits::quiet_NaN()), + BOOST_MATH_CHECK_THROW(find_location(0., 0.5, numeric_limits::quiet_NaN()), std::domain_error); // scale not finite BOOST_MATH_CHECK_THROW(find_location(complement(0., -1., 0.)), std::domain_error); // p below 0 to 1. BOOST_MATH_CHECK_THROW(find_location(complement(0., 2., 0.)), std::domain_error); // p above 0 to 1. BOOST_MATH_CHECK_THROW(find_location(complement(numeric_limits::infinity(), 0.5, 0.)), std::domain_error); // z not finite. - BOOST_MATH_CHECK_THROW(find_location(complement(numeric_limits::quiet_NaN(), -1., 0.)), + BOOST_MATH_CHECK_THROW(find_location(complement(numeric_limits::quiet_NaN(), 0.5, 0.)), std::domain_error); // z not finite - BOOST_MATH_CHECK_THROW(find_location(complement(0., -1., numeric_limits::quiet_NaN())), + BOOST_MATH_CHECK_THROW(find_location(complement(0., 0.5, numeric_limits::quiet_NaN())), std::domain_error); // scale not finite //// Check for ab-use with unsuitable distribution(s) when concept check implemented. @@ -131,12 +131,40 @@ void test_spots(RealType) // Check that can use the complement version. RealType q = 1 - p; // complement. // cout << "find_location >(complement(z, q, sd)) = " << endl; + l = find_location >(complement(z, q, sd, policy<>())); l = find_location >(complement(z, q, sd)); normal_distribution np95pc(l, sd); // Same standard_deviation (scale) but with mean(location) shifted // cout << "Normal distribution with mean = " << l << " has " << "fraction <= " << z << " = " << cdf(np95pc, z) << endl; BOOST_CHECK_CLOSE_FRACTION(q, cdf(np95pc, z), tolerance); + // Test NaN handling + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + if (std::numeric_limits::has_infinity) + { + RealType inf = std::numeric_limits::infinity(); + BOOST_CHECK((boost::math::isnan)(find_location >(static_cast(0), static_cast(-1), static_cast(1), ignore_all_policy()))); + BOOST_CHECK((boost::math::isnan)(find_location >(inf, static_cast(0.5), static_cast(1), ignore_all_policy()))); + BOOST_CHECK((boost::math::isnan)(find_location >(static_cast(0), static_cast(0.5), inf, ignore_all_policy()))); + // Negative scale does not throw. Seems weird... waiting to see if this should be fixed + // BOOST_CHECK((boost::math::isnan)(find_location >(static_cast(0), static_cast(0.5), static_cast(-1), ignore_all_policy()))); + + // Complement + BOOST_CHECK((boost::math::isnan)(find_location >(complement(static_cast(0), static_cast(-1), static_cast(1), ignore_all_policy())))); + BOOST_CHECK((boost::math::isnan)(find_location >(complement(inf, static_cast(0.5), static_cast(1), ignore_all_policy())))); + BOOST_CHECK((boost::math::isnan)(find_location >(complement(static_cast(0), static_cast(0.5), inf, ignore_all_policy())))); + } + + + } // template void test_spots(RealType) BOOST_AUTO_TEST_CASE( test_main ) From 8d57c0b2aee1985f6a7df9940a1e15ff91afc561 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Thu, 28 May 2026 10:12:50 -0700 Subject: [PATCH 78/96] Added find_scale --- test/test_find_scale.cpp | 32 ++++++++++++++++++++++++++++---- 1 file changed, 28 insertions(+), 4 deletions(-) diff --git a/test/test_find_scale.cpp b/test/test_find_scale.cpp index f31e47969d..a3c62f39e8 100644 --- a/test/test_find_scale.cpp +++ b/test/test_find_scale.cpp @@ -84,9 +84,9 @@ void test_spots(RealType) BOOST_MATH_CHECK_THROW(find_scale(complement(0., 2., 0.)), std::domain_error); // p above 0 to 1. BOOST_MATH_CHECK_THROW(find_scale(complement(numeric_limits::infinity(), 0.5, 0.)), std::domain_error); // z not finite. - BOOST_MATH_CHECK_THROW(find_scale(complement(numeric_limits::quiet_NaN(), -1., 0.)), + BOOST_MATH_CHECK_THROW(find_scale(complement(numeric_limits::quiet_NaN(), 0.5, 0.)), std::domain_error); // z not finite - BOOST_MATH_CHECK_THROW(find_scale(complement(0., -1., numeric_limits::quiet_NaN())), + BOOST_MATH_CHECK_THROW(find_scale(complement(0., 0.5, numeric_limits::quiet_NaN())), std::domain_error); // scale not finite BOOST_MATH_CHECK_THROW(find_scale(complement(0., -1., 0.)), std::domain_error); // p below 0 to 1. @@ -110,12 +110,36 @@ void test_spots(RealType) #ifndef BOOST_NO_EXCEPTIONS BOOST_CHECK_NO_THROW(find_scale(0, -1, 1, ignore_domain_policy())); // probability outside [0, 1] - BOOST_CHECK_NO_THROW(find_scale(numeric_limits::infinity(), -1, 1, + BOOST_CHECK_NO_THROW(find_scale(numeric_limits::infinity(), 0.5, 1, ignore_domain_policy())); // z not finite. BOOST_CHECK_NO_THROW(find_scale(complement(0, -1, 1, ignore_domain_policy()))); // probability outside [0, 1] - BOOST_CHECK_NO_THROW(find_scale(complement(numeric_limits::infinity(), -1, 1, + BOOST_CHECK_NO_THROW(find_scale(complement(numeric_limits::infinity(), 0.5, 1, ignore_domain_policy()))); // z not finite. #endif + + // Test NaN handling + typedef policy< + boost::math::policies::domain_error, + boost::math::policies::overflow_error, + boost::math::policies::underflow_error, + boost::math::policies::denorm_error, + boost::math::policies::pole_error, + boost::math::policies::evaluation_error + > ignore_all_policy; + + if (std::numeric_limits::has_infinity) + { + RealType inf = std::numeric_limits::infinity(); + BOOST_CHECK((boost::math::isnan)(find_scale >(static_cast(0), static_cast(-1), static_cast(1), ignore_all_policy()))); + BOOST_CHECK((boost::math::isnan)(find_scale >(inf, static_cast(0.5), static_cast(1), ignore_all_policy()))); + BOOST_CHECK((boost::math::isnan)(find_scale >(static_cast(0), static_cast(0.5), inf, ignore_all_policy()))); + + // Complement + BOOST_CHECK((boost::math::isnan)(find_scale >(complement(static_cast(0), static_cast(-1), static_cast(1), ignore_all_policy())))); + BOOST_CHECK((boost::math::isnan)(find_scale >(complement(inf, static_cast(0.5), static_cast(1), ignore_all_policy())))); + BOOST_CHECK((boost::math::isnan)(find_scale >(complement(static_cast(0), static_cast(0.5), inf, ignore_all_policy())))); + } + RealType l = 0.; // standard normal distribution. RealType sd = static_cast(1); // normal default standard deviation = 1. normal_distribution n01(l, sd); // mean(location) = 0, standard_deviation (scale) = 1. From 7a0975efc91e28c8c1db49c8f2c5e917eda59707 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Thu, 28 May 2026 13:53:18 -0700 Subject: [PATCH 79/96] Fixed find_location negative scale --- .../math/distributions/find_location.hpp | 19 ++++++++++--------- test/test_find_location.cpp | 2 +- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/include/boost/math/distributions/find_location.hpp b/include/boost/math/distributions/find_location.hpp index 33c2e64b62..c23b41eb10 100644 --- a/include/boost/math/distributions/find_location.hpp +++ b/include/boost/math/distributions/find_location.hpp @@ -10,6 +10,7 @@ #include // for all distribution signatures. #include +#include #include #include #include @@ -51,10 +52,10 @@ namespace boost return policies::raise_domain_error( function, "z parameter was %1%, but must be finite!", z, pol); } - if(!(boost::math::isfinite)(scale)) + typename Dist::value_type result; + if(!(boost::math::detail::check_scale)(function, scale, &result, pol)) { - return policies::raise_domain_error( - function, "scale parameter was %1%, but must be finite!", scale, pol); + return result; } //cout << "z " << z << ", p " << p << ", quantile(Dist(), p) " @@ -94,11 +95,11 @@ namespace boost return policies::raise_domain_error( function, "z parameter was %1%, but must be finite!", z, policies::policy<>()); } + typename Dist::value_type result; typename Dist::value_type scale = c.param2; - if(!(boost::math::isfinite)(scale)) + if(!(boost::math::detail::check_scale)(function, scale, &result, policies::policy<>())) { - return policies::raise_domain_error( - function, "scale parameter was %1%, but must be finite!", scale, policies::policy<>()); + return result; } // cout << "z " << c.dist << ", quantile (Dist(), " << c.param1 << ") * scale " << c.param2 << endl; return z - quantile(Dist(), p) * scale; @@ -123,11 +124,11 @@ namespace boost return policies::raise_domain_error( function, "z parameter was %1%, but must be finite!", z, c.param3); } + typename Dist::value_type result; typename Dist::value_type scale = c.param2; - if(!(boost::math::isfinite)(scale)) + if(!(boost::math::detail::check_scale)(function, scale, &result, c.param3)) { - return policies::raise_domain_error( - function, "scale parameter was %1%, but must be finite!", scale, c.param3); + return result; } // cout << "z " << c.dist << ", quantile (Dist(), " << c.param1 << ") * scale " << c.param2 << endl; return z - quantile(Dist(), p) * scale; diff --git a/test/test_find_location.cpp b/test/test_find_location.cpp index e16aaf5ae3..4d376dc254 100644 --- a/test/test_find_location.cpp +++ b/test/test_find_location.cpp @@ -155,7 +155,7 @@ void test_spots(RealType) BOOST_CHECK((boost::math::isnan)(find_location >(inf, static_cast(0.5), static_cast(1), ignore_all_policy()))); BOOST_CHECK((boost::math::isnan)(find_location >(static_cast(0), static_cast(0.5), inf, ignore_all_policy()))); // Negative scale does not throw. Seems weird... waiting to see if this should be fixed - // BOOST_CHECK((boost::math::isnan)(find_location >(static_cast(0), static_cast(0.5), static_cast(-1), ignore_all_policy()))); + BOOST_CHECK((boost::math::isnan)(find_location >(static_cast(0), static_cast(0.5), static_cast(-1), ignore_all_policy()))); // Complement BOOST_CHECK((boost::math::isnan)(find_location >(complement(static_cast(0), static_cast(-1), static_cast(1), ignore_all_policy())))); From bdce34ddf77b22dd4cefe632b7fc28e13f756757 Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Tue, 2 Jun 2026 10:05:01 -0700 Subject: [PATCH 80/96] Added chi squared approximation --- .../math/distributions/non_central_f.hpp | 23 ++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index ccfa18f65c..fcde87b2d2 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -107,6 +108,16 @@ namespace boost bool comp; }; + template + RealType large_v2_approximation(RealType x, RealType v1, RealType p, RealType q) + { // For v2 -> inf approximate f_degreese_of_freedome_finder with chi squared distribution + // with degrees of freedom v1 at the cdf at x * v1 + bool comp = p < q ? false : true; + RealType pval = p < q ? p : q; + chi_squared_distribution d(v1); + return comp ? pval - cdf(complement(d, x*v1)) : cdf(d, x*v1) - pval; + } + template inline RealType find_degrees_of_freedom_f( const RealType x, const RealType v, const RealType nc, const bool find_v1, const RealType p, const RealType q, const Policy& pol) @@ -137,7 +148,17 @@ namespace boost RealType vLarge = sqrt(boost::math::tools::max_value()); RealType vSmall = 1 / vLarge; - if ((f(vLarge) < 0) == (f(vSmall) < 0)){ + // As v2 -> infinity, noncentral f converges to chi-squared distribution. + // Rather than evaluating f(vLarge), we can use the more stable chi-squared approximation + RealType large_difference; + if (find_v1) + { + large_difference = large_v2_approximation(x, v, p, q); + } + else + large_difference = f(vLarge); + + if ((large_difference < 0) == (f(vSmall) < 0)){ return policies::raise_evaluation_error(function, "Can't find degrees of freedom because two degrees of freedom can be found using the given parameters", RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } From a0bd79c0e2c9c29ce3ec29434c06bc36b4bc6edd Mon Sep 17 00:00:00 2001 From: Jacob Hass Date: Tue, 2 Jun 2026 13:18:04 -0700 Subject: [PATCH 81/96] Changed to non-central chi-squared dist --- include/boost/math/distributions/non_central_f.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/boost/math/distributions/non_central_f.hpp b/include/boost/math/distributions/non_central_f.hpp index fcde87b2d2..9bae5c41ff 100644 --- a/include/boost/math/distributions/non_central_f.hpp +++ b/include/boost/math/distributions/non_central_f.hpp @@ -14,7 +14,7 @@ #include #include #include -#include +#include #include #include #include @@ -109,12 +109,12 @@ namespace boost }; template - RealType large_v2_approximation(RealType x, RealType v1, RealType p, RealType q) + RealType large_v2_approximation(RealType x, RealType v1, RealType p, RealType q, RealType nc) { // For v2 -> inf approximate f_degreese_of_freedome_finder with chi squared distribution // with degrees of freedom v1 at the cdf at x * v1 bool comp = p < q ? false : true; RealType pval = p < q ? p : q; - chi_squared_distribution d(v1); + non_central_chi_squared_distribution d(v1, nc); return comp ? pval - cdf(complement(d, x*v1)) : cdf(d, x*v1) - pval; } @@ -153,7 +153,7 @@ namespace boost RealType large_difference; if (find_v1) { - large_difference = large_v2_approximation(x, v, p, q); + large_difference = large_v2_approximation(x, v, p, q, nc); } else large_difference = f(vLarge); From 38ebc9a6cfa84e86fcaec62d0e261cec70f2041e Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Mon, 8 Jun 2026 21:16:39 +0200 Subject: [PATCH 82/96] ENH: add warm start for root finder --- .../boost/math/distributions/students_t.hpp | 29 +++++++++++++++++-- 1 file changed, 27 insertions(+), 2 deletions(-) diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index bc1c7139eb..50f591b690 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -25,6 +25,7 @@ #include #include #include +#include #ifdef _MSC_VER # pragma warning(push) @@ -58,7 +59,7 @@ class students_t_distribution RealType alpha, RealType beta, RealType sd, - RealType hint = 100); + RealType hint = 0); BOOST_MATH_GPU_ENABLED static RealType find_degrees_of_freedom( RealType t, @@ -413,6 +414,30 @@ BOOST_MATH_GPU_ENABLED bool approximate_df_with_edgeworth_expansion(RealType x, return false; } +template +BOOST_MATH_GPU_ENABLED RealType calculate_hill_df_guess( + RealType difference_from_mean, + RealType alpha, + RealType beta, + RealType sd) + + // Hill-corrected normal approximation used as an initial guess + // for find_degrees_of_freedom(difference_from_mean, alpha, beta, sd). + // t_p,nu approx z_p * (1 + (z_p^2 + 1) / (4*nu)) + // nu + 1 = (sd/diff)^2 * (t_alpha + t_beta)^2 + +{ + BOOST_MATH_STD_USING + normal_distribution n(0, 1); + RealType za = quantile(complement(n, alpha)); + RealType zb = quantile(complement(n, beta)); + + RealType v_norm = pow((za + zb) * sd / difference_from_mean, 2) - 1; + RealType hint = v_norm + (za * za - za * zb + zb * zb + 1) / 2; + + return (hint <= 0) ? RealType(1) : hint; +} + } // namespace detail template @@ -434,7 +459,7 @@ BOOST_MATH_GPU_ENABLED RealType students_t_distribution::find_ return error_result; if(hint <= 0) - hint = 1; + hint = detail::calculate_hill_df_guess(difference_from_mean, alpha, beta, sd); detail::sample_size_func f(alpha, beta, sd, difference_from_mean); return detail::solve_for_degrees_of_freedom(f, hint, false, function, Policy()); From ce299579d420ce37f563d561bdb8af7715deb1a7 Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Mon, 8 Jun 2026 21:24:00 +0200 Subject: [PATCH 83/96] ENH: add initial guess for power equation inversion --- include/boost/math/distributions/students_t.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index 50f591b690..94fe686c2d 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -25,7 +25,6 @@ #include #include #include -#include #ifdef _MSC_VER # pragma warning(push) From 4426160fb2afab7a08aae6aca9d62de2d9c3c927 Mon Sep 17 00:00:00 2001 From: Jacob Hass <72838561+JacobHass8@users.noreply.github.com> Date: Mon, 8 Jun 2026 17:21:48 -0700 Subject: [PATCH 84/96] Fix type for testing small x --- test/test_nc_f.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_nc_f.cpp b/test/test_nc_f.cpp index c836ddc5fb..7ad45c1ba3 100644 --- a/test/test_nc_f.cpp +++ b/test/test_nc_f.cpp @@ -399,7 +399,7 @@ void test_spots(RealType, const char* name = nullptr) BOOST_MATH_CHECK_THROW(dist.find_v1(boost::math::complement(x, b, nc, 1)), std::domain_error); BOOST_MATH_CHECK_THROW(dist.find_v1(boost::math::complement(x, b, nc, 0)), std::domain_error); // Check very small values of x an evaluation error is thrown - x = boost::math::tools::epsilon() / 10; + x = boost::math::tools::epsilon() / 10; BOOST_MATH_CHECK_THROW(dist.find_v1(boost::math::complement(x, b, nc, 0.5)), boost::math::evaluation_error); BOOST_MATH_CHECK_THROW(dist.find_v1(x, b, nc, 0.5), boost::math::evaluation_error); From 22da3792e592c3a67c5cba998e88df89807e6fcb Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Tue, 9 Jun 2026 17:51:20 +0100 Subject: [PATCH 85/96] students_t: Add exception handling guards. The degree of freedom finders can call the quantile with arguments which result in an infinite result, as far as our root finders are concerned, any large value as a result will do, but we need to not allow an exception to escape or the root finder terminates. --- .../boost/math/distributions/students_t.hpp | 27 +++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index 94fe686c2d..559b155984 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -315,8 +315,31 @@ struct sample_size_func return 1; } students_t_distribution t(df); - RealType qa = quantile(complement(t, alpha)); - RealType qb = quantile(complement(t, beta)); + RealType qa, qb; +#ifndef BOOST_MATH_NO_EXCEPTIONS + try { +#endif + qa = quantile(complement(t, alpha)); +#ifndef BOOST_MATH_NO_EXCEPTIONS + } + catch (const std::overflow_error&) + { + // Any arbitrary large value will do, must be positive since we calculate qa * qa below: + return tools::max_value(); + } +#endif +#ifndef BOOST_MATH_NO_EXCEPTIONS + try { +#endif + qb = quantile(complement(t, beta)); +#ifndef BOOST_MATH_NO_EXCEPTIONS + } + catch (const std::overflow_error&) + { + // Any arbitrary large value will do, will be negative if beta < 0: + return beta < 0 ? -tools::max_value() : tools::max_value(); + } +#endif qa += qb; qa *= qa; qa *= ratio; From fecc83bd81b71850309c19a5636250f8139c5e99 Mon Sep 17 00:00:00 2001 From: a-leontyev Date: Tue, 9 Jun 2026 21:17:58 +0300 Subject: [PATCH 86/96] Guard against intermediate x*x overflow in Student's t pdf/cdf (#1402) * Guard against intermediate x*x overflow in Student's t pdf/cdf * Return exact tail values instead of raising overflow_error (review feedback) --------- Co-authored-by: LeantionX Co-authored-by: Matt Borland --- .../boost/math/distributions/students_t.hpp | 12 ++ test/Jamfile.v2 | 1 + test/test_students_t_overflow.cpp | 133 ++++++++++++++++++ 3 files changed, 146 insertions(+) create mode 100644 test/test_students_t_overflow.cpp diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index bc1c7139eb..69e5dde43c 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -129,6 +129,11 @@ BOOST_MATH_GPU_ENABLED inline RealType pdf(const students_t_distribution sqrt(tools::max_value())) + { + // Exact limit: the pdf underflows to zero long before x * x can overflow. + return 0; + } RealType basem1 = x * x / df; if(basem1 < 0.125) { @@ -146,6 +151,8 @@ BOOST_MATH_GPU_ENABLED inline RealType pdf(const students_t_distribution BOOST_MATH_GPU_ENABLED inline RealType cdf(const students_t_distribution& dist, const RealType& x) { + BOOST_MATH_STD_USING // for ADL of std functions + RealType error_result; // degrees_of_freedom > 0 or infinity check: RealType df = dist.degrees_of_freedom(); @@ -201,6 +208,11 @@ BOOST_MATH_GPU_ENABLED inline RealType cdf(const students_t_distribution sqrt(tools::max_value())) + { + // Exact tail values: avoids a spurious intermediate overflow in x * x below. + return (x < 0) ? static_cast(0) : static_cast(1); + } RealType x2 = x * x; RealType probability; if(df > 2 * x2) diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index c5bde01635..5b67e474df 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -895,6 +895,7 @@ test-suite distribution_tests : : test_poisson_real_concept ] [ run test_rayleigh.cpp /boost/test//boost_unit_test_framework ] [ run test_students_t.cpp /boost/test//boost_unit_test_framework ] + [ run test_students_t_overflow.cpp /boost/test//boost_unit_test_framework ] [ run test_skew_normal.cpp /boost/test//boost_unit_test_framework ] [ run test_triangular.cpp pch /boost/test//boost_unit_test_framework ] [ run test_uniform.cpp pch /boost/test//boost_unit_test_framework ] diff --git a/test/test_students_t_overflow.cpp b/test/test_students_t_overflow.cpp new file mode 100644 index 0000000000..d4020d3bf4 --- /dev/null +++ b/test/test_students_t_overflow.cpp @@ -0,0 +1,133 @@ +// Copyright Anton Leontev 2026. +// 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) +// +// Regression test for spurious intermediate overflow of x*x in the +// Student's t pdf() and cdf(): for |x| > sqrt(max_value) the functions +// must return the exact tail values (0 for the pdf, 0/1 for the cdf and +// its complement) without raising FE_OVERFLOW and without invoking any +// error handler. + +#define BOOST_TEST_MAIN +#include +#include +#include + +#include +#include +#include + +#ifndef BOOST_MATH_STANDALONE +#include // for real_concept +#include // for cpp_bin_float_50 +#endif + +#include +#include +#include + +#pragma STDC FENV_ACCESS ON + +template +RealType overflowing_deviate() +{ + BOOST_MATH_STD_USING // ADL sqrt for built-in and UDT types + using boost::math::tools::max_value; + return 2 * sqrt(max_value()); +} + +// For |x| > sqrt(max_value) the mathematically exact results have already +// saturated: pdf == 0, cdf == 0 or 1. These must be returned directly. +// A throwing policy is used so that any error-handler invocation would +// fail the test with an unexpected exception. +template +void test_tail_values(RealType) +{ + using namespace boost::math; + using namespace boost::math::policies; + + typedef policy > strict_policy; + typedef students_t_distribution dist_t; + + dist_t dist(static_cast(5)); + const RealType big = overflowing_deviate(); + + BOOST_CHECK((boost::math::isfinite)(big)); + + BOOST_CHECK_EQUAL(pdf(dist, big), static_cast(0)); + BOOST_CHECK_EQUAL(pdf(dist, -big), static_cast(0)); + + BOOST_CHECK_EQUAL(cdf(dist, big), static_cast(1)); + BOOST_CHECK_EQUAL(cdf(dist, -big), static_cast(0)); + + BOOST_CHECK_EQUAL(cdf(complement(dist, big)), static_cast(0)); + BOOST_CHECK_EQUAL(cdf(complement(dist, -big)), static_cast(1)); +} + +// Built-in types: the calls must not leave FE_OVERFLOW set -- avoiding the +// spurious intermediate overflow is the whole point of the guard. +template +void test_no_fp_exceptions(RealType) +{ + using namespace boost::math; + + students_t_distribution dist(static_cast(5)); + const RealType big = overflowing_deviate(); + + std::feclearexcept(FE_ALL_EXCEPT); + RealType p = pdf(dist, big); + RealType c1 = cdf(dist, big); + RealType c2 = cdf(dist, -big); + BOOST_CHECK_EQUAL(std::fetestexcept(FE_OVERFLOW), 0); + BOOST_CHECK_EQUAL(p, static_cast(0)); + BOOST_CHECK_EQUAL(c1, static_cast(1)); + BOOST_CHECK_EQUAL(c2, static_cast(0)); +} + +// Ordinary arguments must be unaffected by the guard. +template +void test_no_regression(RealType) +{ + using namespace boost::math; + + students_t_distribution dist(static_cast(10)); + const RealType tol = boost::math::tools::epsilon() * 100; + + BOOST_CHECK_CLOSE_FRACTION(cdf(dist, static_cast(0)), + static_cast(0.5), tol); + + RealType p = pdf(dist, static_cast(1.5)); + BOOST_CHECK((boost::math::isfinite)(p)); + BOOST_CHECK(p > 0); + + RealType big_but_safe = static_cast(1e6); + BOOST_CHECK_CLOSE_FRACTION(cdf(dist, big_but_safe), + static_cast(1), tol); +} + +BOOST_AUTO_TEST_CASE(students_t_overflow_test) +{ + test_tail_values(0.0F); + test_tail_values(0.0); + test_tail_values(0.0L); + + test_no_fp_exceptions(0.0F); + test_no_fp_exceptions(0.0); + test_no_fp_exceptions(0.0L); + + test_no_regression(0.0F); + test_no_regression(0.0); + test_no_regression(0.0L); + + // real_concept and multiprecision compute distribution constants via + // lexical_cast, which standalone mode disables, so both are skipped there. +#ifndef BOOST_MATH_STANDALONE + test_tail_values(boost::math::concepts::real_concept(0)); + test_no_regression(boost::math::concepts::real_concept(0)); + + test_tail_values(boost::multiprecision::cpp_bin_float_50(0)); + test_no_regression(boost::multiprecision::cpp_bin_float_50(0)); +#endif +} From b1f70c8f988be14a07423947cb3e491bdaa7cae1 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Thu, 11 Jun 2026 17:53:02 +0100 Subject: [PATCH 87/96] pFq: prevent spurious overflow. Make sure we throw the correct exceptions when the result is infinite. Fixes #1404. --- .../hypergeometric_pFq_checked_series.hpp | 33 ++++++++++++------- .../special_functions/hypergeometric_pFq.hpp | 16 +++++++-- test/test_pFq.hpp | 22 +++++++++++++ 3 files changed, 56 insertions(+), 15 deletions(-) diff --git a/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp b/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp index 3adc86f74f..58c0d6b1ea 100644 --- a/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp +++ b/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp @@ -23,6 +23,7 @@ unsigned set_crossover_locations(const Seq& aj, const Seq& bj, const Real& z, unsigned int* crossover_locations) { BOOST_MATH_STD_USING + using nothrow_policy = typename boost::math::policies::normalise, boost::math::policies::rounding_error>::type; unsigned N_terms = 0; if(aj.size() == 1 && bj.size() == 1) @@ -55,13 +56,13 @@ Real t = (-sqrt(sq) - b + z) / 2; if (t >= 0) { - crossover_locations[N_terms] = itrunc(t); + crossover_locations[N_terms] = itrunc(t, nothrow_policy()); ++N_terms; } t = (sqrt(sq) - b + z) / 2; if (t >= 0) { - crossover_locations[N_terms] = itrunc(t); + crossover_locations[N_terms] = itrunc(t, nothrow_policy()); ++N_terms; } } @@ -71,13 +72,13 @@ Real t = (-sqrt(sq) - b - z) / 2; if (t >= 0) { - crossover_locations[N_terms] = itrunc(t); + crossover_locations[N_terms] = itrunc(t, nothrow_policy()); ++N_terms; } t = (sqrt(sq) - b - z) / 2; if (t >= 0) { - crossover_locations[N_terms] = itrunc(t); + crossover_locations[N_terms] = itrunc(t, nothrow_policy()); ++N_terms; } } @@ -110,7 +111,7 @@ unsigned n = 0; for (auto bi = bj.begin(); bi != bj.end(); ++bi, ++n) { - crossover_locations[n] = *bi >= 0 ? 0 : itrunc(-*bi) + 1; + crossover_locations[n] = *bi >= 0 ? 0 : itrunc(-*bi, nothrow_policy()) + 1; } std::sort(crossover_locations, crossover_locations + bj.size(), std::less()); N_terms = (unsigned)bj.size(); @@ -122,6 +123,7 @@ std::pair hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, const Terminal& termination, long long& log_scale) { BOOST_MATH_STD_USING + using nothrow_policy = typename boost::math::policies::normalise>::type; Real result = 1; Real abs_result = 1; Real term = 1; @@ -129,8 +131,15 @@ Real tol = boost::math::policies::get_epsilon(); std::uintmax_t k = 0; Real upper_limit(sqrt(boost::math::tools::max_value())), diff; + if ((tools::max_value() / fabs(z) < upper_limit)) + upper_limit = tools::max_value() / fabs(z); + for (auto pa = aj.begin(); pa != aj.end(); ++pa) + { + if (tools::max_value() / fabs(*pa) < upper_limit) + upper_limit = tools::max_value() / fabs(*pa); + } Real lower_limit(1 / upper_limit); - long long log_scaling_factor = lltrunc(boost::math::tools::log_max_value()) - 2; + long long log_scaling_factor = lltrunc(boost::math::tools::log_max_value(), nothrow_policy()) - 2; Real scaling_factor = exp(Real(log_scaling_factor)); Real term_m1; long long local_scaling = 0; @@ -234,7 +243,7 @@ log_scale += log_scaling_factor; local_scaling += log_scaling_factor; } - if (fabs(abs_result) < lower_limit) + else if ((fabs(abs_result) < lower_limit) && (fabs(abs_result) * scaling_factor < upper_limit)) { abs_result *= scaling_factor; result *= scaling_factor; @@ -321,7 +330,7 @@ else { int ls = 1; - Real p = log_pochhammer(*ai, s, pol, &ls); + Real p = log_pochhammer(static_cast(*ai), s, pol, &ls); s1 *= ls; term += p; loop_error_scale = (std::max)(p, loop_error_scale); @@ -334,7 +343,7 @@ for (auto bi = bj.begin(); bi != bj.end(); ++bi) { int ls = 1; - Real p = log_pochhammer(*bi, s, pol, &ls); + Real p = log_pochhammer(static_cast(*bi), s, pol, &ls); s2 *= ls; term -= p; loop_error_scale = (std::max)(p, loop_error_scale); @@ -371,13 +380,13 @@ if (term <= tools::log_min_value()) { // rescale if we can: - long long scale = lltrunc(floor(term - tools::log_min_value()) - 2); + long long scale = lltrunc(floor(term - tools::log_min_value()) - 2, nothrow_policy()); term -= scale; loop_scale += scale; } if (term > 10) { - int scale = itrunc(floor(term)); + long long scale = lltrunc(floor(term), nothrow_policy()); term -= scale; loop_scale += scale; } @@ -402,7 +411,7 @@ term /= scaling_factor; loop_scale += log_scaling_factor; } - if (fabs(loop_result) < lower_limit) + else if ((fabs(loop_result) < lower_limit) && (fabs(loop_result) * scaling_factor < upper_limit)) { loop_result *= scaling_factor; loop_abs_result *= scaling_factor; diff --git a/include/boost/math/special_functions/hypergeometric_pFq.hpp b/include/boost/math/special_functions/hypergeometric_pFq.hpp index 070b852dd2..96f68e510e 100644 --- a/include/boost/math/special_functions/hypergeometric_pFq.hpp +++ b/include/boost/math/special_functions/hypergeometric_pFq.hpp @@ -58,12 +58,22 @@ namespace boost { BOOST_MATH_STD_USING long long scale = 0; + static const char* function = "boost::math::hypergeometric_pFq<%1%>(%1%,%1%,%1%)"; std::pair r = boost::math::detail::hypergeometric_pFq_checked_series_impl(aj, bj, value_type(z), pol, boost::math::detail::iteration_terminator(boost::math::policies::get_max_series_iterations()), scale); - r.first *= exp(Real(scale)); - r.second *= exp(Real(scale)); + // + // Overflow check: + // + if (static_cast(scale) > tools::log_max_value()) + return (r.first < 0 ? -1 : 1) * policies::raise_overflow_error(function, nullptr, pol); + Real mul = exp(Real(scale)); + if(fabs(r.first) > 1) + if(tools::max_value() / fabs(r.first) < mul) + return (r.first < 0 ? -1 : 1) * policies::raise_overflow_error(function, nullptr, pol); + r.first *= mul; + r.second *= mul; if (p_abs_error) *p_abs_error = static_cast(r.second) * boost::math::tools::epsilon(); - return policies::checked_narrowing_cast(r.first, "boost::math::hypergeometric_pFq<%1%>(%1%,%1%,%1%)"); + return policies::checked_narrowing_cast(r.first, function); } template diff --git a/test/test_pFq.hpp b/test/test_pFq.hpp index 0e7f3a0c54..e9e5bac67e 100644 --- a/test/test_pFq.hpp +++ b/test/test_pFq.hpp @@ -298,6 +298,27 @@ void test_spots_2F2(T, const char*) } } +template +void test_special_cases(T, const char*) +{ + typedef boost::math::policies::policy> nothrow_policy; + typedef boost::math::policies::policy> throw_policy; + T error_rate; + BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ 2 }, { 3 }, static_cast(1e10), &error_rate, throw_policy()), std::overflow_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ static_cast(2) }, { static_cast(3) }, static_cast(1e10), &error_rate, throw_policy()), std::overflow_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ 2 }, { 3 }, boost::math::tools::max_value(), &error_rate, throw_policy()), std::overflow_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ static_cast(2) }, { static_cast(3) }, boost::math::tools::max_value(), &error_rate, throw_policy()), std::overflow_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ boost::math::tools::max_value() }, { static_cast(3) }, static_cast(0.5f), &error_rate, throw_policy()), std::overflow_error); + BOOST_MATH_IF_CONSTEXPR(std::numeric_limits::has_infinity) + { + BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ 2 }, { 3 }, static_cast(1e10), &error_rate, nothrow_policy()), std::numeric_limits::infinity()); + BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ static_cast(2) }, { static_cast(3) }, static_cast(1e10), &error_rate, nothrow_policy()), std::numeric_limits::infinity()); + BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ 2 }, { 3 }, boost::math::tools::max_value(), &error_rate, nothrow_policy()), std::numeric_limits::infinity()); + BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ static_cast(2) }, { static_cast(3) }, boost::math::tools::max_value(), &error_rate, nothrow_policy()), std::numeric_limits::infinity()); + BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ boost::math::tools::max_value() }, { static_cast(3) }, static_cast(0.5f), &error_rate, nothrow_policy()), std::numeric_limits::infinity()); + } +} + template void test_spots(T z, const char* type_name) { @@ -309,4 +330,5 @@ void test_spots(T z, const char* type_name) test_spots_1F2(z, type_name); test_spots_2F2(z, type_name); test_spots_2F1(z, type_name); + test_special_cases(z, type_name); } From 31b1bc1b1fcc658ad969f058b31e93a08b37bb10 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Thu, 11 Jun 2026 19:10:43 +0100 Subject: [PATCH 88/96] pFq: disable one new test for real_concept. --- test/test_pFq.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/test/test_pFq.hpp b/test/test_pFq.hpp index e9e5bac67e..9cfbac30ba 100644 --- a/test/test_pFq.hpp +++ b/test/test_pFq.hpp @@ -308,7 +308,11 @@ void test_special_cases(T, const char*) BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ static_cast(2) }, { static_cast(3) }, static_cast(1e10), &error_rate, throw_policy()), std::overflow_error); BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ 2 }, { 3 }, boost::math::tools::max_value(), &error_rate, throw_policy()), std::overflow_error); BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ static_cast(2) }, { static_cast(3) }, boost::math::tools::max_value(), &error_rate, throw_policy()), std::overflow_error); - BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ boost::math::tools::max_value() }, { static_cast(3) }, static_cast(0.5f), &error_rate, throw_policy()), std::overflow_error); + BOOST_MATH_IF_CONSTEXPR(std::numeric_limits::is_specialized) + { + // Can't figure out how to make this work with real_concept (we get an evaluation_error): + BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ boost::math::tools::max_value() }, { static_cast(3) }, static_cast(0.5f), &error_rate, throw_policy()), std::overflow_error); + } BOOST_MATH_IF_CONSTEXPR(std::numeric_limits::has_infinity) { BOOST_CHECK_EQUAL(boost::math::hypergeometric_pFq({ 2 }, { 3 }, static_cast(1e10), &error_rate, nothrow_policy()), std::numeric_limits::infinity()); From 66d86812abe67da0331d8576c137e370ba1ed852 Mon Sep 17 00:00:00 2001 From: Warren Weckesser Date: Thu, 18 Jun 2026 19:57:29 -0400 Subject: [PATCH 89/96] DOC: Include the use of the bootstrap script to build b2 in README.md [ci skip] --- README.md | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 39b0bddd63..b44755ec55 100644 --- a/README.md +++ b/README.md @@ -138,8 +138,12 @@ The Boost Math Library is located in `libs/math/`. ### Running tests -First, make sure you are in `libs/math/test`. -You can either run all the tests listed in `Jamfile.v2` or run a single test: +The Boost build program `b2` is used to run the tests. To build `b2`, from the +top-level Boost directory, run either `./bootstrap.sh` (Linux/MacOS) or +`bootstrap.bat` (Windows). + +Then move to `libs/math/test`. You can either run all the tests listed in +`Jamfile.v2` or run a single test: test$ ../../../b2 <- run all tests test$ ../../../b2 static_assert_test <- single test From ad99070a188fceb702c0c64c2ce07fe26e9fe4dc Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Fri, 19 Jun 2026 19:12:08 +0100 Subject: [PATCH 90/96] Comment out dead code, tweak testing for more coverage. --- .../detail/hypergeometric_pFq_checked_series.hpp | 5 ++++- test/test_pFq.hpp | 4 ++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp b/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp index 58c0d6b1ea..e37b0ffd6e 100644 --- a/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp +++ b/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp @@ -243,6 +243,9 @@ log_scale += log_scaling_factor; local_scaling += log_scaling_factor; } + /* + * rescaling to avoid underflow is pointless here, given that the first term is 1. + * else if ((fabs(abs_result) < lower_limit) && (fabs(abs_result) * scaling_factor < upper_limit)) { abs_result *= scaling_factor; @@ -251,7 +254,7 @@ log_scale -= log_scaling_factor; local_scaling -= log_scaling_factor; } - + */ if ((abs(result * tol) > abs(term)) && (abs(term0) > abs(term))) break; if (abs_result * tol > abs(result)) diff --git a/test/test_pFq.hpp b/test/test_pFq.hpp index 9cfbac30ba..df054caad3 100644 --- a/test/test_pFq.hpp +++ b/test/test_pFq.hpp @@ -301,8 +301,8 @@ void test_spots_2F2(T, const char*) template void test_special_cases(T, const char*) { - typedef boost::math::policies::policy> nothrow_policy; - typedef boost::math::policies::policy> throw_policy; + typedef boost::math::policies::policy, boost::math::policies::promote_double > nothrow_policy; + typedef boost::math::policies::policy, boost::math::policies::promote_double > throw_policy; T error_rate; BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ 2 }, { 3 }, static_cast(1e10), &error_rate, throw_policy()), std::overflow_error); BOOST_CHECK_THROW(boost::math::hypergeometric_pFq({ static_cast(2) }, { static_cast(3) }, static_cast(1e10), &error_rate, throw_policy()), std::overflow_error); From 45fe9991ca5fc536814641bd0293b6abfb6c0fb5 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sun, 21 Jun 2026 12:01:53 +0100 Subject: [PATCH 91/96] Try and get gcov to register lines we know are hit. --- .../detail/hypergeometric_pFq_checked_series.hpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp b/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp index e37b0ffd6e..895c586482 100644 --- a/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp +++ b/include/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp @@ -132,11 +132,15 @@ std::uintmax_t k = 0; Real upper_limit(sqrt(boost::math::tools::max_value())), diff; if ((tools::max_value() / fabs(z) < upper_limit)) - upper_limit = tools::max_value() / fabs(z); + { + upper_limit = tools::max_value() / fabs(z); + } for (auto pa = aj.begin(); pa != aj.end(); ++pa) { if (tools::max_value() / fabs(*pa) < upper_limit) - upper_limit = tools::max_value() / fabs(*pa); + { + upper_limit = tools::max_value() / fabs(*pa); + } } Real lower_limit(1 / upper_limit); long long log_scaling_factor = lltrunc(boost::math::tools::log_max_value(), nothrow_policy()) - 2; From 135e16b8e6c3750dc95dfd06bb411b8b4eadea4d Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Tue, 30 Jun 2026 17:31:30 -0400 Subject: [PATCH 92/96] Add reproducer test set --- test/Jamfile.v2 | 1 + test/git_issue_1412.cpp | 23 +++++++++++++++++++++++ 2 files changed, 24 insertions(+) create mode 100644 test/git_issue_1412.cpp diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 5b67e474df..3ca1f4eb5e 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -202,6 +202,7 @@ test-suite special_fun : [ run git_issue_1255.cpp ] [ run git_issue_1247.cpp ] [ run git_issue_1308.cpp /boost/test//boost_unit_test_framework ] + [ run git_issue_1412.cpp ] [ run special_functions_test.cpp /boost/test//boost_unit_test_framework ] [ run test_airy.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ] [ run test_bessel_j.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ] diff --git a/test/git_issue_1412.cpp b/test/git_issue_1412.cpp new file mode 100644 index 0000000000..74735e695a --- /dev/null +++ b/test/git_issue_1412.cpp @@ -0,0 +1,23 @@ +// Copyright Matt Borland 2026 +// 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) +// +// Regression test for https://github.com/boostorg/math/issues/1412 + +#define BOOST_MATH_NO_EXCEPTIONS + +// Claim the include guard so any later #include expands to nothing and std::setprecision is unavailable. +// The original reproducer with _SSTREAM broke which iso outside our purview +#define _LIBCPP_IOMANIP // libc++ +#define _GLIBCXX_IOMANIP 1 // libstdc++ + +#include +#include + +int main() +{ + // Exercise a special function so the no-exceptions header path is + // instantiated, not just parsed. + return static_cast(boost::math::laguerre(1u, 0u, 0.5)); +} From bad351c738732cdce7129555f2eeb754ff8adc26 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Tue, 30 Jun 2026 17:31:48 -0400 Subject: [PATCH 93/96] Add patch from issue --- include/boost/math/policies/error_handling.hpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/include/boost/math/policies/error_handling.hpp b/include/boost/math/policies/error_handling.hpp index d7ba694564..adcb56f3e5 100644 --- a/include/boost/math/policies/error_handling.hpp +++ b/include/boost/math/policies/error_handling.hpp @@ -26,9 +26,9 @@ #include #include #include +#ifndef BOOST_MATH_NO_EXCEPTIONS #include #include -#ifndef BOOST_MATH_NO_EXCEPTIONS #include #endif #include @@ -95,6 +95,8 @@ T user_indeterminate_result_error(const char* function, const char* message, con namespace detail { +#ifndef BOOST_MATH_NO_EXCEPTIONS +// Only used to build exception messages; relies on template inline std::string prec_format(const T& val) { @@ -124,7 +126,9 @@ inline std::string prec_format(const std::float128_t& val) return std::string(buffer, r.ptr); } -#endif +#endif // BOOST_MATH_USE_CHARCONV_FOR_CONVERSION + +#endif // BOOST_MATH_NO_EXCEPTIONS inline void replace_all_in_string(std::string& result, const char* what, const char* with) { From 2d00b8cf4aa4c7733e0c8059b798395120d45624 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Wed, 1 Jul 2026 08:48:06 -0400 Subject: [PATCH 94/96] Check all of special functions can be included under same conditions --- test/Jamfile.v2 | 1 + test/git_issue_1412_pt_2.cpp | 23 +++++++++++++++++++++++ 2 files changed, 24 insertions(+) create mode 100644 test/git_issue_1412_pt_2.cpp diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 3ca1f4eb5e..18d5d9ae01 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -203,6 +203,7 @@ test-suite special_fun : [ run git_issue_1247.cpp ] [ run git_issue_1308.cpp /boost/test//boost_unit_test_framework ] [ run git_issue_1412.cpp ] + [ run git_issue_1412_pt_2.cpp ] [ run special_functions_test.cpp /boost/test//boost_unit_test_framework ] [ run test_airy.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ] [ run test_bessel_j.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ] diff --git a/test/git_issue_1412_pt_2.cpp b/test/git_issue_1412_pt_2.cpp new file mode 100644 index 0000000000..5a0519fe3b --- /dev/null +++ b/test/git_issue_1412_pt_2.cpp @@ -0,0 +1,23 @@ +// Copyright Matt Borland 2026 +// 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) +// +// Regression test for https://github.com/boostorg/math/issues/1412 + +#define BOOST_MATH_NO_EXCEPTIONS + +// Claim the include guard so any later #include expands to nothing and std::setprecision is unavailable. +// The original reproducer with _SSTREAM broke which iso outside our purview +#define _LIBCPP_IOMANIP // libc++ +#define _GLIBCXX_IOMANIP 1 // libstdc++ + +#include +#include + +int main() +{ + // Exercise a special function so the no-exceptions header path is + // instantiated, not just parsed. + return static_cast(boost::math::laguerre(1u, 0u, 0.5)); +} From e06a6a657cd37603a370f77bef39745e55fd84c5 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Wed, 1 Jul 2026 09:41:49 -0400 Subject: [PATCH 95/96] Fix UB in sorting results with the cubic root finder --- include/boost/math/tools/cubic_roots.hpp | 26 ++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/include/boost/math/tools/cubic_roots.hpp b/include/boost/math/tools/cubic_roots.hpp index 451c5de813..1097ad5ce4 100644 --- a/include/boost/math/tools/cubic_roots.hpp +++ b/include/boost/math/tools/cubic_roots.hpp @@ -6,11 +6,33 @@ #define BOOST_MATH_TOOLS_CUBIC_ROOTS_HPP #include #include +#include #include #include namespace boost::math::tools { +namespace detail { + +// Orders finite roots ascending and moves NaN entries to the back. +// Sorting NaNs is UB when using vanilla operator< (comparisons to NaN are always false), +// so we need to provide our own sort function that pushes NaN to the end +template +bool roots_less(const Real& lhs, const Real& rhs) +{ + if ((boost::math::isnan)(lhs)) + { + return false; + } + if ((boost::math::isnan)(rhs)) + { + return true; + } + return lhs < rhs; +} + +} // namespace detail + // Solves ax^3 + bx^2 + cx + d = 0. // Only returns the real roots, as types get weird for real coefficients and // complex roots. Follows Numerical Recipes, Chapter 5, section 6. NB: A better @@ -55,7 +77,7 @@ std::array cubic_roots(Real a, Real b, Real c, Real d) { roots[0] = x0; roots[1] = x1; roots[2] = 0; - std::sort(roots.begin(), roots.end()); + std::sort(roots.begin(), roots.end(), detail::roots_less); return roots; } Real p = b / a; @@ -115,7 +137,7 @@ std::array cubic_roots(Real a, Real b, Real c, Real d) { } } } - std::sort(roots.begin(), roots.end()); + std::sort(roots.begin(), roots.end(), detail::roots_less); return roots; } From 9fbea69ab36cdc49b4e3ed4858561652d3ffeb07 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Wed, 1 Jul 2026 10:34:14 -0400 Subject: [PATCH 96/96] Use the instantiate driver for better coverage --- test/git_issue_1412_pt_2.cpp | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/test/git_issue_1412_pt_2.cpp b/test/git_issue_1412_pt_2.cpp index 5a0519fe3b..a591fdf346 100644 --- a/test/git_issue_1412_pt_2.cpp +++ b/test/git_issue_1412_pt_2.cpp @@ -4,20 +4,23 @@ // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) // // Regression test for https://github.com/boostorg/math/issues/1412 +// +// Part 2: Test all of #include as that is what libc++ is going to use #define BOOST_MATH_NO_EXCEPTIONS // Claim the include guard so any later #include expands to nothing and std::setprecision is unavailable. -// The original reproducer with _SSTREAM broke which iso outside our purview +// The original reproducer with _SSTREAM broke which is outside our purview. #define _LIBCPP_IOMANIP // libc++ #define _GLIBCXX_IOMANIP 1 // libstdc++ -#include -#include +#include "compile_test/instantiate.hpp" -int main() +int main(int argc, char* []) { - // Exercise a special function so the no-exceptions header path is - // instantiated, not just parsed. - return static_cast(boost::math::laguerre(1u, 0u, 0.5)); + if (argc > 10000) + { + instantiate(double(0)); + instantiate_mixed(double(0)); + } }