math/special_functions/digamma.hpp

0001 //  (C) Copyright John Maddock 2006.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_SF_DIGAMMA_HPP
0007 #define BOOST_MATH_SF_DIGAMMA_HPP
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #pragma warning(push)
0012 #pragma warning(disable:4702) // Unreachable code (release mode only warning)
0013 #endif
0014
0015 #include <boost/math/special_functions/math_fwd.hpp>
0016 #include <boost/math/tools/rational.hpp>
0017 #include <boost/math/tools/series.hpp>
0018 #include <boost/math/tools/promotion.hpp>
0019 #include <boost/math/policies/error_handling.hpp>
0020 #include <boost/math/constants/constants.hpp>
0021 #include <boost/math/tools/big_constant.hpp>
0022
0023 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
0024 //
0025 // This is the only way we can avoid
0026 // warning: non-standard suffix on floating constant [-Wpedantic]
0027 // when building with -Wall -pedantic.  Neither __extension__
0028 // nor #pragma diagnostic ignored work :(
0029 //
0030 #pragma GCC system_header
0031 #endif
0032
0033 namespace boost{
0034 namespace math{
0035 namespace detail{
0036 //
0037 // Begin by defining the smallest value for which it is safe to
0038 // use the asymptotic expansion for digamma:
0039 //
0040 inline unsigned digamma_large_lim(const std::integral_constant<int, 0>*)
0041 {  return 20;  }
0042 inline unsigned digamma_large_lim(const std::integral_constant<int, 113>*)
0043 {  return 20;  }
0044 inline unsigned digamma_large_lim(const void*)
0045 {  return 10;  }
0046 //
0047 // Implementations of the asymptotic expansion come next,
0048 // the coefficients of the series have been evaluated
0049 // in advance at high precision, and the series truncated
0050 // at the first term that's too small to effect the result.
0051 // Note that the series becomes divergent after a while
0052 // so truncation is very important.
0053 //
0054 // This first one gives 34-digit precision for x >= 20:
0055 //
0056 template <class T>
0057 inline T digamma_imp_large(T x, const std::integral_constant<int, 113>*)
0058 {
0059    BOOST_MATH_STD_USING // ADL of std functions.
0060    static const T P[] = {
0061       BOOST_MATH_BIG_CONSTANT(T, 113, 0.083333333333333333333333333333333333333333333333333),
0062       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0083333333333333333333333333333333333333333333333333),
0063       BOOST_MATH_BIG_CONSTANT(T, 113, 0.003968253968253968253968253968253968253968253968254),
0064       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0041666666666666666666666666666666666666666666666667),
0065       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0075757575757575757575757575757575757575757575757576),
0066       BOOST_MATH_BIG_CONSTANT(T, 113, -0.021092796092796092796092796092796092796092796092796),
0067       BOOST_MATH_BIG_CONSTANT(T, 113, 0.083333333333333333333333333333333333333333333333333),
0068       BOOST_MATH_BIG_CONSTANT(T, 113, -0.44325980392156862745098039215686274509803921568627),
0069       BOOST_MATH_BIG_CONSTANT(T, 113, 3.0539543302701197438039543302701197438039543302701),
0070       BOOST_MATH_BIG_CONSTANT(T, 113, -26.456212121212121212121212121212121212121212121212),
0071       BOOST_MATH_BIG_CONSTANT(T, 113, 281.4601449275362318840579710144927536231884057971),
0072       BOOST_MATH_BIG_CONSTANT(T, 113, -3607.510546398046398046398046398046398046398046398),
0073       BOOST_MATH_BIG_CONSTANT(T, 113, 54827.583333333333333333333333333333333333333333333),
0074       BOOST_MATH_BIG_CONSTANT(T, 113, -974936.82385057471264367816091954022988505747126437),
0075       BOOST_MATH_BIG_CONSTANT(T, 113, 20052695.796688078946143462272494530559046688078946),
0076       BOOST_MATH_BIG_CONSTANT(T, 113, -472384867.72162990196078431372549019607843137254902),
0077       BOOST_MATH_BIG_CONSTANT(T, 113, 12635724795.916666666666666666666666666666666666667)
0078    };
0079    x -= 1;
0080    T result = log(x);
0081    result += 1 / (2 * x);
0082    T z = 1 / (x*x);
0083    result -= z * tools::evaluate_polynomial(P, z);
0084    return result;
0085 }
0086 //
0087 // 19-digit precision for x >= 10:
0088 //
0089 template <class T>
0090 inline T digamma_imp_large(T x, const std::integral_constant<int, 64>*)
0091 {
0092    BOOST_MATH_STD_USING // ADL of std functions.
0093    static const T P[] = {
0094       BOOST_MATH_BIG_CONSTANT(T, 64, 0.083333333333333333333333333333333333333333333333333),
0095       BOOST_MATH_BIG_CONSTANT(T, 64, -0.0083333333333333333333333333333333333333333333333333),
0096       BOOST_MATH_BIG_CONSTANT(T, 64, 0.003968253968253968253968253968253968253968253968254),
0097       BOOST_MATH_BIG_CONSTANT(T, 64, -0.0041666666666666666666666666666666666666666666666667),
0098       BOOST_MATH_BIG_CONSTANT(T, 64, 0.0075757575757575757575757575757575757575757575757576),
0099       BOOST_MATH_BIG_CONSTANT(T, 64, -0.021092796092796092796092796092796092796092796092796),
0100       BOOST_MATH_BIG_CONSTANT(T, 64, 0.083333333333333333333333333333333333333333333333333),
0101       BOOST_MATH_BIG_CONSTANT(T, 64, -0.44325980392156862745098039215686274509803921568627),
0102       BOOST_MATH_BIG_CONSTANT(T, 64, 3.0539543302701197438039543302701197438039543302701),
0103       BOOST_MATH_BIG_CONSTANT(T, 64, -26.456212121212121212121212121212121212121212121212),
0104       BOOST_MATH_BIG_CONSTANT(T, 64, 281.4601449275362318840579710144927536231884057971),
0105    };
0106    x -= 1;
0107    T result = log(x);
0108    result += 1 / (2 * x);
0109    T z = 1 / (x*x);
0110    result -= z * tools::evaluate_polynomial(P, z);
0111    return result;
0112 }
0113 //
0114 // 17-digit precision for x >= 10:
0115 //
0116 template <class T>
0117 inline T digamma_imp_large(T x, const std::integral_constant<int, 53>*)
0118 {
0119    BOOST_MATH_STD_USING // ADL of std functions.
0120    static const T P[] = {
0121       BOOST_MATH_BIG_CONSTANT(T, 53, 0.083333333333333333333333333333333333333333333333333),
0122       BOOST_MATH_BIG_CONSTANT(T, 53, -0.0083333333333333333333333333333333333333333333333333),
0123       BOOST_MATH_BIG_CONSTANT(T, 53, 0.003968253968253968253968253968253968253968253968254),
0124       BOOST_MATH_BIG_CONSTANT(T, 53, -0.0041666666666666666666666666666666666666666666666667),
0125       BOOST_MATH_BIG_CONSTANT(T, 53, 0.0075757575757575757575757575757575757575757575757576),
0126       BOOST_MATH_BIG_CONSTANT(T, 53, -0.021092796092796092796092796092796092796092796092796),
0127       BOOST_MATH_BIG_CONSTANT(T, 53, 0.083333333333333333333333333333333333333333333333333),
0128       BOOST_MATH_BIG_CONSTANT(T, 53, -0.44325980392156862745098039215686274509803921568627)
0129    };
0130    x -= 1;
0131    T result = log(x);
0132    result += 1 / (2 * x);
0133    T z = 1 / (x*x);
0134    result -= z * tools::evaluate_polynomial(P, z);
0135    return result;
0136 }
0137 //
0138 // 9-digit precision for x >= 10:
0139 //
0140 template <class T>
0141 inline T digamma_imp_large(T x, const std::integral_constant<int, 24>*)
0142 {
0143    BOOST_MATH_STD_USING // ADL of std functions.
0144    static const T P[] = {
0145       BOOST_MATH_BIG_CONSTANT(T, 24, 0.083333333333333333333333333333333333333333333333333),
0146       BOOST_MATH_BIG_CONSTANT(T, 24, -0.0083333333333333333333333333333333333333333333333333),
0147       BOOST_MATH_BIG_CONSTANT(T, 24, 0.003968253968253968253968253968253968253968253968254)
0148    };
0149    x -= 1;
0150    T result = log(x);
0151    result += 1 / (2 * x);
0152    T z = 1 / (x*x);
0153    result -= z * tools::evaluate_polynomial(P, z);
0154    return result;
0155 }
0156 //
0157 // Fully generic asymptotic expansion in terms of Bernoulli numbers, see:
0158 // http://functions.wolfram.com/06.14.06.0012.01
0159 //
0160 template <class T>
0161 struct digamma_series_func
0162 {
0163 private:
0164    int k;
0165    T xx;
0166    T term;
0167 public:
0168    digamma_series_func(T x) : k(1), xx(x * x), term(1 / (x * x)) {}
0169    T operator()()
0170    {
0171       T result = term * boost::math::bernoulli_b2n<T>(k) / (2 * k);
0172       term /= xx;
0173       ++k;
0174       return result;
0175    }
0176    typedef T result_type;
0177 };
0178
0179 template <class T, class Policy>
0180 inline T digamma_imp_large(T x, const Policy& pol, const std::integral_constant<int, 0>*)
0181 {
0182    BOOST_MATH_STD_USING
0183    digamma_series_func<T> s(x);
0184    T result = log(x) - 1 / (2 * x);
0185    std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0186    result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, -result);
0187    result = -result;
0188    policies::check_series_iterations<T>("boost::math::digamma<%1%>(%1%)", max_iter, pol);
0189    return result;
0190 }
0191 //
0192 // Now follow rational approximations over the range [1,2].
0193 //
0194 // 35-digit precision:
0195 //
0196 template <class T>
0197 T digamma_imp_1_2(T x, const std::integral_constant<int, 113>*)
0198 {
0199    //
0200    // Now the approximation, we use the form:
0201    //
0202    // digamma(x) = (x - root) * (Y + R(x-1))
0203    //
0204    // Where root is the location of the positive root of digamma,
0205    // Y is a constant, and R is optimised for low absolute error
0206    // compared to Y.
0207    //
0208    // Max error found at 128-bit long double precision:  5.541e-35
0209    // Maximum Deviation Found (approximation error):     1.965e-35
0210    //
0211    static const float Y = 0.99558162689208984375F;
0212
0213    static const T root1 = T(1569415565) / 1073741824uL;
0214    static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
0215    static const T root3 = ((T(111616537) / 1073741824uL) / 1073741824uL) / 1073741824uL;
0216    static const T root4 = (((T(503992070) / 1073741824uL) / 1073741824uL) / 1073741824uL) / 1073741824uL;
0217    static const T root5 = BOOST_MATH_BIG_CONSTANT(T, 113, 0.52112228569249997894452490385577338504019838794544e-36);
0218
0219    static const T P[] = {
0220       BOOST_MATH_BIG_CONSTANT(T, 113, 0.25479851061131551526977464225335883769),
0221       BOOST_MATH_BIG_CONSTANT(T, 113, -0.18684290534374944114622235683619897417),
0222       BOOST_MATH_BIG_CONSTANT(T, 113, -0.80360876047931768958995775910991929922),
0223       BOOST_MATH_BIG_CONSTANT(T, 113, -0.67227342794829064330498117008564270136),
0224       BOOST_MATH_BIG_CONSTANT(T, 113, -0.26569010991230617151285010695543858005),
0225       BOOST_MATH_BIG_CONSTANT(T, 113, -0.05775672694575986971640757748003553385),
0226       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0071432147823164975485922555833274240665),
0227       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00048740753910766168912364555706064993274),
0228       BOOST_MATH_BIG_CONSTANT(T, 113, -0.16454996865214115723416538844975174761e-4),
0229       BOOST_MATH_BIG_CONSTANT(T, 113, -0.20327832297631728077731148515093164955e-6)
0230    };
0231    static const T Q[] = {
0232       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0233       BOOST_MATH_BIG_CONSTANT(T, 113, 2.6210924610812025425088411043163287646),
0234       BOOST_MATH_BIG_CONSTANT(T, 113, 2.6850757078559596612621337395886392594),
0235       BOOST_MATH_BIG_CONSTANT(T, 113, 1.4320913706209965531250495490639289418),
0236       BOOST_MATH_BIG_CONSTANT(T, 113, 0.4410872083455009362557012239501953402),
0237       BOOST_MATH_BIG_CONSTANT(T, 113, 0.081385727399251729505165509278152487225),
0238       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0089478633066857163432104815183858149496),
0239       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00055861622855066424871506755481997374154),
0240       BOOST_MATH_BIG_CONSTANT(T, 113, 0.1760168552357342401304462967950178554e-4),
0241       BOOST_MATH_BIG_CONSTANT(T, 113, 0.20585454493572473724556649516040874384e-6),
0242       BOOST_MATH_BIG_CONSTANT(T, 113, -0.90745971844439990284514121823069162795e-11),
0243       BOOST_MATH_BIG_CONSTANT(T, 113, 0.48857673606545846774761343500033283272e-13),
0244    };
0245    T g = x - root1;
0246    g -= root2;
0247    g -= root3;
0248    g -= root4;
0249    g -= root5;
0250    T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
0251    T result = g * Y + g * r;
0252
0253    return result;
0254 }
0255 //
0256 // 19-digit precision:
0257 //
0258 template <class T>
0259 T digamma_imp_1_2(T x, const std::integral_constant<int, 64>*)
0260 {
0261    //
0262    // Now the approximation, we use the form:
0263    //
0264    // digamma(x) = (x - root) * (Y + R(x-1))
0265    //
0266    // Where root is the location of the positive root of digamma,
0267    // Y is a constant, and R is optimised for low absolute error
0268    // compared to Y.
0269    //
0270    // Max error found at 80-bit long double precision:   5.016e-20
0271    // Maximum Deviation Found (approximation error):     3.575e-20
0272    //
0273    static const float Y = 0.99558162689208984375F;
0274
0275    static const T root1 = T(1569415565) / 1073741824uL;
0276    static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
0277    static const T root3 = BOOST_MATH_BIG_CONSTANT(T, 64, 0.9016312093258695918615325266959189453125e-19);
0278
0279    static const T P[] = {
0280       BOOST_MATH_BIG_CONSTANT(T, 64, 0.254798510611315515235),
0281       BOOST_MATH_BIG_CONSTANT(T, 64, -0.314628554532916496608),
0282       BOOST_MATH_BIG_CONSTANT(T, 64, -0.665836341559876230295),
0283       BOOST_MATH_BIG_CONSTANT(T, 64, -0.314767657147375752913),
0284       BOOST_MATH_BIG_CONSTANT(T, 64, -0.0541156266153505273939),
0285       BOOST_MATH_BIG_CONSTANT(T, 64, -0.00289268368333918761452)
0286    };
0287    static const T Q[] = {
0288       BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0289       BOOST_MATH_BIG_CONSTANT(T, 64, 2.1195759927055347547),
0290       BOOST_MATH_BIG_CONSTANT(T, 64, 1.54350554664961128724),
0291       BOOST_MATH_BIG_CONSTANT(T, 64, 0.486986018231042975162),
0292       BOOST_MATH_BIG_CONSTANT(T, 64, 0.0660481487173569812846),
0293       BOOST_MATH_BIG_CONSTANT(T, 64, 0.00298999662592323990972),
0294       BOOST_MATH_BIG_CONSTANT(T, 64, -0.165079794012604905639e-5),
0295       BOOST_MATH_BIG_CONSTANT(T, 64, 0.317940243105952177571e-7)
0296    };
0297    T g = x - root1;
0298    g -= root2;
0299    g -= root3;
0300    T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
0301    T result = g * Y + g * r;
0302
0303    return result;
0304 }
0305 //
0306 // 18-digit precision:
0307 //
0308 template <class T>
0309 T digamma_imp_1_2(T x, const std::integral_constant<int, 53>*)
0310 {
0311    //
0312    // Now the approximation, we use the form:
0313    //
0314    // digamma(x) = (x - root) * (Y + R(x-1))
0315    //
0316    // Where root is the location of the positive root of digamma,
0317    // Y is a constant, and R is optimised for low absolute error
0318    // compared to Y.
0319    //
0320    // Maximum Deviation Found:               1.466e-18
0321    // At double precision, max error found:  2.452e-17
0322    //
0323    static const float Y = 0.99558162689208984F;
0324
0325    static const T root1 = T(1569415565) / 1073741824uL;
0326    static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
0327    static const T root3 = BOOST_MATH_BIG_CONSTANT(T, 53, 0.9016312093258695918615325266959189453125e-19);
0328
0329    static const T P[] = {
0330       BOOST_MATH_BIG_CONSTANT(T, 53, 0.25479851061131551),
0331       BOOST_MATH_BIG_CONSTANT(T, 53, -0.32555031186804491),
0332       BOOST_MATH_BIG_CONSTANT(T, 53, -0.65031853770896507),
0333       BOOST_MATH_BIG_CONSTANT(T, 53, -0.28919126444774784),
0334       BOOST_MATH_BIG_CONSTANT(T, 53, -0.045251321448739056),
0335       BOOST_MATH_BIG_CONSTANT(T, 53, -0.0020713321167745952)
0336    };
0337    static const T Q[] = {
0338       BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
0339       BOOST_MATH_BIG_CONSTANT(T, 53, 2.0767117023730469),
0340       BOOST_MATH_BIG_CONSTANT(T, 53, 1.4606242909763515),
0341       BOOST_MATH_BIG_CONSTANT(T, 53, 0.43593529692665969),
0342       BOOST_MATH_BIG_CONSTANT(T, 53, 0.054151797245674225),
0343       BOOST_MATH_BIG_CONSTANT(T, 53, 0.0021284987017821144),
0344       BOOST_MATH_BIG_CONSTANT(T, 53, -0.55789841321675513e-6)
0345    };
0346    T g = x - root1;
0347    g -= root2;
0348    g -= root3;
0349    T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
0350    T result = g * Y + g * r;
0351
0352    return result;
0353 }
0354 //
0355 // 9-digit precision:
0356 //
0357 template <class T>
0358 inline T digamma_imp_1_2(T x, const std::integral_constant<int, 24>*)
0359 {
0360    //
0361    // Now the approximation, we use the form:
0362    //
0363    // digamma(x) = (x - root) * (Y + R(x-1))
0364    //
0365    // Where root is the location of the positive root of digamma,
0366    // Y is a constant, and R is optimised for low absolute error
0367    // compared to Y.
0368    //
0369    // Maximum Deviation Found:              3.388e-010
0370    // At float precision, max error found:  2.008725e-008
0371    //
0372    static const float Y = 0.99558162689208984f;
0373    static const T root = 1532632.0f / 1048576;
0374    static const T root_minor = static_cast<T>(0.3700660185912626595423257213284682051735604e-6L);
0375    static const T P[] = {
0376       0.25479851023250261e0f,
0377       -0.44981331915268368e0f,
0378       -0.43916936919946835e0f,
0379       -0.61041765350579073e-1f
0380    };
0381    static const T Q[] = {
0382       0.1e1f,
0383       0.15890202430554952e1f,
0384       0.65341249856146947e0f,
0385       0.63851690523355715e-1f
0386    };
0387    T g = x - root;
0388    g -= root_minor;
0389    T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
0390    T result = g * Y + g * r;
0391
0392    return result;
0393 }
0394
0395 template <class T, class Tag, class Policy>
0396 T digamma_imp(T x, const Tag* t, const Policy& pol)
0397 {
0398    //
0399    // This handles reflection of negative arguments, and all our
0400    // error handling, then forwards to the T-specific approximation.
0401    //
0402    BOOST_MATH_STD_USING // ADL of std functions.
0403
0404    T result = 0;
0405    //
0406    // Check for negative arguments and use reflection:
0407    //
0408    if(x <= -1)
0409    {
0410       // Reflect:
0411       x = 1 - x;
0412       // Argument reduction for tan:
0413       T remainder = x - floor(x);
0414       // Shift to negative if > 0.5:
0415       if(remainder > T(0.5))
0416       {
0417          remainder -= 1;
0418       }
0419       //
0420       // check for evaluation at a negative pole:
0421       //
0422       if(remainder == 0)
0423       {
0424          return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", nullptr, (1-x), pol);
0425       }
0426       result = constants::pi<T>() / tan(constants::pi<T>() * remainder);
0427    }
0428    if(x == 0)
0429       return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", nullptr, x, pol);
0430    //
0431    // If we're above the lower-limit for the
0432    // asymptotic expansion then use it:
0433    //
0434    if(x >= digamma_large_lim(t))
0435    {
0436       result += digamma_imp_large(x, t);
0437    }
0438    else
0439    {
0440       //
0441       // If x > 2 reduce to the interval [1,2]:
0442       //
0443       while(x > 2)
0444       {
0445          x -= 1;
0446          result += 1/x;
0447       }
0448       //
0449       // If x < 1 use recurrence to shift to > 1:
0450       //
0451       while(x < 1)
0452       {
0453          result -= 1/x;
0454          x += 1;
0455       }
0456       result += digamma_imp_1_2(x, t);
0457    }
0458    return result;
0459 }
0460
0461 template <class T, class Policy>
0462 T digamma_imp(T x, const std::integral_constant<int, 0>* t, const Policy& pol)
0463 {
0464    //
0465    // This handles reflection of negative arguments, and all our
0466    // error handling, then forwards to the T-specific approximation.
0467    //
0468    BOOST_MATH_STD_USING // ADL of std functions.
0469
0470    T result = 0;
0471    //
0472    // Check for negative arguments and use reflection:
0473    //
0474    if(x <= -1)
0475    {
0476       // Reflect:
0477       x = 1 - x;
0478       // Argument reduction for tan:
0479       T remainder = x - floor(x);
0480       // Shift to negative if > 0.5:
0481       if(remainder > T(0.5))
0482       {
0483          remainder -= 1;
0484       }
0485       //
0486       // check for evaluation at a negative pole:
0487       //
0488       if(remainder == 0)
0489       {
0490          return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", nullptr, (1 - x), pol);
0491       }
0492       result = constants::pi<T>() / tan(constants::pi<T>() * remainder);
0493    }
0494    if(x == 0)
0495       return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", nullptr, x, pol);
0496    //
0497    // If we're above the lower-limit for the
0498    // asymptotic expansion then use it, the
0499    // limit is a linear interpolation with
0500    // limit = 10 at 50 bit precision and
0501    // limit = 250 at 1000 bit precision.
0502    //
0503    int lim = 10 + ((tools::digits<T>() - 50) * 240L) / 950;
0504    T two_x = ldexp(x, 1);
0505    if(x >= lim)
0506    {
0507       result += digamma_imp_large(x, pol, t);
0508    }
0509    else if(floor(x) == x)
0510    {
0511       //
0512       // Special case for integer arguments, see
0513       // http://functions.wolfram.com/06.14.03.0001.01
0514       //
0515       result = -constants::euler<T, Policy>();
0516       T val = 1;
0517       while(val < x)
0518       {
0519          result += 1 / val;
0520          val += 1;
0521       }
0522    }
0523    else if(floor(two_x) == two_x)
0524    {
0525       //
0526       // Special case for half integer arguments, see:
0527       // http://functions.wolfram.com/06.14.03.0007.01
0528       //
0529       result = -2 * constants::ln_two<T, Policy>() - constants::euler<T, Policy>();
0530       int n = itrunc(x);
0531       if(n)
0532       {
0533          for(int k = 1; k < n; ++k)
0534             result += 1 / T(k);
0535          for(int k = n; k <= 2 * n - 1; ++k)
0536             result += 2 / T(k);
0537       }
0538    }
0539    else
0540    {
0541       //
0542       // Rescale so we can use the asymptotic expansion:
0543       //
0544       while(x < lim)
0545       {
0546          result -= 1 / x;
0547          x += 1;
0548       }
0549       result += digamma_imp_large(x, pol, t);
0550    }
0551    return result;
0552 }
0553 //
0554 // Initializer: ensure all our constants are initialized prior to the first call of main:
0555 //
0556 template <class T, class Policy>
0557 struct digamma_initializer
0558 {
0559    struct init
0560    {
0561       init()
0562       {
0563          typedef typename policies::precision<T, Policy>::type precision_type;
0564          do_init(std::integral_constant<bool, precision_type::value && (precision_type::value <= 113)>());
0565       }
0566       void do_init(const std::true_type&)
0567       {
0568          boost::math::digamma(T(1.5), Policy());
0569          boost::math::digamma(T(500), Policy());
0570       }
0571       void do_init(const std::false_type&){}
0572       void force_instantiate()const{}
0573    };
0574    static const init initializer;
0575    static void force_instantiate()
0576    {
0577       initializer.force_instantiate();
0578    }
0579 };
0580
0581 template <class T, class Policy>
0582 const typename digamma_initializer<T, Policy>::init digamma_initializer<T, Policy>::initializer;
0583
0584 } // namespace detail
0585
0586 template <class T, class Policy>
0587 inline typename tools::promote_args<T>::type
0588    digamma(T x, const Policy&)
0589 {
0590    typedef typename tools::promote_args<T>::type result_type;
0591    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0592    typedef typename policies::precision<T, Policy>::type precision_type;
0593    typedef std::integral_constant<int,
0594       (precision_type::value <= 0) || (precision_type::value > 113) ? 0 :
0595       precision_type::value <= 24 ? 24 :
0596       precision_type::value <= 53 ? 53 :
0597       precision_type::value <= 64 ? 64 :
0598       precision_type::value <= 113 ? 113 : 0 > tag_type;
0599    typedef typename policies::normalise<
0600       Policy,
0601       policies::promote_float<false>,
0602       policies::promote_double<false>,
0603       policies::discrete_quantile<>,
0604       policies::assert_undefined<> >::type forwarding_policy;
0605
0606    // Force initialization of constants:
0607    detail::digamma_initializer<value_type, forwarding_policy>::force_instantiate();
0608
0609    return policies::checked_narrowing_cast<result_type, Policy>(detail::digamma_imp(
0610       static_cast<value_type>(x),
0611       static_cast<const tag_type*>(nullptr), forwarding_policy()), "boost::math::digamma<%1%>(%1%)");
0612 }
0613
0614 template <class T>
0615 inline typename tools::promote_args<T>::type
0616    digamma(T x)
0617 {
0618    return digamma(x, policies::policy<>());
0619 }
0620
0621 } // namespace math
0622 } // namespace boost
0623
0624 #ifdef _MSC_VER
0625 #pragma warning(pop)
0626 #endif
0627
0628 #endif
0629