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