1 // Special functions -*- C++ -*-
3 // Copyright (C) 2006, 2007, 2008, 2009, 2010
4 // Free Software Foundation, Inc.
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 3, or (at your option)
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.
17 // Under Section 7 of GPL version 3, you are granted additional
18 // permissions described in the GCC Runtime Library Exception, version
19 // 3.1, as published by the Free Software Foundation.
21 // You should have received a copy of the GNU General Public License and
22 // a copy of the GCC Runtime Library Exception along with this program;
23 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
24 // <http://www.gnu.org/licenses/>.
26 /** @file tr1/gamma.tcc
27 * This is an internal header file, included by other library headers.
28 * Do not attempt to use it directly. @headername{tr1/cmath}
32 // ISO C++ 14882 TR1: 5.2 Special functions
35 // Written by Edward Smith-Rowland based on:
36 // (1) Handbook of Mathematical Functions,
37 // ed. Milton Abramowitz and Irene A. Stegun,
38 // Dover Publications,
39 // Section 6, pp. 253-266
40 // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl
41 // (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky,
42 // W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992),
43 // 2nd ed, pp. 213-216
44 // (4) Gamma, Exploring Euler's Constant, Julian Havil,
47 #ifndef _GLIBCXX_TR1_GAMMA_TCC
48 #define _GLIBCXX_TR1_GAMMA_TCC 1
50 #include "special_function_util.h"
52 namespace std _GLIBCXX_VISIBILITY(default)
56 // Implementation-space details.
59 _GLIBCXX_BEGIN_NAMESPACE_VERSION
62 * @brief This returns Bernoulli numbers from a table or by summation
65 * Recursion is unstable.
67 * @param __n the order n of the Bernoulli number.
68 * @return The Bernoulli number of order n.
70 template <typename _Tp>
71 _Tp __bernoulli_series(unsigned int __n)
74 static const _Tp __num[28] = {
75 _Tp(1UL), -_Tp(1UL) / _Tp(2UL),
76 _Tp(1UL) / _Tp(6UL), _Tp(0UL),
77 -_Tp(1UL) / _Tp(30UL), _Tp(0UL),
78 _Tp(1UL) / _Tp(42UL), _Tp(0UL),
79 -_Tp(1UL) / _Tp(30UL), _Tp(0UL),
80 _Tp(5UL) / _Tp(66UL), _Tp(0UL),
81 -_Tp(691UL) / _Tp(2730UL), _Tp(0UL),
82 _Tp(7UL) / _Tp(6UL), _Tp(0UL),
83 -_Tp(3617UL) / _Tp(510UL), _Tp(0UL),
84 _Tp(43867UL) / _Tp(798UL), _Tp(0UL),
85 -_Tp(174611) / _Tp(330UL), _Tp(0UL),
86 _Tp(854513UL) / _Tp(138UL), _Tp(0UL),
87 -_Tp(236364091UL) / _Tp(2730UL), _Tp(0UL),
88 _Tp(8553103UL) / _Tp(6UL), _Tp(0UL)
95 return -_Tp(1) / _Tp(2);
97 // Take care of the rest of the odd ones.
101 // Take care of some small evens that are painful for the series.
107 if ((__n / 2) % 2 == 0)
109 for (unsigned int __k = 1; __k <= __n; ++__k)
110 __fact *= __k / (_Tp(2) * __numeric_constants<_Tp>::__pi());
114 for (unsigned int __i = 1; __i < 1000; ++__i)
116 _Tp __term = std::pow(_Tp(__i), -_Tp(__n));
117 if (__term < std::numeric_limits<_Tp>::epsilon())
122 return __fact * __sum;
127 * @brief This returns Bernoulli number \f$B_n\f$.
129 * @param __n the order n of the Bernoulli number.
130 * @return The Bernoulli number of order n.
132 template<typename _Tp>
134 __bernoulli(const int __n)
136 return __bernoulli_series<_Tp>(__n);
141 * @brief Return \f$log(\Gamma(x))\f$ by asymptotic expansion
142 * with Bernoulli number coefficients. This is like
143 * Sterling's approximation.
145 * @param __x The argument of the log of the gamma function.
146 * @return The logarithm of the gamma function.
148 template<typename _Tp>
150 __log_gamma_bernoulli(const _Tp __x)
152 _Tp __lg = (__x - _Tp(0.5L)) * std::log(__x) - __x
153 + _Tp(0.5L) * std::log(_Tp(2)
154 * __numeric_constants<_Tp>::__pi());
156 const _Tp __xx = __x * __x;
157 _Tp __help = _Tp(1) / __x;
158 for ( unsigned int __i = 1; __i < 20; ++__i )
160 const _Tp __2i = _Tp(2 * __i);
161 __help /= __2i * (__2i - _Tp(1)) * __xx;
162 __lg += __bernoulli<_Tp>(2 * __i) * __help;
170 * @brief Return \f$log(\Gamma(x))\f$ by the Lanczos method.
171 * This method dominates all others on the positive axis I think.
173 * @param __x The argument of the log of the gamma function.
174 * @return The logarithm of the gamma function.
176 template<typename _Tp>
178 __log_gamma_lanczos(const _Tp __x)
180 const _Tp __xm1 = __x - _Tp(1);
182 static const _Tp __lanczos_cheb_7[9] = {
183 _Tp( 0.99999999999980993227684700473478L),
184 _Tp( 676.520368121885098567009190444019L),
185 _Tp(-1259.13921672240287047156078755283L),
186 _Tp( 771.3234287776530788486528258894L),
187 _Tp(-176.61502916214059906584551354L),
188 _Tp( 12.507343278686904814458936853L),
189 _Tp(-0.13857109526572011689554707L),
190 _Tp( 9.984369578019570859563e-6L),
191 _Tp( 1.50563273514931155834e-7L)
194 static const _Tp __LOGROOT2PI
195 = _Tp(0.9189385332046727417803297364056176L);
197 _Tp __sum = __lanczos_cheb_7[0];
198 for(unsigned int __k = 1; __k < 9; ++__k)
199 __sum += __lanczos_cheb_7[__k] / (__xm1 + __k);
201 const _Tp __term1 = (__xm1 + _Tp(0.5L))
202 * std::log((__xm1 + _Tp(7.5L))
203 / __numeric_constants<_Tp>::__euler());
204 const _Tp __term2 = __LOGROOT2PI + std::log(__sum);
205 const _Tp __result = __term1 + (__term2 - _Tp(7));
212 * @brief Return \f$ log(|\Gamma(x)|) \f$.
213 * This will return values even for \f$ x < 0 \f$.
214 * To recover the sign of \f$ \Gamma(x) \f$ for
215 * any argument use @a __log_gamma_sign.
217 * @param __x The argument of the log of the gamma function.
218 * @return The logarithm of the gamma function.
220 template<typename _Tp>
222 __log_gamma(const _Tp __x)
225 return __log_gamma_lanczos(__x);
229 = std::abs(std::sin(__numeric_constants<_Tp>::__pi() * __x));
230 if (__sin_fact == _Tp(0))
231 std::__throw_domain_error(__N("Argument is nonpositive integer "
233 return __numeric_constants<_Tp>::__lnpi()
234 - std::log(__sin_fact)
235 - __log_gamma_lanczos(_Tp(1) - __x);
241 * @brief Return the sign of \f$ \Gamma(x) \f$.
242 * At nonpositive integers zero is returned.
244 * @param __x The argument of the gamma function.
245 * @return The sign of the gamma function.
247 template<typename _Tp>
249 __log_gamma_sign(const _Tp __x)
256 = std::sin(__numeric_constants<_Tp>::__pi() * __x);
257 if (__sin_fact > _Tp(0))
259 else if (__sin_fact < _Tp(0))
268 * @brief Return the logarithm of the binomial coefficient.
269 * The binomial coefficient is given by:
271 * \left( \right) = \frac{n!}{(n-k)! k!}
274 * @param __n The first argument of the binomial coefficient.
275 * @param __k The second argument of the binomial coefficient.
276 * @return The binomial coefficient.
278 template<typename _Tp>
280 __log_bincoef(const unsigned int __n, const unsigned int __k)
282 // Max e exponent before overflow.
283 static const _Tp __max_bincoeff
284 = std::numeric_limits<_Tp>::max_exponent10
285 * std::log(_Tp(10)) - _Tp(1);
286 #if _GLIBCXX_USE_C99_MATH_TR1
287 _Tp __coeff = std::tr1::lgamma(_Tp(1 + __n))
288 - std::tr1::lgamma(_Tp(1 + __k))
289 - std::tr1::lgamma(_Tp(1 + __n - __k));
291 _Tp __coeff = __log_gamma(_Tp(1 + __n))
292 - __log_gamma(_Tp(1 + __k))
293 - __log_gamma(_Tp(1 + __n - __k));
299 * @brief Return the binomial coefficient.
300 * The binomial coefficient is given by:
302 * \left( \right) = \frac{n!}{(n-k)! k!}
305 * @param __n The first argument of the binomial coefficient.
306 * @param __k The second argument of the binomial coefficient.
307 * @return The binomial coefficient.
309 template<typename _Tp>
311 __bincoef(const unsigned int __n, const unsigned int __k)
313 // Max e exponent before overflow.
314 static const _Tp __max_bincoeff
315 = std::numeric_limits<_Tp>::max_exponent10
316 * std::log(_Tp(10)) - _Tp(1);
318 const _Tp __log_coeff = __log_bincoef<_Tp>(__n, __k);
319 if (__log_coeff > __max_bincoeff)
320 return std::numeric_limits<_Tp>::quiet_NaN();
322 return std::exp(__log_coeff);
327 * @brief Return \f$ \Gamma(x) \f$.
329 * @param __x The argument of the gamma function.
330 * @return The gamma function.
332 template<typename _Tp>
334 __gamma(const _Tp __x)
336 return std::exp(__log_gamma(__x));
341 * @brief Return the digamma function by series expansion.
342 * The digamma or @f$ \psi(x) @f$ function is defined by
344 * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}
347 * The series is given by:
349 * \psi(x) = -\gamma_E - \frac{1}{x}
350 * \sum_{k=1}^{\infty} \frac{x}{k(x + k)}
353 template<typename _Tp>
355 __psi_series(const _Tp __x)
357 _Tp __sum = -__numeric_constants<_Tp>::__gamma_e() - _Tp(1) / __x;
358 const unsigned int __max_iter = 100000;
359 for (unsigned int __k = 1; __k < __max_iter; ++__k)
361 const _Tp __term = __x / (__k * (__k + __x));
363 if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon())
371 * @brief Return the digamma function for large argument.
372 * The digamma or @f$ \psi(x) @f$ function is defined by
374 * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}
377 * The asymptotic series is given by:
379 * \psi(x) = \ln(x) - \frac{1}{2x}
380 * - \sum_{n=1}^{\infty} \frac{B_{2n}}{2 n x^{2n}}
383 template<typename _Tp>
385 __psi_asymp(const _Tp __x)
387 _Tp __sum = std::log(__x) - _Tp(0.5L) / __x;
388 const _Tp __xx = __x * __x;
390 const unsigned int __max_iter = 100;
391 for (unsigned int __k = 1; __k < __max_iter; ++__k)
393 const _Tp __term = __bernoulli<_Tp>(2 * __k) / (2 * __k * __xp);
395 if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon())
404 * @brief Return the digamma function.
405 * The digamma or @f$ \psi(x) @f$ function is defined by
407 * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}
409 * For negative argument the reflection formula is used:
411 * \psi(x) = \psi(1-x) - \pi \cot(\pi x)
414 template<typename _Tp>
418 const int __n = static_cast<int>(__x + 0.5L);
419 const _Tp __eps = _Tp(4) * std::numeric_limits<_Tp>::epsilon();
420 if (__n <= 0 && std::abs(__x - _Tp(__n)) < __eps)
421 return std::numeric_limits<_Tp>::quiet_NaN();
422 else if (__x < _Tp(0))
424 const _Tp __pi = __numeric_constants<_Tp>::__pi();
425 return __psi(_Tp(1) - __x)
426 - __pi * std::cos(__pi * __x) / std::sin(__pi * __x);
428 else if (__x > _Tp(100))
429 return __psi_asymp(__x);
431 return __psi_series(__x);
436 * @brief Return the polygamma function @f$ \psi^{(n)}(x) @f$.
438 * The polygamma function is related to the Hurwitz zeta function:
440 * \psi^{(n)}(x) = (-1)^{n+1} m! \zeta(m+1,x)
443 template<typename _Tp>
445 __psi(const unsigned int __n, const _Tp __x)
448 std::__throw_domain_error(__N("Argument out of range "
454 const _Tp __hzeta = __hurwitz_zeta(_Tp(__n + 1), __x);
455 #if _GLIBCXX_USE_C99_MATH_TR1
456 const _Tp __ln_nfact = std::tr1::lgamma(_Tp(__n + 1));
458 const _Tp __ln_nfact = __log_gamma(_Tp(__n + 1));
460 _Tp __result = std::exp(__ln_nfact) * __hzeta;
462 __result = -__result;
467 _GLIBCXX_END_NAMESPACE_VERSION
468 } // namespace std::tr1::__detail
472 #endif // _GLIBCXX_TR1_GAMMA_TCC