From a57fb0a894238ecb2d0653b9df178c550e8c2c66 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Mon, 5 Feb 2024 19:16:33 +0000 Subject: [PATCH 1/2] Reinstate root finding protection against huge jumps. Apply error handling more rigorously to any root finding client. Mark evaluation_error's as not reachable for code coverage. Fixes https://github.com/boostorg/math/issues/184. --- include/boost/math/distributions/binomial.hpp | 3 +- .../boost/math/distributions/chi_squared.hpp | 5 ++- .../distributions/detail/generic_mode.hpp | 19 ++++------ .../distributions/detail/generic_quantile.hpp | 5 ++- .../detail/inv_discrete_quantile.hpp | 9 +++-- .../boost/math/distributions/find_scale.hpp | 13 +++---- .../math/distributions/inverse_gaussian.hpp | 23 ++++++++---- .../math/distributions/kolmogorov_smirnov.hpp | 24 +++++++++---- .../math/distributions/non_central_beta.hpp | 34 +++++++----------- .../distributions/non_central_chi_squared.hpp | 34 +++++++----------- .../math/distributions/non_central_t.hpp | 34 +++++++----------- .../boost/math/distributions/skew_normal.hpp | 18 +++++++--- .../boost/math/distributions/students_t.hpp | 5 ++- .../boost/math/special_functions/legendre.hpp | 5 +++ include/boost/math/tools/roots.hpp | 8 ++++- include/boost/math/tools/toms748_solve.hpp | 6 ++-- test/Jamfile.v2 | 1 + test/git_issue_184.cpp | 35 +++++++++++++++++++ 18 files changed, 160 insertions(+), 121 deletions(-) create mode 100644 test/git_issue_184.cpp diff --git a/include/boost/math/distributions/binomial.hpp b/include/boost/math/distributions/binomial.hpp index 05ef01ab8e..cf7451104b 100644 --- a/include/boost/math/distributions/binomial.hpp +++ b/include/boost/math/distributions/binomial.hpp @@ -258,7 +258,7 @@ namespace boost typedef typename Policy::discrete_quantile_type discrete_quantile_type; std::uintmax_t max_iter = policies::get_max_root_iterations(); - return detail::inverse_discrete_quantile( + result = detail::inverse_discrete_quantile( dist, comp ? q : p, comp, @@ -267,6 +267,7 @@ namespace boost RealType(1), discrete_quantile_type(), max_iter); + return result; } // quantile } diff --git a/include/boost/math/distributions/chi_squared.hpp b/include/boost/math/distributions/chi_squared.hpp index 7b65a0da68..f5daddc0ad 100644 --- a/include/boost/math/distributions/chi_squared.hpp +++ b/include/boost/math/distributions/chi_squared.hpp @@ -327,9 +327,8 @@ RealType chi_squared_distribution::find_degrees_of_freedom( RealType result = r.first + (r.second - r.first) / 2; if(max_iter >= policies::get_max_root_iterations()) { - 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" - " or the answer is infinite. Current best guess is %1%", result, Policy()); + policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE + " either there is no answer to how many degrees of freedom are required or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE } return result; } diff --git a/include/boost/math/distributions/detail/generic_mode.hpp b/include/boost/math/distributions/detail/generic_mode.hpp index 96d21b9b0b..19c8b2af01 100644 --- a/include/boost/math/distributions/detail/generic_mode.hpp +++ b/include/boost/math/distributions/detail/generic_mode.hpp @@ -47,9 +47,7 @@ typename Dist::value_type generic_find_mode(const Dist& dist, typename Dist::val // Oops we don't know how to handle this, or even in which // direction we should move in, treat as an evaluation error: // - return policies::raise_evaluation_error( - function, - "Could not locate a starting location for the search for the mode, original guess was %1%", guess, policy_type()); + return policies::raise_evaluation_error(function, "Could not locate a starting location for the search for the mode, original guess was %1%", guess, policy_type()); // LCOV_EXCL_LINE } do { @@ -82,11 +80,9 @@ typename Dist::value_type generic_find_mode(const Dist& dist, typename Dist::val max_iter).first; if(max_iter >= policies::get_max_root_iterations()) { - return policies::raise_evaluation_error( - function, - "Unable to locate solution in a reasonable time:" - " either there is no answer to the mode of the distribution" - " or the answer is infinite. Current best guess is %1%", result, policy_type()); + return policies::raise_evaluation_error(function, // LCOV_EXCL_LINE + "Unable to locate solution in a reasonable time: either there is no answer to the mode of the distribution" // LCOV_EXCL_LINE + " or the answer is infinite. Current best guess is %1%", result, policy_type()); // LCOV_EXCL_LINE } return result; } @@ -135,11 +131,8 @@ typename Dist::value_type generic_find_mode_01(const Dist& dist, typename Dist:: max_iter).first; if(max_iter >= policies::get_max_root_iterations()) { - return policies::raise_evaluation_error( - function, - "Unable to locate solution in a reasonable time:" - " either there is no answer to the mode of the distribution" - " or the answer is infinite. Current best guess is %1%", result, policy_type()); + return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE + " either there is no answer to the mode of the distribution or the answer is infinite. Current best guess is %1%", result, policy_type()); // LCOV_EXCL_LINE } return result; } diff --git a/include/boost/math/distributions/detail/generic_quantile.hpp b/include/boost/math/distributions/detail/generic_quantile.hpp index 67d3f04f88..438ac952f0 100644 --- a/include/boost/math/distributions/detail/generic_quantile.hpp +++ b/include/boost/math/distributions/detail/generic_quantile.hpp @@ -84,9 +84,8 @@ typename Dist::value_type generic_quantile(const Dist& dist, const typename Dist value_type 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:" - " either there is no answer to quantile" - " or the answer is infinite. Current best guess is %1%", result, forwarding_policy()); + return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE + " either there is no answer to quantile or the answer is infinite. Current best guess is %1%", result, forwarding_policy()); // LCOV_EXCL_LINE } return result; } diff --git a/include/boost/math/distributions/detail/inv_discrete_quantile.hpp b/include/boost/math/distributions/detail/inv_discrete_quantile.hpp index 2f5cae7432..739a866660 100644 --- a/include/boost/math/distributions/detail/inv_discrete_quantile.hpp +++ b/include/boost/math/distributions/detail/inv_discrete_quantile.hpp @@ -215,7 +215,7 @@ typename Dist::value_type while(((boost::math::sign)(fb) == (boost::math::sign)(fa)) && (a != b)) { if(count == 0) - return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, policy_type()); + return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, policy_type()); // LCOV_EXCL_LINE a = b; fa = fb; b *= multiplier; @@ -242,7 +242,7 @@ typename Dist::value_type return 0; } if(count == 0) - return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, policy_type()); + return policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, policy_type()); // LCOV_EXCL_LINE b = a; fb = fa; a /= multiplier; @@ -276,6 +276,11 @@ typename Dist::value_type // std::pair r = toms748_solve(f, a, b, fa, fb, tol, count, policy_type()); max_iter += count; + 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 + " either there is no answer to quantile or the answer is infinite. Current best guess is %1%", r.first, policy_type()); // LCOV_EXCL_LINE + } BOOST_MATH_INSTRUMENT_CODE("max_iter = " << max_iter << " count = " << count); return (r.first + r.second) / 2; } diff --git a/include/boost/math/distributions/find_scale.hpp b/include/boost/math/distributions/find_scale.hpp index 25551b87da..cd08ce95a2 100644 --- a/include/boost/math/distributions/find_scale.hpp +++ b/include/boost/math/distributions/find_scale.hpp @@ -80,9 +80,7 @@ namespace boost if (result <= 0) { // If policy isn't to throw, return the scale <= 0. - policies::raise_evaluation_error(function, - "Computed scale (%1%) is <= 0!" " Was the complement intended?", - result, Policy()); + policies::raise_evaluation_error(function, "Computed scale (%1%) is <= 0!" " Was the complement intended?", result, Policy()); // LCOV_EXCL_LINE } return result; } // template find_scale @@ -140,9 +138,7 @@ namespace boost // ( z - location) / (quantile(complement(Dist(), q)) if (result <= 0) { // If policy isn't to throw, return the scale <= 0. - policies::raise_evaluation_error(function, - "Computed scale (%1%) is <= 0!" " Was the complement intended?", - result, Policy()); + policies::raise_evaluation_error(function, "Computed scale (%1%) is <= 0!" " Was the complement intended?", result, Policy()); // LCOV_EXCL_LINE } return result; } // template typename Dist::value_type find_scale @@ -192,9 +188,8 @@ namespace boost // ( z - location) / (quantile(complement(Dist(), q)) if (result <= 0) { // If policy isn't to throw, return the scale <= 0. - policies::raise_evaluation_error(function, - "Computed scale (%1%) is <= 0!" " Was the complement intended?", - result, policies::policy<>()); // This is only the default policy - also Want a version with Policy here. + policies::raise_evaluation_error(function, "Computed scale (%1%) is <= 0!" " Was the complement intended?", // LCOV_EXCL_LINE + result, policies::policy<>()); // This is only the default policy - also Want a version with Policy here. LCOV_EXCL_LINE } return result; } // template typename Dist::value_type find_scale diff --git a/include/boost/math/distributions/inverse_gaussian.hpp b/include/boost/math/distributions/inverse_gaussian.hpp index c8e46913af..b31d1c9257 100644 --- a/include/boost/math/distributions/inverse_gaussian.hpp +++ b/include/boost/math/distributions/inverse_gaussian.hpp @@ -388,11 +388,16 @@ inline RealType quantile(const inverse_gaussian_distribution& // digits used to control how accurate to try to make the result. // To allow user to control accuracy versus speed, int get_digits = policies::digits();// get digits from policy, - std::uintmax_t m = policies::get_max_root_iterations(); // and max iterations. + std::uintmax_t max_iter = policies::get_max_root_iterations(); // and max iterations. using boost::math::tools::newton_raphson_iterate; result = - newton_raphson_iterate(inverse_gaussian_quantile_functor(dist, p), guess, min, max, get_digits, m); - return result; + newton_raphson_iterate(inverse_gaussian_quantile_functor(dist, p), 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 + " either there is no answer to quantile or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE + } + return result; } // quantile template @@ -459,11 +464,15 @@ inline RealType quantile(const complemented2_type::digits; // Maximum possible binary digits accuracy for type T. // digits used to control how accurate to try to make the result. int get_digits = policies::digits(); - std::uintmax_t m = policies::get_max_root_iterations(); + std::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, m); - return result; + result = newton_raphson_iterate(inverse_gaussian_quantile_complement_functor(c.dist, q), 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 + " either there is no answer to quantile or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE + } + return result; } // quantile template diff --git a/include/boost/math/distributions/kolmogorov_smirnov.hpp b/include/boost/math/distributions/kolmogorov_smirnov.hpp index 6ae80632d3..0934344134 100644 --- a/include/boost/math/distributions/kolmogorov_smirnov.hpp +++ b/include/boost/math/distributions/kolmogorov_smirnov.hpp @@ -379,10 +379,16 @@ inline RealType quantile(const kolmogorov_smirnov_distribution RealType k = detail::kolmogorov_smirnov_quantile_guess(p) / sqrt(n); const int get_digits = policies::digits();// get digits from policy, - std::uintmax_t m = policies::get_max_root_iterations(); // and max iterations. + std::uintmax_t max_iter = policies::get_max_root_iterations(); // and max iterations. - return tools::newton_raphson_iterate(detail::kolmogorov_smirnov_quantile_functor(dist, p), - k, RealType(0), RealType(1), get_digits, m); + RealType result = tools::newton_raphson_iterate(detail::kolmogorov_smirnov_quantile_functor(dist, p), + k, RealType(0), RealType(1), 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 + " either there is no answer to quantile or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE + } + return result; } // quantile template @@ -403,11 +409,17 @@ inline RealType quantile(const complemented2_type();// get digits from policy, - std::uintmax_t m = policies::get_max_root_iterations(); // and max iterations. + std::uintmax_t max_iter = policies::get_max_root_iterations(); // and max iterations. - return tools::newton_raphson_iterate( + RealType result = tools::newton_raphson_iterate( detail::kolmogorov_smirnov_complementary_quantile_functor(dist, p), - k, RealType(0), RealType(1), get_digits, m); + k, RealType(0), RealType(1), 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 + " either there is no answer to quantile or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE + } + return result; } // quantile (complemented) template diff --git a/include/boost/math/distributions/non_central_beta.hpp b/include/boost/math/distributions/non_central_beta.hpp index 944f1dc242..b32a605f21 100644 --- a/include/boost/math/distributions/non_central_beta.hpp +++ b/include/boost/math/distributions/non_central_beta.hpp @@ -109,9 +109,7 @@ namespace boost } if(static_cast(count + i - k) > max_iter) { - return policies::raise_evaluation_error( - "cdf(non_central_beta_distribution<%1%>, %1%)", - "Series did not converge, closest value was %1%", sum, pol); + return policies::raise_evaluation_error("cdf(non_central_beta_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE } } return sum; @@ -190,9 +188,7 @@ namespace boost } if(static_cast(i - k) > max_iter) { - return policies::raise_evaluation_error( - "cdf(non_central_beta_distribution<%1%>, %1%)", - "Series did not converge, closest value was %1%", sum, pol); + return policies::raise_evaluation_error("cdf(non_central_beta_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE } last_term = term; } @@ -206,9 +202,7 @@ namespace boost } if(static_cast(count + k - i) > max_iter) { - return policies::raise_evaluation_error( - "cdf(non_central_beta_distribution<%1%>, %1%)", - "Series did not converge, closest value was %1%", sum, pol); + return policies::raise_evaluation_error("cdf(non_central_beta_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE } pois *= i / l2; beta -= xterm; @@ -324,7 +318,7 @@ namespace boost { if(count == 0) { - b = policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, pol); + b = policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, pol); // LCOV_EXCL_LINE return std::make_pair(a, b); } // @@ -362,7 +356,7 @@ namespace boost } if(count == 0) { - a = policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol); + a = policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol); // LCOV_EXCL_LINE return std::make_pair(a, b); } // @@ -508,12 +502,14 @@ namespace boost if(max_iter >= policies::get_max_root_iterations()) { + // LCOV_EXCL_START return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time:" " either there is no answer to quantile of the non central beta distribution" " or the answer is infinite. Current best guess is %1%", policies::checked_narrowing_cast( result, function), Policy()); + // LCOV_EXCL_STOP } return policies::checked_narrowing_cast( result, @@ -583,9 +579,7 @@ namespace boost } if(static_cast(count + i - k) > max_iter) { - return policies::raise_evaluation_error( - "pdf(non_central_beta_distribution<%1%>, %1%)", - "Series did not converge, closest value was %1%", sum, pol); + return policies::raise_evaluation_error("pdf(non_central_beta_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE } } return sum; @@ -817,10 +811,8 @@ namespace boost typedef typename Policy::assert_undefined_type assert_type; static_assert(assert_type::value == 0, "Assert type is undefined."); - return policies::raise_evaluation_error( - function, - "This function is not yet implemented, the only sensible result is %1%.", - std::numeric_limits::quiet_NaN(), Policy()); // infinity? + return policies::raise_evaluation_error(function, "This function is not yet implemented, the only sensible result is %1%.", // LCOV_EXCL_LINE + std::numeric_limits::quiet_NaN(), Policy()); // infinity? LCOV_EXCL_LINE } template @@ -830,10 +822,8 @@ namespace boost typedef typename Policy::assert_undefined_type assert_type; static_assert(assert_type::value == 0, "Assert type is undefined."); - return policies::raise_evaluation_error( - function, - "This function is not yet implemented, the only sensible result is %1%.", - std::numeric_limits::quiet_NaN(), Policy()); // infinity? + return policies::raise_evaluation_error(function, "This function is not yet implemented, the only sensible result is %1%.", // LCOV_EXCL_LINE + std::numeric_limits::quiet_NaN(), Policy()); // infinity? LCOV_EXCL_LINE } // kurtosis_excess template diff --git a/include/boost/math/distributions/non_central_chi_squared.hpp b/include/boost/math/distributions/non_central_chi_squared.hpp index c85c1d08d1..f59be9932c 100644 --- a/include/boost/math/distributions/non_central_chi_squared.hpp +++ b/include/boost/math/distributions/non_central_chi_squared.hpp @@ -101,9 +101,7 @@ namespace boost } //Error check: if(static_cast(i-k) >= max_iter) - return policies::raise_evaluation_error( - "cdf(non_central_chi_squared_distribution<%1%>, %1%)", - "Series did not converge, closest value was %1%", sum, pol); + return policies::raise_evaluation_error("cdf(non_central_chi_squared_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE // // Now backwards iteration: the gamma // function recurrences are unstable in this @@ -175,9 +173,7 @@ namespace boost } //Error check: if(static_cast(i) >= max_iter) - return policies::raise_evaluation_error( - "cdf(non_central_chi_squared_distribution<%1%>, %1%)", - "Series did not converge, closest value was %1%", sum, pol); + return policies::raise_evaluation_error("cdf(non_central_chi_squared_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE return sum; } @@ -274,9 +270,7 @@ namespace boost //Error check: if(static_cast(i) >= max_iter) - return policies::raise_evaluation_error( - "cdf(non_central_chi_squared_distribution<%1%>, %1%)", - "Series did not converge, closest value was %1%", sum, pol); + return policies::raise_evaluation_error("cdf(non_central_chi_squared_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE return sum; } @@ -305,9 +299,7 @@ namespace boost if(pois / sum < errtol) break; if(static_cast(i - k) >= max_iter) - return policies::raise_evaluation_error( - "pdf(non_central_chi_squared_distribution<%1%>, %1%)", - "Series did not converge, closest value was %1%", sum, pol); + return policies::raise_evaluation_error("pdf(non_central_chi_squared_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE pois *= l2 * x2 / ((i + 1) * (n2 + i)); } for(long long i = k - 1; i >= 0; --i) @@ -581,9 +573,8 @@ namespace boost // // Can't 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%", - RealType(std::numeric_limits::quiet_NaN()), Policy()); + return policies::raise_evaluation_error(function, "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%", // LCOV_EXCL_LINE + RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } degrees_of_freedom_finder f(lam, x, p < q ? p : q, p < q ? false : true); tools::eps_tolerance tol(policies::digits()); @@ -600,8 +591,8 @@ namespace boost 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:" - " or there is no answer to problem. Current best guess is %1%", result, Policy()); + 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; } @@ -637,9 +628,8 @@ namespace boost // // Can't do a thing if one of p and q is zero: // - return policies::raise_evaluation_error(function, - "Can't find non centrality parameter when the probability is 0 or 1, only possible answer is %1%", - RealType(std::numeric_limits::quiet_NaN()), Policy()); + return policies::raise_evaluation_error(function, "Can't find non centrality parameter when the probability is 0 or 1, only possible answer is %1%", // LCOV_EXCL_LINE + RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } non_centrality_finder f(v, x, p < q ? p : q, p < q ? false : true); tools::eps_tolerance tol(policies::digits()); @@ -656,8 +646,8 @@ namespace boost 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:" - " or there is no answer to problem. Current best guess is %1%", result, Policy()); + 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; } diff --git a/include/boost/math/distributions/non_central_t.hpp b/include/boost/math/distributions/non_central_t.hpp index 48483f0cae..4651268d59 100644 --- a/include/boost/math/distributions/non_central_t.hpp +++ b/include/boost/math/distributions/non_central_t.hpp @@ -97,9 +97,7 @@ namespace boost ++count; if(count > max_iter) { - return policies::raise_evaluation_error( - "cdf(non_central_t_distribution<%1%>, %1%)", - "Series did not converge, closest value was %1%", sum, pol); + return policies::raise_evaluation_error("cdf(non_central_t_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE } } return sum; @@ -195,9 +193,7 @@ namespace boost last_term = term; if(count > max_iter) { - return policies::raise_evaluation_error( - "cdf(non_central_t_distribution<%1%>, %1%)", - "Series did not converge, closest value was %1%", sum, pol); + return policies::raise_evaluation_error("cdf(non_central_t_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE } ++count; } @@ -427,9 +423,7 @@ namespace boost ++count; if(count > max_iter) { - return policies::raise_evaluation_error( - "pdf(non_central_t_distribution<%1%>, %1%)", - "Series did not converge, closest value was %1%", sum, pol); + return policies::raise_evaluation_error("pdf(non_central_t_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE } } BOOST_MATH_INSTRUMENT_VARIABLE(sum); @@ -444,9 +438,7 @@ namespace boost ++count; if(count > max_iter) { - return policies::raise_evaluation_error( - "pdf(non_central_t_distribution<%1%>, %1%)", - "Series did not converge, closest value was %1%", sum, pol); + return policies::raise_evaluation_error("pdf(non_central_t_distribution<%1%>, %1%)", "Series did not converge, closest value was %1%", sum, pol); // LCOV_EXCL_LINE } } BOOST_MATH_INSTRUMENT_VARIABLE(sum); @@ -645,9 +637,8 @@ namespace boost // // Can't 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%", - RealType(std::numeric_limits::quiet_NaN()), Policy()); + return policies::raise_evaluation_error(function, "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%", // LCOV_EXCL_LINE + RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } t_degrees_of_freedom_finder f(delta, x, p < q ? p : q, p < q ? false : true); tools::eps_tolerance tol(policies::digits()); @@ -661,8 +652,8 @@ namespace boost 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:" - " or there is no answer to problem. Current best guess is %1%", result, Policy()); + 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; } @@ -698,9 +689,8 @@ namespace boost // // Can't do a thing if one of p and q is zero: // - return policies::raise_evaluation_error(function, - "Can't find non-centrality parameter when the probability is 0 or 1, only possible answer is %1%", - RealType(std::numeric_limits::quiet_NaN()), Policy()); + return policies::raise_evaluation_error(function, "Can't find non-centrality parameter when the probability is 0 or 1, only possible answer is %1%", // LCOV_EXCL_LINE + RealType(std::numeric_limits::quiet_NaN()), Policy()); // LCOV_EXCL_LINE } t_non_centrality_finder f(v, x, p < q ? p : q, p < q ? false : true); tools::eps_tolerance tol(policies::digits()); @@ -719,8 +709,8 @@ namespace boost 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:" - " or there is no answer to problem. Current best guess is %1%", result, Policy()); + 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; } diff --git a/include/boost/math/distributions/skew_normal.hpp b/include/boost/math/distributions/skew_normal.hpp index 0f0ffea8a8..f3cc5a6b59 100644 --- a/include/boost/math/distributions/skew_normal.hpp +++ b/include/boost/math/distributions/skew_normal.hpp @@ -561,12 +561,17 @@ namespace boost{ namespace math{ } const int get_digits = policies::digits();// get digits from policy, - std::uintmax_t m = policies::get_max_root_iterations(); // and max iterations. + std::uintmax_t max_iter = policies::get_max_root_iterations(); // and max iterations. skew_normal_distribution helper(0, 1, shape); result = tools::newton_raphson_iterate(detail::skew_normal_mode_functor(helper), result, - search_min, search_max, get_digits, m); + search_min, search_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 + " either there is no answer to quantile or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE + } result = result*scale + location; @@ -680,10 +685,15 @@ namespace boost{ namespace math{ const RealType search_max = support(dist).second; const int get_digits = policies::digits();// get digits from policy, - std::uintmax_t m = policies::get_max_root_iterations(); // and max iterations. + std::uintmax_t max_iter = policies::get_max_root_iterations(); // and max iterations. result = tools::newton_raphson_iterate(detail::skew_normal_quantile_functor(dist, p), result, - search_min, search_max, get_digits, m); + search_min, search_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: either there is no answer to quantile" // LCOV_EXCL_LINE + " or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE + } return result; } // quantile diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index 321c4751b5..b01b8aa0fc 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -335,9 +335,8 @@ RealType students_t_distribution::find_degrees_of_freedom( RealType result = r.first + (r.second - r.first) / 2; if(max_iter >= policies::get_max_root_iterations()) { - 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" - " or the answer is infinite. Current best guess is %1%", result, Policy()); + 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 } return result; } diff --git a/include/boost/math/special_functions/legendre.hpp b/include/boost/math/special_functions/legendre.hpp index 6af798e0d9..5e7f0029a6 100644 --- a/include/boost/math/special_functions/legendre.hpp +++ b/include/boost/math/special_functions/legendre.hpp @@ -206,6 +206,11 @@ std::vector legendre_p_zeros_imp(int n, const Policy& pol) lower_bound, upper_bound, policies::digits(), number_of_iterations); + if (number_of_iterations >= policies::get_max_root_iterations()) + { + policies::raise_evaluation_error("legendre_p_zeros<%1%>", "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE + " either there is no answer or the answer is infinite. Current best guess is %1%", x_nk, Policy()); // LCOV_EXCL_LINE + } BOOST_MATH_ASSERT(lower_bound < x_nk); BOOST_MATH_ASSERT(upper_bound > x_nk); diff --git a/include/boost/math/tools/roots.hpp b/include/boost/math/tools/roots.hpp index af46e2fcdd..040daee811 100644 --- a/include/boost/math/tools/roots.hpp +++ b/include/boost/math/tools/roots.hpp @@ -275,7 +275,13 @@ T newton_raphson_iterate(F f, T guess, T min, T max, int digits, std::uintmax_t& if (fabs(delta * 2) > fabs(delta2)) { // Last two steps haven't converged. - delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2; + T shift = (delta > 0) ? (result - min) / 2 : (result - max) / 2; + if ((result != 0) && (fabs(shift) > fabs(result))) + { + delta = sign(delta) * fabs(result); // protect against huge jumps! + } + else + delta = shift; // reset delta1/2 so we don't take this branch next time round: delta1 = 3 * delta; delta2 = 3 * delta; diff --git a/include/boost/math/tools/toms748_solve.hpp b/include/boost/math/tools/toms748_solve.hpp index 2de491a951..ea93713224 100644 --- a/include/boost/math/tools/toms748_solve.hpp +++ b/include/boost/math/tools/toms748_solve.hpp @@ -57,7 +57,7 @@ struct equal_floor bool operator()(const T& a, const T& b) { BOOST_MATH_STD_USING - return floor(a) == floor(b); + return (floor(a) == floor(b)) || (fabs((b-a)/b) < boost::math::tools::epsilon() * 2); } }; @@ -68,7 +68,7 @@ struct equal_ceil bool operator()(const T& a, const T& b) { BOOST_MATH_STD_USING - return ceil(a) == ceil(b); + return (ceil(a) == ceil(b)) || (fabs((b - a) / b) < boost::math::tools::epsilon() * 2); } }; @@ -79,7 +79,7 @@ struct equal_nearest_integer bool operator()(const T& a, const T& b) { BOOST_MATH_STD_USING - return floor(a + 0.5f) == floor(b + 0.5f); + return (floor(a + 0.5f) == floor(b + 0.5f)) || (fabs((b - a) / b) < boost::math::tools::epsilon() * 2); } }; diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 532b875212..07a7247660 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -170,6 +170,7 @@ test-suite special_fun : [ run git_issue_826.cpp ../../test/build//boost_unit_test_framework ] [ run git_issue_961.cpp ] [ run git_issue_1006.cpp ] + [ run git_issue_184.cpp ] [ run special_functions_test.cpp ../../test/build//boost_unit_test_framework ] [ run test_airy.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] [ run test_bessel_j.cpp test_instances//test_instances pch_light ../../test/build//boost_unit_test_framework ] diff --git a/test/git_issue_184.cpp b/test/git_issue_184.cpp new file mode 100644 index 0000000000..d82790588c --- /dev/null +++ b/test/git_issue_184.cpp @@ -0,0 +1,35 @@ +// Copyright John Maddock 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +#include +#include + +template +void test() +{ + boost::math::skew_normal_distribution skew(573.39724735636185, 77.0, 4.0); + const T x = boost::math::quantile(skew, 0.00285612015554148); + const T y = boost::math::quantile(skew, 0.00285612015554149); + const T z = boost::math::quantile(skew, 0.00285612015554150); + + BOOST_ASSERT(boost::math::isnormal(x)); + BOOST_ASSERT(boost::math::isnormal(y)); + BOOST_ASSERT(boost::math::isnormal(z)); + + BOOST_ASSERT(x <= y); + BOOST_ASSERT(y <= z); +} + + +int main() +{ + test(); + test(); + test(); + + return 0; +} From 04c2c248dfc5e35eeb7638152d5bd7c2985feef2 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Tue, 6 Feb 2024 09:54:49 +0000 Subject: [PATCH 2/2] BOOST_ASSERT->BOOST_MATH_ASSERT. --- test/git_issue_184.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/git_issue_184.cpp b/test/git_issue_184.cpp index d82790588c..cb925c912d 100644 --- a/test/git_issue_184.cpp +++ b/test/git_issue_184.cpp @@ -16,12 +16,12 @@ void test() const T y = boost::math::quantile(skew, 0.00285612015554149); const T z = boost::math::quantile(skew, 0.00285612015554150); - BOOST_ASSERT(boost::math::isnormal(x)); - BOOST_ASSERT(boost::math::isnormal(y)); - BOOST_ASSERT(boost::math::isnormal(z)); + BOOST_MATH_ASSERT(boost::math::isnormal(x)); + BOOST_MATH_ASSERT(boost::math::isnormal(y)); + BOOST_MATH_ASSERT(boost::math::isnormal(z)); - BOOST_ASSERT(x <= y); - BOOST_ASSERT(y <= z); + BOOST_MATH_ASSERT(x <= y); + BOOST_MATH_ASSERT(y <= z); }