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

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