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