contrib/gcc-4.4/libstdc++-v3/include/tr1/gamma.tcc

   1 // Special functions -*- C++ -*-
   2
   3 // Copyright (C) 2006, 2007, 2008, 2009
   4 // Free Software Foundation, Inc.
   5 //
   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)
  10 // any later version.
  11 //
  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.
  16 //
  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.
  20
  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/>.
  25
  26 /** @file tr1/gamma.tcc
  27  *  This is an internal header file, included by other library headers.
  28  *  You should not attempt to use it directly.
  29  */
  30
  31 //
  32 // ISO C++ 14882 TR1: 5.2  Special functions
  33 //
  34
  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,
  45 //       Princeton, 2003.
  46
  47 #ifndef _TR1_GAMMA_TCC
  48 #define _TR1_GAMMA_TCC 1
  49
  50 #include "special_function_util.h"
  51
  52 namespace std
  53 {
  54 namespace tr1
  55 {
  56   // Implementation-space details.
  57   namespace __detail
  58   {
  59
  60     /**
  61      *   @brief This returns Bernoulli numbers from a table or by summation
  62      *          for larger values.
  63      *
  64      *   Recursion is unstable.
  65      *
  66      *   @param __n the order n of the Bernoulli number.
  67      *   @return  The Bernoulli number of order n.
  68      */
  69     template <typename _Tp>
  70     _Tp __bernoulli_series(unsigned int __n)
  71     {
  72
  73       static const _Tp __num[28] = {
  74         _Tp(1UL),                        -_Tp(1UL) / _Tp(2UL),
  75         _Tp(1UL) / _Tp(6UL),             _Tp(0UL),
  76         -_Tp(1UL) / _Tp(30UL),           _Tp(0UL),
  77         _Tp(1UL) / _Tp(42UL),            _Tp(0UL),
  78         -_Tp(1UL) / _Tp(30UL),           _Tp(0UL),
  79         _Tp(5UL) / _Tp(66UL),            _Tp(0UL),
  80         -_Tp(691UL) / _Tp(2730UL),       _Tp(0UL),
  81         _Tp(7UL) / _Tp(6UL),             _Tp(0UL),
  82         -_Tp(3617UL) / _Tp(510UL),       _Tp(0UL),
  83         _Tp(43867UL) / _Tp(798UL),       _Tp(0UL),
  84         -_Tp(174611) / _Tp(330UL),       _Tp(0UL),
  85         _Tp(854513UL) / _Tp(138UL),      _Tp(0UL),
  86         -_Tp(236364091UL) / _Tp(2730UL), _Tp(0UL),
  87         _Tp(8553103UL) / _Tp(6UL),       _Tp(0UL)
  88       };
  89
  90       if (__n == 0)
  91         return _Tp(1);
  92
  93       if (__n == 1)
  94         return -_Tp(1) / _Tp(2);
  95
  96       //  Take care of the rest of the odd ones.
  97       if (__n % 2 == 1)
  98         return _Tp(0);
  99
 100       //  Take care of some small evens that are painful for the series.
 101       if (__n < 28)
 102         return __num[__n];
 103
 104
 105       _Tp __fact = _Tp(1);
 106       if ((__n / 2) % 2 == 0)
 107         __fact *= _Tp(-1);
 108       for (unsigned int __k = 1; __k <= __n; ++__k)
 109         __fact *= __k / (_Tp(2) * __numeric_constants<_Tp>::__pi());
 110       __fact *= _Tp(2);
 111
 112       _Tp __sum = _Tp(0);
 113       for (unsigned int __i = 1; __i < 1000; ++__i)
 114         {
 115           _Tp __term = std::pow(_Tp(__i), -_Tp(__n));
 116           if (__term < std::numeric_limits<_Tp>::epsilon())
 117             break;
 118           __sum += __term;
 119         }
 120
 121       return __fact * __sum;
 122     }
 123
 124
 125     /**
 126      *   @brief This returns Bernoulli number \f$B_n\f$.
 127      *
 128      *   @param __n the order n of the Bernoulli number.
 129      *   @return  The Bernoulli number of order n.
 130      */
 131     template<typename _Tp>
 132     inline _Tp
 133     __bernoulli(const int __n)
 134     {
 135       return __bernoulli_series<_Tp>(__n);
 136     }
 137
 138
 139     /**
 140      *   @brief Return \f$log(\Gamma(x))\f$ by asymptotic expansion
 141      *          with Bernoulli number coefficients.  This is like
 142      *          Sterling's approximation.
 143      *
 144      *   @param __x The argument of the log of the gamma function.
 145      *   @return  The logarithm of the gamma function.
 146      */
 147     template<typename _Tp>
 148     _Tp
 149     __log_gamma_bernoulli(const _Tp __x)
 150     {
 151       _Tp __lg = (__x - _Tp(0.5L)) * std::log(__x) - __x
 152                + _Tp(0.5L) * std::log(_Tp(2)
 153                * __numeric_constants<_Tp>::__pi());
 154
 155       const _Tp __xx = __x * __x;
 156       _Tp __help = _Tp(1) / __x;
 157       for ( unsigned int __i = 1; __i < 20; ++__i )
 158         {
 159           const _Tp __2i = _Tp(2 * __i);
 160           __help /= __2i * (__2i - _Tp(1)) * __xx;
 161           __lg += __bernoulli<_Tp>(2 * __i) * __help;
 162         }
 163
 164       return __lg;
 165     }
 166
 167
 168     /**
 169      *   @brief Return \f$log(\Gamma(x))\f$ by the Lanczos method.
 170      *          This method dominates all others on the positive axis I think.
 171      *
 172      *   @param __x The argument of the log of the gamma function.
 173      *   @return  The logarithm of the gamma function.
 174      */
 175     template<typename _Tp>
 176     _Tp
 177     __log_gamma_lanczos(const _Tp __x)
 178     {
 179       const _Tp __xm1 = __x - _Tp(1);
 180
 181       static const _Tp __lanczos_cheb_7[9] = {
 182        _Tp( 0.99999999999980993227684700473478L),
 183        _Tp( 676.520368121885098567009190444019L),
 184        _Tp(-1259.13921672240287047156078755283L),
 185        _Tp( 771.3234287776530788486528258894L),
 186        _Tp(-176.61502916214059906584551354L),
 187        _Tp( 12.507343278686904814458936853L),
 188        _Tp(-0.13857109526572011689554707L),
 189        _Tp( 9.984369578019570859563e-6L),
 190        _Tp( 1.50563273514931155834e-7L)
 191       };
 192
 193       static const _Tp __LOGROOT2PI
 194           = _Tp(0.9189385332046727417803297364056176L);
 195
 196       _Tp __sum = __lanczos_cheb_7[0];
 197       for(unsigned int __k = 1; __k < 9; ++__k)
 198         __sum += __lanczos_cheb_7[__k] / (__xm1 + __k);
 199
 200       const _Tp __term1 = (__xm1 + _Tp(0.5L))
 201                         * std::log((__xm1 + _Tp(7.5L))
 202                        / __numeric_constants<_Tp>::__euler());
 203       const _Tp __term2 = __LOGROOT2PI + std::log(__sum);
 204       const _Tp __result = __term1 + (__term2 - _Tp(7));
 205
 206       return __result;
 207     }
 208
 209
 210     /**
 211      *   @brief Return \f$ log(|\Gamma(x)|) \f$.
 212      *          This will return values even for \f$ x < 0 \f$.
 213      *          To recover the sign of \f$ \Gamma(x) \f$ for
 214      *          any argument use @a __log_gamma_sign.
 215      *
 216      *   @param __x The argument of the log of the gamma function.
 217      *   @return  The logarithm of the gamma function.
 218      */
 219     template<typename _Tp>
 220     _Tp
 221     __log_gamma(const _Tp __x)
 222     {
 223       if (__x > _Tp(0.5L))
 224         return __log_gamma_lanczos(__x);
 225       else
 226         {
 227           const _Tp __sin_fact
 228                  = std::abs(std::sin(__numeric_constants<_Tp>::__pi() * __x));
 229           if (__sin_fact == _Tp(0))
 230             std::__throw_domain_error(__N("Argument is nonpositive integer "
 231                                           "in __log_gamma"));
 232           return __numeric_constants<_Tp>::__lnpi()
 233                      - std::log(__sin_fact)
 234                      - __log_gamma_lanczos(_Tp(1) - __x);
 235         }
 236     }
 237
 238
 239     /**
 240      *   @brief Return the sign of \f$ \Gamma(x) \f$.
 241      *          At nonpositive integers zero is returned.
 242      *
 243      *   @param __x The argument of the gamma function.
 244      *   @return  The sign of the gamma function.
 245      */
 246     template<typename _Tp>
 247     _Tp
 248     __log_gamma_sign(const _Tp __x)
 249     {
 250       if (__x > _Tp(0))
 251         return _Tp(1);
 252       else
 253         {
 254           const _Tp __sin_fact
 255                   = std::sin(__numeric_constants<_Tp>::__pi() * __x);
 256           if (__sin_fact > _Tp(0))
 257             return (1);
 258           else if (__sin_fact < _Tp(0))
 259             return -_Tp(1);
 260           else
 261             return _Tp(0);
 262         }
 263     }
 264
 265
 266     /**
 267      *   @brief Return the logarithm of the binomial coefficient.
 268      *   The binomial coefficient is given by:
 269      *   @f[
 270      *   \left(  \right) = \frac{n!}{(n-k)! k!}
 271      *   @f]
 272      *
 273      *   @param __n The first argument of the binomial coefficient.
 274      *   @param __k The second argument of the binomial coefficient.
 275      *   @return  The binomial coefficient.
 276      */
 277     template<typename _Tp>
 278     _Tp
 279     __log_bincoef(const unsigned int __n, const unsigned int __k)
 280     {
 281       //  Max e exponent before overflow.
 282       static const _Tp __max_bincoeff
 283                       = std::numeric_limits<_Tp>::max_exponent10
 284                       * std::log(_Tp(10)) - _Tp(1);
 285 #if _GLIBCXX_USE_C99_MATH_TR1
 286       _Tp __coeff =  std::tr1::lgamma(_Tp(1 + __n))
 287                   - std::tr1::lgamma(_Tp(1 + __k))
 288                   - std::tr1::lgamma(_Tp(1 + __n - __k));
 289 #else
 290       _Tp __coeff =  __log_gamma(_Tp(1 + __n))
 291                   - __log_gamma(_Tp(1 + __k))
 292                   - __log_gamma(_Tp(1 + __n - __k));
 293 #endif
 294     }
 295
 296
 297     /**
 298      *   @brief Return the binomial coefficient.
 299      *   The binomial coefficient is given by:
 300      *   @f[
 301      *   \left(  \right) = \frac{n!}{(n-k)! k!}
 302      *   @f]
 303      *
 304      *   @param __n The first argument of the binomial coefficient.
 305      *   @param __k The second argument of the binomial coefficient.
 306      *   @return  The binomial coefficient.
 307      */
 308     template<typename _Tp>
 309     _Tp
 310     __bincoef(const unsigned int __n, const unsigned int __k)
 311     {
 312       //  Max e exponent before overflow.
 313       static const _Tp __max_bincoeff
 314                       = std::numeric_limits<_Tp>::max_exponent10
 315                       * std::log(_Tp(10)) - _Tp(1);
 316
 317       const _Tp __log_coeff = __log_bincoef<_Tp>(__n, __k);
 318       if (__log_coeff > __max_bincoeff)
 319         return std::numeric_limits<_Tp>::quiet_NaN();
 320       else
 321         return std::exp(__log_coeff);
 322     }
 323
 324
 325     /**
 326      *   @brief Return \f$ \Gamma(x) \f$.
 327      *
 328      *   @param __x The argument of the gamma function.
 329      *   @return  The gamma function.
 330      */
 331     template<typename _Tp>
 332     inline _Tp
 333     __gamma(const _Tp __x)
 334     {
 335       return std::exp(__log_gamma(__x));
 336     }
 337
 338
 339     /**
 340      *   @brief  Return the digamma function by series expansion.
 341      *   The digamma or @f$ \psi(x) @f$ function is defined by
 342      *   @f[
 343      *     \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}
 344      *   @f]
 345      *
 346      *   The series is given by:
 347      *   @f[
 348      *     \psi(x) = -\gamma_E - \frac{1}{x}
 349      *              \sum_{k=1}^{\infty} \frac{x}{k(x + k)}
 350      *   @f]
 351      */
 352     template<typename _Tp>
 353     _Tp
 354     __psi_series(const _Tp __x)
 355     {
 356       _Tp __sum = -__numeric_constants<_Tp>::__gamma_e() - _Tp(1) / __x;
 357       const unsigned int __max_iter = 100000;
 358       for (unsigned int __k = 1; __k < __max_iter; ++__k)
 359         {
 360           const _Tp __term = __x / (__k * (__k + __x));
 361           __sum += __term;
 362           if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon())
 363             break;
 364         }
 365       return __sum;
 366     }
 367
 368
 369     /**
 370      *   @brief  Return the digamma function for large argument.
 371      *   The digamma or @f$ \psi(x) @f$ function is defined by
 372      *   @f[
 373      *     \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}
 374      *   @f]
 375      *
 376      *   The asymptotic series is given by:
 377      *   @f[
 378      *     \psi(x) = \ln(x) - \frac{1}{2x}
 379      *             - \sum_{n=1}^{\infty} \frac{B_{2n}}{2 n x^{2n}}
 380      *   @f]
 381      */
 382     template<typename _Tp>
 383     _Tp
 384     __psi_asymp(const _Tp __x)
 385     {
 386       _Tp __sum = std::log(__x) - _Tp(0.5L) / __x;
 387       const _Tp __xx = __x * __x;
 388       _Tp __xp = __xx;
 389       const unsigned int __max_iter = 100;
 390       for (unsigned int __k = 1; __k < __max_iter; ++__k)
 391         {
 392           const _Tp __term = __bernoulli<_Tp>(2 * __k) / (2 * __k * __xp);
 393           __sum -= __term;
 394           if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon())
 395             break;
 396           __xp *= __xx;
 397         }
 398       return __sum;
 399     }
 400
 401
 402     /**
 403      *   @brief  Return the digamma function.
 404      *   The digamma or @f$ \psi(x) @f$ function is defined by
 405      *   @f[
 406      *     \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}
 407      *   @f]
 408      *   For negative argument the reflection formula is used:
 409      *   @f[
 410      *     \psi(x) = \psi(1-x) - \pi \cot(\pi x)
 411      *   @f]
 412      */
 413     template<typename _Tp>
 414     _Tp
 415     __psi(const _Tp __x)
 416     {
 417       const int __n = static_cast<int>(__x + 0.5L);
 418       const _Tp __eps = _Tp(4) * std::numeric_limits<_Tp>::epsilon();
 419       if (__n <= 0 && std::abs(__x - _Tp(__n)) < __eps)
 420         return std::numeric_limits<_Tp>::quiet_NaN();
 421       else if (__x < _Tp(0))
 422         {
 423           const _Tp __pi = __numeric_constants<_Tp>::__pi();
 424           return __psi(_Tp(1) - __x)
 425                - __pi * std::cos(__pi * __x) / std::sin(__pi * __x);
 426         }
 427       else if (__x > _Tp(100))
 428         return __psi_asymp(__x);
 429       else
 430         return __psi_series(__x);
 431     }
 432
 433
 434     /**
 435      *   @brief  Return the polygamma function @f$ \psi^{(n)}(x) @f$.
 436      *
 437      *   The polygamma function is related to the Hurwitz zeta function:
 438      *   @f[
 439      *     \psi^{(n)}(x) = (-1)^{n+1} m! \zeta(m+1,x)
 440      *   @f]
 441      */
 442     template<typename _Tp>
 443     _Tp
 444     __psi(const unsigned int __n, const _Tp __x)
 445     {
 446       if (__x <= _Tp(0))
 447         std::__throw_domain_error(__N("Argument out of range "
 448                                       "in __psi"));
 449       else if (__n == 0)
 450         return __psi(__x);
 451       else
 452         {
 453           const _Tp __hzeta = __hurwitz_zeta(_Tp(__n + 1), __x);
 454 #if _GLIBCXX_USE_C99_MATH_TR1
 455           const _Tp __ln_nfact = std::tr1::lgamma(_Tp(__n + 1));
 456 #else
 457           const _Tp __ln_nfact = __log_gamma(_Tp(__n + 1));
 458 #endif
 459           _Tp __result = std::exp(__ln_nfact) * __hzeta;
 460           if (__n % 2 == 1)
 461             __result = -__result;
 462           return __result;
 463         }
 464     }
 465
 466   } // namespace std::tr1::__detail
 467 }
 468 }
 469
 470 #endif // _TR1_GAMMA_TCC
 471