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 dianostic ignored work :(
30 //
31 #pragma GCC system_header
32 #endif
33 
34 namespace boost{
35 namespace math{
36 namespace detail{
37 //
38 // Begin by defining the smallest value for which it is safe to
39 // use the asymptotic expansion for digamma:
40 //
digamma_large_lim(const mpl::int_<0> *)41 inline unsigned digamma_large_lim(const mpl::int_<0>*)
42 {  return 20;  }
digamma_large_lim(const mpl::int_<113> *)43 inline unsigned digamma_large_lim(const mpl::int_<113>*)
44 {  return 20;  }
digamma_large_lim(const void *)45 inline 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 //
57 template <class T>
digamma_imp_large(T x,const mpl::int_<113> *)58 inline T digamma_imp_large(T x, const mpl::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 //
90 template <class T>
digamma_imp_large(T x,const mpl::int_<64> *)91 inline T digamma_imp_large(T x, const mpl::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 //
117 template <class T>
digamma_imp_large(T x,const mpl::int_<53> *)118 inline T digamma_imp_large(T x, const mpl::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 //
141 template <class T>
digamma_imp_large(T x,const mpl::int_<24> *)142 inline T digamma_imp_large(T x, const mpl::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 //
161 template <class T>
162 struct digamma_series_func
163 {
164 private:
165    int k;
166    T xx;
167    T term;
168 public:
digamma_series_funcboost::math::detail::digamma_series_func169    digamma_series_func(T x) : k(1), xx(x * x), term(1 / (x * x)) {}
operator ()boost::math::detail::digamma_series_func170    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 
180 template <class T, class Policy>
digamma_imp_large(T x,const Policy & pol,const mpl::int_<0> *)181 inline T digamma_imp_large(T x, const Policy& pol, const mpl::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 //
197 template <class T>
digamma_imp_1_2(T x,const mpl::int_<113> *)198 T digamma_imp_1_2(T x, const mpl::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 //
259 template <class T>
digamma_imp_1_2(T x,const mpl::int_<64> *)260 T digamma_imp_1_2(T x, const mpl::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 //
309 template <class T>
digamma_imp_1_2(T x,const mpl::int_<53> *)310 T digamma_imp_1_2(T x, const mpl::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 //
358 template <class T>
digamma_imp_1_2(T x,const mpl::int_<24> *)359 inline T digamma_imp_1_2(T x, const mpl::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 
396 template <class T, class Tag, class Policy>
397 T 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 recurrance 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 
462 template <class T, class Policy>
463 T digamma_imp(T x, const mpl::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 //
557 template <class T, class Policy>
558 struct digamma_initializer
559 {
560    struct init
561    {
initboost::math::detail::digamma_initializer::init562       init()
563       {
564          typedef typename policies::precision<T, Policy>::type precision_type;
565          do_init(mpl::bool_<precision_type::value && (precision_type::value <= 113)>());
566       }
do_initboost::math::detail::digamma_initializer::init567       void do_init(const mpl::true_&)
568       {
569          boost::math::digamma(T(1.5), Policy());
570          boost::math::digamma(T(500), Policy());
571       }
do_initboost::math::detail::digamma_initializer::init572       void do_init(const mpl::false_&){}
force_instantiateboost::math::detail::digamma_initializer::init573       void force_instantiate()const{}
574    };
575    static const init initializer;
force_instantiateboost::math::detail::digamma_initializer576    static void force_instantiate()
577    {
578       initializer.force_instantiate();
579    }
580 };
581 
582 template <class T, class Policy>
583 const typename digamma_initializer<T, Policy>::init digamma_initializer<T, Policy>::initializer;
584 
585 } // namespace detail
586 
587 template <class T, class Policy>
588 inline typename tools::promote_args<T>::type
digamma(T x,const Policy &)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 typename mpl::if_<
595       mpl::or_<
596          mpl::less_equal<precision_type, mpl::int_<0> >,
597          mpl::greater<precision_type, mpl::int_<114> >
598       >,
599       mpl::int_<0>,
600       typename mpl::if_<
601          mpl::less<precision_type, mpl::int_<25> >,
602          mpl::int_<24>,
603          typename mpl::if_<
604             mpl::less<precision_type, mpl::int_<54> >,
605             mpl::int_<53>,
606             typename mpl::if_<
607                mpl::less<precision_type, mpl::int_<65> >,
608                mpl::int_<64>,
609                mpl::int_<113>
610             >::type
611          >::type
612       >::type
613    >::type tag_type;
614 
615    typedef typename policies::normalise<
616       Policy,
617       policies::promote_float<false>,
618       policies::promote_double<false>,
619       policies::discrete_quantile<>,
620       policies::assert_undefined<> >::type forwarding_policy;
621 
622    // Force initialization of constants:
623    detail::digamma_initializer<value_type, forwarding_policy>::force_instantiate();
624 
625    return policies::checked_narrowing_cast<result_type, Policy>(detail::digamma_imp(
626       static_cast<value_type>(x),
627       static_cast<const tag_type*>(0), forwarding_policy()), "boost::math::digamma<%1%>(%1%)");
628 }
629 
630 template <class T>
631 inline typename tools::promote_args<T>::type
digamma(T x)632    digamma(T x)
633 {
634    return digamma(x, policies::policy<>());
635 }
636 
637 } // namespace math
638 } // namespace boost
639 
640 #ifdef _MSC_VER
641 #pragma warning(pop)
642 #endif
643 
644 #endif
645 
646