]>
Commit | Line | Data |
---|---|---|
7c62b943 BK |
1 | // Special functions -*- C++ -*- |
2 | ||
3 | // Copyright (C) 2006-2007 | |
4 | // Free Software Foundation, Inc. | |
5 | // | |
6 | // This file is part of the GNU ISO C++ Library. This library is free | |
7 | // software; you can redistribute it and/or modify it under the | |
8 | // terms of the GNU General Public License as published by the | |
9 | // Free Software Foundation; either version 2, or (at your option) | |
10 | // any later version. | |
11 | // | |
12 | // This library is distributed in the hope that it will be useful, | |
13 | // but WITHOUT ANY WARRANTY; without even the implied warranty of | |
14 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
15 | // GNU General Public License for more details. | |
16 | // | |
17 | // You should have received a copy of the GNU General Public License along | |
18 | // with this library; see the file COPYING. If not, write to the Free | |
19 | // Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, | |
20 | // USA. | |
21 | // | |
22 | // As a special exception, you may use this file as part of a free software | |
23 | // library without restriction. Specifically, if other files instantiate | |
24 | // templates or use macros or inline functions from this file, or you compile | |
25 | // this file and link it with other files to produce an executable, this | |
26 | // file does not by itself cause the resulting executable to be covered by | |
27 | // the GNU General Public License. This exception does not however | |
28 | // invalidate any other reasons why the executable file might be covered by | |
29 | // the GNU General Public License. | |
30 | ||
31 | /** @file tr1/gamma.tcc | |
32 | * This is an internal header file, included by other library headers. | |
33 | * You should not attempt to use it directly. | |
34 | */ | |
35 | ||
36 | // | |
37 | // ISO C++ 14882 TR1: 5.2 Special functions | |
38 | // | |
39 | ||
40 | // Written by Edward Smith-Rowland based on: | |
41 | // (1) Handbook of Mathematical Functions, | |
42 | // ed. Milton Abramowitz and Irene A. Stegun, | |
43 | // Dover Publications, | |
44 | // Section 6, pp. 253-266 | |
45 | // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl | |
46 | // (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky, | |
47 | // W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992), | |
48 | // 2nd ed, pp. 213-216 | |
49 | // (4) Gamma, Exploring Euler's Constant, Julian Havil, | |
50 | // Princeton, 2003. | |
51 | ||
52 | #ifndef _TR1_GAMMA_TCC | |
53 | #define _TR1_GAMMA_TCC 1 | |
54 | ||
55 | #include "special_function_util.h" | |
56 | ||
57 | namespace std | |
58 | { | |
e133ace8 PC |
59 | namespace tr1 |
60 | { | |
7c62b943 BK |
61 | |
62 | /** | |
63 | * @ingroup tr1_math_spec_func | |
64 | * @{ | |
65 | */ | |
66 | ||
67 | // | |
68 | // Implementation-space details. | |
69 | // | |
70 | namespace __detail | |
71 | { | |
72 | ||
73 | /** | |
74 | * @brief This returns Bernoulli numbers from a table or by summation | |
75 | * for larger values. | |
76 | * | |
77 | * Recursion is unstable. | |
78 | * | |
79 | * @param __n the order n of the Bernoulli number. | |
80 | * @return The Bernoulli number of order n. | |
81 | */ | |
82 | template <typename _Tp> | |
83 | _Tp __bernoulli_series(unsigned int __n) | |
84 | { | |
85 | ||
86 | static const _Tp __num[28] = { | |
87 | _Tp(1UL), -_Tp(1UL) / _Tp(2UL), | |
88 | _Tp(1UL) / _Tp(6UL), _Tp(0UL), | |
89 | -_Tp(1UL) / _Tp(30UL), _Tp(0UL), | |
90 | _Tp(1UL) / _Tp(42UL), _Tp(0UL), | |
91 | -_Tp(1UL) / _Tp(30UL), _Tp(0UL), | |
92 | _Tp(5UL) / _Tp(66UL), _Tp(0UL), | |
93 | -_Tp(691UL) / _Tp(2730UL), _Tp(0UL), | |
94 | _Tp(7UL) / _Tp(6UL), _Tp(0UL), | |
95 | -_Tp(3617UL) / _Tp(510UL), _Tp(0UL), | |
96 | _Tp(43867UL) / _Tp(798UL), _Tp(0UL), | |
97 | -_Tp(174611) / _Tp(330UL), _Tp(0UL), | |
98 | _Tp(854513UL) / _Tp(138UL), _Tp(0UL), | |
99 | -_Tp(236364091UL) / _Tp(2730UL), _Tp(0UL), | |
100 | _Tp(8553103UL) / _Tp(6UL), _Tp(0UL) | |
101 | }; | |
102 | ||
103 | if (__n == 0) | |
104 | return _Tp(1); | |
105 | ||
106 | if (__n == 1) | |
107 | return -_Tp(1) / _Tp(2); | |
108 | ||
109 | // Take care of the rest of the odd ones. | |
110 | if (__n % 2 == 1) | |
111 | return _Tp(0); | |
112 | ||
113 | // Take care of some small evens that are painful for the series. | |
114 | if (__n < 28) | |
115 | return __num[__n]; | |
116 | ||
117 | ||
118 | _Tp __fact = _Tp(1); | |
119 | if ((__n / 2) % 2 == 0) | |
120 | __fact *= _Tp(-1); | |
121 | for (unsigned int __k = 1; __k <= __n; ++__k) | |
122 | __fact *= __k / (_Tp(2) * __numeric_constants<_Tp>::__pi()); | |
123 | __fact *= _Tp(2); | |
124 | ||
125 | _Tp __sum = _Tp(0); | |
126 | for (unsigned int __i = 1; __i < 1000; ++__i) | |
127 | { | |
128 | _Tp __term = std::pow(_Tp(__i), -_Tp(__n)); | |
129 | if (__term < std::numeric_limits<_Tp>::epsilon()) | |
130 | break; | |
131 | __sum += __term; | |
132 | } | |
133 | ||
134 | return __fact * __sum; | |
135 | } | |
136 | ||
137 | ||
138 | /** | |
139 | * @brief This returns Bernoulli number \f$B_n\f$. | |
140 | * | |
141 | * @param __n the order n of the Bernoulli number. | |
142 | * @return The Bernoulli number of order n. | |
143 | */ | |
144 | template<typename _Tp> | |
145 | inline _Tp | |
146 | __bernoulli(const int __n) | |
147 | { | |
148 | return __bernoulli_series<_Tp>(__n); | |
149 | } | |
150 | ||
151 | ||
152 | /** | |
153 | * @brief Return \f$log(\Gamma(x))\f$ by asymptotic expansion | |
154 | * with Bernoulli number coefficients. This is like | |
155 | * Sterling's approximation. | |
156 | * | |
157 | * @param __x The argument of the log of the gamma function. | |
158 | * @return The logarithm of the gamma function. | |
159 | */ | |
160 | template<typename _Tp> | |
161 | _Tp | |
162 | __log_gamma_bernoulli(const _Tp __x) | |
163 | { | |
164 | _Tp __lg = (__x - _Tp(0.5L)) * std::log(__x) - __x | |
165 | + _Tp(0.5L) * std::log(_Tp(2) | |
166 | * __numeric_constants<_Tp>::__pi()); | |
167 | ||
168 | const _Tp __xx = __x * __x; | |
169 | _Tp __help = _Tp(1) / __x; | |
170 | for ( unsigned int __i = 1; __i < 20; ++__i ) | |
171 | { | |
172 | const _Tp __2i = _Tp(2 * __i); | |
173 | __help /= __2i * (__2i - _Tp(1)) * __xx; | |
174 | __lg += __bernoulli<_Tp>(2 * __i) * __help; | |
175 | } | |
176 | ||
177 | return __lg; | |
178 | } | |
179 | ||
180 | ||
181 | /** | |
182 | * @brief Return \f$log(\Gamma(x))\f$ by the Lanczos method. | |
183 | * This method dominates all others on the positive axis I think. | |
184 | * | |
185 | * @param __x The argument of the log of the gamma function. | |
186 | * @return The logarithm of the gamma function. | |
187 | */ | |
188 | template<typename _Tp> | |
189 | _Tp | |
190 | __log_gamma_lanczos(const _Tp __x) | |
191 | { | |
192 | const _Tp __xm1 = __x - _Tp(1); | |
193 | ||
194 | static const _Tp __lanczos_cheb_7[9] = { | |
195 | _Tp( 0.99999999999980993227684700473478L), | |
196 | _Tp( 676.520368121885098567009190444019L), | |
197 | _Tp(-1259.13921672240287047156078755283L), | |
198 | _Tp( 771.3234287776530788486528258894L), | |
199 | _Tp(-176.61502916214059906584551354L), | |
200 | _Tp( 12.507343278686904814458936853L), | |
201 | _Tp(-0.13857109526572011689554707L), | |
202 | _Tp( 9.984369578019570859563e-6L), | |
203 | _Tp( 1.50563273514931155834e-7L) | |
204 | }; | |
205 | ||
206 | static const _Tp __LOGROOT2PI | |
207 | = _Tp(0.9189385332046727417803297364056176L); | |
208 | ||
209 | _Tp __sum = __lanczos_cheb_7[0]; | |
210 | for(unsigned int __k = 1; __k < 9; ++__k) | |
211 | __sum += __lanczos_cheb_7[__k] / (__xm1 + __k); | |
212 | ||
213 | const _Tp __term1 = (__xm1 + _Tp(0.5L)) | |
214 | * std::log((__xm1 + _Tp(7.5L)) | |
215 | / __numeric_constants<_Tp>::__euler()); | |
216 | const _Tp __term2 = __LOGROOT2PI + std::log(__sum); | |
217 | const _Tp __result = __term1 + (__term2 - _Tp(7)); | |
218 | ||
219 | return __result; | |
220 | } | |
221 | ||
222 | ||
223 | /** | |
224 | * @brief Return \f$ log(|\Gamma(x)|) \f$. | |
225 | * This will return values even for \f$ x < 0 \f$. | |
226 | * To recover the sign of \f$ \Gamma(x) \f$ for | |
227 | * any argument use @a __log_gamma_sign. | |
228 | * | |
229 | * @param __x The argument of the log of the gamma function. | |
230 | * @return The logarithm of the gamma function. | |
231 | */ | |
232 | template<typename _Tp> | |
233 | _Tp | |
234 | __log_gamma(const _Tp __x) | |
235 | { | |
236 | if (__x > _Tp(0.5L)) | |
237 | return __log_gamma_lanczos(__x); | |
238 | else | |
239 | { | |
240 | const _Tp __sin_fact | |
241 | = std::abs(std::sin(__numeric_constants<_Tp>::__pi() * __x)); | |
242 | if (__sin_fact == _Tp(0)) | |
243 | std::__throw_domain_error(__N("Argument is nonpositive integer " | |
244 | "in __log_gamma")); | |
245 | return __numeric_constants<_Tp>::__lnpi() | |
246 | - std::log(__sin_fact) | |
247 | - __log_gamma_lanczos(_Tp(1) - __x); | |
248 | } | |
249 | } | |
250 | ||
251 | ||
252 | /** | |
253 | * @brief Return the sign of \f$ \Gamma(x) \f$. | |
254 | * At nonpositive integers zero is returned. | |
255 | * | |
256 | * @param __x The argument of the gamma function. | |
257 | * @return The sign of the gamma function. | |
258 | */ | |
259 | template<typename _Tp> | |
260 | _Tp | |
261 | __log_gamma_sign(const _Tp __x) | |
262 | { | |
263 | if (__x > _Tp(0)) | |
264 | return _Tp(1); | |
265 | else | |
266 | { | |
267 | const _Tp __sin_fact | |
268 | = std::sin(__numeric_constants<_Tp>::__pi() * __x); | |
269 | if (__sin_fact > _Tp(0)) | |
270 | return (1); | |
271 | else if (__sin_fact < _Tp(0)) | |
272 | return -_Tp(1); | |
273 | else | |
274 | return _Tp(0); | |
275 | } | |
276 | } | |
277 | ||
278 | ||
279 | /** | |
280 | * @brief Return the logarithm of the binomial coefficient. | |
281 | * The binomial coefficient is given by: | |
282 | * @f[ | |
283 | * \left( \right) = \frac{n!}{(n-k)! k!} | |
284 | * @f] | |
285 | * | |
286 | * @param __n The first argument of the binomial coefficient. | |
287 | * @param __k The second argument of the binomial coefficient. | |
288 | * @return The binomial coefficient. | |
289 | */ | |
290 | template<typename _Tp> | |
291 | _Tp | |
292 | __log_bincoef(const unsigned int __n, const unsigned int __k) | |
293 | { | |
294 | // Max e exponent before overflow. | |
295 | static const _Tp __max_bincoeff | |
296 | = std::numeric_limits<_Tp>::max_exponent10 | |
297 | * std::log(_Tp(10)) - _Tp(1); | |
298 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
e133ace8 PC |
299 | _Tp __coeff = std::tr1::lgamma(_Tp(1 + __n)) |
300 | - std::tr1::lgamma(_Tp(1 + __k)) | |
301 | - std::tr1::lgamma(_Tp(1 + __n - __k)); | |
7c62b943 BK |
302 | #else |
303 | _Tp __coeff = __log_gamma(_Tp(1 + __n)) | |
304 | - __log_gamma(_Tp(1 + __k)) | |
305 | - __log_gamma(_Tp(1 + __n - __k)); | |
306 | #endif | |
307 | } | |
308 | ||
309 | ||
310 | /** | |
311 | * @brief Return the binomial coefficient. | |
312 | * The binomial coefficient is given by: | |
313 | * @f[ | |
314 | * \left( \right) = \frac{n!}{(n-k)! k!} | |
315 | * @f] | |
316 | * | |
317 | * @param __n The first argument of the binomial coefficient. | |
318 | * @param __k The second argument of the binomial coefficient. | |
319 | * @return The binomial coefficient. | |
320 | */ | |
321 | template<typename _Tp> | |
322 | _Tp | |
323 | __bincoef(const unsigned int __n, const unsigned int __k) | |
324 | { | |
325 | // Max e exponent before overflow. | |
326 | static const _Tp __max_bincoeff | |
327 | = std::numeric_limits<_Tp>::max_exponent10 | |
328 | * std::log(_Tp(10)) - _Tp(1); | |
329 | ||
330 | const _Tp __log_coeff = __log_bincoef<_Tp>(__n, __k); | |
331 | if (__log_coeff > __max_bincoeff) | |
332 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
333 | else | |
334 | return std::exp(__log_coeff); | |
335 | } | |
336 | ||
337 | ||
338 | /** | |
339 | * @brief Return \f$ \Gamma(x) \f$. | |
340 | * | |
341 | * @param __x The argument of the gamma function. | |
342 | * @return The gamma function. | |
343 | */ | |
344 | template<typename _Tp> | |
345 | inline _Tp | |
346 | __gamma(const _Tp __x) | |
347 | { | |
348 | return std::exp(__log_gamma(__x)); | |
349 | } | |
350 | ||
351 | ||
352 | /** | |
353 | * @brief Return the digamma function by series expansion. | |
354 | * The digamma or @f$ \psi(x) @f$ function is defined by | |
355 | * @f[ | |
356 | * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)} | |
357 | * @f] | |
358 | * | |
359 | * The series is given by: | |
360 | * @f[ | |
361 | * \psi(x) = -\gamma_E - \frac{1}{x} | |
362 | * \sum_{k=1}^{\infty} \frac{x}{k(x + k)} | |
363 | * @f] | |
364 | */ | |
365 | template<typename _Tp> | |
366 | _Tp | |
367 | __psi_series(const _Tp __x) | |
368 | { | |
369 | _Tp __sum = -__numeric_constants<_Tp>::__gamma_e() - _Tp(1) / __x; | |
370 | const unsigned int __max_iter = 100000; | |
371 | for (unsigned int __k = 1; __k < __max_iter; ++__k) | |
372 | { | |
373 | const _Tp __term = __x / (__k * (__k + __x)); | |
374 | __sum += __term; | |
375 | if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon()) | |
376 | break; | |
377 | } | |
378 | return __sum; | |
379 | } | |
380 | ||
381 | ||
382 | /** | |
383 | * @brief Return the digamma function for large argument. | |
384 | * The digamma or @f$ \psi(x) @f$ function is defined by | |
385 | * @f[ | |
386 | * \psi(x) = \frac{Gamma'(x)}{\Gamma(x)} | |
387 | * @f] | |
388 | * | |
389 | * The asymptotic series is given by: | |
390 | * @f[ | |
391 | * \psi(x) = \ln(x) - \frac{1}{2x} | |
392 | * - \sum_{n=1}^{\infty} \frac{B_{2n}}{2 n x^{2n}} | |
393 | * @f] | |
394 | */ | |
395 | template<typename _Tp> | |
396 | _Tp | |
397 | __psi_asymp(const _Tp __x) | |
398 | { | |
399 | _Tp __sum = std::log(__x) - _Tp(0.5L) / __x; | |
400 | const _Tp __xx = __x * __x; | |
401 | _Tp __xp = __xx; | |
402 | const unsigned int __max_iter = 100; | |
403 | for (unsigned int __k = 1; __k < __max_iter; ++__k) | |
404 | { | |
405 | const _Tp __term = __bernoulli<_Tp>(2 * __k) / (2 * __k * __xp); | |
406 | __sum -= __term; | |
407 | if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon()) | |
408 | break; | |
409 | __xp *= __xx; | |
410 | } | |
411 | return __sum; | |
412 | } | |
413 | ||
414 | ||
415 | /** | |
416 | * @brief Return the digamma function. | |
417 | * The digamma or @f$ \psi(x) @f$ function is defined by | |
418 | * @f[ | |
419 | * \psi(x) = \frac{Gamma'(x)}{\Gamma(x)} | |
420 | * @f] | |
421 | * For negative argument the reflection formula is used: | |
422 | * @f[ | |
423 | * \psi(x) = \psi(1-x) - \pi \cot(\pi x) | |
424 | * @f] | |
425 | */ | |
426 | template<typename _Tp> | |
427 | _Tp | |
428 | __psi(const _Tp __x) | |
429 | { | |
430 | const int __n = static_cast<int>(__x + 0.5L); | |
431 | const _Tp __eps = _Tp(4) * std::numeric_limits<_Tp>::epsilon(); | |
432 | if (__n <= 0 && std::abs(__x - _Tp(__n)) < __eps) | |
433 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
434 | else if (__x < _Tp(0)) | |
435 | { | |
436 | const _Tp __pi = __numeric_constants<_Tp>::__pi(); | |
437 | return __psi(_Tp(1) - __x) | |
438 | - __pi * std::cos(__pi * __x) / std::sin(__pi * __x); | |
439 | } | |
440 | else if (__x > _Tp(100)) | |
441 | return __psi_asymp(__x); | |
442 | else | |
443 | return __psi_series(__x); | |
444 | } | |
445 | ||
446 | ||
447 | /** | |
448 | * @brief Return the polygamma function @f$ \psi^{(n)}(x) @f$. | |
449 | * | |
450 | * The polygamma function is related to the Hurwitz zeta function: | |
451 | * @f[ | |
452 | * \psi^{(n)}(x) = (-1)^{n+1} m! \zeta(m+1,x) | |
453 | * @f] | |
454 | */ | |
455 | template<typename _Tp> | |
456 | _Tp | |
457 | __psi(const unsigned int __n, const _Tp __x) | |
458 | { | |
459 | if (__x <= _Tp(0)) | |
460 | std::__throw_domain_error(__N("Argument out of range " | |
461 | "in __psi")); | |
462 | else if (__n == 0) | |
463 | return __psi(__x); | |
464 | else | |
465 | { | |
466 | const _Tp __hzeta = __hurwitz_zeta(_Tp(__n + 1), __x); | |
467 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
e133ace8 | 468 | const _Tp __ln_nfact = std::tr1::lgamma(_Tp(__n + 1)); |
7c62b943 BK |
469 | #else |
470 | const _Tp __ln_nfact = __log_gamma(_Tp(__n + 1)); | |
471 | #endif | |
472 | _Tp __result = std::exp(__ln_nfact) * __hzeta; | |
473 | if (__n % 2 == 1) | |
474 | __result = -__result; | |
475 | return __result; | |
476 | } | |
477 | } | |
478 | ||
479 | } // namespace std::tr1::__detail | |
480 | ||
481 | /* @} */ // group tr1_math_spec_func | |
482 | ||
e133ace8 | 483 | } |
7c62b943 BK |
484 | } |
485 | ||
486 | #endif // _TR1_GAMMA_TCC | |
487 |