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), unchecked_bernoulli_imp(std::size_t n,const mpl::int_<0> &)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. unchecked_bernoulli_imp(std::size_t n,const mpl::int_<1> &)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; unchecked_bernoulli_imp(std::size_t n,const mpl::int_<2> &)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); unchecked_bernoulli_imp(std::size_t n,const mpl::int_<3> &)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 636