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

Improve accuracy of tgamma_ratio when one argument is very small, thanks to ideas from Rocco Romeo.

Browse Source 
[SVN r83250]
This commit is contained in:
John Maddock
2013-03-02 18:59:50 +00:00
parent cd55f94984
commit 51dd944d85
2 changed files with 76 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -1207,6 +1207,68 @@ T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol)
return tgamma_delta_ratio_imp_lanczos(z, delta, pol, lanczos_type());
}
template <class T, class Policy>
T tgamma_ratio_imp(T x, T y, const Policy& pol)
{
BOOST_MATH_STD_USING
if(x <= 0)
policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got a=%1%).", x, pol);
if(y <= 0)
policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got b=%1%).", y, pol);
if((x < max_factorial<T>::value) && (y < max_factorial<T>::value))
{
// Rather than subtracting values, lets just call the gamma functions directly:
return boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
}
T prefix = 1;
if(x < 1)
{
if(y < 2 * max_factorial<T>::value)
{
// We need to sidestep on x as well, otherwise we'll underflow
// before we get to factor in the prefix term:
prefix /= x;
x += 1;
while(y >= max_factorial<T>::value)
{
y -= 1;
prefix /= y;
}
return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
}
//
// result is almost certainly going to underflow to zero, try logs just in case:
//
return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
}
if(y < 1)
{
if(x < 2 * max_factorial<T>::value)
{
// We need to sidestep on y as well, otherwise we'll overflow
// before we get to factor in the prefix term:
prefix *= y;
y += 1;
while(x >= max_factorial<T>::value)
{
x -= 1;
prefix *= x;
}
return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
}
//
// Result will almost certainly overflow, try logs just in case:
//
return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
}
//
// Regular case, x and y both large and similar in magnitude:
//
return boost::math::tgamma_delta_ratio(x, y - x, pol);
}
template <class T, class Policy>
T gamma_p_derivative_imp(T a, T x, const Policy& pol)
{
@@ -1609,7 +1671,7 @@ inline typename tools::promote_args<T1, T2>::type
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(static_cast<value_type>(b) - static_cast<value_type>(a)), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(b), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
}
template <class T1, class T2>
inline typename tools::promote_args<T1, T2>::type

13
test/test_tgamma_ratio.hpp
View File

@@ -116,5 +116,18 @@ void test_tgamma_ratio(T, const char* name)
do_test_tgamma_ratio<T>(tgamma_ratio_data, name, "tgamma ratios");
//
// A few special spot tests:
//
BOOST_MATH_STD_USING
T tol = boost::math::tools::epsilon<T>() * 20;
if(std::numeric_limits<T>::max_exponent > 200)
{
BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_ratio(ldexp(T(1), -500), T(180.25)), 8.0113754557649679470816892372669519037339812035512e-178L, tol);
BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_ratio(ldexp(T(1), -525), T(192.25)), 1.5966560279353205461166489184101261541784867035063e-197L, tol);
BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_ratio(T(182.25), ldexp(T(1), -500)), 4.077990437521002194346763299159975185747917450788e+181L, tol);
BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_ratio(T(193.25), ldexp(T(1), -525)), 1.2040790040958522422697601672703926839178050326148e+199L, tol);
BOOST_CHECK_CLOSE_FRACTION(boost::math::tgamma_ratio(T(193.25), T(194.75)), 0.00037151765099653237632823607820104961270831942138159L, tol);
}
}