libstdc++
|
00001 // Special functions -*- C++ -*- 00002 00003 // Copyright (C) 2006, 2007, 2008, 2009 00004 // Free Software Foundation, Inc. 00005 // 00006 // This file is part of the GNU ISO C++ Library. This library is free 00007 // software; you can redistribute it and/or modify it under the 00008 // terms of the GNU General Public License as published by the 00009 // Free Software Foundation; either version 3, or (at your option) 00010 // any later version. 00011 // 00012 // This library is distributed in the hope that it will be useful, 00013 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 // GNU General Public License for more details. 00016 // 00017 // Under Section 7 of GPL version 3, you are granted additional 00018 // permissions described in the GCC Runtime Library Exception, version 00019 // 3.1, as published by the Free Software Foundation. 00020 00021 // You should have received a copy of the GNU General Public License and 00022 // a copy of the GCC Runtime Library Exception along with this program; 00023 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 00024 // <http://www.gnu.org/licenses/>. 00025 00026 /** @file tr1/gamma.tcc 00027 * This is an internal header file, included by other library headers. 00028 * You should not attempt to use it directly. 00029 */ 00030 00031 // 00032 // ISO C++ 14882 TR1: 5.2 Special functions 00033 // 00034 00035 // Written by Edward Smith-Rowland based on: 00036 // (1) Handbook of Mathematical Functions, 00037 // ed. Milton Abramowitz and Irene A. Stegun, 00038 // Dover Publications, 00039 // Section 6, pp. 253-266 00040 // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl 00041 // (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky, 00042 // W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992), 00043 // 2nd ed, pp. 213-216 00044 // (4) Gamma, Exploring Euler's Constant, Julian Havil, 00045 // Princeton, 2003. 00046 00047 #ifndef _TR1_GAMMA_TCC 00048 #define _TR1_GAMMA_TCC 1 00049 00050 #include "special_function_util.h" 00051 00052 namespace std 00053 { 00054 namespace tr1 00055 { 00056 // Implementation-space details. 00057 namespace __detail 00058 { 00059 00060 /** 00061 * @brief This returns Bernoulli numbers from a table or by summation 00062 * for larger values. 00063 * 00064 * Recursion is unstable. 00065 * 00066 * @param __n the order n of the Bernoulli number. 00067 * @return The Bernoulli number of order n. 00068 */ 00069 template <typename _Tp> 00070 _Tp __bernoulli_series(unsigned int __n) 00071 { 00072 00073 static const _Tp __num[28] = { 00074 _Tp(1UL), -_Tp(1UL) / _Tp(2UL), 00075 _Tp(1UL) / _Tp(6UL), _Tp(0UL), 00076 -_Tp(1UL) / _Tp(30UL), _Tp(0UL), 00077 _Tp(1UL) / _Tp(42UL), _Tp(0UL), 00078 -_Tp(1UL) / _Tp(30UL), _Tp(0UL), 00079 _Tp(5UL) / _Tp(66UL), _Tp(0UL), 00080 -_Tp(691UL) / _Tp(2730UL), _Tp(0UL), 00081 _Tp(7UL) / _Tp(6UL), _Tp(0UL), 00082 -_Tp(3617UL) / _Tp(510UL), _Tp(0UL), 00083 _Tp(43867UL) / _Tp(798UL), _Tp(0UL), 00084 -_Tp(174611) / _Tp(330UL), _Tp(0UL), 00085 _Tp(854513UL) / _Tp(138UL), _Tp(0UL), 00086 -_Tp(236364091UL) / _Tp(2730UL), _Tp(0UL), 00087 _Tp(8553103UL) / _Tp(6UL), _Tp(0UL) 00088 }; 00089 00090 if (__n == 0) 00091 return _Tp(1); 00092 00093 if (__n == 1) 00094 return -_Tp(1) / _Tp(2); 00095 00096 // Take care of the rest of the odd ones. 00097 if (__n % 2 == 1) 00098 return _Tp(0); 00099 00100 // Take care of some small evens that are painful for the series. 00101 if (__n < 28) 00102 return __num[__n]; 00103 00104 00105 _Tp __fact = _Tp(1); 00106 if ((__n / 2) % 2 == 0) 00107 __fact *= _Tp(-1); 00108 for (unsigned int __k = 1; __k <= __n; ++__k) 00109 __fact *= __k / (_Tp(2) * __numeric_constants<_Tp>::__pi()); 00110 __fact *= _Tp(2); 00111 00112 _Tp __sum = _Tp(0); 00113 for (unsigned int __i = 1; __i < 1000; ++__i) 00114 { 00115 _Tp __term = std::pow(_Tp(__i), -_Tp(__n)); 00116 if (__term < std::numeric_limits<_Tp>::epsilon()) 00117 break; 00118 __sum += __term; 00119 } 00120 00121 return __fact * __sum; 00122 } 00123 00124 00125 /** 00126 * @brief This returns Bernoulli number \f$B_n\f$. 00127 * 00128 * @param __n the order n of the Bernoulli number. 00129 * @return The Bernoulli number of order n. 00130 */ 00131 template<typename _Tp> 00132 inline _Tp 00133 __bernoulli(const int __n) 00134 { 00135 return __bernoulli_series<_Tp>(__n); 00136 } 00137 00138 00139 /** 00140 * @brief Return \f$log(\Gamma(x))\f$ by asymptotic expansion 00141 * with Bernoulli number coefficients. This is like 00142 * Sterling's approximation. 00143 * 00144 * @param __x The argument of the log of the gamma function. 00145 * @return The logarithm of the gamma function. 00146 */ 00147 template<typename _Tp> 00148 _Tp 00149 __log_gamma_bernoulli(const _Tp __x) 00150 { 00151 _Tp __lg = (__x - _Tp(0.5L)) * std::log(__x) - __x 00152 + _Tp(0.5L) * std::log(_Tp(2) 00153 * __numeric_constants<_Tp>::__pi()); 00154 00155 const _Tp __xx = __x * __x; 00156 _Tp __help = _Tp(1) / __x; 00157 for ( unsigned int __i = 1; __i < 20; ++__i ) 00158 { 00159 const _Tp __2i = _Tp(2 * __i); 00160 __help /= __2i * (__2i - _Tp(1)) * __xx; 00161 __lg += __bernoulli<_Tp>(2 * __i) * __help; 00162 } 00163 00164 return __lg; 00165 } 00166 00167 00168 /** 00169 * @brief Return \f$log(\Gamma(x))\f$ by the Lanczos method. 00170 * This method dominates all others on the positive axis I think. 00171 * 00172 * @param __x The argument of the log of the gamma function. 00173 * @return The logarithm of the gamma function. 00174 */ 00175 template<typename _Tp> 00176 _Tp 00177 __log_gamma_lanczos(const _Tp __x) 00178 { 00179 const _Tp __xm1 = __x - _Tp(1); 00180 00181 static const _Tp __lanczos_cheb_7[9] = { 00182 _Tp( 0.99999999999980993227684700473478L), 00183 _Tp( 676.520368121885098567009190444019L), 00184 _Tp(-1259.13921672240287047156078755283L), 00185 _Tp( 771.3234287776530788486528258894L), 00186 _Tp(-176.61502916214059906584551354L), 00187 _Tp( 12.507343278686904814458936853L), 00188 _Tp(-0.13857109526572011689554707L), 00189 _Tp( 9.984369578019570859563e-6L), 00190 _Tp( 1.50563273514931155834e-7L) 00191 }; 00192 00193 static const _Tp __LOGROOT2PI 00194 = _Tp(0.9189385332046727417803297364056176L); 00195 00196 _Tp __sum = __lanczos_cheb_7[0]; 00197 for(unsigned int __k = 1; __k < 9; ++__k) 00198 __sum += __lanczos_cheb_7[__k] / (__xm1 + __k); 00199 00200 const _Tp __term1 = (__xm1 + _Tp(0.5L)) 00201 * std::log((__xm1 + _Tp(7.5L)) 00202 / __numeric_constants<_Tp>::__euler()); 00203 const _Tp __term2 = __LOGROOT2PI + std::log(__sum); 00204 const _Tp __result = __term1 + (__term2 - _Tp(7)); 00205 00206 return __result; 00207 } 00208 00209 00210 /** 00211 * @brief Return \f$ log(|\Gamma(x)|) \f$. 00212 * This will return values even for \f$ x < 0 \f$. 00213 * To recover the sign of \f$ \Gamma(x) \f$ for 00214 * any argument use @a __log_gamma_sign. 00215 * 00216 * @param __x The argument of the log of the gamma function. 00217 * @return The logarithm of the gamma function. 00218 */ 00219 template<typename _Tp> 00220 _Tp 00221 __log_gamma(const _Tp __x) 00222 { 00223 if (__x > _Tp(0.5L)) 00224 return __log_gamma_lanczos(__x); 00225 else 00226 { 00227 const _Tp __sin_fact 00228 = std::abs(std::sin(__numeric_constants<_Tp>::__pi() * __x)); 00229 if (__sin_fact == _Tp(0)) 00230 std::__throw_domain_error(__N("Argument is nonpositive integer " 00231 "in __log_gamma")); 00232 return __numeric_constants<_Tp>::__lnpi() 00233 - std::log(__sin_fact) 00234 - __log_gamma_lanczos(_Tp(1) - __x); 00235 } 00236 } 00237 00238 00239 /** 00240 * @brief Return the sign of \f$ \Gamma(x) \f$. 00241 * At nonpositive integers zero is returned. 00242 * 00243 * @param __x The argument of the gamma function. 00244 * @return The sign of the gamma function. 00245 */ 00246 template<typename _Tp> 00247 _Tp 00248 __log_gamma_sign(const _Tp __x) 00249 { 00250 if (__x > _Tp(0)) 00251 return _Tp(1); 00252 else 00253 { 00254 const _Tp __sin_fact 00255 = std::sin(__numeric_constants<_Tp>::__pi() * __x); 00256 if (__sin_fact > _Tp(0)) 00257 return (1); 00258 else if (__sin_fact < _Tp(0)) 00259 return -_Tp(1); 00260 else 00261 return _Tp(0); 00262 } 00263 } 00264 00265 00266 /** 00267 * @brief Return the logarithm of the binomial coefficient. 00268 * The binomial coefficient is given by: 00269 * @f[ 00270 * \left( \right) = \frac{n!}{(n-k)! k!} 00271 * @f] 00272 * 00273 * @param __n The first argument of the binomial coefficient. 00274 * @param __k The second argument of the binomial coefficient. 00275 * @return The binomial coefficient. 00276 */ 00277 template<typename _Tp> 00278 _Tp 00279 __log_bincoef(const unsigned int __n, const unsigned int __k) 00280 { 00281 // Max e exponent before overflow. 00282 static const _Tp __max_bincoeff 00283 = std::numeric_limits<_Tp>::max_exponent10 00284 * std::log(_Tp(10)) - _Tp(1); 00285 #if _GLIBCXX_USE_C99_MATH_TR1 00286 _Tp __coeff = std::tr1::lgamma(_Tp(1 + __n)) 00287 - std::tr1::lgamma(_Tp(1 + __k)) 00288 - std::tr1::lgamma(_Tp(1 + __n - __k)); 00289 #else 00290 _Tp __coeff = __log_gamma(_Tp(1 + __n)) 00291 - __log_gamma(_Tp(1 + __k)) 00292 - __log_gamma(_Tp(1 + __n - __k)); 00293 #endif 00294 } 00295 00296 00297 /** 00298 * @brief Return the binomial coefficient. 00299 * The binomial coefficient is given by: 00300 * @f[ 00301 * \left( \right) = \frac{n!}{(n-k)! k!} 00302 * @f] 00303 * 00304 * @param __n The first argument of the binomial coefficient. 00305 * @param __k The second argument of the binomial coefficient. 00306 * @return The binomial coefficient. 00307 */ 00308 template<typename _Tp> 00309 _Tp 00310 __bincoef(const unsigned int __n, const unsigned int __k) 00311 { 00312 // Max e exponent before overflow. 00313 static const _Tp __max_bincoeff 00314 = std::numeric_limits<_Tp>::max_exponent10 00315 * std::log(_Tp(10)) - _Tp(1); 00316 00317 const _Tp __log_coeff = __log_bincoef<_Tp>(__n, __k); 00318 if (__log_coeff > __max_bincoeff) 00319 return std::numeric_limits<_Tp>::quiet_NaN(); 00320 else 00321 return std::exp(__log_coeff); 00322 } 00323 00324 00325 /** 00326 * @brief Return \f$ \Gamma(x) \f$. 00327 * 00328 * @param __x The argument of the gamma function. 00329 * @return The gamma function. 00330 */ 00331 template<typename _Tp> 00332 inline _Tp 00333 __gamma(const _Tp __x) 00334 { 00335 return std::exp(__log_gamma(__x)); 00336 } 00337 00338 00339 /** 00340 * @brief Return the digamma function by series expansion. 00341 * The digamma or @f$ \psi(x) @f$ function is defined by 00342 * @f[ 00343 * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)} 00344 * @f] 00345 * 00346 * The series is given by: 00347 * @f[ 00348 * \psi(x) = -\gamma_E - \frac{1}{x} 00349 * \sum_{k=1}^{\infty} \frac{x}{k(x + k)} 00350 * @f] 00351 */ 00352 template<typename _Tp> 00353 _Tp 00354 __psi_series(const _Tp __x) 00355 { 00356 _Tp __sum = -__numeric_constants<_Tp>::__gamma_e() - _Tp(1) / __x; 00357 const unsigned int __max_iter = 100000; 00358 for (unsigned int __k = 1; __k < __max_iter; ++__k) 00359 { 00360 const _Tp __term = __x / (__k * (__k + __x)); 00361 __sum += __term; 00362 if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon()) 00363 break; 00364 } 00365 return __sum; 00366 } 00367 00368 00369 /** 00370 * @brief Return the digamma function for large argument. 00371 * The digamma or @f$ \psi(x) @f$ function is defined by 00372 * @f[ 00373 * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)} 00374 * @f] 00375 * 00376 * The asymptotic series is given by: 00377 * @f[ 00378 * \psi(x) = \ln(x) - \frac{1}{2x} 00379 * - \sum_{n=1}^{\infty} \frac{B_{2n}}{2 n x^{2n}} 00380 * @f] 00381 */ 00382 template<typename _Tp> 00383 _Tp 00384 __psi_asymp(const _Tp __x) 00385 { 00386 _Tp __sum = std::log(__x) - _Tp(0.5L) / __x; 00387 const _Tp __xx = __x * __x; 00388 _Tp __xp = __xx; 00389 const unsigned int __max_iter = 100; 00390 for (unsigned int __k = 1; __k < __max_iter; ++__k) 00391 { 00392 const _Tp __term = __bernoulli<_Tp>(2 * __k) / (2 * __k * __xp); 00393 __sum -= __term; 00394 if (std::abs(__term / __sum) < std::numeric_limits<_Tp>::epsilon()) 00395 break; 00396 __xp *= __xx; 00397 } 00398 return __sum; 00399 } 00400 00401 00402 /** 00403 * @brief Return the digamma function. 00404 * The digamma or @f$ \psi(x) @f$ function is defined by 00405 * @f[ 00406 * \psi(x) = \frac{\Gamma'(x)}{\Gamma(x)} 00407 * @f] 00408 * For negative argument the reflection formula is used: 00409 * @f[ 00410 * \psi(x) = \psi(1-x) - \pi \cot(\pi x) 00411 * @f] 00412 */ 00413 template<typename _Tp> 00414 _Tp 00415 __psi(const _Tp __x) 00416 { 00417 const int __n = static_cast<int>(__x + 0.5L); 00418 const _Tp __eps = _Tp(4) * std::numeric_limits<_Tp>::epsilon(); 00419 if (__n <= 0 && std::abs(__x - _Tp(__n)) < __eps) 00420 return std::numeric_limits<_Tp>::quiet_NaN(); 00421 else if (__x < _Tp(0)) 00422 { 00423 const _Tp __pi = __numeric_constants<_Tp>::__pi(); 00424 return __psi(_Tp(1) - __x) 00425 - __pi * std::cos(__pi * __x) / std::sin(__pi * __x); 00426 } 00427 else if (__x > _Tp(100)) 00428 return __psi_asymp(__x); 00429 else 00430 return __psi_series(__x); 00431 } 00432 00433 00434 /** 00435 * @brief Return the polygamma function @f$ \psi^{(n)}(x) @f$. 00436 * 00437 * The polygamma function is related to the Hurwitz zeta function: 00438 * @f[ 00439 * \psi^{(n)}(x) = (-1)^{n+1} m! \zeta(m+1,x) 00440 * @f] 00441 */ 00442 template<typename _Tp> 00443 _Tp 00444 __psi(const unsigned int __n, const _Tp __x) 00445 { 00446 if (__x <= _Tp(0)) 00447 std::__throw_domain_error(__N("Argument out of range " 00448 "in __psi")); 00449 else if (__n == 0) 00450 return __psi(__x); 00451 else 00452 { 00453 const _Tp __hzeta = __hurwitz_zeta(_Tp(__n + 1), __x); 00454 #if _GLIBCXX_USE_C99_MATH_TR1 00455 const _Tp __ln_nfact = std::tr1::lgamma(_Tp(__n + 1)); 00456 #else 00457 const _Tp __ln_nfact = __log_gamma(_Tp(__n + 1)); 00458 #endif 00459 _Tp __result = std::exp(__ln_nfact) * __hzeta; 00460 if (__n % 2 == 1) 00461 __result = -__result; 00462 return __result; 00463 } 00464 } 00465 00466 } // namespace std::tr1::__detail 00467 } 00468 } 00469 00470 #endif // _TR1_GAMMA_TCC 00471