summaryrefslogtreecommitdiff
path: root/boost/math/special_functions/gamma.hpp
diff options
context:
space:
mode:
authorDongHun Kwak <dh0128.kwak@samsung.com>2016-10-06 10:41:18 +0900
committerDongHun Kwak <dh0128.kwak@samsung.com>2016-10-06 10:43:11 +0900
commitf763a99a501650eff2c60288aa6f10ef916d769e (patch)
tree02af7e13f9a38c888ebf340fe764cbe7dae99da9 /boost/math/special_functions/gamma.hpp
parent5cde13f21d36c7224b0e13d11c4b49379ae5210d (diff)
downloadboost-f763a99a501650eff2c60288aa6f10ef916d769e.tar.gz
boost-f763a99a501650eff2c60288aa6f10ef916d769e.tar.bz2
boost-f763a99a501650eff2c60288aa6f10ef916d769e.zip
Imported Upstream version 1.62.0upstream/1.62.0
Change-Id: I9d4c1ddb7b7d8f0069217ecc582700f9fda6dd4c Signed-off-by: DongHun Kwak <dh0128.kwak@samsung.com>
Diffstat (limited to 'boost/math/special_functions/gamma.hpp')
-rw-r--r--boost/math/special_functions/gamma.hpp101
1 files changed, 71 insertions, 30 deletions
diff --git a/boost/math/special_functions/gamma.hpp b/boost/math/special_functions/gamma.hpp
index 3a3191a807..a2b30f8551 100644
--- a/boost/math/special_functions/gamma.hpp
+++ b/boost/math/special_functions/gamma.hpp
@@ -33,6 +33,7 @@
#include <boost/math/special_functions/detail/unchecked_factorial.hpp>
#include <boost/math/special_functions/detail/lgamma_small.hpp>
#include <boost/math/special_functions/bernoulli.hpp>
+#include <boost/math/special_functions/zeta.hpp>
#include <boost/type_traits/is_convertible.hpp>
#include <boost/assert.hpp>
#include <boost/mpl/greater.hpp>
@@ -495,6 +496,37 @@ T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&)
}
template <class T, class Policy>
+inline T log_gamma_near_1(const T& z, Policy const& pol)
+{
+ //
+ // This is for the multiprecision case where there is
+ // no lanczos support...
+ //
+ BOOST_MATH_STD_USING // ADL of std names
+
+ BOOST_ASSERT(fabs(z) < 1);
+
+ T result = -constants::euler<T>() * z;
+
+ T power_term = z * z;
+ T term;
+ unsigned j = 0;
+
+ do
+ {
+ term = boost::math::zeta<T>(j + 2, pol) * power_term / (j + 2);
+ if(j & 1)
+ result -= term;
+ else
+ result += term;
+ power_term *= z;
+ ++j;
+ } while(fabs(result) * tools::epsilon<T>() < fabs(term));
+
+ return result;
+}
+
+template <class T, class Policy>
T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign)
{
BOOST_MATH_STD_USING
@@ -529,7 +561,15 @@ T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sig
// Here we simply take the logarithm of tgamma(). This is somewhat
// inefficient, but simple. The rationale is that the argument here
// is relatively small and overflow is not expected to be likely.
- if (z > -tools::root_epsilon<T>())
+ if(fabs(z - 1) < 0.25)
+ {
+ return log_gamma_near_1(T(zz - 1), pol);
+ }
+ else if(fabs(z - 2) < 0.25)
+ {
+ return log_gamma_near_1(T(zz - 2), pol) + log(zz - 1);
+ }
+ else if (z > -tools::root_epsilon<T>())
{
// Reflection formula may fail if z is very close to zero, let the series
// expansion for tgamma close to zero do the work:
@@ -634,7 +674,7 @@ T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& l)
mpl::greater<precision_type, mpl::int_<113> >
>,
typename mpl::if_<
- is_same<Lanczos, lanczos::lanczos24m113>,
+ mpl::and_<is_same<Lanczos, lanczos::lanczos24m113>, mpl::greater<precision_type, mpl::int_<0> > >,
mpl::int_<113>,
mpl::int_<0>
>::type,
@@ -680,32 +720,16 @@ T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& l)
}
template <class T, class Policy>
-inline T tgammap1m1_imp(T dz, Policy const& pol,
- const ::boost::math::lanczos::undefined_lanczos& l)
+inline T tgammap1m1_imp(T z, Policy const& pol,
+ const ::boost::math::lanczos::undefined_lanczos&)
{
BOOST_MATH_STD_USING // ADL of std names
- //
- // There should be a better solution than this, but the
- // algebra isn't easy for the general case....
- // Start by subracting 1 from tgamma:
- //
- T result = gamma_imp(T(1 + dz), pol, l) - 1;
- BOOST_MATH_INSTRUMENT_CODE(result);
- //
- // Test the level of cancellation error observed: we loose one bit
- // for each power of 2 the result is less than 1. If we would get
- // more bits from our most precise lgamma rational approximation,
- // then use that instead:
- //
- BOOST_MATH_INSTRUMENT_CODE((dz > -0.5));
- BOOST_MATH_INSTRUMENT_CODE((dz < 2));
- BOOST_MATH_INSTRUMENT_CODE((ldexp(1.0, boost::math::policies::digits<T, Policy>()) * fabs(result) < 1e34));
- if((dz > -0.5) && (dz < 2) && (ldexp(1.0, boost::math::policies::digits<T, Policy>()) * fabs(result) < 1e34))
+
+ if(fabs(z) < 0.55)
{
- result = tgammap1m1_imp(dz, pol, boost::math::lanczos::lanczos24m113());
- BOOST_MATH_INSTRUMENT_CODE(result);
+ return boost::math::expm1(log_gamma_near_1(z, pol));
}
- return result;
+ return boost::math::expm1(boost::math::lgamma(1 + z, pol));
}
//
@@ -1396,16 +1420,26 @@ T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos&
}
T zgh = static_cast<T>(z + Lanczos::g() - constants::half<T>());
T result;
- if(fabs(delta) < 10)
+ if(z + delta == z)
{
- result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
+ if(fabs(delta) < 10)
+ result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
+ else
+ result = 1;
}
else
{
- result = pow(zgh / (zgh + delta), z - constants::half<T>());
+ if(fabs(delta) < 10)
+ {
+ result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
+ }
+ else
+ {
+ result = pow(zgh / (zgh + delta), z - constants::half<T>());
+ }
+ // Split the calculation up to avoid spurious overflow:
+ result *= Lanczos::lanczos_sum(z) / Lanczos::lanczos_sum(T(z + delta));
}
- // Split the calculation up to avoid spurious overflow:
- result *= Lanczos::lanczos_sum(z) / Lanczos::lanczos_sum(T(z + delta));
result *= pow(constants::e<T>() / (zgh + delta), delta);
return result;
}
@@ -1663,7 +1697,14 @@ struct igamma_initializer
template <int N>
static void do_init(const mpl::int_<N>&)
{
- boost::math::gamma_p(static_cast<T>(400), static_cast<T>(400), Policy());
+ // If std::numeric_limits<T>::digits is zero, we must not call
+ // our inituialization code here as the precision presumably
+ // varies at runtime, and will not have been set yet. Plus the
+ // code requiring initialization isn't called when digits == 0.
+ if(std::numeric_limits<T>::digits)
+ {
+ boost::math::gamma_p(static_cast<T>(400), static_cast<T>(400), Policy());
+ }
}
static void do_init(const mpl::int_<53>&){}
void force_instantiate()const{}