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

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