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

Subject: [Boost-commit] svn:boost r56455 - in trunk: boost/math/special_functions boost/math/tools libs/math/test
From: john_at_[hidden]
Date: 2009-09-28 13:06:40

Author: johnmaddock
Date: 2009-09-28 13:06:39 EDT (Mon, 28 Sep 2009)
New Revision: 56455
URL: http://svn.boost.org/trac/boost/changeset/56455

Log:
A few more minor performance tweaks for issue #3407.
Text files modified: 
   trunk/boost/math/special_functions/gamma.hpp |   218 +++++++++++++++++++++++---------------- 
   trunk/boost/math/tools/series.hpp            |    12 +-                                      
   trunk/libs/math/test/test_igamma_inv.cpp     |     4                                         
   3 files changed, 137 insertions(+), 97 deletions(-)

Modified: trunk/boost/math/special_functions/gamma.hpp
==============================================================================
--- trunk/boost/math/special_functions/gamma.hpp	(original)
+++ trunk/boost/math/special_functions/gamma.hpp	2009-09-28 13:06:39 EDT (Mon, 28 Sep 2009)
@@ -324,7 +324,7 @@
 };
 
 template <class T, class Policy>
-inline T lower_gamma_series(T a, T z, const Policy& pol)
+inline T lower_gamma_series(T a, T z, const Policy& pol, T init_value = 0)
 {
    // Multiply result by ((z^a) * (e^-z) / a) to get the full
    // lower incomplete integral. Then divide by tgamma(a)
@@ -332,12 +332,7 @@
    lower_incomplete_gamma_series<T> s(a, z);
    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
    int bits = policies::digits<T, Policy>();
-#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
-   T zero = 0;
-   T result = boost::math::tools::sum_series(s, bits, max_iter, zero);
-#else
-   T result = boost::math::tools::sum_series(s, bits, max_iter);
-#endif
+   T result = boost::math::tools::sum_series(s, bits, max_iter, init_value);
    policies::check_series_iterations("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
    return result;
 }
@@ -733,7 +728,7 @@
 // Upper gamma fraction for very small a:
 //
 template <class T, class Policy>
-inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0)
+inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0, bool invert = false, T* pderivative = 0)
 {
    BOOST_MATH_STD_USING  // ADL of std functions.
    //
@@ -747,27 +742,29 @@
    result -= p;
    result /= a;
    detail::small_gamma2_series<T> s(a, x);
-   boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
-#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
-   T zero = 0;
-   result -= (p + 1) * tools::sum_series(s, boost::math::policies::digits<T, Policy>(), max_iter, zero);
-#else
-   result -= (p + 1) * tools::sum_series(s, boost::math::policies::digits<T, Policy>(), max_iter);
-#endif
+   boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>() - 10;
+   p += 1;
+   if(pderivative)
+      *pderivative = p / (*pgam * exp(x));
+   T init_value = invert ? *pgam : 0;
+   result = -p * tools::sum_series(s, boost::math::policies::digits<T, Policy>(), max_iter, (init_value - result) / p);
    policies::check_series_iterations("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
+   if(invert)
+      result = -result;
    return result;
 }
 //
 // Upper gamma fraction for integer a:
 //
-template <class T>
-inline T finite_gamma_q(T a, T x)
+template <class T, class Policy>
+inline T finite_gamma_q(T a, T x, Policy const& pol, T* pderivative = 0)
 {
    //
    // Calculates normalised Q when a is an integer:
    //
    BOOST_MATH_STD_USING
-   T sum = exp(-x);
+   T e = exp(-x);
+   T sum = e;
    if(sum != 0)
    {
       T term = sum;
@@ -778,6 +775,10 @@
          sum += term;
       }
    }
+   if(pderivative)
+   {
+      *pderivative = e * pow(x, a) / boost::math::unchecked_factorial<T>(itrunc(T(a - 1), pol));
+   }
    return sum;
 }
 //
@@ -838,30 +839,32 @@
 
    BOOST_ASSERT((p_derivative == 0) || (normalised == true));
 
-   bool is_int = floor(a) == a;
-   bool is_half_int = (floor(2 * a) == 2 * a) && !is_int;
+   bool is_int, is_half_int;
    bool is_small_a = (a < 30) && (a <= x + 1);
+   if(is_small_a)
+   {
+      T fa = floor(a);
+      is_int = (fa == a);
+      is_half_int = is_int ? false : (fabs(fa - a) == 0.5f);
+   }
+   else
+   {
+      is_int = is_half_int = false;
+   }
+
+   int eval_method;
    
