]>
Commit | Line | Data |
---|---|---|
b34f60ac | 1 | // Special functions -*- C++ -*- |
2 | ||
fbd26352 | 3 | // Copyright (C) 2006-2019 Free Software Foundation, Inc. |
b34f60ac | 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 | |
6bc9506f | 8 | // Free Software Foundation; either version 3, or (at your option) |
b34f60ac | 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 | // | |
6bc9506f | 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/>. | |
b34f60ac | 24 | |
25 | /** @file tr1/riemann_zeta.tcc | |
26 | * This is an internal header file, included by other library headers. | |
5846aeac | 27 | * Do not attempt to use it directly. @headername{tr1/cmath} |
b34f60ac | 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. by Milton Abramowitz and Irene A. Stegun, | |
37 | // Dover Publications, New-York, Section 5, pp. 807-808. | |
38 | // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl | |
39 | // (3) Gamma, Exploring Euler's Constant, Julian Havil, | |
40 | // Princeton, 2003. | |
41 | ||
c17b0a1c | 42 | #ifndef _GLIBCXX_TR1_RIEMANN_ZETA_TCC |
43 | #define _GLIBCXX_TR1_RIEMANN_ZETA_TCC 1 | |
b34f60ac | 44 | |
45 | #include "special_function_util.h" | |
46 | ||
2948dd21 | 47 | namespace std _GLIBCXX_VISIBILITY(default) |
b34f60ac | 48 | { |
ae6a4ce9 | 49 | _GLIBCXX_BEGIN_NAMESPACE_VERSION |
50 | ||
dc043b87 | 51 | #if _GLIBCXX_USE_STD_SPEC_FUNCS |
ca5b12a4 | 52 | # define _GLIBCXX_MATH_NS ::std |
53 | #elif defined(_GLIBCXX_TR1_CMATH) | |
c17b0a1c | 54 | namespace tr1 |
55 | { | |
ca5b12a4 | 56 | # define _GLIBCXX_MATH_NS ::std::tr1 |
57 | #else | |
58 | # error do not include this header directly, use <cmath> or <tr1/cmath> | |
59 | #endif | |
b34f60ac | 60 | // [5.2] Special functions |
61 | ||
b34f60ac | 62 | // Implementation-space details. |
b34f60ac | 63 | namespace __detail |
64 | { | |
b34f60ac | 65 | /** |
66 | * @brief Compute the Riemann zeta function @f$ \zeta(s) @f$ | |
67 | * by summation for s > 1. | |
68 | * | |
69 | * The Riemann zeta function is defined by: | |
70 | * \f[ | |
71 | * \zeta(s) = \sum_{k=1}^{\infty} \frac{1}{k^{s}} for s > 1 | |
72 | * \f] | |
73 | * For s < 1 use the reflection formula: | |
74 | * \f[ | |
75 | * \zeta(s) = 2^s \pi^{s-1} \Gamma(1-s) \zeta(1-s) | |
76 | * \f] | |
77 | */ | |
78 | template<typename _Tp> | |
79 | _Tp | |
cd7f5f45 | 80 | __riemann_zeta_sum(_Tp __s) |
b34f60ac | 81 | { |
82 | // A user shouldn't get to this. | |
83 | if (__s < _Tp(1)) | |
84 | std::__throw_domain_error(__N("Bad argument in zeta sum.")); | |
85 | ||
86 | const unsigned int max_iter = 10000; | |
87 | _Tp __zeta = _Tp(0); | |
88 | for (unsigned int __k = 1; __k < max_iter; ++__k) | |
89 | { | |
90 | _Tp __term = std::pow(static_cast<_Tp>(__k), -__s); | |
91 | if (__term < std::numeric_limits<_Tp>::epsilon()) | |
92 | { | |
93 | break; | |
94 | } | |
95 | __zeta += __term; | |
96 | } | |
97 | ||
98 | return __zeta; | |
99 | } | |
100 | ||
101 | ||
102 | /** | |
103 | * @brief Evaluate the Riemann zeta function @f$ \zeta(s) @f$ | |
104 | * by an alternate series for s > 0. | |
105 | * | |
106 | * The Riemann zeta function is defined by: | |
107 | * \f[ | |
108 | * \zeta(s) = \sum_{k=1}^{\infty} \frac{1}{k^{s}} for s > 1 | |
109 | * \f] | |
110 | * For s < 1 use the reflection formula: | |
111 | * \f[ | |
112 | * \zeta(s) = 2^s \pi^{s-1} \Gamma(1-s) \zeta(1-s) | |
113 | * \f] | |
114 | */ | |
115 | template<typename _Tp> | |
116 | _Tp | |
cd7f5f45 | 117 | __riemann_zeta_alt(_Tp __s) |
b34f60ac | 118 | { |
119 | _Tp __sgn = _Tp(1); | |
120 | _Tp __zeta = _Tp(0); | |
121 | for (unsigned int __i = 1; __i < 10000000; ++__i) | |
122 | { | |
123 | _Tp __term = __sgn / std::pow(__i, __s); | |
124 | if (std::abs(__term) < std::numeric_limits<_Tp>::epsilon()) | |
125 | break; | |
126 | __zeta += __term; | |
127 | __sgn *= _Tp(-1); | |
128 | } | |
129 | __zeta /= _Tp(1) - std::pow(_Tp(2), _Tp(1) - __s); | |
130 | ||
131 | return __zeta; | |
132 | } | |
133 | ||
134 | ||
135 | /** | |
136 | * @brief Evaluate the Riemann zeta function by series for all s != 1. | |
137 | * Convergence is great until largish negative numbers. | |
138 | * Then the convergence of the > 0 sum gets better. | |
139 | * | |
140 | * The series is: | |
141 | * \f[ | |
142 | * \zeta(s) = \frac{1}{1-2^{1-s}} | |
143 | * \sum_{n=0}^{\infty} \frac{1}{2^{n+1}} | |
144 | * \sum_{k=0}^{n} (-1)^k \frac{n!}{(n-k)!k!} (k+1)^{-s} | |
145 | * \f] | |
146 | * Havil 2003, p. 206. | |
147 | * | |
148 | * The Riemann zeta function is defined by: | |
149 | * \f[ | |
150 | * \zeta(s) = \sum_{k=1}^{\infty} \frac{1}{k^{s}} for s > 1 | |
151 | * \f] | |
152 | * For s < 1 use the reflection formula: | |
153 | * \f[ | |
154 | * \zeta(s) = 2^s \pi^{s-1} \Gamma(1-s) \zeta(1-s) | |
155 | * \f] | |
156 | */ | |
157 | template<typename _Tp> | |
158 | _Tp | |
cd7f5f45 | 159 | __riemann_zeta_glob(_Tp __s) |
b34f60ac | 160 | { |
161 | _Tp __zeta = _Tp(0); | |
162 | ||
163 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); | |
164 | // Max e exponent before overflow. | |
165 | const _Tp __max_bincoeff = std::numeric_limits<_Tp>::max_exponent10 | |
166 | * std::log(_Tp(10)) - _Tp(1); | |
167 | ||
168 | // This series works until the binomial coefficient blows up | |
169 | // so use reflection. | |
170 | if (__s < _Tp(0)) | |
171 | { | |
172 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
ca5b12a4 | 173 | if (_GLIBCXX_MATH_NS::fmod(__s,_Tp(2)) == _Tp(0)) |
b34f60ac | 174 | return _Tp(0); |
175 | else | |
176 | #endif | |
177 | { | |
178 | _Tp __zeta = __riemann_zeta_glob(_Tp(1) - __s); | |
179 | __zeta *= std::pow(_Tp(2) | |
180 | * __numeric_constants<_Tp>::__pi(), __s) | |
181 | * std::sin(__numeric_constants<_Tp>::__pi_2() * __s) | |
182 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
ca5b12a4 | 183 | * std::exp(_GLIBCXX_MATH_NS::lgamma(_Tp(1) - __s)) |
b34f60ac | 184 | #else |
185 | * std::exp(__log_gamma(_Tp(1) - __s)) | |
186 | #endif | |
187 | / __numeric_constants<_Tp>::__pi(); | |
188 | return __zeta; | |
189 | } | |
190 | } | |
191 | ||
192 | _Tp __num = _Tp(0.5L); | |
193 | const unsigned int __maxit = 10000; | |
194 | for (unsigned int __i = 0; __i < __maxit; ++__i) | |
195 | { | |
196 | bool __punt = false; | |
197 | _Tp __sgn = _Tp(1); | |
198 | _Tp __term = _Tp(0); | |
199 | for (unsigned int __j = 0; __j <= __i; ++__j) | |
200 | { | |
201 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
ca5b12a4 | 202 | _Tp __bincoeff = _GLIBCXX_MATH_NS::lgamma(_Tp(1 + __i)) |
203 | - _GLIBCXX_MATH_NS::lgamma(_Tp(1 + __j)) | |
204 | - _GLIBCXX_MATH_NS::lgamma(_Tp(1 + __i - __j)); | |
b34f60ac | 205 | #else |
206 | _Tp __bincoeff = __log_gamma(_Tp(1 + __i)) | |
207 | - __log_gamma(_Tp(1 + __j)) | |
208 | - __log_gamma(_Tp(1 + __i - __j)); | |
209 | #endif | |
210 | if (__bincoeff > __max_bincoeff) | |
211 | { | |
212 | // This only gets hit for x << 0. | |
213 | __punt = true; | |
214 | break; | |
215 | } | |
216 | __bincoeff = std::exp(__bincoeff); | |
217 | __term += __sgn * __bincoeff * std::pow(_Tp(1 + __j), -__s); | |
218 | __sgn *= _Tp(-1); | |
219 | } | |
220 | if (__punt) | |
221 | break; | |
222 | __term *= __num; | |
223 | __zeta += __term; | |
224 | if (std::abs(__term/__zeta) < __eps) | |
225 | break; | |
226 | __num *= _Tp(0.5L); | |
227 | } | |
228 | ||
229 | __zeta /= _Tp(1) - std::pow(_Tp(2), _Tp(1) - __s); | |
230 | ||
231 | return __zeta; | |
232 | } | |
233 | ||
234 | ||
235 | /** | |
236 | * @brief Compute the Riemann zeta function @f$ \zeta(s) @f$ | |
237 | * using the product over prime factors. | |
238 | * \f[ | |
239 | * \zeta(s) = \Pi_{i=1}^\infty \frac{1}{1 - p_i^{-s}} | |
240 | * \f] | |
241 | * where @f$ {p_i} @f$ are the prime numbers. | |
242 | * | |
243 | * The Riemann zeta function is defined by: | |
244 | * \f[ | |
245 | * \zeta(s) = \sum_{k=1}^{\infty} \frac{1}{k^{s}} for s > 1 | |
246 | * \f] | |
247 | * For s < 1 use the reflection formula: | |
248 | * \f[ | |
249 | * \zeta(s) = 2^s \pi^{s-1} \Gamma(1-s) \zeta(1-s) | |
250 | * \f] | |
251 | */ | |
252 | template<typename _Tp> | |
253 | _Tp | |
cd7f5f45 | 254 | __riemann_zeta_product(_Tp __s) |
b34f60ac | 255 | { |
256 | static const _Tp __prime[] = { | |
257 | _Tp(2), _Tp(3), _Tp(5), _Tp(7), _Tp(11), _Tp(13), _Tp(17), _Tp(19), | |
258 | _Tp(23), _Tp(29), _Tp(31), _Tp(37), _Tp(41), _Tp(43), _Tp(47), | |
259 | _Tp(53), _Tp(59), _Tp(61), _Tp(67), _Tp(71), _Tp(73), _Tp(79), | |
260 | _Tp(83), _Tp(89), _Tp(97), _Tp(101), _Tp(103), _Tp(107), _Tp(109) | |
261 | }; | |
262 | static const unsigned int __num_primes = sizeof(__prime) / sizeof(_Tp); | |
263 | ||
264 | _Tp __zeta = _Tp(1); | |
265 | for (unsigned int __i = 0; __i < __num_primes; ++__i) | |
266 | { | |
267 | const _Tp __fact = _Tp(1) - std::pow(__prime[__i], -__s); | |
268 | __zeta *= __fact; | |
269 | if (_Tp(1) - __fact < std::numeric_limits<_Tp>::epsilon()) | |
270 | break; | |
271 | } | |
272 | ||
273 | __zeta = _Tp(1) / __zeta; | |
274 | ||
275 | return __zeta; | |
276 | } | |
277 | ||
278 | ||
279 | /** | |
280 | * @brief Return the Riemann zeta function @f$ \zeta(s) @f$. | |
281 | * | |
282 | * The Riemann zeta function is defined by: | |
283 | * \f[ | |
284 | * \zeta(s) = \sum_{k=1}^{\infty} k^{-s} for s > 1 | |
285 | * \frac{(2\pi)^s}{pi} sin(\frac{\pi s}{2}) | |
286 | * \Gamma (1 - s) \zeta (1 - s) for s < 1 | |
287 | * \f] | |
288 | * For s < 1 use the reflection formula: | |
289 | * \f[ | |
290 | * \zeta(s) = 2^s \pi^{s-1} \Gamma(1-s) \zeta(1-s) | |
291 | * \f] | |
292 | */ | |
293 | template<typename _Tp> | |
294 | _Tp | |
cd7f5f45 | 295 | __riemann_zeta(_Tp __s) |
b34f60ac | 296 | { |
297 | if (__isnan(__s)) | |
298 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
299 | else if (__s == _Tp(1)) | |
300 | return std::numeric_limits<_Tp>::infinity(); | |
301 | else if (__s < -_Tp(19)) | |
302 | { | |
303 | _Tp __zeta = __riemann_zeta_product(_Tp(1) - __s); | |
304 | __zeta *= std::pow(_Tp(2) * __numeric_constants<_Tp>::__pi(), __s) | |
305 | * std::sin(__numeric_constants<_Tp>::__pi_2() * __s) | |
306 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
ca5b12a4 | 307 | * std::exp(_GLIBCXX_MATH_NS::lgamma(_Tp(1) - __s)) |
b34f60ac | 308 | #else |
309 | * std::exp(__log_gamma(_Tp(1) - __s)) | |
310 | #endif | |
311 | / __numeric_constants<_Tp>::__pi(); | |
312 | return __zeta; | |
313 | } | |
314 | else if (__s < _Tp(20)) | |
315 | { | |
316 | // Global double sum or McLaurin? | |
317 | bool __glob = true; | |
318 | if (__glob) | |
319 | return __riemann_zeta_glob(__s); | |
320 | else | |
321 | { | |
322 | if (__s > _Tp(1)) | |
323 | return __riemann_zeta_sum(__s); | |
324 | else | |
325 | { | |
326 | _Tp __zeta = std::pow(_Tp(2) | |
327 | * __numeric_constants<_Tp>::__pi(), __s) | |
328 | * std::sin(__numeric_constants<_Tp>::__pi_2() * __s) | |
329 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
ca5b12a4 | 330 | * _GLIBCXX_MATH_NS::tgamma(_Tp(1) - __s) |
b34f60ac | 331 | #else |
332 | * std::exp(__log_gamma(_Tp(1) - __s)) | |
333 | #endif | |
334 | * __riemann_zeta_sum(_Tp(1) - __s); | |
335 | return __zeta; | |
336 | } | |
337 | } | |
338 | } | |
339 | else | |
340 | return __riemann_zeta_product(__s); | |
341 | } | |
342 | ||
343 | ||
344 | /** | |
345 | * @brief Return the Hurwitz zeta function @f$ \zeta(x,s) @f$ | |
346 | * for all s != 1 and x > -1. | |
347 | * | |
348 | * The Hurwitz zeta function is defined by: | |
349 | * @f[ | |
350 | * \zeta(x,s) = \sum_{n=0}^{\infty} \frac{1}{(n + x)^s} | |
351 | * @f] | |
352 | * The Riemann zeta function is a special case: | |
353 | * @f[ | |
354 | * \zeta(s) = \zeta(1,s) | |
355 | * @f] | |
356 | * | |
357 | * This functions uses the double sum that converges for s != 1 | |
358 | * and x > -1: | |
359 | * @f[ | |
360 | * \zeta(x,s) = \frac{1}{s-1} | |
361 | * \sum_{n=0}^{\infty} \frac{1}{n + 1} | |
362 | * \sum_{k=0}^{n} (-1)^k \frac{n!}{(n-k)!k!} (x+k)^{-s} | |
363 | * @f] | |
364 | */ | |
365 | template<typename _Tp> | |
366 | _Tp | |
cd7f5f45 | 367 | __hurwitz_zeta_glob(_Tp __a, _Tp __s) |
b34f60ac | 368 | { |
369 | _Tp __zeta = _Tp(0); | |
370 | ||
371 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); | |
372 | // Max e exponent before overflow. | |
373 | const _Tp __max_bincoeff = std::numeric_limits<_Tp>::max_exponent10 | |
374 | * std::log(_Tp(10)) - _Tp(1); | |
375 | ||
376 | const unsigned int __maxit = 10000; | |
377 | for (unsigned int __i = 0; __i < __maxit; ++__i) | |
378 | { | |
379 | bool __punt = false; | |
380 | _Tp __sgn = _Tp(1); | |
381 | _Tp __term = _Tp(0); | |
382 | for (unsigned int __j = 0; __j <= __i; ++__j) | |
383 | { | |
384 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
ca5b12a4 | 385 | _Tp __bincoeff = _GLIBCXX_MATH_NS::lgamma(_Tp(1 + __i)) |
386 | - _GLIBCXX_MATH_NS::lgamma(_Tp(1 + __j)) | |
387 | - _GLIBCXX_MATH_NS::lgamma(_Tp(1 + __i - __j)); | |
b34f60ac | 388 | #else |
389 | _Tp __bincoeff = __log_gamma(_Tp(1 + __i)) | |
390 | - __log_gamma(_Tp(1 + __j)) | |
391 | - __log_gamma(_Tp(1 + __i - __j)); | |
392 | #endif | |
393 | if (__bincoeff > __max_bincoeff) | |
394 | { | |
395 | // This only gets hit for x << 0. | |
396 | __punt = true; | |
397 | break; | |
398 | } | |
399 | __bincoeff = std::exp(__bincoeff); | |
400 | __term += __sgn * __bincoeff * std::pow(_Tp(__a + __j), -__s); | |
401 | __sgn *= _Tp(-1); | |
402 | } | |
403 | if (__punt) | |
404 | break; | |
405 | __term /= _Tp(__i + 1); | |
406 | if (std::abs(__term / __zeta) < __eps) | |
407 | break; | |
408 | __zeta += __term; | |
409 | } | |
410 | ||
411 | __zeta /= __s - _Tp(1); | |
412 | ||
413 | return __zeta; | |
414 | } | |
415 | ||
416 | ||
417 | /** | |
418 | * @brief Return the Hurwitz zeta function @f$ \zeta(x,s) @f$ | |
419 | * for all s != 1 and x > -1. | |
420 | * | |
421 | * The Hurwitz zeta function is defined by: | |
422 | * @f[ | |
423 | * \zeta(x,s) = \sum_{n=0}^{\infty} \frac{1}{(n + x)^s} | |
424 | * @f] | |
425 | * The Riemann zeta function is a special case: | |
426 | * @f[ | |
427 | * \zeta(s) = \zeta(1,s) | |
428 | * @f] | |
429 | */ | |
430 | template<typename _Tp> | |
431 | inline _Tp | |
cd7f5f45 | 432 | __hurwitz_zeta(_Tp __a, _Tp __s) |
433 | { return __hurwitz_zeta_glob(__a, __s); } | |
ca5b12a4 | 434 | } // namespace __detail |
435 | #undef _GLIBCXX_MATH_NS | |
dc043b87 | 436 | #if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH) |
ca5b12a4 | 437 | } // namespace tr1 |
438 | #endif | |
ae6a4ce9 | 439 | |
440 | _GLIBCXX_END_NAMESPACE_VERSION | |
b34f60ac | 441 | } |
442 | ||
c17b0a1c | 443 | #endif // _GLIBCXX_TR1_RIEMANN_ZETA_TCC |