diff options
author | DongHun Kwak <dh0128.kwak@samsung.com> | 2016-10-06 10:41:18 +0900 |
---|---|---|
committer | DongHun Kwak <dh0128.kwak@samsung.com> | 2016-10-06 10:43:11 +0900 |
commit | f763a99a501650eff2c60288aa6f10ef916d769e (patch) | |
tree | 02af7e13f9a38c888ebf340fe764cbe7dae99da9 /boost/math/special_functions/gamma.hpp | |
parent | 5cde13f21d36c7224b0e13d11c4b49379ae5210d (diff) | |
download | boost-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.hpp | 101 |
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{} |