]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | ||
92f5a8d4 TL |
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__ | |
f67539c2 | 29 | // nor #pragma diagnostic ignored work :( |
92f5a8d4 TL |
30 | // |
31 | #pragma GCC system_header | |
32 | #endif | |
33 | ||
7c673cae FG |
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 | // | |
f67539c2 | 41 | inline unsigned digamma_large_lim(const boost::integral_constant<int, 0>*) |
7c673cae | 42 | { return 20; } |
f67539c2 | 43 | inline unsigned digamma_large_lim(const boost::integral_constant<int, 113>*) |
7c673cae FG |
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> | |
f67539c2 | 58 | inline T digamma_imp_large(T x, const boost::integral_constant<int, 113>*) |
7c673cae FG |
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> | |
f67539c2 | 91 | inline T digamma_imp_large(T x, const boost::integral_constant<int, 64>*) |
7c673cae FG |
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> | |
f67539c2 | 118 | inline T digamma_imp_large(T x, const boost::integral_constant<int, 53>*) |
7c673cae FG |
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> | |
f67539c2 | 142 | inline T digamma_imp_large(T x, const boost::integral_constant<int, 24>*) |
7c673cae FG |
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> | |
f67539c2 | 181 | inline T digamma_imp_large(T x, const Policy& pol, const boost::integral_constant<int, 0>*) |
7c673cae FG |
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> | |
f67539c2 | 198 | T digamma_imp_1_2(T x, const boost::integral_constant<int, 113>*) |
7c673cae FG |
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> | |
f67539c2 | 260 | T digamma_imp_1_2(T x, const boost::integral_constant<int, 64>*) |
7c673cae FG |
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> | |
f67539c2 | 310 | T digamma_imp_1_2(T x, const boost::integral_constant<int, 53>*) |
7c673cae FG |
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> | |
f67539c2 | 359 | inline T digamma_imp_1_2(T x, const boost::integral_constant<int, 24>*) |
7c673cae FG |
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 | // | |
f67539c2 | 450 | // If x < 1 use recurrence to shift to > 1: |
7c673cae FG |
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> | |
f67539c2 | 463 | T digamma_imp(T x, const boost::integral_constant<int, 0>* t, const Policy& pol) |
7c673cae FG |
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; | |
f67539c2 | 565 | do_init(boost::integral_constant<bool, precision_type::value && (precision_type::value <= 113)>()); |
7c673cae | 566 | } |
f67539c2 | 567 | void do_init(const boost::true_type&) |
7c673cae FG |
568 | { |
569 | boost::math::digamma(T(1.5), Policy()); | |
570 | boost::math::digamma(T(500), Policy()); | |
571 | } | |
f67539c2 | 572 | void do_init(const false_type&){} |
7c673cae FG |
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; | |
f67539c2 TL |
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; | |
7c673cae FG |
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 |