-   if(is_int && is_small_a && (x > 0.6))
+   if(is_int && (x > 0.6))
    {
       // calculate Q via finite sum:
       invert = !invert;
-      result = finite_gamma_q(a, x);
-      if(normalised == false)
-         result *= boost::math::tgamma(a, pol);
-      // TODO: calculate derivative inside sum:
-      if(p_derivative)
-         *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
+      eval_method = 0;
    }
-   else if(is_half_int && is_small_a && (x > 0.2))
+   else if(is_half_int && (x > 0.2))
    {
       // calculate Q via finite sum for half integer a:
       invert = !invert;
-      result = finite_half_gamma_q(a, x, p_derivative, pol);
-      if(normalised == false)
-         result *= boost::math::tgamma(a, pol);
-      if(p_derivative && (*p_derivative == 0))
-         *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
+      eval_method = 1;
    }
    else if(x < 0.5)
    {
@@ -870,23 +873,11 @@
       //
       if(-0.4 / log(x) < a)
       {
-         // Compute P:
-         result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
-         if(p_derivative)
-            *p_derivative = result;
-         if(result != 0)
-            result *= detail::lower_gamma_series(a, x, pol) / a;
+         eval_method = 2;
       }
       else
       {
-         // Compute Q:
-         invert = !invert;
-         T g;
-         result = tgamma_small_upper_part(a, x, pol, &g);
-         if(normalised)
-            result /= g;
-         if(p_derivative)
-            *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
+         eval_method = 3;
       }
    }
    else if(x < 1.1)
@@ -896,23 +887,11 @@
       //
       if(x * 0.75f < a)
       {
-         // Compute P:
-         result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
-         if(p_derivative)
-            *p_derivative = result;
-         if(result != 0)
-            result *= detail::lower_gamma_series(a, x, pol) / a;
+         eval_method = 2;
       }
       else
       {
-         // Compute Q:
-         invert = !invert;
-         T g;
-         result = tgamma_small_upper_part(a, x, pol, &g);
-         if(normalised)
-            result /= g;
-         if(p_derivative)
-            *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
+         eval_method = 3;
       }
    }
    else
@@ -950,6 +929,93 @@
       }
       if(use_temme)
       {
+         eval_method = 5;
+      }
+      else
+      {
+         //
+         // Regular case where the result will not be too close to 0.5.
+         //
+         // Changeover here occurs at P ~ Q ~ 0.5
+         // Note that series computation of P is about x2 faster than continued fraction
+         // calculation of Q, so try and use the CF only when really necessary, especially
+         // for small x.
+         //
+         if(x - (1 / (3 * x)) < a)
+         {
+            eval_method = 2;
+         }
+         else
+         {
+            eval_method = 4;
+            invert = !invert;
+         }
+      }
+   }
+
+   switch(eval_method)
+   {
+   case 0:
+      {
+         result = finite_gamma_q(a, x, pol, p_derivative);
+         if(normalised == false)
+            result *= boost::math::tgamma(a, pol);
+         break;
+      }
+   case 1:
+      {
+         result = finite_half_gamma_q(a, x, p_derivative, pol);
+         if(normalised == false)
+            result *= boost::math::tgamma(a, pol);
+         if(p_derivative && (*p_derivative == 0))
+            *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
+         break;
+      }
+   case 2:
+      {
+         // Compute P:
+         result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
+         if(p_derivative)
+            *p_derivative = result;
+         if(result != 0)
+         {
+            T init_value = 0;
+            if(invert)
+            {
+               init_value = -a * (normalised ? 1 : boost::math::tgamma(a, pol)) / result;
+            }
+            result *= detail::lower_gamma_series(a, x, pol, init_value) / a;
+            if(invert)
+            {
+               invert = false;
+               result = -result;
+            }
+         }
+         break;
+      }
+   case 3:
+      {
+         // Compute Q:
+         invert = !invert;
+         T g;
+         result = tgamma_small_upper_part(a, x, pol, &g, invert, p_derivative);
+         invert = false;
+         if(normalised)
+            result /= g;
+         break;
+      }
+   case 4:
+      {
+         // Compute Q:
+         result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
+         if(p_derivative)
+            *p_derivative = result;
+         if(result != 0)
+            result *= upper_gamma_fraction(a, x, policies::digits<T, Policy>());
+         break;
+      }
+   case 5:
+      {
          //
          // Use compile time dispatch to the appropriate
          // Temme asymptotic expansion.  This may be dead code
@@ -980,33 +1046,7 @@
             invert = !invert;
          if(p_derivative)
             *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
-      }
-      else
-      {
-         //
-         // Regular case where the result will not be too close to 0.5.
-         //
-         // Changeover here occurs at P ~ Q ~ 0.5
-         // Note that series computation of P is about x2 faster than continued fraction
-         // calculation of Q, so try and use the CF only when really necessary, especially
-         // for small x.
-         //
-         result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
-         if(p_derivative)
-            *p_derivative = result;
-         if(x - (1 / (3 * x)) < a)
-         {
-            // Compute P:
-            if(result != 0)
-               result *= detail::lower_gamma_series(a, x, pol) / a;
-         }
-         else
-         {
-            // Compute Q:
-            invert = !invert;
-            if(result != 0)
-               result *= upper_gamma_fraction(a, x, policies::digits<T, Policy>());
-         }
+         break;
       }
    }
 
