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