[Boost-commit] svn:boost r83250 - in trunk: boost/math/special_functions libs/math/test

Subject: [Boost-commit] svn:boost r83250 - in trunk: boost/math/special_functions libs/math/test
From: john_at_[hidden]
Date: 2013-03-02 13:59:50

Author: johnmaddock
Date: 2013-03-02 13:59:50 EST (Sat, 02 Mar 2013)
New Revision: 83250
URL: http://svn.boost.org/trac/boost/changeset/83250

Log:
Improve accuracy of tgamma_ratio when one argument is very small, thanks to ideas from Rocco Romeo.
Text files modified: 
   trunk/boost/math/special_functions/gamma.hpp |    64 +++++++++++++++++++++++++++++++++++++++ 
   trunk/libs/math/test/test_tgamma_ratio.hpp   |    13 ++++++++                                
   2 files changed, 76 insertions(+), 1 deletions(-)

Modified: trunk/boost/math/special_functions/gamma.hpp
==============================================================================
--- trunk/boost/math/special_functions/gamma.hpp	(original)
+++ trunk/boost/math/special_functions/gamma.hpp	2013-03-02 13:59:50 EST (Sat, 02 Mar 2013)
@@ -1208,6 +1208,68 @@
 }
 
 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 @@
       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 

Modified: trunk/libs/math/test/test_tgamma_ratio.hpp
==============================================================================
--- trunk/libs/math/test/test_tgamma_ratio.hpp	(original)
+++ trunk/libs/math/test/test_tgamma_ratio.hpp	2013-03-02 13:59:50 EST (Sat, 02 Mar 2013)
@@ -116,5 +116,18 @@
 
    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);
+   }
 }