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