diff options
Diffstat (limited to 'gcc-4.4.3/libstdc++-v3/include/tr1/gamma.tcc')
| -rw-r--r-- | gcc-4.4.3/libstdc++-v3/include/tr1/gamma.tcc | 471 |
1 files changed, 471 insertions, 0 deletions
diff --git a/gcc-4.4.3/libstdc++-v3/include/tr1/gamma.tcc b/gcc-4.4.3/libstdc++-v3/include/tr1/gamma.tcc new file mode 100644 index 000000000..f456da32b --- /dev/null +++ b/gcc-4.4.3/libstdc++-v3/include/tr1/gamma.tcc @@ -0,0 +1,471 @@ +// Special functions -*- C++ -*- + +// Copyright (C) 2006, 2007, 2008, 2009 +// Free Software Foundation, Inc. +// +// This file is part of the GNU ISO C++ Library. This library is free +// software; you can redistribute it and/or modify it under the +// terms of the GNU General Public License as published by the +// Free Software Foundation; either version 3, or (at your option) +// any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// Under Section 7 of GPL version 3, you are granted additional +// permissions described in the GCC Runtime Library Exception, version +// 3.1, as published by the Free Software Foundation. + +// You should have received a copy of the GNU General Public License and +// a copy of the GCC Runtime Library Exception along with this program; +// see the files COPYING3 and COPYING.RUNTIME respectively. If not, see +// <http://www.gnu.org/licenses/>. + +/** @file tr1/gamma.tcc + * This is an internal header file, included by other library headers. + * You should not attempt to use it directly. + */ + +// +// ISO C++ 14882 TR1: 5.2 Special functions +// + +// Written by Edward Smith-Rowland based on: +// (1) Handbook of Mathematical Functions, +// ed. Milton Abramowitz and Irene A. Stegun, +// Dover Publications, +// Section 6, pp. 253-266 +// (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl +// (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky, +// W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992), +// 2nd ed, pp. 213-216 +// (4) Gamma, Exploring Euler's Constant, Julian Havil, +// Princeton, 2003. + +#ifndef _TR1_GAMMA_TCC +#define _TR1_GAMMA_TCC 1 + +#include "special_function_util.h" + +namespace std +{ +namespace tr1 +{ + // Implementation-space details. + namespace __detail + { + + /** + * @brief This returns Bernoulli numbers from a table or by summation + * for larger values. + * + * Recursion is unstable. + * + * @param __n the order n of the Bernoulli number. + * @return The Bernoulli number of order n. + */ + template <typename _Tp> + _Tp __bernoulli_series(unsigned int __n) + { + + static const _Tp __num[28] = { + _Tp(1UL), -_Tp(1UL) / _Tp(2UL), + _Tp(1UL) / _Tp(6UL), _Tp(0UL), + -_Tp(1UL) / _Tp(30UL), _Tp(0UL), + _Tp(1UL) / _Tp(42UL), _Tp(0UL), + -_Tp(1UL) / _Tp(30UL), _Tp(0UL), + _Tp(5UL) / _Tp(66UL), _Tp(0UL), + -_Tp(691UL) / _Tp(2730UL), _Tp(0UL), + _Tp(7UL) / _Tp(6UL), _Tp(0UL), + -_Tp(3617UL) / _Tp(510UL), _Tp(0UL), + _Tp(43867UL) / _Tp(798UL), _Tp(0UL), + -_Tp(174611) / _Tp(330UL), _Tp(0UL), + _Tp(854513UL) / _Tp(138UL), _Tp(0UL), + -_Tp(236364091UL) / _Tp(2730UL), _Tp(0UL), + _Tp(8553103UL) / _Tp(6UL), _Tp(0UL) + }; + + if (__n == 0) + return _Tp(1); + + if (__n == 1) + return -_Tp(1) / _Tp(2); + + // Take care of the rest of the odd ones. + if (__n % 2 == 1) + return _Tp(0); + + // Take care of some small evens that are painful for the series. + if (__n < 28) + return __num[__n]; + + + _Tp __fact = _Tp(1); + if ((__n / 2) % 2 == 0) + __fact *= _Tp(-1); + for (unsigned int __k = 1; __k <= __n; ++__k) + __fact *= __k / (_Tp(2) * __numeric_constants<_Tp>::__pi()); + __fact *= _Tp(2); + + _Tp __sum = _Tp(0); + for (unsigned int __i = 1; __i < 1000; ++__i) + { + _Tp __term = std::pow(_Tp(__i), -_Tp(__n)); + if (__term < std::numeric_limits<_Tp>::epsilon()) + break; + __sum += __term; + } + + return __fact * __sum; + } + + + /** + * @brief This returns Bernoulli number \f$B_n\f$. + * + * @param __n the order n of the Bernoulli number. + * @return The Bernoulli number of order n. + */ + template<typename _Tp> + inline _Tp + __bernoulli(const int __n) + { + return __bernoulli_series<_Tp>(__n); + } + + + /** + * @brief Return \f$log(\Gamma(x))\f$ by asymptotic expansion + * with Bernoulli number coefficients. This is like + * Sterling's approximation. + * + * @param __x The argument of the log of the gamma function. + * @return The logarithm of the gamma function. + */ + template<typename _Tp> + _Tp + __log_gamma_bernoulli(const _Tp __x) + { + _Tp __lg = (__x - _Tp(0.5L)) * std::log(__x) - __x + + _Tp(0.5L) * std::log(_Tp(2) + * __numeric_constants<_Tp>::__pi()); + + const _Tp __xx = __x * __x; + _Tp __help = _Tp(1) / __x; + for ( unsigned int __i = 1; __i < 20; ++__i ) + { + const _Tp __2i = _Tp(2 * __i); + __help /= __2i * (__2i - _Tp(1)) * __xx; + __lg += __bernoulli<_Tp>(2 * __i) * __help; + } + + return __lg; + } + + + /** + * @brief Return \f$log(\Gamma(x))\f$ by the Lanczos method. + * This method dominates all others on the positive axis I think. + * + * @param __x The argument of the log of the gamma function. + * @return The logarithm of the gamma function. + */ + template<typename _Tp> + _Tp + __log_gamma_lanczos(const _Tp __x) + { + const _Tp __xm1 = __x - _Tp(1); + + static const _Tp __lanczos_cheb_7[9] = { + _Tp( 0.99999999999980993227684700473478L), + _Tp( 676.520368121885098567009190444019L), + _Tp(-1259.13921672240287047156078755283L), + _Tp( 771.3234287776530788486528258894L), + _Tp(-176.61502916214059906584551354L), + _Tp( 12.507343278686904814458936853L), + _Tp(-0.13857109526572011689554707L), + _Tp( 9.984369578019570859563e-6L), + _Tp( 1.50563273514931155834e-7L) + }; + + static const _Tp __LOGROOT2PI + = _Tp(0.9189385332046727417803297364056176L); + + _Tp __sum = __lanczos_cheb_7[0]; + for(unsigned int __k = 1; __k < 9; ++__k) + __sum += __lanczos_cheb_7[__k] / (__xm1 + __k); + + const _Tp __term1 = (__xm1 + _Tp(0.5L)) + * std::log((__xm1 + _Tp(7.5L)) + / __numeric_constants<_Tp>::__euler()); + const _Tp __term2 = __LOGROOT2PI + std::log(__sum); + const _Tp __result = __term1 + (__term2 - _Tp(7)); + + return __result; + } + + + /** + * @brief Return \f$ log(|\Gamma(x)|) \f$. + * This will return values even for \f$ x < 0 \f$. + * To recover the sign of \f$ \Gamma(x) \f$ for + * any argument use @a __log_gamma_sign. + * + * @param __x The argument of the log of the gamma function. + * @return The logarithm of the gamma function. + */ + template<typename _Tp> + _Tp + __log_gamma(const _Tp __x) + { + if (__x > _Tp(0.5L)) + return __log_gamma_lanczos(__x); + else + { + const _Tp __sin_fact + = std::abs(std::sin(__numeric_constants<_Tp>::__pi() * __x)); + if (__sin_fact == _Tp(0)) + std::__throw_domain_error(__N("Argument is nonpositive integer " + "in __log_gamma")); + return __numeric_constants<_Tp>::__lnpi() + - std::log(__sin_fact) + - __log_gamma_lanczos(_Tp(1) - __x); + } + } + + + /** + * @brief Return the sign of \f$ \Gamma(x) \f$. + * At nonpositive integers zero is returned. + * + * @param __x The argument of the gamma function. + * @return The sign of the gamma function. + */ + template<typename _Tp> + _Tp + __log_gamma_sign(const _Tp __x) + { + if (__x > _Tp(0)) + return _Tp(1); + else + { + const _Tp __sin_fact + = std::sin(__numeric_constants<_Tp>::__pi() * __x); + if (__sin_fact > _Tp(0)) + return (1); + else if (__sin_fact < _Tp(0)) + return -_Tp(1); + else + return _Tp(0); + } + } + + + /** + * @brief Return the logarithm of the binomial coefficient. + * The binomial coefficient is given by: + * @f[ + * \left( \right) = \frac{n!}{(n-k)! k!} + * @f] + * + * @param __n The first argument of the binomial coefficient. + * @param __k The second argument of the binomial coefficient. + * @return The binomial coefficient. + */ + template<typename _Tp> + _Tp + __log_bincoef(const unsigned int __n, const unsigned int __k) + { + // Max e exponent before overflow. + static const _Tp __max_bincoeff + = std::numeric_limits<_Tp>::max_exponent10 + * std::log(_Tp(10)) - _Tp(1); +#if _GLIBCXX_USE_C99_MATH_TR1 + _Tp __coeff = std::tr1::lgamma(_Tp(1 + __n)) + - std::tr1::lgamma(_Tp(1 + __k)) + - std::tr1::lgamma(_Tp(1 + __n - __k)); +#else + _Tp __coeff = __log_gamma(_Tp(1 + __n)) + - __log_gamma(_Tp(1 + __k)) + - __log_gamma(_Tp(1 + __n - __k)); +#endif + } + + + /** + * @brief Return the binomial coefficient. + * The binomial coefficient is given by: + * @f[ + * \left( \right) = \frac{n!}{(n-k)! k!} + * @f] + * + * @param __n The first argument of the binomial coefficient. + * @param __k The second argument of the binomial coefficient. + * @return The binomial coefficient. + */ + template<typename _Tp> + _Tp + __bincoef(const unsigned int __n, const unsigned int __k) + { + // Max e exponent before overflow. + static const _Tp __max_bincoeff + = std::numeric_limits<_Tp>::max_exponent10 + * std::log(_Tp(10)) - _Tp(1); + + const _Tp __log_coeff = __log_bincoef<_Tp>(__n, __k); + if (__log_coeff > __max_bincoeff) + return std::numeric_limits<_Tp>::quiet_NaN(); + else + return std::exp(__log_coeff); + } + + + /** + * @brief Return \f$ \Gamma(x) \f$. + * + * @param __x The argument of the gamma function. + * @return The gamma function. + */ + template<typename _Tp> + inline _Tp + __gamma(const _Tp __x) + { + return std::exp(__log_gamma(__x)); + } + + + /** + * @brief Return the digamma function by series expansion. + * The digamma or @f$ \psi(x) @f$ function is defined by + * @f[ + * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)} + * @f] + * + * The series is given by: + * @f[ + * \psi(x) = -\gamma_E - \frac{1}{x} + * \sum_{k=1}^{\infty} \frac{x}{k(x + k)} + * @f] + */ + template<typename _Tp> + _Tp + __psi_series(const _Tp __x) + { + _Tp __sum = -__numeric_constants<_Tp>::__gamma_e() - _Tp(1) / __x; + const unsigned int __max_iter = 100000; + for (unsigned int __k = 1; __k < __max_iter; ++__k) + { + const _Tp __term = __x / (__k * (__k + __x)); + __sum += __term; + if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon()) + break; + } + return __sum; + } + + + /** + * @brief Return the digamma function for large argument. + * The digamma or @f$ \psi(x) @f$ function is defined by + * @f[ + * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)} + * @f] + * + * The asymptotic series is given by: + * @f[ + * \psi(x) = \ln(x) - \frac{1}{2x} + * - \sum_{n=1}^{\infty} \frac{B_{2n}}{2 n x^{2n}} + * @f] + */ + template<typename _Tp> + _Tp + __psi_asymp(const _Tp __x) + { + _Tp __sum = std::log(__x) - _Tp(0.5L) / __x; + const _Tp __xx = __x * __x; + _Tp __xp = __xx; + const unsigned int __max_iter = 100; + for (unsigned int __k = 1; __k < __max_iter; ++__k) + { + const _Tp __term = __bernoulli<_Tp>(2 * __k) / (2 * __k * __xp); + __sum -= __term; + if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon()) + break; + __xp *= __xx; + } + return __sum; + } + + + /** + * @brief Return the digamma function. + * The digamma or @f$ \psi(x) @f$ function is defined by + * @f[ + * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)} + * @f] + * For negative argument the reflection formula is used: + * @f[ + * \psi(x) = \psi(1-x) - \pi \cot(\pi x) + * @f] + */ + template<typename _Tp> + _Tp + __psi(const _Tp __x) + { + const int __n = static_cast<int>(__x + 0.5L); + const _Tp __eps = _Tp(4) * std::numeric_limits<_Tp>::epsilon(); + if (__n <= 0 && std::abs(__x - _Tp(__n)) < __eps) + return std::numeric_limits<_Tp>::quiet_NaN(); + else if (__x < _Tp(0)) + { + const _Tp __pi = __numeric_constants<_Tp>::__pi(); + return __psi(_Tp(1) - __x) + - __pi * std::cos(__pi * __x) / std::sin(__pi * __x); + } + else if (__x > _Tp(100)) + return __psi_asymp(__x); + else + return __psi_series(__x); + } + + + /** + * @brief Return the polygamma function @f$ \psi^{(n)}(x) @f$. + * + * The polygamma function is related to the Hurwitz zeta function: + * @f[ + * \psi^{(n)}(x) = (-1)^{n+1} m! \zeta(m+1,x) + * @f] + */ + template<typename _Tp> + _Tp + __psi(const unsigned int __n, const _Tp __x) + { + if (__x <= _Tp(0)) + std::__throw_domain_error(__N("Argument out of range " + "in __psi")); + else if (__n == 0) + return __psi(__x); + else + { + const _Tp __hzeta = __hurwitz_zeta(_Tp(__n + 1), __x); +#if _GLIBCXX_USE_C99_MATH_TR1 + const _Tp __ln_nfact = std::tr1::lgamma(_Tp(__n + 1)); +#else + const _Tp __ln_nfact = __log_gamma(_Tp(__n + 1)); +#endif + _Tp __result = std::exp(__ln_nfact) * __hzeta; + if (__n % 2 == 1) + __result = -__result; + return __result; + } + } + + } // namespace std::tr1::__detail +} +} + +#endif // _TR1_GAMMA_TCC + |