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 | |
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 | // |
41 | inline unsigned digamma_large_lim(const boost::integral_constant<int, 0>*) |
42 | { return 20; } |
43 | inline unsigned digamma_large_lim(const boost::integral_constant<int, 113>*) |
44 | { return 20; } |
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> |
58 | inline 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 | // |
90 | template <class T> |
91 | inline 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 | // |
117 | template <class T> |
118 | inline 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 | // |
141 | template <class T> |
142 | inline 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 | // |
161 | template <class T> |
162 | struct digamma_series_func |
163 | { |
164 | private: |
165 | int k; |
166 | T xx; |
167 | T term; |
168 | public: |
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 | |
180 | template <class T, class Policy> |
181 | inline 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 | // |
197 | template <class T> |
198 | T 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 | // |
259 | template <class T> |
260 | T 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 | // |
309 | template <class T> |
310 | T 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 | // |
358 | template <class T> |
359 | inline 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 | |
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 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 | |
462 | template <class T, class Policy> |
463 | T 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 | // |
557 | template <class T, class Policy> |
558 | struct 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 | |
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 |
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 | |
615 | template <class T> |
616 | inline 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 | |