1// Special functions -*- C++ -*-
2
3// Copyright (C) 2006-2022 Free Software Foundation, Inc.
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
8// Free Software Foundation; either version 3, or (at your option)
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//
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/>.
24
25/** @file tr1/gamma.tcc
26 * This is an internal header file, included by other library headers.
27 * Do not attempt to use it directly. @headername{tr1/cmath}
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
46#ifndef _GLIBCXX_TR1_GAMMA_TCC
47#define _GLIBCXX_TR1_GAMMA_TCC 1
48
49#include <tr1/special_function_util.h>
50
51namespace std _GLIBCXX_VISIBILITY(default)
52{
53_GLIBCXX_BEGIN_NAMESPACE_VERSION
54
55#if _GLIBCXX_USE_STD_SPEC_FUNCS
56# define _GLIBCXX_MATH_NS ::std
57#elif defined(_GLIBCXX_TR1_CMATH)
58namespace tr1
59{
60# define _GLIBCXX_MATH_NS ::std::tr1
61#else
62# error do not include this header directly, use <cmath> or <tr1/cmath>
63#endif
64 // Implementation-space details.
65 namespace __detail
66 {
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>
77 _Tp
78 __bernoulli_series(unsigned int __n)
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
141 __bernoulli(int __n)
142 { return __bernoulli_series<_Tp>(__n); }
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
155 __log_gamma_bernoulli(_Tp __x)
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
183 __log_gamma_lanczos(_Tp __x)
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
227 __log_gamma(_Tp __x)
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
254 __log_gamma_sign(_Tp __x)
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
285 __log_bincoef(unsigned int __n, unsigned int __k)
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
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));
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
316 __bincoef(unsigned int __n, unsigned int __k)
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
339 __gamma(_Tp __x)
340 { return std::exp(__log_gamma(__x)); }
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
358 __psi_series(_Tp __x)
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[
377 * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}
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
388 __psi_asymp(_Tp __x)
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[
410 * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}
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
419 __psi(_Tp __x)
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
448 __psi(unsigned int __n, _Tp __x)
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
459 const _Tp __ln_nfact = _GLIBCXX_MATH_NS::lgamma(_Tp(__n + 1));
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 }
469 } // namespace __detail
470#undef _GLIBCXX_MATH_NS
471#if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH)
472} // namespace tr1
473#endif
474
475_GLIBCXX_END_NAMESPACE_VERSION
476} // namespace std
477
478#endif // _GLIBCXX_TR1_GAMMA_TCC
479
480