gamma.tcc

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

Generated on Thu Jul 23 21:16:06 2009 for libstdc++ by  doxygen 1.5.8