diff options
Diffstat (limited to 'boost/math/special_functions/digamma.hpp')
| -rw-r--r-- | boost/math/special_functions/digamma.hpp | 176 |
1 files changed, 161 insertions, 15 deletions
diff --git a/boost/math/special_functions/digamma.hpp b/boost/math/special_functions/digamma.hpp index 785cd75c5e..718eaf9278 100644 --- a/boost/math/special_functions/digamma.hpp +++ b/boost/math/special_functions/digamma.hpp @@ -12,6 +12,7 @@ #include <boost/math/special_functions/math_fwd.hpp> #include <boost/math/tools/rational.hpp> +#include <boost/math/tools/series.hpp> #include <boost/math/tools/promotion.hpp> #include <boost/math/policies/error_handling.hpp> #include <boost/math/constants/constants.hpp> @@ -27,7 +28,8 @@ namespace detail{ // inline unsigned digamma_large_lim(const mpl::int_<0>*) { return 20; } - +inline unsigned digamma_large_lim(const mpl::int_<113>*) +{ return 20; } inline unsigned digamma_large_lim(const void*) { return 10; } // @@ -41,7 +43,7 @@ inline unsigned digamma_large_lim(const void*) // This first one gives 34-digit precision for x >= 20: // template <class T> -inline T digamma_imp_large(T x, const mpl::int_<0>*) +inline T digamma_imp_large(T x, const mpl::int_<113>*) { BOOST_MATH_STD_USING // ADL of std functions. static const T P[] = { @@ -141,12 +143,47 @@ inline T digamma_imp_large(T x, const mpl::int_<24>*) return result; } // +// Fully generic asymptotic expansion in terms of Bernoulli numbers, see: +// http://functions.wolfram.com/06.14.06.0012.01 +// +template <class T> +struct digamma_series_func +{ +private: + int k; + T xx; + T term; +public: + digamma_series_func(T x) : k(1), xx(x * x), term(1 / (x * x)) {} + T operator()() + { + T result = term * boost::math::bernoulli_b2n<T>(k) / (2 * k); + term /= xx; + ++k; + return result; + } + typedef T result_type; +}; + +template <class T, class Policy> +inline T digamma_imp_large(T x, const Policy& pol, const mpl::int_<0>*) +{ + BOOST_MATH_STD_USING + digamma_series_func<T> s(x); + T result = log(x) - 1 / (2 * x); + boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); + result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, -result); + result = -result; + policies::check_series_iterations<T>("boost::math::digamma<%1%>(%1%)", max_iter, pol); + return result; +} +// // Now follow rational approximations over the range [1,2]. // // 35-digit precision: // template <class T> -T digamma_imp_1_2(T x, const mpl::int_<0>*) +T digamma_imp_1_2(T x, const mpl::int_<113>*) { // // Now the approximation, we use the form: @@ -325,16 +362,16 @@ inline T digamma_imp_1_2(T x, const mpl::int_<24>*) static const T root = 1532632.0f / 1048576; static const T root_minor = static_cast<T>(0.3700660185912626595423257213284682051735604e-6L); static const T P[] = { - 0.25479851023250261e0, - -0.44981331915268368e0, - -0.43916936919946835e0, - -0.61041765350579073e-1 + 0.25479851023250261e0f, + -0.44981331915268368e0f, + -0.43916936919946835e0f, + -0.61041765350579073e-1f }; static const T Q[] = { 0.1e1, - 0.15890202430554952e1, - 0.65341249856146947e0, - 0.63851690523355715e-1 + 0.15890202430554952e1f, + 0.65341249856146947e0f, + 0.63851690523355715e-1f }; T g = x - root; g -= root_minor; @@ -410,6 +447,98 @@ T digamma_imp(T x, const Tag* t, const Policy& pol) return result; } +template <class T, class Policy> +T digamma_imp(T x, const mpl::int_<0>* t, const Policy& pol) +{ + // + // This handles reflection of negative arguments, and all our + // error handling, then forwards to the T-specific approximation. + // + BOOST_MATH_STD_USING // ADL of std functions. + + T result = 0; + // + // Check for negative arguments and use reflection: + // + if(x <= -1) + { + // Reflect: + x = 1 - x; + // Argument reduction for tan: + T remainder = x - floor(x); + // Shift to negative if > 0.5: + if(remainder > 0.5) + { + remainder -= 1; + } + // + // check for evaluation at a negative pole: + // + if(remainder == 0) + { + return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, (1 - x), pol); + } + result = constants::pi<T>() / tan(constants::pi<T>() * remainder); + } + if(x == 0) + return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, x, pol); + // + // If we're above the lower-limit for the + // asymptotic expansion then use it, the + // limit is a linear interpolation with + // limit = 10 at 50 bit precision and + // limit = 250 at 1000 bit precision. + // + T lim = 10 + (tools::digits<T>() - 50) * 240 / 950; + T two_x = ldexp(x, 1); + if(x >= lim) + { + result += digamma_imp_large(x, pol, t); + } + else if(floor(x) == x) + { + // + // Special case for integer arguments, see + // http://functions.wolfram.com/06.14.03.0001.01 + // + result = -constants::euler<T, Policy>(); + T val = 1; + while(val < x) + { + result += 1 / val; + val += 1; + } + } + else if(floor(two_x) == two_x) + { + // + // Special case for half integer arguments, see: + // http://functions.wolfram.com/06.14.03.0007.01 + // + result = -2 * constants::ln_two<T, Policy>() - constants::euler<T, Policy>(); + int n = itrunc(x); + if(n) + { + for(int k = 1; k < n; ++k) + result += 1 / T(k); + for(int k = n; k <= 2 * n - 1; ++k) + result += 2 / T(k); + } + } + else + { + // + // Rescale so we can use the asymptotic expansion: + // + while(x < lim) + { + result -= 1 / x; + x += 1; + } + result += digamma_imp_large(x, pol, t); + } + return result; +} // // Initializer: ensure all our constants are initialized prior to the first call of main: // @@ -420,9 +549,15 @@ struct digamma_initializer { init() { + typedef typename policies::precision<T, Policy>::type precision_type; + do_init(mpl::bool_<precision_type::value && (precision_type::value <= 113)>()); + } + void do_init(const mpl::true_&) + { boost::math::digamma(T(1.5), Policy()); boost::math::digamma(T(500), Policy()); } + void do_init(const mpl::false_&){} void force_instantiate()const{} }; static const init initializer; @@ -439,7 +574,7 @@ const typename digamma_initializer<T, Policy>::init digamma_initializer<T, Polic template <class T, class Policy> inline typename tools::promote_args<T>::type - digamma(T x, const Policy& pol) + digamma(T x, const Policy&) { typedef typename tools::promote_args<T>::type result_type; typedef typename policies::evaluation<result_type, Policy>::type value_type; @@ -447,7 +582,7 @@ inline typename tools::promote_args<T>::type typedef typename mpl::if_< mpl::or_< mpl::less_equal<precision_type, mpl::int_<0> >, - mpl::greater<precision_type, mpl::int_<64> > + mpl::greater<precision_type, mpl::int_<114> > >, mpl::int_<0>, typename mpl::if_< @@ -456,17 +591,28 @@ inline typename tools::promote_args<T>::type typename mpl::if_< mpl::less<precision_type, mpl::int_<54> >, mpl::int_<53>, - mpl::int_<64> + typename mpl::if_< + mpl::less<precision_type, mpl::int_<65> >, + mpl::int_<64>, + mpl::int_<113> + >::type >::type >::type >::type tag_type; + typedef typename policies::normalise< + Policy, + policies::promote_float<false>, + policies::promote_double<false>, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + // Force initialization of constants: - detail::digamma_initializer<result_type, Policy>::force_instantiate(); + detail::digamma_initializer<value_type, forwarding_policy>::force_instantiate(); return policies::checked_narrowing_cast<result_type, Policy>(detail::digamma_imp( static_cast<value_type>(x), - static_cast<const tag_type*>(0), pol), "boost::math::digamma<%1%>(%1%)"); + static_cast<const tag_type*>(0), forwarding_policy()), "boost::math::digamma<%1%>(%1%)"); } template <class T> |