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

Fix ibeta_inv for very small p. (#962)

Browse Source 
* Fix ibeta_inv for very small p.
Change assert's in temme_method_1_ibeta_inverse to corrections when guess goes out of range.
Change handling of non-convergence in second_order_root_finder to use bracketing when the end points are many orders of magnitude apart.
Fixes: https://github.com/boostorg/math/issues/961.

* Add missing copyright.
[CI SKIP]
This commit is contained in:
jzmaddock
2023-03-05 13:18:27 +00:00
committed by GitHub
parent 82ccb85907
commit bf3bc2e6c2
4 changed files with 110 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

19
include/boost/math/special_functions/detail/ibeta_inverse.hpp
View File

@@ -118,9 +118,17 @@ T temme_method_1_ibeta_inverse(T a, T b, T z, const Policy& pol)
x = 0.5;
else
x = (1 + eta * sqrt((1 + c) / eta_2)) / 2;
BOOST_MATH_ASSERT(x >= 0);
BOOST_MATH_ASSERT(x <= 1);
//
// These are post-conditions of the method, but the addition above
// may result in us being out by 1ulp either side of the boundary,
// so just check that we're in bounds and adjust as needed.
// See https://github.com/boostorg/math/issues/961
//
if (x < 0)
x = 0;
else if (x > 1)
x = 1;
BOOST_MATH_ASSERT(eta * (x - 0.5) >= 0);
#ifdef BOOST_INSTRUMENT
std::cout << "Estimating x with Temme method 1: " << x << std::endl;
@@ -901,6 +909,8 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
if(x < lower)
x = lower;
}
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
std::uintmax_t max_iter_used = 0;
//
// Figure out how many digits to iterate towards:
//
@@ -923,10 +933,9 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
// Now iterate, we can use either p or q as the target here
// depending on which is smaller:
//
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
x = boost::math::tools::halley_iterate(
boost::math::detail::ibeta_roots<T, Policy>(a, b, (p < q ? p : q), (p < q ? false : true)), x, lower, upper, digits, max_iter);
policies::check_root_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%)", max_iter, pol);
policies::check_root_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%)", max_iter + max_iter_used, pol);
//
// We don't really want these asserts here, but they are useful for sanity
// checking that we have the limits right, uncomment if you suspect bugs *only*.

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

@@ -586,8 +586,8 @@ namespace detail {
// we can jump way out of bounds if we're not careful.
// See https://svn.boost.org/trac/boost/ticket/8314.
delta = f0 / f1;
if (fabs(delta) > 2 * fabs(guess))
delta = (delta < 0 ? -1 : 1) * 2 * fabs(guess);
if (fabs(delta) > 2 * fabs(result))
delta = (delta < 0 ? -1 : 1) * 2 * fabs(result);
}
}
else
@@ -600,9 +600,19 @@ namespace detail {
if ((convergence > 0.8) && (convergence < 2))
{
// last two steps haven't converged.
delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
if ((result != 0) && (fabs(delta) > result))
delta = sign(delta) * fabs(result) * 0.9f; // protect against huge jumps!
if (fabs(min) < 1 ? fabs(1000 * min) < fabs(max) : fabs(max / min) > 1000)
{
if(delta > 0)
delta = bracket_root_towards_min(f, result, f0, min, max, count);
else
delta = bracket_root_towards_max(f, result, f0, min, max, count);
}
else
{
delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
if ((result != 0) && (fabs(delta) > result))
delta = sign(delta) * fabs(result) * 0.9f; // protect against huge jumps!
}
// reset delta2 so that this branch will *not* be taken on the
// next iteration:
delta2 = delta * 3;
@@ -641,6 +651,8 @@ namespace detail {
result = guess - delta;
if (result <= min)
result = float_next(min);
if (result >= max)
result = float_prior(max);
guess = min;
continue;
}
@@ -669,6 +681,8 @@ namespace detail {
result = guess - delta;
if (result >= max)
result = float_prior(max);
if (result <= min)
result = float_next(min);
guess = min;
continue;
}