Skip to content

More and better handling of ibeta small args. #884

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Nov 28, 2022
Merged

More and better handling of ibeta small args. #884

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions include/boost/math/special_functions/beta.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -570,6 +570,9 @@ T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_deriva
T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b));

if (!(boost::math::isfinite)(result))
result = 0;

T l1 = log(cgh / bgh) * (b - 0.5f);
T l2 = log(x * cgh / agh) * a;
//
Expand Down
40 changes: 37 additions & 3 deletions include/boost/math/special_functions/detail/ibeta_inverse.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -638,6 +638,12 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
#endif
{
bet = boost::math::beta(a, b, pol);

typedef typename Policy::overflow_error_type overflow_type;

BOOST_IF_CONSTEXPR(overflow_type::value != boost::math::policies::throw_on_error)
if(bet > tools::max_value<T>())
bet = tools::max_value<T>();
}
#ifndef BOOST_NO_EXCEPTIONS
catch (const std::overflow_error&)
Expand Down Expand Up @@ -686,7 +692,7 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
invert = !invert;
xs = 1 - xs;
}
if (a < tools::min_value<T>())
if ((a < tools::min_value<T>()) && (b > tools::min_value<T>()))
{
if (py)
{
Expand Down Expand Up @@ -715,6 +721,8 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
if (overflow || !(boost::math::isfinite)(bet))
{
xg = exp((boost::math::lgamma(a + 1, pol) + boost::math::lgamma(b, pol) - boost::math::lgamma(a + b, pol) + log(p)) / a);
if (xg > 2 / tools::epsilon<T>())
xg = 2 / tools::epsilon<T>();
}
else
xg = pow(a * p * bet, 1/a);
Expand Down Expand Up @@ -814,15 +822,41 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
}
if(pow(p, 1/a) < 0.5)
{
x = pow(p * a * boost::math::beta(a, b, pol), 1 / a);
#ifndef BOOST_NO_EXCEPTIONS
try
{
#endif
x = pow(p * a * boost::math::beta(a, b, pol), 1 / a);
if ((x > 1) || !(boost::math::isfinite)(x))
x = 1;
#ifndef BOOST_NO_EXCEPTIONS
}
catch (const std::overflow_error&)
{
x = 1;
}
#endif
if(x == 0)
x = boost::math::tools::min_value<T>();
y = 1 - x;
}
else /*if(pow(q, 1/b) < 0.1)*/
{
// model a distorted quarter circle:
y = pow(1 - pow(p, b * boost::math::beta(a, b, pol)), 1/b);
#ifndef BOOST_NO_EXCEPTIONS
try
{
#endif
y = pow(1 - pow(p, b * boost::math::beta(a, b, pol)), 1/b);
if ((y > 1) || !(boost::math::isfinite)(y))
y = 1;
#ifndef BOOST_NO_EXCEPTIONS
}
catch (const std::overflow_error&)
{
y = 1;
}
#endif
if(y == 0)
y = boost::math::tools::min_value<T>();
x = 1 - y;
Expand Down
2 changes: 1 addition & 1 deletion include/boost/math/special_functions/gamma.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -882,7 +882,7 @@ T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l)
//
// TODO: is this still required? Lanczos approx should be better now?
//
if(z <= tools::log_min_value<T>())
if((z <= tools::log_min_value<T>()) || (a < 1 / tools::max_value<T>()))
{
// Oh dear, have to use logs, should be free of cancellation errors though:
return exp(a * log(z) - z - lgamma_imp(a, pol, l));
Expand Down
11 changes: 11 additions & 0 deletions test/test_ibeta_inv.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -305,5 +305,16 @@ void test_spots(T)
BOOST_MATH_CHECK_THROW(::boost::math::ibeta_inv(static_cast<T>(2.125), -n, static_cast<T>(0.125)), std::domain_error);
BOOST_MATH_CHECK_THROW(::boost::math::ibeta_inv(static_cast<T>(2.125), static_cast<T>(1.125), -n), std::domain_error);
}
if (std::numeric_limits<T>::has_denorm)
{
T m = std::numeric_limits<T>::denorm_min();
T small = 2 * (std::numeric_limits<T>::min)();
BOOST_CHECK((boost::math::isfinite)(boost::math::ibeta_inv(m, static_cast<T>(2.125), static_cast<T>(0.125))));
BOOST_CHECK((boost::math::isfinite)(boost::math::ibeta_inv(m, m, static_cast<T>(0.125))));
BOOST_CHECK_LT(boost::math::ibeta_inv(m, static_cast<T>(12.125), static_cast<T>(0.125)), small);
BOOST_CHECK((boost::math::isfinite)(boost::math::ibeta_inv(static_cast<T>(2.125), m, static_cast<T>(0.125))));
BOOST_CHECK((boost::math::isfinite)(boost::math::ibeta_inv(static_cast<T>(12.125), m, static_cast<T>(0.125))));
BOOST_CHECK((boost::math::isfinite)(boost::math::ibeta_inv(m, m, static_cast<T>(0.125))));
}
}