gamma.tcc
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
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
00057 namespace __detail
00058 {
00059
00060
00061
00062
00063
00064
00065
00066
00067
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
00097 if (__n % 2 == 1)
00098 return _Tp(0);
00099
00100
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
00127
00128
00129
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
00141
00142
00143
00144
00145
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
00170
00171
00172
00173
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
00212
00213
00214
00215
00216
00217
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
00241
00242
00243
00244
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
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277 template<typename _Tp>
00278 _Tp
00279 __log_bincoef(const unsigned int __n, const unsigned int __k)
00280 {
00281
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
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308 template<typename _Tp>
00309 _Tp
00310 __bincoef(const unsigned int __n, const unsigned int __k)
00311 {
00312
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
00327
00328
00329
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
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
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
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
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
00404
00405
00406
00407
00408
00409
00410
00411
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
00436
00437
00438
00439
00440
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 }
00467 }
00468 }
00469
00470 #endif // _TR1_GAMMA_TCC
00471