math/special_functions/digamma.hpp

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