libstdc++-v3/include/tr1/gamma.tcc

   1 // Special functions -*- C++ -*-
   2
   3 // Copyright (C) 2006, 2007, 2008, 2009, 2010
   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  *  Do not attempt to use it directly. @headername{tr1/cmath}
  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 _GLIBCXX_TR1_GAMMA_TCC
  48 #define _GLIBCXX_TR1_GAMMA_TCC 1
  49
  50 #include "special_function_util.h"
  51
  52 namespace std _GLIBCXX_VISIBILITY(default)
  53 {
  54 namespace tr1
  55 {
  56   // Implementation-space details.
  57   namespace __detail
  58   {
  59   _GLIBCXX_BEGIN_NAMESPACE_VERSION
  60
  61     /**
  62      *   @brief This returns Bernoulli numbers from a table or by summation
  63      *          for larger values.
  64      *
  65      *   Recursion is unstable.
  66      *
  67      *   @param __n the order n of the Bernoulli number.
  68      *   @return  The Bernoulli number of order n.
  69      */
  70     template <typename _Tp>
  71     _Tp __bernoulli_series(unsigned int __n)
  72     {
  73
  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)
  89       };
  90
  91       if (__n == 0)
  92         return _Tp(1);
  93
  94       if (__n == 1)
  95         return -_Tp(1) / _Tp(2);
  96
  97       //  Take care of the rest of the odd ones.
  98       if (__n % 2 == 1)
  99         return _Tp(0);
 100
 101       //  Take care of some small evens that are painful for the series.
 102       if (__n < 28)
 103         return __num[__n];
 104
 105
 106       _Tp __fact = _Tp(1);
 107       if ((__n / 2) % 2 == 0)
 108         __fact *= _Tp(-1);
 109       for (unsigned int __k = 1; __k <= __n; ++__k)
 110         __fact *= __k / (_Tp(2) * __numeric_constants<_Tp>::__pi());
 111       __fact *= _Tp(2);
 112
 113       _Tp __sum = _Tp(0);
 114       for (unsigned int __i = 1; __i < 1000; ++__i)
 115         {
 116           _Tp __term = std::pow(_Tp(__i), -_Tp(__n));
 117           if (__term < std::numeric_limits<_Tp>::epsilon())
 118             break;
 119           __sum += __term;
 120         }
 121
 122       return __fact * __sum;
 123     }
 124
 125
 126     /**
 127      *   @brief This returns Bernoulli number \f$B_n\f$.
 128      *
 129      *   @param __n the order n of the Bernoulli number.
 130      *   @return  The Bernoulli number of order n.
 131      */
 132     template<typename _Tp>
 133     inline _Tp
 134     __bernoulli(const int __n)
 135     {
 136       return __bernoulli_series<_Tp>(__n);
 137     }
 138
 139
 140     /**
 141      *   @brief Return \f$log(\Gamma(x))\f$ by asymptotic expansion
 142      *          with Bernoulli number coefficients.  This is like
 143      *          Sterling's approximation.
 144      *
 145      *   @param __x The argument of the log of the gamma function.
 146      *   @return  The logarithm of the gamma function.
 147      */
 148     template<typename _Tp>
 149     _Tp
 150     __log_gamma_bernoulli(const _Tp __x)
 151     {
 152       _Tp __lg = (__x - _Tp(0.5L)) * std::log(__x) - __x
 153                + _Tp(0.5L) * std::log(_Tp(2)
 154                * __numeric_constants<_Tp>::__pi());
 155
 156       const _Tp __xx = __x * __x;
 157       _Tp __help = _Tp(1) / __x;
 158       for ( unsigned int __i = 1; __i < 20; ++__i )
 159         {
 160           const _Tp __2i = _Tp(2 * __i);
 161           __help /= __2i * (__2i - _Tp(1)) * __xx;
 162           __lg += __bernoulli<_Tp>(2 * __i) * __help;
 163         }
 164
 165       return __lg;
 166     }
 167
 168
 169     /**
 170      *   @brief Return \f$log(\Gamma(x))\f$ by the Lanczos method.
 171      *          This method dominates all others on the positive axis I think.
 172      *
 173      *   @param __x The argument of the log of the gamma function.
 174      *   @return  The logarithm of the gamma function.
 175      */
 176     template<typename _Tp>
 177     _Tp
 178     __log_gamma_lanczos(const _Tp __x)
 179     {
 180       const _Tp __xm1 = __x - _Tp(1);
 181
 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)
 192       };
 193
 194       static const _Tp __LOGROOT2PI
 195           = _Tp(0.9189385332046727417803297364056176L);
 196
 197       _Tp __sum = __lanczos_cheb_7[0];
 198       for(unsigned int __k = 1; __k < 9; ++__k)
 199         __sum += __lanczos_cheb_7[__k] / (__xm1 + __k);
 200
 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));
 206
 207       return __result;
 208     }
 209
 210
 211     /**
 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.
 216      *
 217      *   @param __x The argument of the log of the gamma function.
 218      *   @return  The logarithm of the gamma function.
 219      */
 220     template<typename _Tp>
 221     _Tp
 222     __log_gamma(const _Tp __x)
 223     {
 224       if (__x > _Tp(0.5L))
 225         return __log_gamma_lanczos(__x);
 226       else
 227         {
 228           const _Tp __sin_fact
 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 "
 232                                           "in __log_gamma"));
 233           return __numeric_constants<_Tp>::__lnpi()
 234                      - std::log(__sin_fact)
 235                      - __log_gamma_lanczos(_Tp(1) - __x);
 236         }
 237     }
 238
 239
 240     /**
 241      *   @brief Return the sign of \f$ \Gamma(x) \f$.
 242      *          At nonpositive integers zero is returned.
 243      *
 244      *   @param __x The argument of the gamma function.
 245      *   @return  The sign of the gamma function.
 246      */
 247     template<typename _Tp>
 248     _Tp
 249     __log_gamma_sign(const _Tp __x)
 250     {
 251       if (__x > _Tp(0))
 252         return _Tp(1);
 253       else
 254         {
 255           const _Tp __sin_fact
 256                   = std::sin(__numeric_constants<_Tp>::__pi() * __x);
 257           if (__sin_fact > _Tp(0))
 258             return (1);
 259           else if (__sin_fact < _Tp(0))
 260             return -_Tp(1);
 261           else
 262             return _Tp(0);
 263         }
 264     }
 265
 266
 267     /**
 268      *   @brief Return the logarithm of the binomial coefficient.
 269      *   The binomial coefficient is given by:
 270      *   @f[
 271      *   \left(  \right) = \frac{n!}{(n-k)! k!}
 272      *   @f]
 273      *
 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.
 277      */
 278     template<typename _Tp>
 279     _Tp
 280     __log_bincoef(const unsigned int __n, const unsigned int __k)
 281     {
 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));
 290 #else
 291       _Tp __coeff =  __log_gamma(_Tp(1 + __n))
 292                   - __log_gamma(_Tp(1 + __k))
 293                   - __log_gamma(_Tp(1 + __n - __k));
 294 #endif
 295     }
 296
 297
 298     /**
 299      *   @brief Return the binomial coefficient.
 300      *   The binomial coefficient is given by:
 301      *   @f[
 302      *   \left(  \right) = \frac{n!}{(n-k)! k!}
 303      *   @f]
 304      *
 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.
 308      */
 309     template<typename _Tp>
 310     _Tp
 311     __bincoef(const unsigned int __n, const unsigned int __k)
 312     {
 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);
 317
 318       const _Tp __log_coeff = __log_bincoef<_Tp>(__n, __k);
 319       if (__log_coeff > __max_bincoeff)
 320         return std::numeric_limits<_Tp>::quiet_NaN();
 321       else
 322         return std::exp(__log_coeff);
 323     }
 324
 325
 326     /**
 327      *   @brief Return \f$ \Gamma(x) \f$.
 328      *
 329      *   @param __x The argument of the gamma function.
 330      *   @return  The gamma function.
 331      */
 332     template<typename _Tp>
 333     inline _Tp
 334     __gamma(const _Tp __x)
 335     {
 336       return std::exp(__log_gamma(__x));
 337     }
 338
 339
 340     /**
 341      *   @brief  Return the digamma function by series expansion.
 342      *   The digamma or @f$ \psi(x) @f$ function is defined by
 343      *   @f[
 344      *     \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}
 345      *   @f]
 346      *
 347      *   The series is given by:
 348      *   @f[
 349      *     \psi(x) = -\gamma_E - \frac{1}{x}
 350      *              \sum_{k=1}^{\infty} \frac{x}{k(x + k)}
 351      *   @f]
 352      */
 353     template<typename _Tp>
 354     _Tp
 355     __psi_series(const _Tp __x)
 356     {
 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)
 360         {
 361           const _Tp __term = __x / (__k * (__k + __x));
 362           __sum += __term;
 363           if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon())
 364             break;
 365         }
 366       return __sum;
 367     }
 368
 369
 370     /**
 371      *   @brief  Return the digamma function for large argument.
 372      *   The digamma or @f$ \psi(x) @f$ function is defined by
 373      *   @f[
 374      *     \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}
 375      *   @f]
 376      *
 377      *   The asymptotic series is given by:
 378      *   @f[
 379      *     \psi(x) = \ln(x) - \frac{1}{2x}
 380      *             - \sum_{n=1}^{\infty} \frac{B_{2n}}{2 n x^{2n}}
 381      *   @f]
 382      */
 383     template<typename _Tp>
 384     _Tp
 385     __psi_asymp(const _Tp __x)
 386     {
 387       _Tp __sum = std::log(__x) - _Tp(0.5L) / __x;
 388       const _Tp __xx = __x * __x;
 389       _Tp __xp = __xx;
 390       const unsigned int __max_iter = 100;
 391       for (unsigned int __k = 1; __k < __max_iter; ++__k)
 392         {
 393           const _Tp __term = __bernoulli<_Tp>(2 * __k) / (2 * __k * __xp);
 394           __sum -= __term;
 395           if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon())
 396             break;
 397           __xp *= __xx;
 398         }
 399       return __sum;
 400     }
 401
 402
 403     /**
 404      *   @brief  Return the digamma function.
 405      *   The digamma or @f$ \psi(x) @f$ function is defined by
 406      *   @f[
 407      *     \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)}
 408      *   @f]
 409      *   For negative argument the reflection formula is used:
 410      *   @f[
 411      *     \psi(x) = \psi(1-x) - \pi \cot(\pi x)
 412      *   @f]
 413      */
 414     template<typename _Tp>
 415     _Tp
 416     __psi(const _Tp __x)
 417     {
 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))
 423         {
 424           const _Tp __pi = __numeric_constants<_Tp>::__pi();
 425           return __psi(_Tp(1) - __x)
 426                - __pi * std::cos(__pi * __x) / std::sin(__pi * __x);
 427         }
 428       else if (__x > _Tp(100))
 429         return __psi_asymp(__x);
 430       else
 431         return __psi_series(__x);
 432     }
 433
 434
 435     /**
 436      *   @brief  Return the polygamma function @f$ \psi^{(n)}(x) @f$.
 437      *
 438      *   The polygamma function is related to the Hurwitz zeta function:
 439      *   @f[
 440      *     \psi^{(n)}(x) = (-1)^{n+1} m! \zeta(m+1,x)
 441      *   @f]
 442      */
 443     template<typename _Tp>
 444     _Tp
 445     __psi(const unsigned int __n, const _Tp __x)
 446     {
 447       if (__x <= _Tp(0))
 448         std::__throw_domain_error(__N("Argument out of range "
 449                                       "in __psi"));
 450       else if (__n == 0)
 451         return __psi(__x);
 452       else
 453         {
 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));
 457 #else
 458           const _Tp __ln_nfact = __log_gamma(_Tp(__n + 1));
 459 #endif
 460           _Tp __result = std::exp(__ln_nfact) * __hzeta;
 461           if (__n % 2 == 1)
 462             __result = -__result;
 463           return __result;
 464         }
 465     }
 466
 467   _GLIBCXX_END_NAMESPACE_VERSION
 468   } // namespace std::tr1::__detail
 469 }
 470 }
 471
 472 #endif // _GLIBCXX_TR1_GAMMA_TCC
 473