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
}