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

More and better handling of ibeta small args.

Browse Source
This commit is contained in:
jzmaddock
2022-11-27 19:14:14 +00:00
parent 0dc6a70caa
commit 20fc2addea
4 changed files with 52 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
include/boost/math/special_functions/beta.hpp
View File

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

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

@@ -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&)
@@ -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)
{
@@ -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);
@@ -814,7 +822,20 @@ 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;
@@ -822,7 +843,20 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
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;

2
include/boost/math/special_functions/gamma.hpp
View File

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