]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | |
2 | /////////////////////////////////////////////////////////////////////////////// | |
3 | // Copyright 2013 Nikhar Agrawal | |
4 | // Copyright 2013 Christopher Kormanyos | |
5 | // Copyright 2014 John Maddock | |
6 | // Copyright 2013 Paul Bristow | |
7 | // Distributed under the Boost | |
8 | // Software License, Version 1.0. (See accompanying file | |
9 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
10 | ||
11 | #ifndef _BOOST_POLYGAMMA_DETAIL_2013_07_30_HPP_ | |
12 | #define _BOOST_POLYGAMMA_DETAIL_2013_07_30_HPP_ | |
13 | ||
14 | #include <cmath> | |
15 | #include <limits> | |
16 | #include <boost/cstdint.hpp> | |
17 | #include <boost/math/policies/policy.hpp> | |
18 | #include <boost/math/special_functions/bernoulli.hpp> | |
19 | #include <boost/math/special_functions/trunc.hpp> | |
20 | #include <boost/math/special_functions/zeta.hpp> | |
21 | #include <boost/math/special_functions/digamma.hpp> | |
22 | #include <boost/math/special_functions/sin_pi.hpp> | |
23 | #include <boost/math/special_functions/cos_pi.hpp> | |
24 | #include <boost/math/special_functions/pow.hpp> | |
25 | #include <boost/mpl/if.hpp> | |
26 | #include <boost/mpl/int.hpp> | |
27 | #include <boost/static_assert.hpp> | |
28 | #include <boost/type_traits/is_convertible.hpp> | |
29 | ||
30 | #ifdef _MSC_VER | |
31 | #pragma once | |
32 | #pragma warning(push) | |
33 | #pragma warning(disable:4702) // Unreachable code (release mode only warning) | |
34 | #endif | |
35 | ||
36 | namespace boost { namespace math { namespace detail{ | |
37 | ||
38 | template<class T, class Policy> | |
39 | T polygamma_atinfinityplus(const int n, const T& x, const Policy& pol, const char* function) // for large values of x such as for x> 400 | |
40 | { | |
41 | // See http://functions.wolfram.com/GammaBetaErf/PolyGamma2/06/02/0001/ | |
42 | BOOST_MATH_STD_USING | |
43 | // | |
44 | // sum == current value of accumulated sum. | |
45 | // term == value of current term to be added to sum. | |
46 | // part_term == value of current term excluding the Bernoulli number part | |
47 | // | |
48 | if(n + x == x) | |
49 | { | |
50 | // x is crazy large, just concentrate on the first part of the expression and use logs: | |
51 | if(n == 1) return 1 / x; | |
52 | T nlx = n * log(x); | |
53 | if((nlx < tools::log_max_value<T>()) && (n < (int)max_factorial<T>::value)) | |
54 | return ((n & 1) ? 1 : -1) * boost::math::factorial<T>(n - 1) * pow(x, -n); | |
55 | else | |
56 | return ((n & 1) ? 1 : -1) * exp(boost::math::lgamma(T(n), pol) - n * log(x)); | |
57 | } | |
58 | T term, sum, part_term; | |
59 | T x_squared = x * x; | |
60 | // | |
61 | // Start by setting part_term to: | |
62 | // | |
63 | // (n-1)! / x^(n+1) | |
64 | // | |
65 | // which is common to both the first term of the series (with k = 1) | |
66 | // and to the leading part. | |
67 | // We can then get to the leading term by: | |
68 | // | |
69 | // part_term * (n + 2 * x) / 2 | |
70 | // | |
71 | // and to the first term in the series | |
72 | // (excluding the Bernoulli number) by: | |
73 | // | |
74 | // part_term n * (n + 1) / (2x) | |
75 | // | |
76 | // If either the factorial would overflow, | |
77 | // or the power term underflows, this just gets set to 0 and then we | |
78 | // know that we have to use logs for the initial terms: | |
79 | // | |
80 | part_term = ((n > (int)boost::math::max_factorial<T>::value) && (T(n) * n > tools::log_max_value<T>())) | |
81 | ? T(0) : static_cast<T>(boost::math::factorial<T>(n - 1, pol) * pow(x, -n - 1)); | |
82 | if(part_term == 0) | |
83 | { | |
84 | // Either n is very large, or the power term underflows, | |
85 | // set the initial values of part_term, term and sum via logs: | |
86 | part_term = static_cast<T>(boost::math::lgamma(n, pol) - (n + 1) * log(x)); | |
87 | sum = exp(part_term + log(n + 2 * x) - boost::math::constants::ln_two<T>()); | |
88 | part_term += log(T(n) * (n + 1)) - boost::math::constants::ln_two<T>() - log(x); | |
89 | part_term = exp(part_term); | |
90 | } | |
91 | else | |
92 | { | |
93 | sum = part_term * (n + 2 * x) / 2; | |
94 | part_term *= (T(n) * (n + 1)) / 2; | |
95 | part_term /= x; | |
96 | } | |
97 | // | |
98 | // If the leading term is 0, so is the result: | |
99 | // | |
100 | if(sum == 0) | |
101 | return sum; | |
102 | ||
103 | for(unsigned k = 1;;) | |
104 | { | |
105 | term = part_term * boost::math::bernoulli_b2n<T>(k, pol); | |
106 | sum += term; | |
107 | // | |
108 | // Normal termination condition: | |
109 | // | |
110 | if(fabs(term / sum) < tools::epsilon<T>()) | |
111 | break; | |
112 | // | |
113 | // Increment our counter, and move part_term on to the next value: | |
114 | // | |
115 | ++k; | |
116 | part_term *= T(n + 2 * k - 2) * (n - 1 + 2 * k); | |
117 | part_term /= (2 * k - 1) * 2 * k; | |
118 | part_term /= x_squared; | |
119 | // | |
120 | // Emergency get out termination condition: | |
121 | // | |
122 | if(k > policies::get_max_series_iterations<Policy>()) | |
123 | { | |
124 | return policies::raise_evaluation_error(function, "Series did not converge, closest value was %1%", sum, pol); | |
125 | } | |
126 | } | |
127 | ||
128 | if((n - 1) & 1) | |
129 | sum = -sum; | |
130 | ||
131 | return sum; | |
132 | } | |
133 | ||
134 | template<class T, class Policy> | |
135 | T polygamma_attransitionplus(const int n, const T& x, const Policy& pol, const char* function) | |
136 | { | |
137 | // See: http://functions.wolfram.com/GammaBetaErf/PolyGamma2/16/01/01/0017/ | |
138 | ||
139 | // Use N = (0.4 * digits) + (4 * n) for target value for x: | |
140 | BOOST_MATH_STD_USING | |
141 | const int d4d = static_cast<int>(0.4F * policies::digits_base10<T, Policy>()); | |
142 | const int N = d4d + (4 * n); | |
143 | const int m = n; | |
144 | const int iter = N - itrunc(x); | |
145 | ||
146 | if(iter > (int)policies::get_max_series_iterations<Policy>()) | |
147 | return policies::raise_evaluation_error<T>(function, ("Exceeded maximum series evaluations evaluating at n = " + boost::lexical_cast<std::string>(n) + " and x = %1%").c_str(), x, pol); | |
148 | ||
149 | const int minus_m_minus_one = -m - 1; | |
150 | ||
151 | T z(x); | |
152 | T sum0(0); | |
153 | T z_plus_k_pow_minus_m_minus_one(0); | |
154 | ||
155 | // Forward recursion to larger x, need to check for overflow first though: | |
156 | if(log(z + iter) * minus_m_minus_one > -tools::log_max_value<T>()) | |
157 | { | |
158 | for(int k = 1; k <= iter; ++k) | |
159 | { | |
160 | z_plus_k_pow_minus_m_minus_one = pow(z, minus_m_minus_one); | |
161 | sum0 += z_plus_k_pow_minus_m_minus_one; | |
162 | z += 1; | |
163 | } | |
164 | sum0 *= boost::math::factorial<T>(n); | |
165 | } | |
166 | else | |
167 | { | |
168 | for(int k = 1; k <= iter; ++k) | |
169 | { | |
170 | T log_term = log(z) * minus_m_minus_one + boost::math::lgamma(T(n + 1), pol); | |
171 | sum0 += exp(log_term); | |
172 | z += 1; | |
173 | } | |
174 | } | |
175 | if((n - 1) & 1) | |
176 | sum0 = -sum0; | |
177 | ||
178 | return sum0 + polygamma_atinfinityplus(n, z, pol, function); | |
179 | } | |
180 | ||
181 | template <class T, class Policy> | |
182 | T polygamma_nearzero(int n, T x, const Policy& pol, const char* function) | |
183 | { | |
184 | BOOST_MATH_STD_USING | |
185 | // | |
186 | // If we take this expansion for polygamma: http://functions.wolfram.com/06.15.06.0003.02 | |
187 | // and substitute in this expression for polygamma(n, 1): http://functions.wolfram.com/06.15.03.0009.01 | |
188 | // we get an alternating series for polygamma when x is small in terms of zeta functions of | |
189 | // integer arguments (which are easy to evaluate, at least when the integer is even). | |
190 | // | |
191 | // In order to avoid spurious overflow, save the n! term for later, and rescale at the end: | |
192 | // | |
193 | T scale = boost::math::factorial<T>(n, pol); | |
194 | // | |
195 | // "factorial_part" contains everything except the zeta function | |
196 | // evaluations in each term: | |
197 | // | |
198 | T factorial_part = 1; | |
199 | // | |
200 | // "prefix" is what we'll be adding the accumulated sum to, it will | |
201 | // be n! / z^(n+1), but since we're scaling by n! it's just | |
202 | // 1 / z^(n+1) for now: | |
203 | // | |
204 | T prefix = pow(x, n + 1); | |
205 | if(prefix == 0) | |
206 | return boost::math::policies::raise_overflow_error<T>(function, 0, pol); | |
207 | prefix = 1 / prefix; | |
208 | // | |
209 | // First term in the series is necessarily < zeta(2) < 2, so | |
210 | // ignore the sum if it will have no effect on the result anyway: | |
211 | // | |
212 | if(prefix > 2 / policies::get_epsilon<T, Policy>()) | |
213 | return ((n & 1) ? 1 : -1) * | |
214 | (tools::max_value<T>() / prefix < scale ? policies::raise_overflow_error<T>(function, 0, pol) : prefix * scale); | |
215 | // | |
216 | // As this is an alternating series we could accelerate it using | |
217 | // "Convergence Acceleration of Alternating Series", | |
218 | // Henri Cohen, Fernando Rodriguez Villegas, and Don Zagier, Experimental Mathematics, 1999. | |
219 | // In practice however, it appears not to make any difference to the number of terms | |
220 | // required except in some edge cases which are filtered out anyway before we get here. | |
221 | // | |
222 | T sum = prefix; | |
223 | for(unsigned k = 0;;) | |
224 | { | |
225 | // Get the k'th term: | |
226 | T term = factorial_part * boost::math::zeta(T(k + n + 1), pol); | |
227 | sum += term; | |
228 | // Termination condition: | |
229 | if(fabs(term) < fabs(sum * boost::math::policies::get_epsilon<T, Policy>())) | |
230 | break; | |
231 | // | |
232 | // Move on k and factorial_part: | |
233 | // | |
234 | ++k; | |
235 | factorial_part *= (-x * (n + k)) / k; | |
236 | // | |
237 | // Last chance exit: | |
238 | // | |
239 | if(k > policies::get_max_series_iterations<Policy>()) | |
240 | return policies::raise_evaluation_error<T>(function, "Series did not converge, best value is %1%", sum, pol); | |
241 | } | |
242 | // | |
243 | // We need to multiply by the scale, at each stage checking for oveflow: | |
244 | // | |
245 | if(boost::math::tools::max_value<T>() / scale < sum) | |
246 | return boost::math::policies::raise_overflow_error<T>(function, 0, pol); | |
247 | sum *= scale; | |
248 | return n & 1 ? sum : T(-sum); | |
249 | } | |
250 | ||
251 | // | |
252 | // Helper function which figures out which slot our coefficient is in | |
253 | // given an angle multiplier for the cosine term of power: | |
254 | // | |
255 | template <class Table> | |
256 | typename Table::value_type::reference dereference_table(Table& table, unsigned row, unsigned power) | |
257 | { | |
258 | return table[row][power / 2]; | |
259 | } | |
260 | ||
261 | ||
262 | ||
263 | template <class T, class Policy> | |
264 | T poly_cot_pi(int n, T x, T xc, const Policy& pol, const char* function) | |
265 | { | |
266 | BOOST_MATH_STD_USING | |
267 | // Return n'th derivative of cot(pi*x) at x, these are simply | |
268 | // tabulated for up to n = 9, beyond that it is possible to | |
269 | // calculate coefficients as follows: | |
270 | // | |
271 | // The general form of each derivative is: | |
272 | // | |
273 | // pi^n * SUM{k=0, n} C[k,n] * cos^k(pi * x) * csc^(n+1)(pi * x) | |
274 | // | |
275 | // With constant C[0,1] = -1 and all other C[k,n] = 0; | |
276 | // Then for each k < n+1: | |
277 | // C[k-1, n+1] -= k * C[k, n]; | |
278 | // C[k+1, n+1] += (k-n-1) * C[k, n]; | |
279 | // | |
280 | // Note that there are many different ways of representing this derivative thanks to | |
281 | // the many trigomonetric identies available. In particular, the sum of powers of | |
282 | // cosines could be replaced by a sum of cosine multiple angles, and indeed if you | |
283 | // plug the derivative into Mathematica this is the form it will give. The two | |
284 | // forms are related via the Chebeshev polynomials of the first kind and | |
285 | // T_n(cos(x)) = cos(n x). The polynomial form has the great advantage that | |
286 | // all the cosine terms are zero at half integer arguments - right where this | |
287 | // function has it's minumum - thus avoiding cancellation error in this region. | |
288 | // | |
289 | // And finally, since every other term in the polynomials is zero, we can save | |
290 | // space by only storing the non-zero terms. This greatly complexifies | |
291 | // subscripting the tables in the calculation, but halves the storage space | |
292 | // (and complexity for that matter). | |
293 | // | |
294 | T s = fabs(x) < fabs(xc) ? boost::math::sin_pi(x, pol) : boost::math::sin_pi(xc, pol); | |
295 | T c = boost::math::cos_pi(x, pol); | |
296 | switch(n) | |
297 | { | |
298 | case 1: | |
299 | return -constants::pi<T, Policy>() / (s * s); | |
300 | case 2: | |
301 | { | |
302 | return 2 * constants::pi<T, Policy>() * constants::pi<T, Policy>() * c / boost::math::pow<3>(s, pol); | |
303 | } | |
304 | case 3: | |
305 | { | |
306 | int P[] = { -2, -4 }; | |
307 | return boost::math::pow<3>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<4>(s, pol); | |
308 | } | |
309 | case 4: | |
310 | { | |
311 | int P[] = { 16, 8 }; | |
312 | return boost::math::pow<4>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<5>(s, pol); | |
313 | } | |
314 | case 5: | |
315 | { | |
316 | int P[] = { -16, -88, -16 }; | |
317 | return boost::math::pow<5>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<6>(s, pol); | |
318 | } | |
319 | case 6: | |
320 | { | |
321 | int P[] = { 272, 416, 32 }; | |
322 | return boost::math::pow<6>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<7>(s, pol); | |
323 | } | |
324 | case 7: | |
325 | { | |
326 | int P[] = { -272, -2880, -1824, -64 }; | |
327 | return boost::math::pow<7>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<8>(s, pol); | |
328 | } | |
329 | case 8: | |
330 | { | |
331 | int P[] = { 7936, 24576, 7680, 128 }; | |
332 | return boost::math::pow<8>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<9>(s, pol); | |
333 | } | |
334 | case 9: | |
335 | { | |
336 | int P[] = { -7936, -137216, -185856, -31616, -256 }; | |
337 | return boost::math::pow<9>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<10>(s, pol); | |
338 | } | |
339 | case 10: | |
340 | { | |
341 | int P[] = { 353792, 1841152, 1304832, 128512, 512 }; | |
342 | return boost::math::pow<10>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<11>(s, pol); | |
343 | } | |
344 | case 11: | |
345 | { | |
346 | int P[] = { -353792, -9061376, -21253376, -8728576, -518656, -1024}; | |
347 | return boost::math::pow<11>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<12>(s, pol); | |
348 | } | |
349 | case 12: | |
350 | { | |
351 | int P[] = { 22368256, 175627264, 222398464, 56520704, 2084864, 2048 }; | |
352 | return boost::math::pow<12>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<13>(s, pol); | |
353 | } | |
354 | #ifndef BOOST_NO_LONG_LONG | |
355 | case 13: | |
356 | { | |
357 | long long P[] = { -22368256LL, -795300864LL, -2868264960LL, -2174832640LL, -357888000LL, -8361984LL, -4096 }; | |
358 | return boost::math::pow<13>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<14>(s, pol); | |
359 | } | |
360 | case 14: | |
361 | { | |
362 | long long P[] = { 1903757312LL, 21016670208LL, 41731645440LL, 20261765120LL, 2230947840LL, 33497088LL, 8192 }; | |
363 | return boost::math::pow<14>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<15>(s, pol); | |
364 | } | |
365 | case 15: | |
366 | { | |
367 | long long P[] = { -1903757312LL, -89702612992LL, -460858269696LL, -559148810240LL, -182172651520LL, -13754155008LL, -134094848LL, -16384 }; | |
368 | return boost::math::pow<15>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<16>(s, pol); | |
369 | } | |
370 | case 16: | |
371 | { | |
372 | long long P[] = { 209865342976LL, 3099269660672LL, 8885192097792LL, 7048869314560LL, 1594922762240LL, 84134068224LL, 536608768LL, 32768 }; | |
373 | return boost::math::pow<16>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<17>(s, pol); | |
374 | } | |
375 | case 17: | |
376 | { | |
377 | long long P[] = { -209865342976LL, -12655654469632LL, -87815735738368LL, -155964390375424LL, -84842998005760LL, -13684856848384LL, -511780323328LL, -2146926592LL, -65536 }; | |
378 | return boost::math::pow<17>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<18>(s, pol); | |
379 | } | |
380 | case 18: | |
381 | { | |
382 | long long P[] = { 29088885112832LL, 553753414467584LL, 2165206642589696LL, 2550316668551168LL, 985278548541440LL, 115620218667008LL, 3100738912256LL, 8588754944LL, 131072 }; | |
383 | return boost::math::pow<18>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<19>(s, pol); | |
384 | } | |
385 | case 19: | |
386 | { | |
387 | long long P[] = { -29088885112832LL, -2184860175433728LL, -19686087844429824LL, -48165109676113920LL, -39471306959486976LL, -11124607890751488LL, -965271355195392LL, -18733264797696LL, -34357248000LL, -262144 }; | |
388 | return boost::math::pow<19>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<20>(s, pol); | |
389 | } | |
390 | case 20: | |
391 | { | |
392 | long long P[] = { 4951498053124096LL, 118071834535526400LL, 603968063567560704LL, 990081991141490688LL, 584901762421358592LL, 122829335169859584LL, 7984436548730880LL, 112949304754176LL, 137433710592LL, 524288 }; | |
393 | return boost::math::pow<20>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<21>(s, pol); | |
394 | } | |
395 | #endif | |
396 | } | |
397 | ||
398 | // | |
399 | // We'll have to compute the coefficients up to n, | |
400 | // complexity is O(n^2) which we don't worry about for now | |
401 | // as the values are computed once and then cached. | |
402 | // However, if the final evaluation would have too many | |
403 | // terms just bail out right away: | |
404 | // | |
405 | if((unsigned)n / 2u > policies::get_max_series_iterations<Policy>()) | |
406 | return policies::raise_evaluation_error<T>(function, "The value of n is so large that we're unable to compute the result in reasonable time, best guess is %1%", 0, pol); | |
407 | #ifdef BOOST_HAS_THREADS | |
408 | static boost::detail::lightweight_mutex m; | |
409 | boost::detail::lightweight_mutex::scoped_lock l(m); | |
410 | #endif | |
411 | static int digits = tools::digits<T>(); | |
412 | static std::vector<std::vector<T> > table(1, std::vector<T>(1, T(-1))); | |
413 | ||
414 | int current_digits = tools::digits<T>(); | |
415 | ||
416 | if(digits != current_digits) | |
417 | { | |
418 | // Oh my... our precision has changed! | |
419 | table = std::vector<std::vector<T> >(1, std::vector<T>(1, T(-1))); | |
420 | digits = current_digits; | |
421 | } | |
422 | ||
423 | int index = n - 1; | |
424 | ||
425 | if(index >= (int)table.size()) | |
426 | { | |
427 | for(int i = (int)table.size() - 1; i < index; ++i) | |
428 | { | |
429 | int offset = i & 1; // 1 if the first cos power is 0, otherwise 0. | |
430 | int sin_order = i + 2; // order of the sin term | |
431 | int max_cos_order = sin_order - 1; // largest order of the polynomial of cos terms | |
432 | int max_columns = (max_cos_order - offset) / 2; // How many entries there are in the current row. | |
433 | int next_offset = offset ? 0 : 1; | |
434 | int next_max_columns = (max_cos_order + 1 - next_offset) / 2; // How many entries there will be in the next row | |
435 | table.push_back(std::vector<T>(next_max_columns + 1, T(0))); | |
436 | ||
437 | for(int column = 0; column <= max_columns; ++column) | |
438 | { | |
439 | int cos_order = 2 * column + offset; // order of the cosine term in entry "column" | |
440 | BOOST_ASSERT(column < (int)table[i].size()); | |
441 | BOOST_ASSERT((cos_order + 1) / 2 < (int)table[i + 1].size()); | |
442 | table[i + 1][(cos_order + 1) / 2] += ((cos_order - sin_order) * table[i][column]) / (sin_order - 1); | |
443 | if(cos_order) | |
444 | table[i + 1][(cos_order - 1) / 2] += (-cos_order * table[i][column]) / (sin_order - 1); | |
445 | } | |
446 | } | |
447 | ||
448 | } | |
449 | T sum = boost::math::tools::evaluate_even_polynomial(&table[index][0], c, table[index].size()); | |
450 | if(index & 1) | |
451 | sum *= c; // First coeffient is order 1, and really an odd polynomial. | |
452 | if(sum == 0) | |
453 | return sum; | |
454 | // | |
455 | // The remaining terms are computed using logs since the powers and factorials | |
456 | // get real large real quick: | |
457 | // | |
458 | T power_terms = n * log(boost::math::constants::pi<T>()); | |
459 | if(s == 0) | |
460 | return sum * boost::math::policies::raise_overflow_error<T>(function, 0, pol); | |
461 | power_terms -= log(fabs(s)) * (n + 1); | |
462 | power_terms += boost::math::lgamma(T(n)); | |
463 | power_terms += log(fabs(sum)); | |
464 | ||
465 | if(power_terms > boost::math::tools::log_max_value<T>()) | |
466 | return sum * boost::math::policies::raise_overflow_error<T>(function, 0, pol); | |
467 | ||
468 | return exp(power_terms) * ((s < 0) && ((n + 1) & 1) ? -1 : 1) * boost::math::sign(sum); | |
469 | } | |
470 | ||
471 | template <class T, class Policy> | |
472 | struct polygamma_initializer | |
473 | { | |
474 | struct init | |
475 | { | |
476 | init() | |
477 | { | |
478 | // Forces initialization of our table of coefficients and mutex: | |
479 | boost::math::polygamma(30, T(-2.5f), Policy()); | |
480 | } | |
481 | void force_instantiate()const{} | |
482 | }; | |
483 | static const init initializer; | |
484 | static void force_instantiate() | |
485 | { | |
486 | initializer.force_instantiate(); | |
487 | } | |
488 | }; | |
489 | ||
490 | template <class T, class Policy> | |
491 | const typename polygamma_initializer<T, Policy>::init polygamma_initializer<T, Policy>::initializer; | |
492 | ||
493 | template<class T, class Policy> | |
494 | inline T polygamma_imp(const int n, T x, const Policy &pol) | |
495 | { | |
496 | BOOST_MATH_STD_USING | |
497 | static const char* function = "boost::math::polygamma<%1%>(int, %1%)"; | |
498 | polygamma_initializer<T, Policy>::initializer.force_instantiate(); | |
499 | if(n < 0) | |
500 | return policies::raise_domain_error<T>(function, "Order must be >= 0, but got %1%", static_cast<T>(n), pol); | |
501 | if(x < 0) | |
502 | { | |
503 | if(floor(x) == x) | |
504 | { | |
505 | // | |
506 | // Result is infinity if x is odd, and a pole error if x is even. | |
507 | // | |
508 | if(lltrunc(x) & 1) | |
509 | return policies::raise_overflow_error<T>(function, 0, pol); | |
510 | else | |
511 | return policies::raise_pole_error<T>(function, "Evaluation at negative integer %1%", x, pol); | |
512 | } | |
513 | T z = 1 - x; | |
514 | T result = polygamma_imp(n, z, pol) + constants::pi<T, Policy>() * poly_cot_pi(n, z, x, pol, function); | |
515 | return n & 1 ? T(-result) : result; | |
516 | } | |
517 | // | |
518 | // Limit for use of small-x-series is chosen | |
519 | // so that the series doesn't go too divergent | |
520 | // in the first few terms. Ordinarily this | |
521 | // would mean setting the limit to ~ 1 / n, | |
522 | // but we can tolerate a small amount of divergence: | |
523 | // | |
524 | T small_x_limit = (std::min)(T(T(5) / n), T(0.25f)); | |
525 | if(x < small_x_limit) | |
526 | { | |
527 | return polygamma_nearzero(n, x, pol, function); | |
528 | } | |
529 | else if(x > 0.4F * policies::digits_base10<T, Policy>() + 4.0f * n) | |
530 | { | |
531 | return polygamma_atinfinityplus(n, x, pol, function); | |
532 | } | |
533 | else if(x == 1) | |
534 | { | |
535 | return (n & 1 ? 1 : -1) * boost::math::factorial<T>(n, pol) * boost::math::zeta(T(n + 1), pol); | |
536 | } | |
537 | else if(x == 0.5f) | |
538 | { | |
539 | T result = (n & 1 ? 1 : -1) * boost::math::factorial<T>(n, pol) * boost::math::zeta(T(n + 1), pol); | |
540 | if(fabs(result) >= ldexp(tools::max_value<T>(), -n - 1)) | |
541 | return boost::math::sign(result) * policies::raise_overflow_error<T>(function, 0, pol); | |
542 | result *= ldexp(T(1), n + 1) - 1; | |
543 | return result; | |
544 | } | |
545 | else | |
546 | { | |
547 | return polygamma_attransitionplus(n, x, pol, function); | |
548 | } | |
549 | } | |
550 | ||
551 | } } } // namespace boost::math::detail | |
552 | ||
553 | #ifdef _MSC_VER | |
554 | #pragma warning(pop) | |
555 | #endif | |
556 | ||
557 | #endif // _BOOST_POLYGAMMA_DETAIL_2013_07_30_HPP_ | |
558 |