boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Merge pull request #1080 from boostorg/issue184_2024

Browse Source 
Reinstate root finding protection against huge jumps.
This commit is contained in:
jzmaddock
2024-02-06 12:34:39 +00:00
committed by GitHub
parent 75dcb3e3d5 04c2c248df
commit 03ea9c8d7d
18 changed files with 160 additions and 121 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
include/boost/math/distributions/binomial.hpp
View File

@@ -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<Policy>();
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
}

5
include/boost/math/distributions/chi_squared.hpp
View File

@@ -327,9 +327,8 @@ RealType chi_squared_distribution<RealType, Policy>::find_degrees_of_freedom(
RealType result = r.first + (r.second - r.first) / 2;
if(max_iter >= policies::get_max_root_iterations<Policy>())
{
policies::raise_evaluation_error<RealType>(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<RealType>(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;
}

19
include/boost/math/distributions/detail/generic_mode.hpp
View File

@@ -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<policy_type>())
{
return policies::raise_evaluation_error<value_type>(
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<value_type>(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<policy_type>())
{
return policies::raise_evaluation_error<value_type>(
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<value_type>(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;
}

5
include/boost/math/distributions/detail/generic_quantile.hpp
View File

@@ -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<forwarding_policy>())
{
return policies::raise_evaluation_error<value_type>(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<value_type>(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;
}

9
include/boost/math/distributions/detail/inv_discrete_quantile.hpp
View File

@@ -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<value_type, value_type> r = toms748_solve(f, a, b, fa, fb, tol, count, policy_type());
max_iter += count;
if (max_iter >= policies::get_max_root_iterations<policy_type>())
{
return policies::raise_evaluation_error<value_type>(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;
}

13
include/boost/math/distributions/find_scale.hpp
View File

@@ -80,9 +80,7 @@ namespace boost
if (result <= 0)
{ // If policy isn't to throw, return the scale <= 0.
policies::raise_evaluation_error<typename Dist::value_type>(function,
"Computed scale (%1%) is <= 0!" " Was the complement intended?",
result, Policy());
policies::raise_evaluation_error<typename Dist::value_type>(function, "Computed scale (%1%) is <= 0!" " Was the complement intended?", result, Policy()); // LCOV_EXCL_LINE
}
return result;
} // template <class Dist, class Policy> 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<typename Dist::value_type>(function,
"Computed scale (%1%) is <= 0!" " Was the complement intended?",
result, Policy());
policies::raise_evaluation_error<typename Dist::value_type>(function, "Computed scale (%1%) is <= 0!" " Was the complement intended?", result, Policy()); // LCOV_EXCL_LINE
}
return result;
} // template <class Dist, class Policy, class Real1, class Real2, class Real3> 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<typename Dist::value_type>(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<typename Dist::value_type>(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 <class Dist, class Real1, class Real2, class Real3> typename Dist::value_type find_scale

23
include/boost/math/distributions/inverse_gaussian.hpp
View File

@@ -388,11 +388,16 @@ inline RealType quantile(const inverse_gaussian_distribution<RealType, Policy>&
// 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<RealType, Policy>();// get digits from policy,
std::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); // and max iterations.
using boost::math::tools::newton_raphson_iterate;
result =
newton_raphson_iterate(inverse_gaussian_quantile_functor<RealType, Policy>(dist, p), guess, min, max, get_digits, m);
return result;
newton_raphson_iterate(inverse_gaussian_quantile_functor<RealType, Policy>(dist, p), guess, min, max, get_digits, max_iter);
if (max_iter >= policies::get_max_root_iterations<Policy>())
{
return policies::raise_evaluation_error<RealType>(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 <class RealType, class Policy>
@@ -459,11 +464,15 @@ inline RealType quantile(const complemented2_type<inverse_gaussian_distribution<
// int digits = std::numeric_limits<RealType>::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<RealType, Policy>();
std::uintmax_t m = policies::get_max_root_iterations<Policy>();
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
using boost::math::tools::newton_raphson_iterate;
result =
newton_raphson_iterate(inverse_gaussian_quantile_complement_functor<RealType, Policy>(c.dist, q), guess, min, max, get_digits, m);
return result;
result = newton_raphson_iterate(inverse_gaussian_quantile_complement_functor<RealType, Policy>(c.dist, q), guess, min, max, get_digits, max_iter);
if (max_iter >= policies::get_max_root_iterations<Policy>())
{
return policies::raise_evaluation_error<RealType>(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 <class RealType, class Policy>

24
include/boost/math/distributions/kolmogorov_smirnov.hpp
View File

@@ -379,10 +379,16 @@ inline RealType quantile(const kolmogorov_smirnov_distribution<RealType, Policy>
RealType k = detail::kolmogorov_smirnov_quantile_guess(p) / sqrt(n);
const int get_digits = policies::digits<RealType, Policy>();// get digits from policy,
std::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); // and max iterations.
return tools::newton_raphson_iterate(detail::kolmogorov_smirnov_quantile_functor<RealType, Policy>(dist, p),
k, RealType(0), RealType(1), get_digits, m);
RealType result = tools::newton_raphson_iterate(detail::kolmogorov_smirnov_quantile_functor<RealType, Policy>(dist, p),
k, RealType(0), RealType(1), get_digits, max_iter);
if (max_iter >= policies::get_max_root_iterations<Policy>())
{
return policies::raise_evaluation_error<RealType>(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 <class RealType, class Policy>
@@ -403,11 +409,17 @@ inline RealType quantile(const complemented2_type<kolmogorov_smirnov_distributio
RealType k = detail::kolmogorov_smirnov_quantile_guess(RealType(1-p)) / sqrt(n);
const int get_digits = policies::digits<RealType, Policy>();// get digits from policy,
std::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); // and max iterations.
return tools::newton_raphson_iterate(
RealType result = tools::newton_raphson_iterate(
detail::kolmogorov_smirnov_complementary_quantile_functor<RealType, Policy>(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<Policy>())
{
return policies::raise_evaluation_error<RealType>(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 <class RealType, class Policy>

34
include/boost/math/distributions/non_central_beta.hpp
View File

@@ -109,9 +109,7 @@ namespace boost
}
if(static_cast<std::uintmax_t>(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<std::uintmax_t>(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<std::uintmax_t>(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<Policy>())
{
// LCOV_EXCL_START
return policies::raise_evaluation_error<RealType>(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<RealType, forwarding_policy>(
result,
function), Policy());
// LCOV_EXCL_STOP
}
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
result,
@@ -583,9 +579,7 @@ namespace boost
}
if(static_cast<std::uintmax_t>(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<RealType>(
function,
"This function is not yet implemented, the only sensible result is %1%.",
std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity?
return policies::raise_evaluation_error<RealType>(function, "This function is not yet implemented, the only sensible result is %1%.", // LCOV_EXCL_LINE
std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity? LCOV_EXCL_LINE
}
template <class RealType, class Policy>
@@ -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<RealType>(
function,
"This function is not yet implemented, the only sensible result is %1%.",
std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity?
return policies::raise_evaluation_error<RealType>(function, "This function is not yet implemented, the only sensible result is %1%.", // LCOV_EXCL_LINE
std::numeric_limits<RealType>::quiet_NaN(), Policy()); // infinity? LCOV_EXCL_LINE
} // kurtosis_excess
template <class RealType, class Policy>

34
include/boost/math/distributions/non_central_chi_squared.hpp
View File

@@ -101,9 +101,7 @@ namespace boost
}
//Error check:
if(static_cast<std::uintmax_t>(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<std::uintmax_t>(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<std::uintmax_t>(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<std::uintmax_t>(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<RealType>(function,
"Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%",
RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
return policies::raise_evaluation_error<RealType>(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<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE
}
degrees_of_freedom_finder<RealType, Policy> f(lam, x, p < q ? p : q, p < q ? false : true);
tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
@@ -600,8 +591,8 @@ namespace boost
RealType result = ir.first + (ir.second - ir.first) / 2;
if(max_iter >= policies::get_max_root_iterations<Policy>())
{
return policies::raise_evaluation_error<RealType>(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<RealType>(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<RealType>(function,
"Can't find non centrality parameter when the probability is 0 or 1, only possible answer is %1%",
RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
return policies::raise_evaluation_error<RealType>(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<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE
}
non_centrality_finder<RealType, Policy> f(v, x, p < q ? p : q, p < q ? false : true);
tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
@@ -656,8 +646,8 @@ namespace boost
RealType result = ir.first + (ir.second - ir.first) / 2;
if(max_iter >= policies::get_max_root_iterations<Policy>())
{
return policies::raise_evaluation_error<RealType>(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<RealType>(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;
}

34
include/boost/math/distributions/non_central_t.hpp
View File

@@ -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<RealType>(function,
"Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%",
RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
return policies::raise_evaluation_error<RealType>(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<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE
}
t_degrees_of_freedom_finder<RealType, Policy> f(delta, x, p < q ? p : q, p < q ? false : true);
tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
@@ -661,8 +652,8 @@ namespace boost
RealType result = ir.first + (ir.second - ir.first) / 2;
if(max_iter >= policies::get_max_root_iterations<Policy>())
{
return policies::raise_evaluation_error<RealType>(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<RealType>(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<RealType>(function,
"Can't find non-centrality parameter when the probability is 0 or 1, only possible answer is %1%",
RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
return policies::raise_evaluation_error<RealType>(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<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE
}
t_non_centrality_finder<RealType, Policy> f(v, x, p < q ? p : q, p < q ? false : true);
tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
@@ -719,8 +709,8 @@ namespace boost
RealType result = ir.first + (ir.second - ir.first) / 2;
if(max_iter >= policies::get_max_root_iterations<Policy>())
{
return policies::raise_evaluation_error<RealType>(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<RealType>(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;
}

18
include/boost/math/distributions/skew_normal.hpp
View File

@@ -561,12 +561,17 @@ namespace boost{ namespace math{
}
const int get_digits = policies::digits<RealType, Policy>();// get digits from policy,
std::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); // and max iterations.
skew_normal_distribution<RealType, Policy> helper(0, 1, shape);
result = tools::newton_raphson_iterate(detail::skew_normal_mode_functor<RealType, Policy>(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<Policy>())
{
return policies::raise_evaluation_error<RealType>(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<RealType, Policy>();// get digits from policy,
std::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); // and max iterations.
result = tools::newton_raphson_iterate(detail::skew_normal_quantile_functor<RealType, Policy>(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<Policy>())
{
return policies::raise_evaluation_error<RealType>(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

5
include/boost/math/distributions/students_t.hpp
View File

@@ -335,9 +335,8 @@ RealType students_t_distribution<RealType, Policy>::find_degrees_of_freedom(
RealType result = r.first + (r.second - r.first) / 2;
if(max_iter >= policies::get_max_root_iterations<Policy>())
{
return policies::raise_evaluation_error<RealType>(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<RealType>(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;
}

5
include/boost/math/special_functions/legendre.hpp
View File

@@ -206,6 +206,11 @@ std::vector<T> legendre_p_zeros_imp(int n, const Policy& pol)
lower_bound, upper_bound,
policies::digits<T, Policy>(),
number_of_iterations);
if (number_of_iterations >= policies::get_max_root_iterations<Policy>())
{
policies::raise_evaluation_error<T>("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);

8
include/boost/math/tools/roots.hpp
View File

@@ -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;

6
include/boost/math/tools/toms748_solve.hpp
View File

@@ -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<T>() * 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<T>() * 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<T>() * 2);
}
};

1
test/Jamfile.v2
View File

@@ -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 ]

35
test/git_issue_184.cpp Normal file
View File

@@ -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 <iostream>
#include <boost/math/distributions/skew_normal.hpp>
template <class T>
void test()
{
boost::math::skew_normal_distribution<T> 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_MATH_ASSERT(boost::math::isnormal(x));
BOOST_MATH_ASSERT(boost::math::isnormal(y));
BOOST_MATH_ASSERT(boost::math::isnormal(z));
BOOST_MATH_ASSERT(x <= y);
BOOST_MATH_ASSERT(y <= z);
}
int main()
{
test<float>();
test<double>();
test<long double>();
return 0;
}