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
33namespace boost{
34namespace math{
35namespace detail{
36//
37// Begin by defining the smallest value for which it is safe to
38// use the asymptotic expansion for digamma:
39//
40inline unsigned digamma_large_lim(const std::integral_constant<int, 0>*)
41{ return 20; }
42inline unsigned digamma_large_lim(const std::integral_constant<int, 113>*)
43{ return 20; }
44inline 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//
56template <class T>
57inline 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//
89template <class T>
90inline 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//
116template <class T>
117inline 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 0.083333333333333333333333333333333333333333333333333,
122 -0.0083333333333333333333333333333333333333333333333333,
123 0.003968253968253968253968253968253968253968253968254,
124 -0.0041666666666666666666666666666666666666666666666667,
125 0.0075757575757575757575757575757575757575757575757576,
126 -0.021092796092796092796092796092796092796092796092796,
127 0.083333333333333333333333333333333333333333333333333,
128 -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//
140template <class T>
141inline 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 0.083333333333333333333333333333333333333333333333333f,
146 -0.0083333333333333333333333333333333333333333333333333f,
147 0.003968253968253968253968253968253968253968253968254f
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//
160template <class T>
161struct digamma_series_func
162{
163private:
164 int k;
165 T xx;
166 T term;
167public:
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
179template <class T, class Policy>
180inline 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//
196template <class T>
197T 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 // LCOV_EXCL_START
212 static const float Y = 0.99558162689208984375F;
213
214 static const T root1 = T(1569415565) / 1073741824uL;
215 static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
216 static const T root3 = ((T(111616537) / 1073741824uL) / 1073741824uL) / 1073741824uL;
217 static const T root4 = (((T(503992070) / 1073741824uL) / 1073741824uL) / 1073741824uL) / 1073741824uL;
218 static const T root5 = BOOST_MATH_BIG_CONSTANT(T, 113, 0.52112228569249997894452490385577338504019838794544e-36);
219
220 static const T P[] = {
221 BOOST_MATH_BIG_CONSTANT(T, 113, 0.25479851061131551526977464225335883769),
222 BOOST_MATH_BIG_CONSTANT(T, 113, -0.18684290534374944114622235683619897417),
223 BOOST_MATH_BIG_CONSTANT(T, 113, -0.80360876047931768958995775910991929922),
224 BOOST_MATH_BIG_CONSTANT(T, 113, -0.67227342794829064330498117008564270136),
225 BOOST_MATH_BIG_CONSTANT(T, 113, -0.26569010991230617151285010695543858005),
226 BOOST_MATH_BIG_CONSTANT(T, 113, -0.05775672694575986971640757748003553385),
227 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0071432147823164975485922555833274240665),
228 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00048740753910766168912364555706064993274),
229 BOOST_MATH_BIG_CONSTANT(T, 113, -0.16454996865214115723416538844975174761e-4),
230 BOOST_MATH_BIG_CONSTANT(T, 113, -0.20327832297631728077731148515093164955e-6)
231 };
232 static const T Q[] = {
233 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
234 BOOST_MATH_BIG_CONSTANT(T, 113, 2.6210924610812025425088411043163287646),
235 BOOST_MATH_BIG_CONSTANT(T, 113, 2.6850757078559596612621337395886392594),
236 BOOST_MATH_BIG_CONSTANT(T, 113, 1.4320913706209965531250495490639289418),
237 BOOST_MATH_BIG_CONSTANT(T, 113, 0.4410872083455009362557012239501953402),
238 BOOST_MATH_BIG_CONSTANT(T, 113, 0.081385727399251729505165509278152487225),
239 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0089478633066857163432104815183858149496),
240 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00055861622855066424871506755481997374154),
241 BOOST_MATH_BIG_CONSTANT(T, 113, 0.1760168552357342401304462967950178554e-4),
242 BOOST_MATH_BIG_CONSTANT(T, 113, 0.20585454493572473724556649516040874384e-6),
243 BOOST_MATH_BIG_CONSTANT(T, 113, -0.90745971844439990284514121823069162795e-11),
244 BOOST_MATH_BIG_CONSTANT(T, 113, 0.48857673606545846774761343500033283272e-13),
245 };
246 // LCOV_EXCL_STOP
247 T g = x - root1;
248 g -= root2;
249 g -= root3;
250 g -= root4;
251 g -= root5;
252 T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
253 T result = g * Y + g * r;
254
255 return result;
256}
257//
258// 19-digit precision:
259//
260template <class T>
261T digamma_imp_1_2(T x, const std::integral_constant<int, 64>*)
262{
263 //
264 // Now the approximation, we use the form:
265 //
266 // digamma(x) = (x - root) * (Y + R(x-1))
267 //
268 // Where root is the location of the positive root of digamma,
269 // Y is a constant, and R is optimised for low absolute error
270 // compared to Y.
271 //
272 // Max error found at 80-bit long double precision: 5.016e-20
273 // Maximum Deviation Found (approximation error): 3.575e-20
274 //
275 // LCOV_EXCL_START
276 static const float Y = 0.99558162689208984375F;
277
278 static const T root1 = T(1569415565) / 1073741824uL;
279 static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
280 static const T root3 = BOOST_MATH_BIG_CONSTANT(T, 64, 0.9016312093258695918615325266959189453125e-19);
281
282 static const T P[] = {
283 BOOST_MATH_BIG_CONSTANT(T, 64, 0.254798510611315515235),
284 BOOST_MATH_BIG_CONSTANT(T, 64, -0.314628554532916496608),
285 BOOST_MATH_BIG_CONSTANT(T, 64, -0.665836341559876230295),
286 BOOST_MATH_BIG_CONSTANT(T, 64, -0.314767657147375752913),
287 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0541156266153505273939),
288 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00289268368333918761452)
289 };
290 static const T Q[] = {
291 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
292 BOOST_MATH_BIG_CONSTANT(T, 64, 2.1195759927055347547),
293 BOOST_MATH_BIG_CONSTANT(T, 64, 1.54350554664961128724),
294 BOOST_MATH_BIG_CONSTANT(T, 64, 0.486986018231042975162),
295 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0660481487173569812846),
296 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00298999662592323990972),
297 BOOST_MATH_BIG_CONSTANT(T, 64, -0.165079794012604905639e-5),
298 BOOST_MATH_BIG_CONSTANT(T, 64, 0.317940243105952177571e-7)
299 };
300 // LCOV_EXCL_STOP
301 T g = x - root1;
302 g -= root2;
303 g -= root3;
304 T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
305 T result = g * Y + g * r;
306
307 return result;
308}
309//
310// 18-digit precision:
311//
312template <class T>
313T digamma_imp_1_2(T x, const std::integral_constant<int, 53>*)
314{
315 //
316 // Now the approximation, we use the form:
317 //
318 // digamma(x) = (x - root) * (Y + R(x-1))
319 //
320 // Where root is the location of the positive root of digamma,
321 // Y is a constant, and R is optimised for low absolute error
322 // compared to Y.
323 //
324 // Maximum Deviation Found: 1.466e-18
325 // At double precision, max error found: 2.452e-17
326 //
327 // LCOV_EXCL_START
328 static const float Y = 0.99558162689208984F;
329
330 static const T root1 = T(1569415565) / 1073741824uL;
331 static const T root2 = (T(381566830) / 1073741824uL) / 1073741824uL;
332 static const T root3 = BOOST_MATH_BIG_CONSTANT(T, 53, 0.9016312093258695918615325266959189453125e-19);
333
334 static const T P[] = {
335 BOOST_MATH_BIG_CONSTANT(T, 53, 0.25479851061131551),
336 BOOST_MATH_BIG_CONSTANT(T, 53, -0.32555031186804491),
337 BOOST_MATH_BIG_CONSTANT(T, 53, -0.65031853770896507),
338 BOOST_MATH_BIG_CONSTANT(T, 53, -0.28919126444774784),
339 BOOST_MATH_BIG_CONSTANT(T, 53, -0.045251321448739056),
340 BOOST_MATH_BIG_CONSTANT(T, 53, -0.0020713321167745952)
341 };
342 static const T Q[] = {
343 BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
344 BOOST_MATH_BIG_CONSTANT(T, 53, 2.0767117023730469),
345 BOOST_MATH_BIG_CONSTANT(T, 53, 1.4606242909763515),
346 BOOST_MATH_BIG_CONSTANT(T, 53, 0.43593529692665969),
347 BOOST_MATH_BIG_CONSTANT(T, 53, 0.054151797245674225),
348 BOOST_MATH_BIG_CONSTANT(T, 53, 0.0021284987017821144),
349 BOOST_MATH_BIG_CONSTANT(T, 53, -0.55789841321675513e-6)
350 };
351 // LCOV_EXCL_STOP
352 T g = x - root1;
353 g -= root2;
354 g -= root3;
355 T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
356 T result = g * Y + g * r;
357
358 return result;
359}
360//
361// 9-digit precision:
362//
363template <class T>
364inline T digamma_imp_1_2(T x, const std::integral_constant<int, 24>*)
365{
366 //
367 // Now the approximation, we use the form:
368 //
369 // digamma(x) = (x - root) * (Y + R(x-1))
370 //
371 // Where root is the location of the positive root of digamma,
372 // Y is a constant, and R is optimised for low absolute error
373 // compared to Y.
374 //
375 // Maximum Deviation Found: 3.388e-010
376 // At float precision, max error found: 2.008725e-008
377 //
378 // LCOV_EXCL_START
379 static const float Y = 0.99558162689208984f;
380 static const T root = 1532632.0f / 1048576;
381 static const T root_minor = static_cast<T>(0.3700660185912626595423257213284682051735604e-6L);
382 static const T P[] = {
383 0.25479851023250261e0f,
384 -0.44981331915268368e0f,
385 -0.43916936919946835e0f,
386 -0.61041765350579073e-1f
387 };
388 static const T Q[] = {
389 0.1e1f,
390 0.15890202430554952e1f,
391 0.65341249856146947e0f,
392 0.63851690523355715e-1f
393 };
394 // LCOV_EXCL_STOP
395 T g = x - root;
396 g -= root_minor;
397 T r = tools::evaluate_polynomial(P, T(x-1)) / tools::evaluate_polynomial(Q, T(x-1));
398 T result = g * Y + g * r;
399
400 return result;
401}
402
403template <class T, class Tag, class Policy>
404T digamma_imp(T x, const Tag* t, const Policy& pol)
405{
406 //
407 // This handles reflection of negative arguments, and all our
408 // error handling, then forwards to the T-specific approximation.
409 //
410 BOOST_MATH_STD_USING // ADL of std functions.
411
412 T result = 0;
413 //
414 // Check for negative arguments and use reflection:
415 //
416 if(x <= -1)
417 {
418 // Reflect:
419 x = 1 - x;
420 // Argument reduction for tan:
421 T remainder = x - floor(x);
422 // Shift to negative if > 0.5:
423 if(remainder > T(0.5))
424 {
425 remainder -= 1;
426 }
427 //
428 // check for evaluation at a negative pole:
429 //
430 if(remainder == 0)
431 {
432 return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", nullptr, (1-x), pol);
433 }
434 result = constants::pi<T>() / tan(constants::pi<T>() * remainder);
435 }
436 if(x == 0)
437 return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", nullptr, x, pol);
438 //
439 // If we're above the lower-limit for the
440 // asymptotic expansion then use it:
441 //
442 if(x >= digamma_large_lim(t))
443 {
444 result += digamma_imp_large(x, t);
445 }
446 else
447 {
448 //
449 // If x > 2 reduce to the interval [1,2]:
450 //
451 while(x > 2)
452 {
453 x -= 1;
454 result += 1/x;
455 }
456 //
457 // If x < 1 use recurrence to shift to > 1:
458 //
459 while(x < 1)
460 {
461 result -= 1/x;
462 x += 1;
463 }
464 result += digamma_imp_1_2(x, t);
465 }
466 return result;
467}
468
469template <class T, class Policy>
470T digamma_imp(T x, const std::integral_constant<int, 0>* t, const Policy& pol)
471{
472 //
473 // This handles reflection of negative arguments, and all our
474 // error handling, then forwards to the T-specific approximation.
475 //
476 // This is covered by our real_concept tests, but these are disabled for
477 // code coverage runs for performance reasons.
478 // LCOV_EXCL_START
479 //
480 BOOST_MATH_STD_USING // ADL of std functions.
481
482 T result = 0;
483 //
484 // Check for negative arguments and use reflection:
485 //
486 if(x <= -1)
487 {
488 // Reflect:
489 x = 1 - x;
490 // Argument reduction for tan:
491 T remainder = x - floor(x);
492 // Shift to negative if > 0.5:
493 if(remainder > T(0.5))
494 {
495 remainder -= 1;
496 }
497 //
498 // check for evaluation at a negative pole:
499 //
500 if(remainder == 0)
501 {
502 return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", nullptr, (1 - x), pol);
503 }
504 result = constants::pi<T>() / tan(constants::pi<T>() * remainder);
505 }
506 if(x == 0)
507 return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", nullptr, x, pol);
508 //
509 // If we're above the lower-limit for the
510 // asymptotic expansion then use it, the
511 // limit is a linear interpolation with
512 // limit = 10 at 50 bit precision and
513 // limit = 250 at 1000 bit precision.
514 //
515 int lim = 10 + ((tools::digits<T>() - 50) * 240L) / 950;
516 T two_x = ldexp(x, 1);
517 if(x >= lim)
518 {
519 result += digamma_imp_large(x, pol, t);
520 }
521 else if(floor(x) == x)
522 {
523 //
524 // Special case for integer arguments, see
525 // http://functions.wolfram.com/06.14.03.0001.01
526 //
527 result = -constants::euler<T, Policy>();
528 T val = 1;
529 while(val < x)
530 {
531 result += 1 / val;
532 val += 1;
533 }
534 }
535 else if(floor(two_x) == two_x)
536 {
537 //
538 // Special case for half integer arguments, see:
539 // http://functions.wolfram.com/06.14.03.0007.01
540 //
541 result = -2 * constants::ln_two<T, Policy>() - constants::euler<T, Policy>();
542 int n = itrunc(x);
543 if(n)
544 {
545 for(int k = 1; k < n; ++k)
546 result += 1 / T(k);
547 for(int k = n; k <= 2 * n - 1; ++k)
548 result += 2 / T(k);
549 }
550 }
551 else
552 {
553 //
554 // Rescale so we can use the asymptotic expansion:
555 //
556 while(x < lim)
557 {
558 result -= 1 / x;
559 x += 1;
560 }
561 result += digamma_imp_large(x, pol, t);
562 }
563 return result;
564 // LCOV_EXCL_STOP
565}
566
567} // namespace detail
568
569template <class T, class Policy>
570inline typename tools::promote_args<T>::type
571 digamma(T x, const Policy&)
572{
573 typedef typename tools::promote_args<T>::type result_type;
574 typedef typename policies::evaluation<result_type, Policy>::type value_type;
575 typedef typename policies::precision<T, Policy>::type precision_type;
576 typedef std::integral_constant<int,
577 (precision_type::value <= 0) || (precision_type::value > 113) ? 0 :
578 precision_type::value <= 24 ? 24 :
579 precision_type::value <= 53 ? 53 :
580 precision_type::value <= 64 ? 64 :
581 precision_type::value <= 113 ? 113 : 0 > tag_type;
582 typedef typename policies::normalise<
583 Policy,
584 policies::promote_float<false>,
585 policies::promote_double<false>,
586 policies::discrete_quantile<>,
587 policies::assert_undefined<> >::type forwarding_policy;
588
589 return policies::checked_narrowing_cast<result_type, Policy>(detail::digamma_imp(
590 static_cast<value_type>(x),
591 static_cast<const tag_type*>(nullptr), forwarding_policy()), "boost::math::digamma<%1%>(%1%)");
592}
593
594template <class T>
595inline typename tools::promote_args<T>::type
596 digamma(T x)
597{
598 return digamma(x, policies::policy<>());
599}
600
601} // namespace math
602} // namespace boost
603
604#ifdef _MSC_VER
605#pragma warning(pop)
606#endif
607
608#endif
609
610

source code of boost/libs/math/include/boost/math/special_functions/digamma.hpp