gamma.tcc

Go to the documentation of this file.
00001 // Special functions -*- C++ -*-
00002 
00003 // Copyright (C) 2006, 2007, 2008
00004 // Free Software Foundation, Inc.
00005 //
00006 // This file is part of the GNU ISO C++ Library.  This library is free
00007 // software; you can redistribute it and/or modify it under the
00008 // terms of the GNU General Public License as published by the
00009 // Free Software Foundation; either version 2, or (at your option)
00010 // any later version.
00011 //
00012 // This library is distributed in the hope that it will be useful,
00013 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00014 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015 // GNU General Public License for more details.
00016 //
00017 // You should have received a copy of the GNU General Public License along
00018 // with this library; see the file COPYING.  If not, write to the Free
00019 // Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
00020 // USA.
00021 //
00022 // As a special exception, you may use this file as part of a free software
00023 // library without restriction.  Specifically, if other files instantiate
00024 // templates or use macros or inline functions from this file, or you compile
00025 // this file and link it with other files to produce an executable, this
00026 // file does not by itself cause the resulting executable to be covered by
00027 // the GNU General Public License.  This exception does not however
00028 // invalidate any other reasons why the executable file might be covered by
00029 // the GNU General Public License.
00030 
00031 /** @file tr1/gamma.tcc
00032  *  This is an internal header file, included by other library headers.
00033  *  You should not attempt to use it directly.
00034  */
00035 
00036 //
00037 // ISO C++ 14882 TR1: 5.2  Special functions
00038 //
00039 
00040 // Written by Edward Smith-Rowland based on:
00041 //   (1) Handbook of Mathematical Functions,
00042 //       ed. Milton Abramowitz and Irene A. Stegun,
00043 //       Dover Publications,
00044 //       Section 6, pp. 253-266
00045 //   (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl
00046 //   (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky,
00047 //       W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992),
00048 //       2nd ed, pp. 213-216
00049 //   (4) Gamma, Exploring Euler's Constant, Julian Havil,
00050 //       Princeton, 2003.
00051 
00052 #ifndef _TR1_GAMMA_TCC
00053 #define _TR1_GAMMA_TCC 1
00054 
00055 #include "special_function_util.h"
00056 
00057 namespace std
00058 {
00059 namespace tr1
00060 {
00061   // Implementation-space details.
00062   namespace __detail
00063   {
00064 
00065     /**
00066      *   @brief This returns Bernoulli numbers from a table or by summation
00067      *          for larger values.
00068      *
00069      *   Recursion is unstable.
00070      *
00071      *   @param __n the order n of the Bernoulli number.
00072      *   @return  The Bernoulli number of order n.
00073      */
00074     template <typename _Tp>
00075     _Tp __bernoulli_series(unsigned int __n)
00076     {
00077 
00078       static const _Tp __num[28] = {
00079         _Tp(1UL),                        -_Tp(1UL) / _Tp(2UL),
00080         _Tp(1UL) / _Tp(6UL),             _Tp(0UL),
00081         -_Tp(1UL) / _Tp(30UL),           _Tp(0UL),
00082         _Tp(1UL) / _Tp(42UL),            _Tp(0UL),
00083         -_Tp(1UL) / _Tp(30UL),           _Tp(0UL),
00084         _Tp(5UL) / _Tp(66UL),            _Tp(0UL),
00085         -_Tp(691UL) / _Tp(2730UL),       _Tp(0UL),
00086         _Tp(7UL) / _Tp(6UL),             _Tp(0UL),
00087         -_Tp(3617UL) / _Tp(510UL),       _Tp(0UL),
00088         _Tp(43867UL) / _Tp(798UL),       _Tp(0UL),
00089         -_Tp(174611) / _Tp(330UL),       _Tp(0UL),
00090         _Tp(854513UL) / _Tp(138UL),      _Tp(0UL),
00091         -_Tp(236364091UL) / _Tp(2730UL), _Tp(0UL),
00092         _Tp(8553103UL) / _Tp(6UL),       _Tp(0UL)
00093       };
00094 
00095       if (__n == 0)
00096         return _Tp(1);
00097 
00098       if (__n == 1)
00099         return -_Tp(1) / _Tp(2);
00100 
00101       //  Take care of the rest of the odd ones.
00102       if (__n % 2 == 1)
00103         return _Tp(0);
00104 
00105       //  Take care of some small evens that are painful for the series.
00106       if (__n < 28)
00107         return __num[__n];
00108 
00109 
00110       _Tp __fact = _Tp(1);
00111       if ((__n / 2) % 2 == 0)
00112         __fact *= _Tp(-1);
00113       for (unsigned int __k = 1; __k <= __n; ++__k)
00114         __fact *= __k / (_Tp(2) * __numeric_constants<_Tp>::__pi());
00115       __fact *= _Tp(2);
00116 
00117       _Tp __sum = _Tp(0);
00118       for (unsigned int __i = 1; __i < 1000; ++__i)
00119         {
00120           _Tp __term = std::pow(_Tp(__i), -_Tp(__n));
00121           if (__term < std::numeric_limits<_Tp>::epsilon())
00122             break;
00123           __sum += __term;
00124         }
00125 
00126       return __fact * __sum;
00127     }
00128 
00129 
00130     /**
00131      *   @brief This returns Bernoulli number \f$B_n\f$.
00132      *
00133      *   @param __n the order n of the Bernoulli number.
00134      *   @return  The Bernoulli number of order n.
00135      */
00136     template<typename _Tp>
00137     inline _Tp
00138     __bernoulli(const int __n)
00139     {
00140       return __bernoulli_series<_Tp>(__n);
00141     }
00142 
00143 
00144     /**
00145      *   @brief Return \f$log(\Gamma(x))\f$ by asymptotic expansion
00146      *          with Bernoulli number coefficients.  This is like
00147      *          Sterling's approximation.
00148      *
00149      *   @param __x The argument of the log of the gamma function.
00150      *   @return  The logarithm of the gamma function.
00151      */
00152     template<typename _Tp>
00153     _Tp
00154     __log_gamma_bernoulli(const _Tp __x)
00155     {
00156       _Tp __lg = (__x - _Tp(0.5L)) * std::log(__x) - __x
00157                + _Tp(0.5L) * std::log(_Tp(2)
00158                * __numeric_constants<_Tp>::__pi());
00159 
00160       const _Tp __xx = __x * __x;
00161       _Tp __help = _Tp(1) / __x;
00162       for ( unsigned int __i = 1; __i < 20; ++__i )
00163         {
00164           const _Tp __2i = _Tp(2 * __i);
00165           __help /= __2i * (__2i - _Tp(1)) * __xx;
00166           __lg += __bernoulli<_Tp>(2 * __i) * __help;
00167         }
00168 
00169       return __lg;
00170     }
00171 
00172 
00173     /**
00174      *   @brief Return \f$log(\Gamma(x))\f$ by the Lanczos method.
00175      *          This method dominates all others on the positive axis I think.
00176      *
00177      *   @param __x The argument of the log of the gamma function.
00178      *   @return  The logarithm of the gamma function.
00179      */
00180     template<typename _Tp>
00181     _Tp
00182     __log_gamma_lanczos(const _Tp __x)
00183     {
00184       const _Tp __xm1 = __x - _Tp(1);
00185 
00186       static const _Tp __lanczos_cheb_7[9] = {
00187        _Tp( 0.99999999999980993227684700473478L),
00188        _Tp( 676.520368121885098567009190444019L),
00189        _Tp(-1259.13921672240287047156078755283L),
00190        _Tp( 771.3234287776530788486528258894L),
00191        _Tp(-176.61502916214059906584551354L),
00192        _Tp( 12.507343278686904814458936853L),
00193        _Tp(-0.13857109526572011689554707L),
00194        _Tp( 9.984369578019570859563e-6L),
00195        _Tp( 1.50563273514931155834e-7L)
00196       };
00197 
00198       static const _Tp __LOGROOT2PI
00199           = _Tp(0.9189385332046727417803297364056176L);
00200 
00201       _Tp __sum = __lanczos_cheb_7[0];
00202       for(unsigned int __k = 1; __k < 9; ++__k)
00203         __sum += __lanczos_cheb_7[__k] / (__xm1 + __k);
00204 
00205       const _Tp __term1 = (__xm1 + _Tp(0.5L))
00206                         * std::log((__xm1 + _Tp(7.5L))
00207                        / __numeric_constants<_Tp>::__euler());
00208       const _Tp __term2 = __LOGROOT2PI + std::log(__sum);
00209       const _Tp __result = __term1 + (__term2 - _Tp(7));
00210 
00211       return __result;
00212     }
00213 
00214 
00215     /**
00216      *   @brief Return \f$ log(|\Gamma(x)|) \f$.
00217      *          This will return values even for \f$ x < 0 \f$.
00218      *          To recover the sign of \f$ \Gamma(x) \f$ for
00219      *          any argument use @a __log_gamma_sign.
00220      *
00221      *   @param __x The argument of the log of the gamma function.
00222      *   @return  The logarithm of the gamma function.
00223      */
00224     template<typename _Tp>
00225     _Tp
00226     __log_gamma(const _Tp __x)
00227     {
00228       if (__x > _Tp(0.5L))
00229         return __log_gamma_lanczos(__x);
00230       else
00231         {
00232           const _Tp __sin_fact
00233                  = std::abs(std::sin(__numeric_constants<_Tp>::__pi() * __x));
00234           if (__sin_fact == _Tp(0))
00235             std::__throw_domain_error(__N("Argument is nonpositive integer "
00236                                           "in __log_gamma"));
00237           return __numeric_constants<_Tp>::__lnpi()
00238                      - std::log(__sin_fact)
00239                      - __log_gamma_lanczos(_Tp(1) - __x);
00240         }
00241     }
00242 
00243 
00244     /**
00245      *   @brief Return the sign of \f$ \Gamma(x) \f$.
00246      *          At nonpositive integers zero is returned.
00247      *
00248      *   @param __x The argument of the gamma function.
00249      *   @return  The sign of the gamma function.
00250      */
00251     template<typename _Tp>
00252     _Tp
00253     __log_gamma_sign(const _Tp __x)
00254     {
00255       if (__x > _Tp(0))
00256         return _Tp(1);
00257       else
00258         {
00259           const _Tp __sin_fact
00260                   = std::sin(__numeric_constants<_Tp>::__pi() * __x);
00261           if (__sin_fact > _Tp(0))
00262             return (1);
00263           else if (__sin_fact < _Tp(0))
00264             return -_Tp(1);
00265           else
00266             return _Tp(0);
00267         }
00268     }
00269 
00270 
00271     /**
00272      *   @brief Return the logarithm of the binomial coefficient.
00273      *   The binomial coefficient is given by:
00274      *   @f[
00275      *   \left(  \right) = \frac{n!}{(n-k)! k!}
00276      *   @f]
00277      *
00278      *   @param __n The first argument of the binomial coefficient.
00279      *   @param __k The second argument of the binomial coefficient.
00280      *   @return  The binomial coefficient.
00281      */
00282     template<typename _Tp>
00283     _Tp
00284     __log_bincoef(const unsigned int __n, const unsigned int __k)
00285     {
00286       //  Max e exponent before overflow.
00287       static const _Tp __max_bincoeff
00288                       = std::numeric_limits<_Tp>::max_exponent10
00289                       * std::log(_Tp(10)) - _Tp(1);
00290 #if _GLIBCXX_USE_C99_MATH_TR1
00291       _Tp __coeff =  std::tr1::lgamma(_Tp(1 + __n))
00292                   - std::tr1::lgamma(_Tp(1 + __k))
00293                   - std::tr1::lgamma(_Tp(1 + __n - __k));
00294 #else
00295       _Tp __coeff =  __log_gamma(_Tp(1 + __n))
00296                   - __log_gamma(_Tp(1 + __k))
00297                   - __log_gamma(_Tp(1 + __n - __k));
00298 #endif
00299     }
00300 
00301 
00302     /**
00303      *   @brief Return the binomial coefficient.
00304      *   The binomial coefficient is given by:
00305      *   @f[
00306      *   \left(  \right) = \frac{n!}{(n-k)! k!}
00307      *   @f]
00308      *
00309      *   @param __n The first argument of the binomial coefficient.
00310      *   @param __k The second argument of the binomial coefficient.
00311      *   @return  The binomial coefficient.
00312      */
00313     template<typename _Tp>
00314     _Tp
00315     __bincoef(const unsigned int __n, const unsigned int __k)
00316     {
00317       //  Max e exponent before overflow.
00318       static const _Tp __max_bincoeff
00319                       = std::numeric_limits<_Tp>::max_exponent10
00320                       * std::log(_Tp(10)) - _Tp(1);
00321 
00322       const _Tp __log_coeff = __log_bincoef<_Tp>(__n, __k);
00323       if (__log_coeff > __max_bincoeff)
00324         return std::numeric_limits<_Tp>::quiet_NaN();
00325       else
00326         return std::exp(__log_coeff);
00327     }
00328 
00329 
00330     /**
00331      *   @brief Return \f$ \Gamma(x) \f$.
00332      *
00333      *   @param __x The argument of the gamma function.
00334      *   @return  The gamma function.
00335      */
00336     template<typename _Tp>
00337     inline _Tp
00338     __gamma(const _Tp __x)
00339     {
00340       return std::exp(__log_gamma(__x));
00341     }
00342 
00343 
00344     /**
00345      *   @brief  Return the digamma function by series expansion.
00346      *   The digamma or @f$ \psi(x) @f$ function is defined by
00347      *   @f[
00348      *     \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}
00349      *   @f]
00350      *
00351      *   The series is given by:
00352      *   @f[
00353      *     \psi(x) = -\gamma_E - \frac{1}{x}
00354      *              \sum_{k=1}^{\infty} \frac{x}{k(x + k)}
00355      *   @f]
00356      */
00357     template<typename _Tp>
00358     _Tp
00359     __psi_series(const _Tp __x)
00360     {
00361       _Tp __sum = -__numeric_constants<_Tp>::__gamma_e() - _Tp(1) / __x;
00362       const unsigned int __max_iter = 100000;
00363       for (unsigned int __k = 1; __k < __max_iter; ++__k)
00364         {
00365           const _Tp __term = __x / (__k * (__k + __x));
00366           __sum += __term;
00367           if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon())
00368             break;
00369         }
00370       return __sum;
00371     }
00372 
00373 
00374     /**
00375      *   @brief  Return the digamma function for large argument.
00376      *   The digamma or @f$ \psi(x) @f$ function is defined by
00377      *   @f[
00378      *     \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}
00379      *   @f]
00380      *
00381      *   The asymptotic series is given by:
00382      *   @f[
00383      *     \psi(x) = \ln(x) - \frac{1}{2x}
00384      *             - \sum_{n=1}^{\infty} \frac{B_{2n}}{2 n x^{2n}}
00385      *   @f]
00386      */
00387     template<typename _Tp>
00388     _Tp
00389     __psi_asymp(const _Tp __x)
00390     {
00391       _Tp __sum = std::log(__x) - _Tp(0.5L) / __x;
00392       const _Tp __xx = __x * __x;
00393       _Tp __xp = __xx;
00394       const unsigned int __max_iter = 100;
00395       for (unsigned int __k = 1; __k < __max_iter; ++__k)
00396         {
00397           const _Tp __term = __bernoulli<_Tp>(2 * __k) / (2 * __k * __xp);
00398           __sum -= __term;
00399           if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon())
00400             break;
00401           __xp *= __xx;
00402         }
00403       return __sum;
00404     }
00405 
00406 
00407     /**
00408      *   @brief  Return the digamma function.
00409      *   The digamma or @f$ \psi(x) @f$ function is defined by
00410      *   @f[
00411      *     \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}
00412      *   @f]
00413      *   For negative argument the reflection formula is used:
00414      *   @f[
00415      *     \psi(x) = \psi(1-x) - \pi \cot(\pi x)
00416      *   @f]
00417      */
00418     template<typename _Tp>
00419     _Tp
00420     __psi(const _Tp __x)
00421     {
00422       const int __n = static_cast<int>(__x + 0.5L);
00423       const _Tp __eps = _Tp(4) * std::numeric_limits<_Tp>::epsilon();
00424       if (__n <= 0 && std::abs(__x - _Tp(__n)) < __eps)
00425         return std::numeric_limits<_Tp>::quiet_NaN();
00426       else if (__x < _Tp(0))
00427         {
00428           const _Tp __pi = __numeric_constants<_Tp>::__pi();
00429           return __psi(_Tp(1) - __x)
00430                - __pi * std::cos(__pi * __x) / std::sin(__pi * __x);
00431         }
00432       else if (__x > _Tp(100))
00433         return __psi_asymp(__x);
00434       else
00435         return __psi_series(__x);
00436     }
00437 
00438 
00439     /**
00440      *   @brief  Return the polygamma function @f$ \psi^{(n)}(x) @f$.
00441      * 
00442      *   The polygamma function is related to the Hurwitz zeta function:
00443      *   @f[
00444      *     \psi^{(n)}(x) = (-1)^{n+1} m! \zeta(m+1,x)
00445      *   @f]
00446      */
00447     template<typename _Tp>
00448     _Tp
00449     __psi(const unsigned int __n, const _Tp __x)
00450     {
00451       if (__x <= _Tp(0))
00452         std::__throw_domain_error(__N("Argument out of range "
00453                                       "in __psi"));
00454       else if (__n == 0)
00455         return __psi(__x);
00456       else
00457         {
00458           const _Tp __hzeta = __hurwitz_zeta(_Tp(__n + 1), __x);
00459 #if _GLIBCXX_USE_C99_MATH_TR1
00460           const _Tp __ln_nfact = std::tr1::lgamma(_Tp(__n + 1));
00461 #else
00462           const _Tp __ln_nfact = __log_gamma(_Tp(__n + 1));
00463 #endif
00464           _Tp __result = std::exp(__ln_nfact) * __hzeta;
00465           if (__n % 2 == 1)
00466             __result = -__result;
00467           return __result;
00468         }
00469     }
00470 
00471   } // namespace std::tr1::__detail
00472 }
00473 }
00474 
00475 #endif // _TR1_GAMMA_TCC
00476 

Generated on Sat Oct 25 05:09:02 2008 for libstdc++ by  doxygen 1.5.6