riemann_zeta.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
00048 #ifndef _GLIBCXX_TR1_RIEMANN_ZETA_TCC
00049 #define _GLIBCXX_TR1_RIEMANN_ZETA_TCC 1
00050
00051 #include "special_function_util.h"
00052
00053 namespace std
00054 {
00055 namespace tr1
00056 {
00057
00058
00059
00060
00061 namespace __detail
00062 {
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077 template<typename _Tp>
00078 _Tp
00079 __riemann_zeta_sum(const _Tp __s)
00080 {
00081
00082 if (__s < _Tp(1))
00083 std::__throw_domain_error(__N("Bad argument in zeta sum."));
00084
00085 const unsigned int max_iter = 10000;
00086 _Tp __zeta = _Tp(0);
00087 for (unsigned int __k = 1; __k < max_iter; ++__k)
00088 {
00089 _Tp __term = std::pow(static_cast<_Tp>(__k), -__s);
00090 if (__term < std::numeric_limits<_Tp>::epsilon())
00091 {
00092 break;
00093 }
00094 __zeta += __term;
00095 }
00096
00097 return __zeta;
00098 }
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114 template<typename _Tp>
00115 _Tp
00116 __riemann_zeta_alt(const _Tp __s)
00117 {
00118 _Tp __sgn = _Tp(1);
00119 _Tp __zeta = _Tp(0);
00120 for (unsigned int __i = 1; __i < 10000000; ++__i)
00121 {
00122 _Tp __term = __sgn / std::pow(__i, __s);
00123 if (std::abs(__term) < std::numeric_limits<_Tp>::epsilon())
00124 break;
00125 __zeta += __term;
00126 __sgn *= _Tp(-1);
00127 }
00128 __zeta /= _Tp(1) - std::pow(_Tp(2), _Tp(1) - __s);
00129
00130 return __zeta;
00131 }
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156 template<typename _Tp>
00157 _Tp
00158 __riemann_zeta_glob(const _Tp __s)
00159 {
00160 _Tp __zeta = _Tp(0);
00161
00162 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
00163
00164 const _Tp __max_bincoeff = std::numeric_limits<_Tp>::max_exponent10
00165 * std::log(_Tp(10)) - _Tp(1);
00166
00167
00168
00169 if (__s < _Tp(0))
00170 {
00171 #if _GLIBCXX_USE_C99_MATH_TR1
00172 if (std::tr1::fmod(__s,_Tp(2)) == _Tp(0))
00173 return _Tp(0);
00174 else
00175 #endif
00176 {
00177 _Tp __zeta = __riemann_zeta_glob(_Tp(1) - __s);
00178 __zeta *= std::pow(_Tp(2)
00179 * __numeric_constants<_Tp>::__pi(), __s)
00180 * std::sin(__numeric_constants<_Tp>::__pi_2() * __s)
00181 #if _GLIBCXX_USE_C99_MATH_TR1
00182 * std::exp(std::tr1::lgamma(_Tp(1) - __s))
00183 #else
00184 * std::exp(__log_gamma(_Tp(1) - __s))
00185 #endif
00186 / __numeric_constants<_Tp>::__pi();
00187 return __zeta;
00188 }
00189 }
00190
00191 _Tp __num = _Tp(0.5L);
00192 const unsigned int __maxit = 10000;
00193 for (unsigned int __i = 0; __i < __maxit; ++__i)
00194 {
00195 bool __punt = false;
00196 _Tp __sgn = _Tp(1);
00197 _Tp __term = _Tp(0);
00198 for (unsigned int __j = 0; __j <= __i; ++__j)
00199 {
00200 #if _GLIBCXX_USE_C99_MATH_TR1
00201 _Tp __bincoeff = std::tr1::lgamma(_Tp(1 + __i))
00202 - std::tr1::lgamma(_Tp(1 + __j))
00203 - std::tr1::lgamma(_Tp(1 + __i - __j));
00204 #else
00205 _Tp __bincoeff = __log_gamma(_Tp(1 + __i))
00206 - __log_gamma(_Tp(1 + __j))
00207 - __log_gamma(_Tp(1 + __i - __j));
00208 #endif
00209 if (__bincoeff > __max_bincoeff)
00210 {
00211
00212 __punt = true;
00213 break;
00214 }
00215 __bincoeff = std::exp(__bincoeff);
00216 __term += __sgn * __bincoeff * std::pow(_Tp(1 + __j), -__s);
00217 __sgn *= _Tp(-1);
00218 }
00219 if (__punt)
00220 break;
00221 __term *= __num;
00222 __zeta += __term;
00223 if (std::abs(__term/__zeta) < __eps)
00224 break;
00225 __num *= _Tp(0.5L);
00226 }
00227
00228 __zeta /= _Tp(1) - std::pow(_Tp(2), _Tp(1) - __s);
00229
00230 return __zeta;
00231 }
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251 template<typename _Tp>
00252 _Tp
00253 __riemann_zeta_product(const _Tp __s)
00254 {
00255 static const _Tp __prime[] = {
00256 _Tp(2), _Tp(3), _Tp(5), _Tp(7), _Tp(11), _Tp(13), _Tp(17), _Tp(19),
00257 _Tp(23), _Tp(29), _Tp(31), _Tp(37), _Tp(41), _Tp(43), _Tp(47),
00258 _Tp(53), _Tp(59), _Tp(61), _Tp(67), _Tp(71), _Tp(73), _Tp(79),
00259 _Tp(83), _Tp(89), _Tp(97), _Tp(101), _Tp(103), _Tp(107), _Tp(109)
00260 };
00261 static const unsigned int __num_primes = sizeof(__prime) / sizeof(_Tp);
00262
00263 _Tp __zeta = _Tp(1);
00264 for (unsigned int __i = 0; __i < __num_primes; ++__i)
00265 {
00266 const _Tp __fact = _Tp(1) - std::pow(__prime[__i], -__s);
00267 __zeta *= __fact;
00268 if (_Tp(1) - __fact < std::numeric_limits<_Tp>::epsilon())
00269 break;
00270 }
00271
00272 __zeta = _Tp(1) / __zeta;
00273
00274 return __zeta;
00275 }
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292 template<typename _Tp>
00293 _Tp
00294 __riemann_zeta(const _Tp __s)
00295 {
00296 if (__isnan(__s))
00297 return std::numeric_limits<_Tp>::quiet_NaN();
00298 else if (__s == _Tp(1))
00299 return std::numeric_limits<_Tp>::infinity();
00300 else if (__s < -_Tp(19))
00301 {
00302 _Tp __zeta = __riemann_zeta_product(_Tp(1) - __s);
00303 __zeta *= std::pow(_Tp(2) * __numeric_constants<_Tp>::__pi(), __s)
00304 * std::sin(__numeric_constants<_Tp>::__pi_2() * __s)
00305 #if _GLIBCXX_USE_C99_MATH_TR1
00306 * std::exp(std::tr1::lgamma(_Tp(1) - __s))
00307 #else
00308 * std::exp(__log_gamma(_Tp(1) - __s))
00309 #endif
00310 / __numeric_constants<_Tp>::__pi();
00311 return __zeta;
00312 }
00313 else if (__s < _Tp(20))
00314 {
00315
00316 bool __glob = true;
00317 if (__glob)
00318 return __riemann_zeta_glob(__s);
00319 else
00320 {
00321 if (__s > _Tp(1))
00322 return __riemann_zeta_sum(__s);
00323 else
00324 {
00325 _Tp __zeta = std::pow(_Tp(2)
00326 * __numeric_constants<_Tp>::__pi(), __s)
00327 * std::sin(__numeric_constants<_Tp>::__pi_2() * __s)
00328 #if _GLIBCXX_USE_C99_MATH_TR1
00329 * std::tr1::tgamma(_Tp(1) - __s)
00330 #else
00331 * std::exp(__log_gamma(_Tp(1) - __s))
00332 #endif
00333 * __riemann_zeta_sum(_Tp(1) - __s);
00334 return __zeta;
00335 }
00336 }
00337 }
00338 else
00339 return __riemann_zeta_product(__s);
00340 }
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364 template<typename _Tp>
00365 _Tp
00366 __hurwitz_zeta_glob(const _Tp __a, const _Tp __s)
00367 {
00368 _Tp __zeta = _Tp(0);
00369
00370 const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
00371
00372 const _Tp __max_bincoeff = std::numeric_limits<_Tp>::max_exponent10
00373 * std::log(_Tp(10)) - _Tp(1);
00374
00375 const unsigned int __maxit = 10000;
00376 for (unsigned int __i = 0; __i < __maxit; ++__i)
00377 {
00378 bool __punt = false;
00379 _Tp __sgn = _Tp(1);
00380 _Tp __term = _Tp(0);
00381 for (unsigned int __j = 0; __j <= __i; ++__j)
00382 {
00383 #if _GLIBCXX_USE_C99_MATH_TR1
00384 _Tp __bincoeff = std::tr1::lgamma(_Tp(1 + __i))
00385 - std::tr1::lgamma(_Tp(1 + __j))
00386 - std::tr1::lgamma(_Tp(1 + __i - __j));
00387 #else
00388 _Tp __bincoeff = __log_gamma(_Tp(1 + __i))
00389 - __log_gamma(_Tp(1 + __j))
00390 - __log_gamma(_Tp(1 + __i - __j));
00391 #endif
00392 if (__bincoeff > __max_bincoeff)
00393 {
00394
00395 __punt = true;
00396 break;
00397 }
00398 __bincoeff = std::exp(__bincoeff);
00399 __term += __sgn * __bincoeff * std::pow(_Tp(__a + __j), -__s);
00400 __sgn *= _Tp(-1);
00401 }
00402 if (__punt)
00403 break;
00404 __term /= _Tp(__i + 1);
00405 if (std::abs(__term / __zeta) < __eps)
00406 break;
00407 __zeta += __term;
00408 }
00409
00410 __zeta /= __s - _Tp(1);
00411
00412 return __zeta;
00413 }
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429 template<typename _Tp>
00430 inline _Tp
00431 __hurwitz_zeta(const _Tp __a, const _Tp __s)
00432 {
00433 return __hurwitz_zeta_glob(__a, __s);
00434 }
00435
00436 }
00437 }
00438 }
00439
00440 #endif // _GLIBCXX_TR1_RIEMANN_ZETA_TCC