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

Merge pull request #875 from boostorg/issue873

Browse Source 
Correct root finders in extreme cases.
This commit is contained in:
jzmaddock
2022-11-15 09:00:39 +00:00
committed by GitHub
parent 03ac231b83 9a9bb9ddc8
commit 2cdabeb30e
2 changed files with 40 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -372,13 +372,19 @@ namespace detail {
T bracket_root_towards_max(F f, T guess, const T& f0, T& min, T& max, std::uintmax_t& count) noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
{
using std::fabs;
using std::ldexp;
using std::abs;
using std::frexp;
if(count < 2)
return guess - (max + min) / 2; // Not enough counts left to do anything!!
//
// Move guess towards max until we bracket the root, updating min and max as we go:
//
int e;
frexp(max / guess, &e);
e = abs(e);
T guess0 = guess;
T multiplier = 2;
T multiplier = e < 64 ? static_cast<T>(2) : static_cast<T>(ldexp(T(1), e / 32));
T f_current = f0;
if (fabs(min) < fabs(max))
{
@@ -392,7 +398,7 @@ namespace detail {
f_current = -f_current; // There must be a change of sign!
break;
}
multiplier *= 2;
multiplier *= e > 1024 ? 8 : 2;
unpack_0(f(guess), f_current);
}
}
@@ -411,7 +417,7 @@ namespace detail {
f_current = -f_current; // There must be a change of sign!
break;
}
multiplier *= 2;
multiplier *= e > 1024 ? 8 : 2;
unpack_0(f(guess), f_current);
}
}
@@ -429,13 +435,19 @@ namespace detail {
T bracket_root_towards_min(F f, T guess, const T& f0, T& min, T& max, std::uintmax_t& count) noexcept(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval<F>()(std::declval<T>())))
{
using std::fabs;
using std::ldexp;
using std::abs;
using std::frexp;
if (count < 2)
return guess - (max + min) / 2; // Not enough counts left to do anything!!
//
// Move guess towards min until we bracket the root, updating min and max as we go:
//
int e;
frexp(guess / min, &e);
e = abs(e);
T guess0 = guess;
T multiplier = 2;
T multiplier = e < 64 ? static_cast<T>(2) : static_cast<T>(ldexp(T(1), e / 32));
T f_current = f0;
if (fabs(min) < fabs(max))
@@ -450,7 +462,7 @@ namespace detail {
f_current = -f_current; // There must be a change of sign!
break;
}
multiplier *= 2;
multiplier *= e > 1024 ? 8 : 2;
unpack_0(f(guess), f_current);
}
}
@@ -469,7 +481,7 @@ namespace detail {
f_current = -f_current; // There must be a change of sign!
break;
}
multiplier *= 2;
multiplier *= e > 1024 ? 8 : 2;
unpack_0(f(guess), f_current);
}
}
@@ -615,6 +627,8 @@ namespace detail {
}
delta = bracket_root_towards_min(f, guess, f0, min, max, count);
result = guess - delta;
if (result <= min)
result = float_next(min);
guess = min;
continue;
}
@@ -641,6 +655,8 @@ namespace detail {
}
delta = bracket_root_towards_max(f, guess, f0, min, max, count);
result = guess - delta;
if (result >= max)
result = float_prior(max);
guess = min;
continue;
}

18
test/test_ibeta_inv.hpp
View File

@@ -270,5 +270,23 @@ void test_spots(T)
static_cast<T>(0.125),
static_cast<T>(0.125)),
static_cast<T>(0.99999994039535522460937500000000000000000000000L), tolerance);
//
// Bug cases, issue 873:
//
if ((std::numeric_limits<T>::max)() > static_cast<T>(1e50))
{
BOOST_CHECK_CLOSE(
::boost::math::ibeta_inv(
static_cast<T>(1e50L),
static_cast<T>(10),
static_cast<T>(1) / static_cast<T>(10)),
static_cast<T>(1), tolerance);
BOOST_CHECK_CLOSE(
::boost::math::ibetac_inv(
static_cast<T>(1e50L),
static_cast<T>(10),
static_cast<T>(1) / static_cast<T>(10)),
static_cast<T>(1), tolerance);
}
}