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 #endif
  12
  13 #include <boost/math/tools/rational.hpp>
  14 #include <boost/math/tools/promotion.hpp>
  15 #include <boost/math/policies/error_handling.hpp>
  16 #include <boost/math/constants/constants.hpp>
  17 #include <boost/mpl/comparison.hpp>
  18 #include <boost/math/tools/big_constant.hpp>
  19
  20 namespace boost{
  21 namespace math{
  22 namespace detail{
  23 //
  24 // Begin by defining the smallest value for which it is safe to
  25 // use the asymptotic expansion for digamma:
  26 //
  27 inline unsigned digamma_large_lim(const mpl::int_<0>*)
  28 {  return 20;  }
  29
  30 inline unsigned digamma_large_lim(const void*)
  31 {  return 10;  }
  32 //
  33 // Implementations of the asymptotic expansion come next,
  34 // the coefficients of the series have been evaluated
  35 // in advance at high precision, and the series truncated
  36 // at the first term that's too small to effect the result.
  37 // Note that the series becomes divergent after a while
  38 // so truncation is very important.
  39 //
  40 // This first one gives 34-digit precision for x >= 20:
  41 //
  42 template <class T>
  43 inline T digamma_imp_large(T x, const mpl::int_<0>*)
  44 {
  45    BOOST_MATH_STD_USING // ADL of std functions.
  46    static const T P[] = {
  47       BOOST_MATH_BIG_CONSTANT(T, 113, 0.083333333333333333333333333333333333333333333333333),
  48       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0083333333333333333333333333333333333333333333333333),
  49       BOOST_MATH_BIG_CONSTANT(T, 113, 0.003968253968253968253968253968253968253968253968254),
  50       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0041666666666666666666666666666666666666666666666667),
  51       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0075757575757575757575757575757575757575757575757576),
  52       BOOST_MATH_BIG_CONSTANT(T, 113, -0.021092796092796092796092796092796092796092796092796),
  53       BOOST_MATH_BIG_CONSTANT(T, 113, 0.083333333333333333333333333333333333333333333333333),
  54       BOOST_MATH_BIG_CONSTANT(T, 113, -0.44325980392156862745098039215686274509803921568627),
  55       BOOST_MATH_BIG_CONSTANT(T, 113, 3.0539543302701197438039543302701197438039543302701),
  56       BOOST_MATH_BIG_CONSTANT(T, 113, -26.456212121212121212121212121212121212121212121212),
  57       BOOST_MATH_BIG_CONSTANT(T, 113, 281.4601449275362318840579710144927536231884057971),
  58       BOOST_MATH_BIG_CONSTANT(T, 113, -3607.510546398046398046398046398046398046398046398),
  59       BOOST_MATH_BIG_CONSTANT(T, 113, 54827.583333333333333333333333333333333333333333333),
  60       BOOST_MATH_BIG_CONSTANT(T, 113, -974936.82385057471264367816091954022988505747126437),
  61       BOOST_MATH_BIG_CONSTANT(T, 113, 20052695.796688078946143462272494530559046688078946),
  62       BOOST_MATH_BIG_CONSTANT(T, 113, -472384867.72162990196078431372549019607843137254902),
  63       BOOST_MATH_BIG_CONSTANT(T, 113, 12635724795.916666666666666666666666666666666666667)
  64    };
  65    x -= 1;
  66    T result = log(x);
  67    result += 1 / (2 * x);
  68    T z = 1 / (x*x);
  69    result -= z * tools::evaluate_polynomial(P, z);
  70    return result;
  71 }
  72 //
  73 // 19-digit precision for x >= 10:
  74 //
  75 template <class T>
  76 inline T digamma_imp_large(T x, const mpl::int_<64>*)
  77 {
  78    BOOST_MATH_STD_USING // ADL of std functions.
  79    static const T P[] = {
  80       BOOST_MATH_BIG_CONSTANT(T, 64, 0.083333333333333333333333333333333333333333333333333),
  81       BOOST_MATH_BIG_CONSTANT(T, 64, -0.0083333333333333333333333333333333333333333333333333),
  82       BOOST_MATH_BIG_CONSTANT(T, 64, 0.003968253968253968253968253968253968253968253968254),
  83       BOOST_MATH_BIG_CONSTANT(T, 64, -0.0041666666666666666666666666666666666666666666666667),
  84       BOOST_MATH_BIG_CONSTANT(T, 64, 0.0075757575757575757575757575757575757575757575757576),
  85       BOOST_MATH_BIG_CONSTANT(T, 64, -0.021092796092796092796092796092796092796092796092796),
  86       BOOST_MATH_BIG_CONSTANT(T, 64, 0.083333333333333333333333333333333333333333333333333),
  87       BOOST_MATH_BIG_CONSTANT(T, 64, -0.44325980392156862745098039215686274509803921568627),
  88       BOOST_MATH_BIG_CONSTANT(T, 64, 3.0539543302701197438039543302701197438039543302701),
  89       BOOST_MATH_BIG_CONSTANT(T, 64, -26.456212121212121212121212121212121212121212121212),
  90       BOOST_MATH_BIG_CONSTANT(T, 64, 281.4601449275362318840579710144927536231884057971),
  91    };
  92    x -= 1;
  93    T result = log(x);
  94    result += 1 / (2 * x);
  95    T z = 1 / (x*x);
  96    result -= z * tools::evaluate_polynomial(P, z);
  97    return result;
  98 }
  99 //
 100 // 17-digit precision for x >= 10:
 101 //
 102 template <class T>
 103 inline T digamma_imp_large(T x, const mpl::int_<53>*)
 104 {
 105    BOOST_MATH_STD_USING // ADL of std functions.
 106    static const T P[] = {
 107       BOOST_MATH_BIG_CONSTANT(T, 53, 0.083333333333333333333333333333333333333333333333333),
 108       BOOST_MATH_BIG_CONSTANT(T, 53, -0.0083333333333333333333333333333333333333333333333333),
 109       BOOST_MATH_BIG_CONSTANT(T, 53, 0.003968253968253968253968253968253968253968253968254),
 110       BOOST_MATH_BIG_CONSTANT(T, 53, -0.0041666666666666666666666666666666666666666666666667),
 111       BOOST_MATH_BIG_CONSTANT(T, 53, 0.0075757575757575757575757575757575757575757575757576),
 112       BOOST_MATH_BIG_CONSTANT(T, 53, -0.021092796092796092796092796092796092796092796092796),
 113       BOOST_MATH_BIG_CONSTANT(T, 53, 0.083333333333333333333333333333333333333333333333333),
 114       BOOST_MATH_BIG_CONSTANT(T, 53, -0.44325980392156862745098039215686274509803921568627)
 115    };
 116    x -= 1;
 117    T result = log(x);
 118    result += 1 / (2 * x);
 119    T z = 1 / (x*x);
 120    result -= z * tools::evaluate_polynomial(P, z);
 121    return result;
 122 }
 123 //
 124 // 9-digit precision for x >= 10:
 125 //
 126 template <class T>
 127 inline T digamma_imp_large(T x, const mpl::int_<24>*)
 128 {
 129    BOOST_MATH_STD_USING // ADL of std functions.
 130    static const T P[] = {
 131       BOOST_MATH_BIG_CONSTANT(T, 24, 0.083333333333333333333333333333333333333333333333333),
 132       BOOST_MATH_BIG_CONSTANT(T, 24, -0.0083333333333333333333333333333333333333333333333333),
 133       BOOST_MATH_BIG_CONSTANT(T, 24, 0.003968253968253968253968253968253968253968253968254)
 134    };
 135    x -= 1;
 136    T result = log(x);
 137    result += 1 / (2 * x);
 138    T z = 1 / (x*x);
 139    result -= z * tools::evaluate_polynomial(P, z);
 140    return result;
 141 }
 142 //
 143 // Now follow rational approximations over the range [1,2].
 144 //
 145 // 35-digit precision:
 146 //
 147 template <class T>
 148 T digamma_imp_1_2(T x, const mpl::int_<0>*)
 149 {
 150    //
 151    // Now the approximation, we use the form:
 152    //
 153    // digamma(x) = (x - root) * (Y + R(x-1))
 154    //
 155    // Where root is the location of the positive root of digamma,
 156    // Y is a constant, and R is optimised for low absolute error
 157    // compared to Y.
 158    //
 159    // Max error found at 128-bit long double precision:  5.541e-35
 160    // Maximum Deviation Found (approximation error):     1.965e-35
 161    //
 162    static const float Y = 0.99558162689208984375F;
 163
 164    static const T root1 = T(1569415565) / 1073741824uL;
 165    static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
 166    static const T root3 = ((T(111616537) / 1073741824uL) / 1073741824uL) / 1073741824uL;
 167    static const T root4 = (((T(503992070) / 1073741824uL) / 1073741824uL) / 1073741824uL) / 1073741824uL;
 168    static const T root5 = BOOST_MATH_BIG_CONSTANT(T, 113, 0.52112228569249997894452490385577338504019838794544e-36);
 169
 170    static const T P[] = {
 171       BOOST_MATH_BIG_CONSTANT(T, 113, 0.25479851061131551526977464225335883769),
 172       BOOST_MATH_BIG_CONSTANT(T, 113, -0.18684290534374944114622235683619897417),
 173       BOOST_MATH_BIG_CONSTANT(T, 113, -0.80360876047931768958995775910991929922),
 174       BOOST_MATH_BIG_CONSTANT(T, 113, -0.67227342794829064330498117008564270136),
 175       BOOST_MATH_BIG_CONSTANT(T, 113, -0.26569010991230617151285010695543858005),
 176       BOOST_MATH_BIG_CONSTANT(T, 113, -0.05775672694575986971640757748003553385),
 177       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0071432147823164975485922555833274240665),
 178       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00048740753910766168912364555706064993274),
 179       BOOST_MATH_BIG_CONSTANT(T, 113, -0.16454996865214115723416538844975174761e-4),
 180       BOOST_MATH_BIG_CONSTANT(T, 113, -0.20327832297631728077731148515093164955e-6)
 181    };
 182    static const T Q[] = {
 183       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
 184       BOOST_MATH_BIG_CONSTANT(T, 113, 2.6210924610812025425088411043163287646),
 185       BOOST_MATH_BIG_CONSTANT(T, 113, 2.6850757078559596612621337395886392594),
 186       BOOST_MATH_BIG_CONSTANT(T, 113, 1.4320913706209965531250495490639289418),
 187       BOOST_MATH_BIG_CONSTANT(T, 113, 0.4410872083455009362557012239501953402),
 188       BOOST_MATH_BIG_CONSTANT(T, 113, 0.081385727399251729505165509278152487225),
 189       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0089478633066857163432104815183858149496),
 190       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00055861622855066424871506755481997374154),
 191       BOOST_MATH_BIG_CONSTANT(T, 113, 0.1760168552357342401304462967950178554e-4),
 192       BOOST_MATH_BIG_CONSTANT(T, 113, 0.20585454493572473724556649516040874384e-6),
 193       BOOST_MATH_BIG_CONSTANT(T, 113, -0.90745971844439990284514121823069162795e-11),
 194       BOOST_MATH_BIG_CONSTANT(T, 113, 0.48857673606545846774761343500033283272e-13),
 195    };
 196    T g = x - root1;
 197    g -= root2;
 198    g -= root3;
 199    g -= root4;
 200    g -= root5;
 201    T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
 202    T result = g * Y + g * r;
 203
 204    return result;
 205 }
 206 //
 207 // 19-digit precision:
 208 //
 209 template <class T>
 210 T digamma_imp_1_2(T x, const mpl::int_<64>*)
 211 {
 212    //
 213    // Now the approximation, we use the form:
 214    //
 215    // digamma(x) = (x - root) * (Y + R(x-1))
 216    //
 217    // Where root is the location of the positive root of digamma,
 218    // Y is a constant, and R is optimised for low absolute error
 219    // compared to Y.
 220    //
 221    // Max error found at 80-bit long double precision:   5.016e-20
 222    // Maximum Deviation Found (approximation error):     3.575e-20
 223    //
 224    static const float Y = 0.99558162689208984375F;
 225
 226    static const T root1 = T(1569415565) / 1073741824uL;
 227    static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
 228    static const T root3 = BOOST_MATH_BIG_CONSTANT(T, 64, 0.9016312093258695918615325266959189453125e-19);
 229
 230    static const T P[] = {
 231       BOOST_MATH_BIG_CONSTANT(T, 64, 0.254798510611315515235),
 232       BOOST_MATH_BIG_CONSTANT(T, 64, -0.314628554532916496608),
 233       BOOST_MATH_BIG_CONSTANT(T, 64, -0.665836341559876230295),
 234       BOOST_MATH_BIG_CONSTANT(T, 64, -0.314767657147375752913),
 235       BOOST_MATH_BIG_CONSTANT(T, 64, -0.0541156266153505273939),
 236       BOOST_MATH_BIG_CONSTANT(T, 64, -0.00289268368333918761452)
 237    };
 238    static const T Q[] = {
 239       BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
 240       BOOST_MATH_BIG_CONSTANT(T, 64, 2.1195759927055347547),
 241       BOOST_MATH_BIG_CONSTANT(T, 64, 1.54350554664961128724),
 242       BOOST_MATH_BIG_CONSTANT(T, 64, 0.486986018231042975162),
 243       BOOST_MATH_BIG_CONSTANT(T, 64, 0.0660481487173569812846),
 244       BOOST_MATH_BIG_CONSTANT(T, 64, 0.00298999662592323990972),
 245       BOOST_MATH_BIG_CONSTANT(T, 64, -0.165079794012604905639e-5),
 246       BOOST_MATH_BIG_CONSTANT(T, 64, 0.317940243105952177571e-7)
 247    };
 248    T g = x - root1;
 249    g -= root2;
 250    g -= root3;
 251    T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
 252    T result = g * Y + g * r;
 253
 254    return result;
 255 }
 256 //
 257 // 18-digit precision:
 258 //
 259 template <class T>
 260 T digamma_imp_1_2(T x, const mpl::int_<53>*)
 261 {
 262    //
 263    // Now the approximation, we use the form:
 264    //
 265    // digamma(x) = (x - root) * (Y + R(x-1))
 266    //
 267    // Where root is the location of the positive root of digamma,
 268    // Y is a constant, and R is optimised for low absolute error
 269    // compared to Y.
 270    //
 271    // Maximum Deviation Found:               1.466e-18
 272    // At double precision, max error found:  2.452e-17
 273    //
 274    static const float Y = 0.99558162689208984F;
 275
 276    static const T root1 = T(1569415565) / 1073741824uL;
 277    static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
 278    static const T root3 = BOOST_MATH_BIG_CONSTANT(T, 53, 0.9016312093258695918615325266959189453125e-19);
 279
 280    static const T P[] = {
 281       BOOST_MATH_BIG_CONSTANT(T, 53, 0.25479851061131551),
 282       BOOST_MATH_BIG_CONSTANT(T, 53, -0.32555031186804491),
 283       BOOST_MATH_BIG_CONSTANT(T, 53, -0.65031853770896507),
 284       BOOST_MATH_BIG_CONSTANT(T, 53, -0.28919126444774784),
 285       BOOST_MATH_BIG_CONSTANT(T, 53, -0.045251321448739056),
 286       BOOST_MATH_BIG_CONSTANT(T, 53, -0.0020713321167745952)
 287    };
 288    static const T Q[] = {
 289       BOOST_MATH_BIG_CONSTANT(T, 53, 1),
 290       BOOST_MATH_BIG_CONSTANT(T, 53, 2.0767117023730469),
 291       BOOST_MATH_BIG_CONSTANT(T, 53, 1.4606242909763515),
 292       BOOST_MATH_BIG_CONSTANT(T, 53, 0.43593529692665969),
 293       BOOST_MATH_BIG_CONSTANT(T, 53, 0.054151797245674225),
 294       BOOST_MATH_BIG_CONSTANT(T, 53, 0.0021284987017821144),
 295       BOOST_MATH_BIG_CONSTANT(T, 53, -0.55789841321675513e-6)
 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 // 9-digit precision:
 307 //
 308 template <class T>
 309 inline T digamma_imp_1_2(T x, const mpl::int_<24>*)
 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:              3.388e-010
 321    // At float precision, max error found:  2.008725e-008
 322    //
 323    static const float Y = 0.99558162689208984f;
 324    static const T root = 1532632.0f / 1048576;
 325    static const T root_minor = static_cast<T>(0.3700660185912626595423257213284682051735604e-6L);
 326    static const T P[] = {
 327       0.25479851023250261e0,
 328       -0.44981331915268368e0,
 329       -0.43916936919946835e0,
 330       -0.61041765350579073e-1
 331    };
 332    static const T Q[] = {
 333       0.1e1,
 334       0.15890202430554952e1,
 335       0.65341249856146947e0,
 336       0.63851690523355715e-1
 337    };
 338    T g = x - root;
 339    g -= root_minor;
 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 template <class T, class Tag, class Policy>
 347 T digamma_imp(T x, const Tag* t, const Policy& pol)
 348 {
 349    //
 350    // This handles reflection of negative arguments, and all our
 351    // error handling, then forwards to the T-specific approximation.
 352    //
 353    BOOST_MATH_STD_USING // ADL of std functions.
 354
 355    T result = 0;
 356    //
 357    // Check for negative arguments and use reflection:
 358    //
 359    if(x < 0)
 360    {
 361       // Reflect:
 362       x = 1 - x;
 363       // Argument reduction for tan:
 364       T remainder = x - floor(x);
 365       // Shift to negative if > 0.5:
 366       if(remainder > 0.5)
 367       {
 368          remainder -= 1;
 369       }
 370       //
 371       // check for evaluation at a negative pole:
 372       //
 373       if(remainder == 0)
 374       {
 375          return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, (1-x), pol);
 376       }
 377       result = constants::pi<T>() / tan(constants::pi<T>() * remainder);
 378    }
 379    //
 380    // If we're above the lower-limit for the
 381    // asymptotic expansion then use it:
 382    //
 383    if(x >= digamma_large_lim(t))
 384    {
 385       result += digamma_imp_large(x, t);
 386    }
 387    else
 388    {
 389       //
 390       // If x > 2 reduce to the interval [1,2]:
 391       //
 392       while(x > 2)
 393       {
 394          x -= 1;
 395          result += 1/x;
 396       }
 397       //
 398       // If x < 1 use recurrance to shift to > 1:
 399       //
 400       if(x < 1)
 401       {
 402          result = -1/x;
 403          x += 1;
 404       }
 405       result += digamma_imp_1_2(x, t);
 406    }
 407    return result;
 408 }
 409
 410 //
 411 // Initializer: ensure all our constants are initialized prior to the first call of main:
 412 //
 413 template <class T, class Policy>
 414 struct digamma_initializer
 415 {
 416    struct init
 417    {
 418       init()
 419       {
 420          boost::math::digamma(T(1.5), Policy());
 421          boost::math::digamma(T(500), Policy());
 422       }
 423       void force_instantiate()const{}
 424    };
 425    static const init initializer;
 426    static void force_instantiate()
 427    {
 428       initializer.force_instantiate();
 429    }
 430 };
 431
 432 template <class T, class Policy>
 433 const typename digamma_initializer<T, Policy>::init digamma_initializer<T, Policy>::initializer;
 434
 435 } // namespace detail
 436
 437 template <class T, class Policy>
 438 inline typename tools::promote_args<T>::type
 439    digamma(T x, const Policy& pol)
 440 {
 441    typedef typename tools::promote_args<T>::type result_type;
 442    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 443    typedef typename policies::precision<T, Policy>::type precision_type;
 444    typedef typename mpl::if_<
 445       mpl::or_<
 446          mpl::less_equal<precision_type, mpl::int_<0> >,
 447          mpl::greater<precision_type, mpl::int_<64> >
 448       >,
 449       mpl::int_<0>,
 450       typename mpl::if_<
 451          mpl::less<precision_type, mpl::int_<25> >,
 452          mpl::int_<24>,
 453          typename mpl::if_<
 454             mpl::less<precision_type, mpl::int_<54> >,
 455             mpl::int_<53>,
 456             mpl::int_<64>
 457          >::type
 458       >::type
 459    >::type tag_type;
 460
 461    // Force initialization of constants:
 462    detail::digamma_initializer<result_type, Policy>::force_instantiate();
 463
 464    return policies::checked_narrowing_cast<result_type, Policy>(detail::digamma_imp(
 465       static_cast<value_type>(x),
 466       static_cast<const tag_type*>(0), pol), "boost::math::digamma<%1%>(%1%)");
 467 }
 468
 469 template <class T>
 470 inline typename tools::promote_args<T>::type
 471    digamma(T x)
 472 {
 473    return digamma(x, policies::policy<>());
 474 }
 475
 476 } // namespace math
 477 } // namespace boost
 478 #endif
 479