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