@@ -1115,7 +1155,7 @@
          //
          if((z <= max_factorial<T>::value) && (z + delta <= max_factorial<T>::value))
          {
-            return unchecked_factorial<T>((unsigned)itrunc(z, pol) - 1) / unchecked_factorial<T>((unsigned)itrunc(T(z + delta)) - 1);
+            return unchecked_factorial<T>((unsigned)itrunc(z, pol) - 1) / unchecked_factorial<T>((unsigned)itrunc(T(z + delta), pol) - 1);
          }
       }
       if(fabs(delta) < 20)

Modified: trunk/boost/math/tools/series.hpp
==============================================================================
--- trunk/boost/math/tools/series.hpp	(original)
+++ trunk/boost/math/tools/series.hpp	2009-09-28 13:06:39 EDT (Mon, 28 Sep 2009)
@@ -20,7 +20,7 @@
 // Simple series summation come first:
 //
 template <class Functor>
-typename Functor::result_type sum_series(Functor& func, int bits)
+inline typename Functor::result_type sum_series(Functor& func, int bits)
 {
    BOOST_MATH_STD_USING
 
@@ -38,7 +38,7 @@
 }
 
 template <class Functor>
-typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms)
+inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms)
 {
    BOOST_MATH_STD_USING
 
@@ -62,7 +62,7 @@
 }
 
 template <class Functor, class U>
-typename Functor::result_type sum_series(Functor& func, int bits, U init_value)
+inline typename Functor::result_type sum_series(Functor& func, int bits, U init_value)
 {
    BOOST_MATH_STD_USING
 
@@ -81,7 +81,7 @@
 }
 
 template <class Functor, class U>
-typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms, U init_value)
+inline typename Functor::result_type sum_series(Functor& func, int bits, boost::uintmax_t& max_terms, U init_value)
 {
    BOOST_MATH_STD_USING
 
@@ -117,7 +117,7 @@
 // in any case the result is still much better than a naive summation.
 //
 template <class Functor>
-typename Functor::result_type kahan_sum_series(Functor& func, int bits)
+inline typename Functor::result_type kahan_sum_series(Functor& func, int bits)
 {
    BOOST_MATH_STD_USING
 
@@ -140,7 +140,7 @@
 }
 
 template <class Functor>
-typename Functor::result_type kahan_sum_series(Functor& func, int bits, boost::uintmax_t& max_terms)
+inline typename Functor::result_type kahan_sum_series(Functor& func, int bits, boost::uintmax_t& max_terms)
 {
    BOOST_MATH_STD_USING
 

Modified: trunk/libs/math/test/test_igamma_inv.cpp
==============================================================================
--- trunk/libs/math/test/test_igamma_inv.cpp	(original)
+++ trunk/libs/math/test/test_igamma_inv.cpp	2009-09-28 13:06:39 EDT (Mon, 28 Sep 2009)
@@ -326,7 +326,7 @@
          data,
          bind_func(funcp, 0, 1),
          extract_result(2));
-      handle_test_result(result, data[result.worst()], result.worst(), type_name, "other::gamma_p_inv", test_name);
+      print_test_result(result, data[result.worst()], result.worst(), type_name, "other::gamma_q");
       //
       // test gamma_q_inv(T, T) against data:
       //
@@ -335,7 +335,7 @@
          data,
          bind_func(funcp, 0, 1),
          extract_result(3));
-      handle_test_result(result, data[result.worst()], result.worst(), type_name, "other::gamma_q_inv", test_name);
+      print_test_result(result, data[result.worst()], result.worst(), type_name, "other::gamma_q");
    }
 #endif
 }