mingw/lib/gcc/mingw32/4.3.3/include/c++/tr1/gamma.tcc

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