1 // Special functions -*- C++ -*- 2 3 // Copyright (C) 2006-2014 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 "special_function_util.h" 50 51 namespace std _GLIBCXX_VISIBILITY(default) 52 { 53 namespace tr1 54 { 55 // Implementation-space details. 56 namespace __detail 57 { 58 _GLIBCXX_BEGIN_NAMESPACE_VERSION 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 __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 __bernoulli(int __n)134 __bernoulli(int __n) 135 { return __bernoulli_series<_Tp>(__n); } 136 137 138 /** 139 * @brief Return \f$log(\Gamma(x))\f$ by asymptotic expansion 140 * with Bernoulli number coefficients. This is like 141 * Sterling's approximation. 142 * 143 * @param __x The argument of the log of the gamma function. 144 * @return The logarithm of the gamma function. 145 */ 146 template<typename _Tp> 147 _Tp __log_gamma_bernoulli(_Tp __x)148 __log_gamma_bernoulli(_Tp __x) 149 { 150 _Tp __lg = (__x - _Tp(0.5L)) * std::log(__x) - __x 151 + _Tp(0.5L) * std::log(_Tp(2) 152 * __numeric_constants<_Tp>::__pi()); 153 154 const _Tp __xx = __x * __x; 155 _Tp __help = _Tp(1) / __x; 156 for ( unsigned int __i = 1; __i < 20; ++__i ) 157 { 158 const _Tp __2i = _Tp(2 * __i); 159 __help /= __2i * (__2i - _Tp(1)) * __xx; 160 __lg += __bernoulli<_Tp>(2 * __i) * __help; 161 } 162 163 return __lg; 164 } 165 166 167 /** 168 * @brief Return \f$log(\Gamma(x))\f$ by the Lanczos method. 169 * This method dominates all others on the positive axis I think. 170 * 171 * @param __x The argument of the log of the gamma function. 172 * @return The logarithm of the gamma function. 173 */ 174 template<typename _Tp> 175 _Tp __log_gamma_lanczos(_Tp __x)176 __log_gamma_lanczos(_Tp __x) 177 { 178 const _Tp __xm1 = __x - _Tp(1); 179 180 static const _Tp __lanczos_cheb_7[9] = { 181 _Tp( 0.99999999999980993227684700473478L), 182 _Tp( 676.520368121885098567009190444019L), 183 _Tp(-1259.13921672240287047156078755283L), 184 _Tp( 771.3234287776530788486528258894L), 185 _Tp(-176.61502916214059906584551354L), 186 _Tp( 12.507343278686904814458936853L), 187 _Tp(-0.13857109526572011689554707L), 188 _Tp( 9.984369578019570859563e-6L), 189 _Tp( 1.50563273514931155834e-7L) 190 }; 191 192 static const _Tp __LOGROOT2PI 193 = _Tp(0.9189385332046727417803297364056176L); 194 195 _Tp __sum = __lanczos_cheb_7[0]; 196 for(unsigned int __k = 1; __k < 9; ++__k) 197 __sum += __lanczos_cheb_7[__k] / (__xm1 + __k); 198 199 const _Tp __term1 = (__xm1 + _Tp(0.5L)) 200 * std::log((__xm1 + _Tp(7.5L)) 201 / __numeric_constants<_Tp>::__euler()); 202 const _Tp __term2 = __LOGROOT2PI + std::log(__sum); 203 const _Tp __result = __term1 + (__term2 - _Tp(7)); 204 205 return __result; 206 } 207 208 209 /** 210 * @brief Return \f$ log(|\Gamma(x)|) \f$. 211 * This will return values even for \f$ x < 0 \f$. 212 * To recover the sign of \f$ \Gamma(x) \f$ for 213 * any argument use @a __log_gamma_sign. 214 * 215 * @param __x The argument of the log of the gamma function. 216 * @return The logarithm of the gamma function. 217 */ 218 template<typename _Tp> 219 _Tp __log_gamma(_Tp __x)220 __log_gamma(_Tp __x) 221 { 222 if (__x > _Tp(0.5L)) 223 return __log_gamma_lanczos(__x); 224 else 225 { 226 const _Tp __sin_fact 227 = std::abs(std::sin(__numeric_constants<_Tp>::__pi() * __x)); 228 if (__sin_fact == _Tp(0)) 229 std::__throw_domain_error(__N("Argument is nonpositive integer " 230 "in __log_gamma")); 231 return __numeric_constants<_Tp>::__lnpi() 232 - std::log(__sin_fact) 233 - __log_gamma_lanczos(_Tp(1) - __x); 234 } 235 } 236 237 238 /** 239 * @brief Return the sign of \f$ \Gamma(x) \f$. 240 * At nonpositive integers zero is returned. 241 * 242 * @param __x The argument of the gamma function. 243 * @return The sign of the gamma function. 244 */ 245 template<typename _Tp> 246 _Tp __log_gamma_sign(_Tp __x)247 __log_gamma_sign(_Tp __x) 248 { 249 if (__x > _Tp(0)) 250 return _Tp(1); 251 else 252 { 253 const _Tp __sin_fact 254 = std::sin(__numeric_constants<_Tp>::__pi() * __x); 255 if (__sin_fact > _Tp(0)) 256 return (1); 257 else if (__sin_fact < _Tp(0)) 258 return -_Tp(1); 259 else 260 return _Tp(0); 261 } 262 } 263 264 265 /** 266 * @brief Return the logarithm of the binomial coefficient. 267 * The binomial coefficient is given by: 268 * @f[ 269 * \left( \right) = \frac{n!}{(n-k)! k!} 270 * @f] 271 * 272 * @param __n The first argument of the binomial coefficient. 273 * @param __k The second argument of the binomial coefficient. 274 * @return The binomial coefficient. 275 */ 276 template<typename _Tp> 277 _Tp __log_bincoef(unsigned int __n,unsigned int __k)278 __log_bincoef(unsigned int __n, unsigned int __k) 279 { 280 // Max e exponent before overflow. 281 static const _Tp __max_bincoeff 282 = std::numeric_limits<_Tp>::max_exponent10 283 * std::log(_Tp(10)) - _Tp(1); 284 #if _GLIBCXX_USE_C99_MATH_TR1 285 _Tp __coeff = std::tr1::lgamma(_Tp(1 + __n)) 286 - std::tr1::lgamma(_Tp(1 + __k)) 287 - std::tr1::lgamma(_Tp(1 + __n - __k)); 288 #else 289 _Tp __coeff = __log_gamma(_Tp(1 + __n)) 290 - __log_gamma(_Tp(1 + __k)) 291 - __log_gamma(_Tp(1 + __n - __k)); 292 #endif 293 } 294 295 296 /** 297 * @brief Return the binomial coefficient. 298 * The binomial coefficient is given by: 299 * @f[ 300 * \left( \right) = \frac{n!}{(n-k)! k!} 301 * @f] 302 * 303 * @param __n The first argument of the binomial coefficient. 304 * @param __k The second argument of the binomial coefficient. 305 * @return The binomial coefficient. 306 */ 307 template<typename _Tp> 308 _Tp __bincoef(unsigned int __n,unsigned int __k)309 __bincoef(unsigned int __n, unsigned int __k) 310 { 311 // Max e exponent before overflow. 312 static const _Tp __max_bincoeff 313 = std::numeric_limits<_Tp>::max_exponent10 314 * std::log(_Tp(10)) - _Tp(1); 315 316 const _Tp __log_coeff = __log_bincoef<_Tp>(__n, __k); 317 if (__log_coeff > __max_bincoeff) 318 return std::numeric_limits<_Tp>::quiet_NaN(); 319 else 320 return std::exp(__log_coeff); 321 } 322 323 324 /** 325 * @brief Return \f$ \Gamma(x) \f$. 326 * 327 * @param __x The argument of the gamma function. 328 * @return The gamma function. 329 */ 330 template<typename _Tp> 331 inline _Tp __gamma(_Tp __x)332 __gamma(_Tp __x) 333 { return std::exp(__log_gamma(__x)); } 334 335 336 /** 337 * @brief Return the digamma function by series expansion. 338 * The digamma or @f$ \psi(x) @f$ function is defined by 339 * @f[ 340 * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)} 341 * @f] 342 * 343 * The series is given by: 344 * @f[ 345 * \psi(x) = -\gamma_E - \frac{1}{x} 346 * \sum_{k=1}^{\infty} \frac{x}{k(x + k)} 347 * @f] 348 */ 349 template<typename _Tp> 350 _Tp __psi_series(_Tp __x)351 __psi_series(_Tp __x) 352 { 353 _Tp __sum = -__numeric_constants<_Tp>::__gamma_e() - _Tp(1) / __x; 354 const unsigned int __max_iter = 100000; 355 for (unsigned int __k = 1; __k < __max_iter; ++__k) 356 { 357 const _Tp __term = __x / (__k * (__k + __x)); 358 __sum += __term; 359 if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon()) 360 break; 361 } 362 return __sum; 363 } 364 365 366 /** 367 * @brief Return the digamma function for large argument. 368 * The digamma or @f$ \psi(x) @f$ function is defined by 369 * @f[ 370 * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)} 371 * @f] 372 * 373 * The asymptotic series is given by: 374 * @f[ 375 * \psi(x) = \ln(x) - \frac{1}{2x} 376 * - \sum_{n=1}^{\infty} \frac{B_{2n}}{2 n x^{2n}} 377 * @f] 378 */ 379 template<typename _Tp> 380 _Tp __psi_asymp(_Tp __x)381 __psi_asymp(_Tp __x) 382 { 383 _Tp __sum = std::log(__x) - _Tp(0.5L) / __x; 384 const _Tp __xx = __x * __x; 385 _Tp __xp = __xx; 386 const unsigned int __max_iter = 100; 387 for (unsigned int __k = 1; __k < __max_iter; ++__k) 388 { 389 const _Tp __term = __bernoulli<_Tp>(2 * __k) / (2 * __k * __xp); 390 __sum -= __term; 391 if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon()) 392 break; 393 __xp *= __xx; 394 } 395 return __sum; 396 } 397 398 399 /** 400 * @brief Return the digamma function. 401 * The digamma or @f$ \psi(x) @f$ function is defined by 402 * @f[ 403 * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)} 404 * @f] 405 * For negative argument the reflection formula is used: 406 * @f[ 407 * \psi(x) = \psi(1-x) - \pi \cot(\pi x) 408 * @f] 409 */ 410 template<typename _Tp> 411 _Tp __psi(_Tp __x)412 __psi(_Tp __x) 413 { 414 const int __n = static_cast<int>(__x + 0.5L); 415 const _Tp __eps = _Tp(4) * std::numeric_limits<_Tp>::epsilon(); 416 if (__n <= 0 && std::abs(__x - _Tp(__n)) < __eps) 417 return std::numeric_limits<_Tp>::quiet_NaN(); 418 else if (__x < _Tp(0)) 419 { 420 const _Tp __pi = __numeric_constants<_Tp>::__pi(); 421 return __psi(_Tp(1) - __x) 422 - __pi * std::cos(__pi * __x) / std::sin(__pi * __x); 423 } 424 else if (__x > _Tp(100)) 425 return __psi_asymp(__x); 426 else 427 return __psi_series(__x); 428 } 429 430 431 /** 432 * @brief Return the polygamma function @f$ \psi^{(n)}(x) @f$. 433 * 434 * The polygamma function is related to the Hurwitz zeta function: 435 * @f[ 436 * \psi^{(n)}(x) = (-1)^{n+1} m! \zeta(m+1,x) 437 * @f] 438 */ 439 template<typename _Tp> 440 _Tp __psi(unsigned int __n,_Tp __x)441 __psi(unsigned int __n, _Tp __x) 442 { 443 if (__x <= _Tp(0)) 444 std::__throw_domain_error(__N("Argument out of range " 445 "in __psi")); 446 else if (__n == 0) 447 return __psi(__x); 448 else 449 { 450 const _Tp __hzeta = __hurwitz_zeta(_Tp(__n + 1), __x); 451 #if _GLIBCXX_USE_C99_MATH_TR1 452 const _Tp __ln_nfact = std::tr1::lgamma(_Tp(__n + 1)); 453 #else 454 const _Tp __ln_nfact = __log_gamma(_Tp(__n + 1)); 455 #endif 456 _Tp __result = std::exp(__ln_nfact) * __hzeta; 457 if (__n % 2 == 1) 458 __result = -__result; 459 return __result; 460 } 461 } 462 463 _GLIBCXX_END_NAMESPACE_VERSION 464 } // namespace std::tr1::__detail 465 } 466 } 467 468 #endif // _GLIBCXX_TR1_GAMMA_TCC 469 470