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