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 474
