ceph/src/boost/boost/math/special_functions/digamma.hpp

   1 //  (C) Copyright John Maddock 2006.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_SF_DIGAMMA_HPP
   7 #define BOOST_MATH_SF_DIGAMMA_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #pragma warning(push)
  12 #pragma warning(disable:4702) // Unreachable code (release mode only warning)
  13 #endif
  14
  15 #include <boost/math/special_functions/math_fwd.hpp>
  16 #include <boost/math/tools/rational.hpp>
  17 #include <boost/math/tools/series.hpp>
  18 #include <boost/math/tools/promotion.hpp>
  19 #include <boost/math/policies/error_handling.hpp>
  20 #include <boost/math/constants/constants.hpp>
  21 #include <boost/math/tools/big_constant.hpp>
  22
  23 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
  24 //
  25 // This is the only way we can avoid
  26 // warning: non-standard suffix on floating constant [-Wpedantic]
  27 // when building with -Wall -pedantic.  Neither __extension__
  28 // nor #pragma diagnostic ignored work :(
  29 //
  30 #pragma GCC system_header
  31 #endif
  32
  33 namespace boost{
  34 namespace math{
  35 namespace detail{
  36 //
  37 // Begin by defining the smallest value for which it is safe to
  38 // use the asymptotic expansion for digamma:
  39 //
  40 inline unsigned digamma_large_lim(const std::integral_constant<int, 0>*)
  41 {  return 20;  }
  42 inline unsigned digamma_large_lim(const std::integral_constant<int, 113>*)
  43 {  return 20;  }
  44 inline unsigned digamma_large_lim(const void*)
  45 {  return 10;  }
  46 //
  47 // Implementations of the asymptotic expansion come next,
  48 // the coefficients of the series have been evaluated
  49 // in advance at high precision, and the series truncated
  50 // at the first term that's too small to effect the result.
  51 // Note that the series becomes divergent after a while
  52 // so truncation is very important.
  53 //
  54 // This first one gives 34-digit precision for x >= 20:
  55 //
  56 template <class T>
  57 inline T digamma_imp_large(T x, const std::integral_constant<int, 113>*)
  58 {
  59    BOOST_MATH_STD_USING // ADL of std functions.
  60    static const T P[] = {
  61       BOOST_MATH_BIG_CONSTANT(T, 113, 0.083333333333333333333333333333333333333333333333333),
  62       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0083333333333333333333333333333333333333333333333333),
  63       BOOST_MATH_BIG_CONSTANT(T, 113, 0.003968253968253968253968253968253968253968253968254),
  64       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0041666666666666666666666666666666666666666666666667),
  65       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0075757575757575757575757575757575757575757575757576),
  66       BOOST_MATH_BIG_CONSTANT(T, 113, -0.021092796092796092796092796092796092796092796092796),
  67       BOOST_MATH_BIG_CONSTANT(T, 113, 0.083333333333333333333333333333333333333333333333333),
  68       BOOST_MATH_BIG_CONSTANT(T, 113, -0.44325980392156862745098039215686274509803921568627),
  69       BOOST_MATH_BIG_CONSTANT(T, 113, 3.0539543302701197438039543302701197438039543302701),
  70       BOOST_MATH_BIG_CONSTANT(T, 113, -26.456212121212121212121212121212121212121212121212),
  71       BOOST_MATH_BIG_CONSTANT(T, 113, 281.4601449275362318840579710144927536231884057971),
  72       BOOST_MATH_BIG_CONSTANT(T, 113, -3607.510546398046398046398046398046398046398046398),
  73       BOOST_MATH_BIG_CONSTANT(T, 113, 54827.583333333333333333333333333333333333333333333),
  74       BOOST_MATH_BIG_CONSTANT(T, 113, -974936.82385057471264367816091954022988505747126437),
  75       BOOST_MATH_BIG_CONSTANT(T, 113, 20052695.796688078946143462272494530559046688078946),
  76       BOOST_MATH_BIG_CONSTANT(T, 113, -472384867.72162990196078431372549019607843137254902),
  77       BOOST_MATH_BIG_CONSTANT(T, 113, 12635724795.916666666666666666666666666666666666667)
  78    };
  79    x -= 1;
  80    T result = log(x);
  81    result += 1 / (2 * x);
  82    T z = 1 / (x*x);
  83    result -= z * tools::evaluate_polynomial(P, z);
  84    return result;
  85 }
  86 //
  87 // 19-digit precision for x >= 10:
  88 //
  89 template <class T>
  90 inline T digamma_imp_large(T x, const std::integral_constant<int, 64>*)
  91 {
  92    BOOST_MATH_STD_USING // ADL of std functions.
  93    static const T P[] = {
  94       BOOST_MATH_BIG_CONSTANT(T, 64, 0.083333333333333333333333333333333333333333333333333),
  95       BOOST_MATH_BIG_CONSTANT(T, 64, -0.0083333333333333333333333333333333333333333333333333),
  96       BOOST_MATH_BIG_CONSTANT(T, 64, 0.003968253968253968253968253968253968253968253968254),
  97       BOOST_MATH_BIG_CONSTANT(T, 64, -0.0041666666666666666666666666666666666666666666666667),
  98       BOOST_MATH_BIG_CONSTANT(T, 64, 0.0075757575757575757575757575757575757575757575757576),
  99       BOOST_MATH_BIG_CONSTANT(T, 64, -0.021092796092796092796092796092796092796092796092796),
 100       BOOST_MATH_BIG_CONSTANT(T, 64, 0.083333333333333333333333333333333333333333333333333),
 101       BOOST_MATH_BIG_CONSTANT(T, 64, -0.44325980392156862745098039215686274509803921568627),
 102       BOOST_MATH_BIG_CONSTANT(T, 64, 3.0539543302701197438039543302701197438039543302701),
 103       BOOST_MATH_BIG_CONSTANT(T, 64, -26.456212121212121212121212121212121212121212121212),
 104       BOOST_MATH_BIG_CONSTANT(T, 64, 281.4601449275362318840579710144927536231884057971),
 105    };
 106    x -= 1;
 107    T result = log(x);
 108    result += 1 / (2 * x);
 109    T z = 1 / (x*x);
 110    result -= z * tools::evaluate_polynomial(P, z);
 111    return result;
 112 }
 113 //
 114 // 17-digit precision for x >= 10:
 115 //
 116 template <class T>
 117 inline T digamma_imp_large(T x, const std::integral_constant<int, 53>*)
 118 {
 119    BOOST_MATH_STD_USING // ADL of std functions.
 120    static const T P[] = {
 121       BOOST_MATH_BIG_CONSTANT(T, 53, 0.083333333333333333333333333333333333333333333333333),
 122       BOOST_MATH_BIG_CONSTANT(T, 53, -0.0083333333333333333333333333333333333333333333333333),
 123       BOOST_MATH_BIG_CONSTANT(T, 53, 0.003968253968253968253968253968253968253968253968254),
 124       BOOST_MATH_BIG_CONSTANT(T, 53, -0.0041666666666666666666666666666666666666666666666667),
 125       BOOST_MATH_BIG_CONSTANT(T, 53, 0.0075757575757575757575757575757575757575757575757576),
 126       BOOST_MATH_BIG_CONSTANT(T, 53, -0.021092796092796092796092796092796092796092796092796),
 127       BOOST_MATH_BIG_CONSTANT(T, 53, 0.083333333333333333333333333333333333333333333333333),
 128       BOOST_MATH_BIG_CONSTANT(T, 53, -0.44325980392156862745098039215686274509803921568627)
 129    };
 130    x -= 1;
 131    T result = log(x);
 132    result += 1 / (2 * x);
 133    T z = 1 / (x*x);
 134    result -= z * tools::evaluate_polynomial(P, z);
 135    return result;
 136 }
 137 //
 138 // 9-digit precision for x >= 10:
 139 //
 140 template <class T>
 141 inline T digamma_imp_large(T x, const std::integral_constant<int, 24>*)
 142 {
 143    BOOST_MATH_STD_USING // ADL of std functions.
 144    static const T P[] = {
 145       BOOST_MATH_BIG_CONSTANT(T, 24, 0.083333333333333333333333333333333333333333333333333),
 146       BOOST_MATH_BIG_CONSTANT(T, 24, -0.0083333333333333333333333333333333333333333333333333),
 147       BOOST_MATH_BIG_CONSTANT(T, 24, 0.003968253968253968253968253968253968253968253968254)
 148    };
 149    x -= 1;
 150    T result = log(x);
 151    result += 1 / (2 * x);
 152    T z = 1 / (x*x);
 153    result -= z * tools::evaluate_polynomial(P, z);
 154    return result;
 155 }
 156 //
 157 // Fully generic asymptotic expansion in terms of Bernoulli numbers, see:
 158 // http://functions.wolfram.com/06.14.06.0012.01
 159 //
 160 template <class T>
 161 struct digamma_series_func
 162 {
 163 private:
 164    int k;
 165    T xx;
 166    T term;
 167 public:
 168    digamma_series_func(T x) : k(1), xx(x * x), term(1 / (x * x)) {}
 169    T operator()()
 170    {
 171       T result = term * boost::math::bernoulli_b2n<T>(k) / (2 * k);
 172       term /= xx;
 173       ++k;
 174       return result;
 175    }
 176    typedef T result_type;
 177 };
 178
 179 template <class T, class Policy>
 180 inline T digamma_imp_large(T x, const Policy& pol, const std::integral_constant<int, 0>*)
 181 {
 182    BOOST_MATH_STD_USING
 183    digamma_series_func<T> s(x);
 184    T result = log(x) - 1 / (2 * x);
 185    std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 186    result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, -result);
 187    result = -result;
 188    policies::check_series_iterations<T>("boost::math::digamma<%1%>(%1%)", max_iter, pol);
 189    return result;
 190 }
 191 //
 192 // Now follow rational approximations over the range [1,2].
 193 //
 194 // 35-digit precision:
 195 //
 196 template <class T>
 197 T digamma_imp_1_2(T x, const std::integral_constant<int, 113>*)
 198 {
 199    //
 200    // Now the approximation, we use the form:
 201    //
 202    // digamma(x) = (x - root) * (Y + R(x-1))
 203    //
 204    // Where root is the location of the positive root of digamma,
 205    // Y is a constant, and R is optimised for low absolute error
 206    // compared to Y.
 207    //
 208    // Max error found at 128-bit long double precision:  5.541e-35
 209    // Maximum Deviation Found (approximation error):     1.965e-35
 210    //
 211    static const float Y = 0.99558162689208984375F;
 212
 213    static const T root1 = T(1569415565) / 1073741824uL;
 214    static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
 215    static const T root3 = ((T(111616537) / 1073741824uL) / 1073741824uL) / 1073741824uL;
 216    static const T root4 = (((T(503992070) / 1073741824uL) / 1073741824uL) / 1073741824uL) / 1073741824uL;
 217    static const T root5 = BOOST_MATH_BIG_CONSTANT(T, 113, 0.52112228569249997894452490385577338504019838794544e-36);
 218
 219    static const T P[] = {
 220       BOOST_MATH_BIG_CONSTANT(T, 113, 0.25479851061131551526977464225335883769),
 221       BOOST_MATH_BIG_CONSTANT(T, 113, -0.18684290534374944114622235683619897417),
 222       BOOST_MATH_BIG_CONSTANT(T, 113, -0.80360876047931768958995775910991929922),
 223       BOOST_MATH_BIG_CONSTANT(T, 113, -0.67227342794829064330498117008564270136),
 224       BOOST_MATH_BIG_CONSTANT(T, 113, -0.26569010991230617151285010695543858005),
 225       BOOST_MATH_BIG_CONSTANT(T, 113, -0.05775672694575986971640757748003553385),
 226       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0071432147823164975485922555833274240665),
 227       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00048740753910766168912364555706064993274),
 228       BOOST_MATH_BIG_CONSTANT(T, 113, -0.16454996865214115723416538844975174761e-4),
 229       BOOST_MATH_BIG_CONSTANT(T, 113, -0.20327832297631728077731148515093164955e-6)
 230    };
 231    static const T Q[] = {
 232       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
 233       BOOST_MATH_BIG_CONSTANT(T, 113, 2.6210924610812025425088411043163287646),
 234       BOOST_MATH_BIG_CONSTANT(T, 113, 2.6850757078559596612621337395886392594),
 235       BOOST_MATH_BIG_CONSTANT(T, 113, 1.4320913706209965531250495490639289418),
 236       BOOST_MATH_BIG_CONSTANT(T, 113, 0.4410872083455009362557012239501953402),
 237       BOOST_MATH_BIG_CONSTANT(T, 113, 0.081385727399251729505165509278152487225),
 238       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0089478633066857163432104815183858149496),
 239       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00055861622855066424871506755481997374154),
 240       BOOST_MATH_BIG_CONSTANT(T, 113, 0.1760168552357342401304462967950178554e-4),
 241       BOOST_MATH_BIG_CONSTANT(T, 113, 0.20585454493572473724556649516040874384e-6),
 242       BOOST_MATH_BIG_CONSTANT(T, 113, -0.90745971844439990284514121823069162795e-11),
 243       BOOST_MATH_BIG_CONSTANT(T, 113, 0.48857673606545846774761343500033283272e-13),
 244    };
 245    T g = x - root1;
 246    g -= root2;
 247    g -= root3;
 248    g -= root4;
 249    g -= root5;
 250    T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
 251    T result = g * Y + g * r;
 252
 253    return result;
 254 }
 255 //
 256 // 19-digit precision:
 257 //
 258 template <class T>
 259 T digamma_imp_1_2(T x, const std::integral_constant<int, 64>*)
 260 {
 261    //
 262    // Now the approximation, we use the form:
 263    //
 264    // digamma(x) = (x - root) * (Y + R(x-1))
 265    //
 266    // Where root is the location of the positive root of digamma,
 267    // Y is a constant, and R is optimised for low absolute error
 268    // compared to Y.
 269    //
 270    // Max error found at 80-bit long double precision:   5.016e-20
 271    // Maximum Deviation Found (approximation error):     3.575e-20
 272    //
 273    static const float Y = 0.99558162689208984375F;
 274
 275    static const T root1 = T(1569415565) / 1073741824uL;
 276    static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
 277    static const T root3 = BOOST_MATH_BIG_CONSTANT(T, 64, 0.9016312093258695918615325266959189453125e-19);
 278
 279    static const T P[] = {
 280       BOOST_MATH_BIG_CONSTANT(T, 64, 0.254798510611315515235),
 281       BOOST_MATH_BIG_CONSTANT(T, 64, -0.314628554532916496608),
 282       BOOST_MATH_BIG_CONSTANT(T, 64, -0.665836341559876230295),
 283       BOOST_MATH_BIG_CONSTANT(T, 64, -0.314767657147375752913),
 284       BOOST_MATH_BIG_CONSTANT(T, 64, -0.0541156266153505273939),
 285       BOOST_MATH_BIG_CONSTANT(T, 64, -0.00289268368333918761452)
 286    };
 287    static const T Q[] = {
 288       BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
 289       BOOST_MATH_BIG_CONSTANT(T, 64, 2.1195759927055347547),
 290       BOOST_MATH_BIG_CONSTANT(T, 64, 1.54350554664961128724),
 291       BOOST_MATH_BIG_CONSTANT(T, 64, 0.486986018231042975162),
 292       BOOST_MATH_BIG_CONSTANT(T, 64, 0.0660481487173569812846),
 293       BOOST_MATH_BIG_CONSTANT(T, 64, 0.00298999662592323990972),
 294       BOOST_MATH_BIG_CONSTANT(T, 64, -0.165079794012604905639e-5),
 295       BOOST_MATH_BIG_CONSTANT(T, 64, 0.317940243105952177571e-7)
 296    };
 297    T g = x - root1;
 298    g -= root2;
 299    g -= root3;
 300    T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
 301    T result = g * Y + g * r;
 302
 303    return result;
 304 }
 305 //
 306 // 18-digit precision:
 307 //
 308 template <class T>
 309 T digamma_imp_1_2(T x, const std::integral_constant<int, 53>*)
 310 {
 311    //
 312    // Now the approximation, we use the form:
 313    //
 314    // digamma(x) = (x - root) * (Y + R(x-1))
 315    //
 316    // Where root is the location of the positive root of digamma,
 317    // Y is a constant, and R is optimised for low absolute error
 318    // compared to Y.
 319    //
 320    // Maximum Deviation Found:               1.466e-18
 321    // At double precision, max error found:  2.452e-17
 322    //
 323    static const float Y = 0.99558162689208984F;
 324
 325    static const T root1 = T(1569415565) / 1073741824uL;
 326    static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
 327    static const T root3 = BOOST_MATH_BIG_CONSTANT(T, 53, 0.9016312093258695918615325266959189453125e-19);
 328
 329    static const T P[] = {
 330       BOOST_MATH_BIG_CONSTANT(T, 53, 0.25479851061131551),
 331       BOOST_MATH_BIG_CONSTANT(T, 53, -0.32555031186804491),
 332       BOOST_MATH_BIG_CONSTANT(T, 53, -0.65031853770896507),
 333       BOOST_MATH_BIG_CONSTANT(T, 53, -0.28919126444774784),
 334       BOOST_MATH_BIG_CONSTANT(T, 53, -0.045251321448739056),
 335       BOOST_MATH_BIG_CONSTANT(T, 53, -0.0020713321167745952)
 336    };
 337    static const T Q[] = {
 338       BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
 339       BOOST_MATH_BIG_CONSTANT(T, 53, 2.0767117023730469),
 340       BOOST_MATH_BIG_CONSTANT(T, 53, 1.4606242909763515),
 341       BOOST_MATH_BIG_CONSTANT(T, 53, 0.43593529692665969),
 342       BOOST_MATH_BIG_CONSTANT(T, 53, 0.054151797245674225),
 343       BOOST_MATH_BIG_CONSTANT(T, 53, 0.0021284987017821144),
 344       BOOST_MATH_BIG_CONSTANT(T, 53, -0.55789841321675513e-6)
 345    };
 346    T g = x - root1;
 347    g -= root2;
 348    g -= root3;
 349    T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
 350    T result = g * Y + g * r;
 351
 352    return result;
 353 }
 354 //
 355 // 9-digit precision:
 356 //
 357 template <class T>
 358 inline T digamma_imp_1_2(T x, const std::integral_constant<int, 24>*)
 359 {
 360    //
 361    // Now the approximation, we use the form:
 362    //
 363    // digamma(x) = (x - root) * (Y + R(x-1))
 364    //
 365    // Where root is the location of the positive root of digamma,
 366    // Y is a constant, and R is optimised for low absolute error
 367    // compared to Y.
 368    //
 369    // Maximum Deviation Found:              3.388e-010
 370    // At float precision, max error found:  2.008725e-008
 371    //
 372    static const float Y = 0.99558162689208984f;
 373    static const T root = 1532632.0f / 1048576;
 374    static const T root_minor = static_cast<T>(0.3700660185912626595423257213284682051735604e-6L);
 375    static const T P[] = {
 376       0.25479851023250261e0f,
 377       -0.44981331915268368e0f,
 378       -0.43916936919946835e0f,
 379       -0.61041765350579073e-1f
 380    };
 381    static const T Q[] = {
 382       0.1e1,
 383       0.15890202430554952e1f,
 384       0.65341249856146947e0f,
 385       0.63851690523355715e-1f
 386    };
 387    T g = x - root;
 388    g -= root_minor;
 389    T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
 390    T result = g * Y + g * r;
 391
 392    return result;
 393 }
 394
 395 template <class T, class Tag, class Policy>
 396 T digamma_imp(T x, const Tag* t, const Policy& pol)
 397 {
 398    //
 399    // This handles reflection of negative arguments, and all our
 400    // error handling, then forwards to the T-specific approximation.
 401    //
 402    BOOST_MATH_STD_USING // ADL of std functions.
 403
 404    T result = 0;
 405    //
 406    // Check for negative arguments and use reflection:
 407    //
 408    if(x <= -1)
 409    {
 410       // Reflect:
 411       x = 1 - x;
 412       // Argument reduction for tan:
 413       T remainder = x - floor(x);
 414       // Shift to negative if > 0.5:
 415       if(remainder > 0.5)
 416       {
 417          remainder -= 1;
 418       }
 419       //
 420       // check for evaluation at a negative pole:
 421       //
 422       if(remainder == 0)
 423       {
 424          return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, (1-x), pol);
 425       }
 426       result = constants::pi<T>() / tan(constants::pi<T>() * remainder);
 427    }
 428    if(x == 0)
 429       return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, x, pol);
 430    //
 431    // If we're above the lower-limit for the
 432    // asymptotic expansion then use it:
 433    //
 434    if(x >= digamma_large_lim(t))
 435    {
 436       result += digamma_imp_large(x, t);
 437    }
 438    else
 439    {
 440       //
 441       // If x > 2 reduce to the interval [1,2]:
 442       //
 443       while(x > 2)
 444       {
 445          x -= 1;
 446          result += 1/x;
 447       }
 448       //
 449       // If x < 1 use recurrence to shift to > 1:
 450       //
 451       while(x < 1)
 452       {
 453          result -= 1/x;
 454          x += 1;
 455       }
 456       result += digamma_imp_1_2(x, t);
 457    }
 458    return result;
 459 }
 460
 461 template <class T, class Policy>
 462 T digamma_imp(T x, const std::integral_constant<int, 0>* t, const Policy& pol)
 463 {
 464    //
 465    // This handles reflection of negative arguments, and all our
 466    // error handling, then forwards to the T-specific approximation.
 467    //
 468    BOOST_MATH_STD_USING // ADL of std functions.
 469
 470    T result = 0;
 471    //
 472    // Check for negative arguments and use reflection:
 473    //
 474    if(x <= -1)
 475    {
 476       // Reflect:
 477       x = 1 - x;
 478       // Argument reduction for tan:
 479       T remainder = x - floor(x);
 480       // Shift to negative if > 0.5:
 481       if(remainder > 0.5)
 482       {
 483          remainder -= 1;
 484       }
 485       //
 486       // check for evaluation at a negative pole:
 487       //
 488       if(remainder == 0)
 489       {
 490          return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, (1 - x), pol);
 491       }
 492       result = constants::pi<T>() / tan(constants::pi<T>() * remainder);
 493    }
 494    if(x == 0)
 495       return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, x, pol);
 496    //
 497    // If we're above the lower-limit for the
 498    // asymptotic expansion then use it, the
 499    // limit is a linear interpolation with
 500    // limit = 10 at 50 bit precision and
 501    // limit = 250 at 1000 bit precision.
 502    //
 503    int lim = 10 + ((tools::digits<T>() - 50) * 240L) / 950;
 504    T two_x = ldexp(x, 1);
 505    if(x >= lim)
 506    {
 507       result += digamma_imp_large(x, pol, t);
 508    }
 509    else if(floor(x) == x)
 510    {
 511       //
 512       // Special case for integer arguments, see
 513       // http://functions.wolfram.com/06.14.03.0001.01
 514       //
 515       result = -constants::euler<T, Policy>();
 516       T val = 1;
 517       while(val < x)
 518       {
 519          result += 1 / val;
 520          val += 1;
 521       }
 522    }
 523    else if(floor(two_x) == two_x)
 524    {
 525       //
 526       // Special case for half integer arguments, see:
 527       // http://functions.wolfram.com/06.14.03.0007.01
 528       //
 529       result = -2 * constants::ln_two<T, Policy>() - constants::euler<T, Policy>();
 530       int n = itrunc(x);
 531       if(n)
 532       {
 533          for(int k = 1; k < n; ++k)
 534             result += 1 / T(k);
 535          for(int k = n; k <= 2 * n - 1; ++k)
 536             result += 2 / T(k);
 537       }
 538    }
 539    else
 540    {
 541       //
 542       // Rescale so we can use the asymptotic expansion:
 543       //
 544       while(x < lim)
 545       {
 546          result -= 1 / x;
 547          x += 1;
 548       }
 549       result += digamma_imp_large(x, pol, t);
 550    }
 551    return result;
 552 }
 553 //
 554 // Initializer: ensure all our constants are initialized prior to the first call of main:
 555 //
 556 template <class T, class Policy>
 557 struct digamma_initializer
 558 {
 559    struct init
 560    {
 561       init()
 562       {
 563          typedef typename policies::precision<T, Policy>::type precision_type;
 564          do_init(std::integral_constant<bool, precision_type::value && (precision_type::value <= 113)>());
 565       }
 566       void do_init(const std::true_type&)
 567       {
 568          boost::math::digamma(T(1.5), Policy());
 569          boost::math::digamma(T(500), Policy());
 570       }
 571       void do_init(const std::false_type&){}
 572       void force_instantiate()const{}
 573    };
 574    static const init initializer;
 575    static void force_instantiate()
 576    {
 577       initializer.force_instantiate();
 578    }
 579 };
 580
 581 template <class T, class Policy>
 582 const typename digamma_initializer<T, Policy>::init digamma_initializer<T, Policy>::initializer;
 583
 584 } // namespace detail
 585
 586 template <class T, class Policy>
 587 inline typename tools::promote_args<T>::type
 588    digamma(T x, const Policy&)
 589 {
 590    typedef typename tools::promote_args<T>::type result_type;
 591    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 592    typedef typename policies::precision<T, Policy>::type precision_type;
 593    typedef std::integral_constant<int,
 594       (precision_type::value <= 0) || (precision_type::value > 113) ? 0 :
 595       precision_type::value <= 24 ? 24 :
 596       precision_type::value <= 53 ? 53 :
 597       precision_type::value <= 64 ? 64 :
 598       precision_type::value <= 113 ? 113 : 0 > tag_type;
 599    typedef typename policies::normalise<
 600       Policy,
 601       policies::promote_float<false>,
 602       policies::promote_double<false>,
 603       policies::discrete_quantile<>,
 604       policies::assert_undefined<> >::type forwarding_policy;
 605
 606    // Force initialization of constants:
 607    detail::digamma_initializer<value_type, forwarding_policy>::force_instantiate();
 608
 609    return policies::checked_narrowing_cast<result_type, Policy>(detail::digamma_imp(
 610       static_cast<value_type>(x),
 611       static_cast<const tag_type*>(0), forwarding_policy()), "boost::math::digamma<%1%>(%1%)");
 612 }
 613
 614 template <class T>
 615 inline typename tools::promote_args<T>::type
 616    digamma(T x)
 617 {
 618    return digamma(x, policies::policy<>());
 619 }
 620
 621 } // namespace math
 622 } // namespace boost
 623
 624 #ifdef _MSC_VER
 625 #pragma warning(pop)
 626 #endif
 627
 628 #endif
 629