Where Online Learning is simpler!
The C and C++ Include Header Files
cat -n /usr/include/c++/15/tr1/gamma.tcc
1 // Special functions -*- C++ -*- 2 3 // Copyright (C) 2006-2025 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
