libstdc++
riemann_zeta.tcc
Go to the documentation of this file.
1 // Special functions -*- C++ -*-
2 
3 // Copyright (C) 2006, 2007, 2008, 2009
4 // Free Software Foundation, Inc.
5 //
6 // This file is part of the GNU ISO C++ Library. This library is free
7 // software; you can redistribute it and/or modify it under the
8 // terms of the GNU General Public License as published by the
9 // Free Software Foundation; either version 3, or (at your option)
10 // any later version.
11 //
12 // This library is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
16 //
17 // Under Section 7 of GPL version 3, you are granted additional
18 // permissions described in the GCC Runtime Library Exception, version
19 // 3.1, as published by the Free Software Foundation.
20 
21 // You should have received a copy of the GNU General Public License and
22 // a copy of the GCC Runtime Library Exception along with this program;
23 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
24 // <http://www.gnu.org/licenses/>.
25 
26 /** @file tr1/riemann_zeta.tcc
27  * This is an internal header file, included by other library headers.
28  * You should not attempt to use it directly.
29  */
30 
31 //
32 // ISO C++ 14882 TR1: 5.2 Special functions
33 //
34 
35 // Written by Edward Smith-Rowland based on:
36 // (1) Handbook of Mathematical Functions,
37 // Ed. by Milton Abramowitz and Irene A. Stegun,
38 // Dover Publications, New-York, Section 5, pp. 807-808.
39 // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl
40 // (3) Gamma, Exploring Euler's Constant, Julian Havil,
41 // Princeton, 2003.
42 
43 #ifndef _GLIBCXX_TR1_RIEMANN_ZETA_TCC
44 #define _GLIBCXX_TR1_RIEMANN_ZETA_TCC 1
45 
46 #include "special_function_util.h"
47 
48 namespace std
49 {
50 namespace tr1
51 {
52 
53  // [5.2] Special functions
54 
55  // Implementation-space details.
56  namespace __detail
57  {
58 
59  /**
60  * @brief Compute the Riemann zeta function @f$ \zeta(s) @f$
61  * by summation for s > 1.
62  *
63  * The Riemann zeta function is defined by:
64  * \f[
65  * \zeta(s) = \sum_{k=1}^{\infty} \frac{1}{k^{s}} for s > 1
66  * \f]
67  * For s < 1 use the reflection formula:
68  * \f[
69  * \zeta(s) = 2^s \pi^{s-1} \Gamma(1-s) \zeta(1-s)
70  * \f]
71  */
72  template<typename _Tp>
73  _Tp
74  __riemann_zeta_sum(const _Tp __s)
75  {
76  // A user shouldn't get to this.
77  if (__s < _Tp(1))
78  std::__throw_domain_error(__N("Bad argument in zeta sum."));
79 
80  const unsigned int max_iter = 10000;
81  _Tp __zeta = _Tp(0);
82  for (unsigned int __k = 1; __k < max_iter; ++__k)
83  {
84  _Tp __term = std::pow(static_cast<_Tp>(__k), -__s);
85  if (__term < std::numeric_limits<_Tp>::epsilon())
86  {
87  break;
88  }
89  __zeta += __term;
90  }
91 
92  return __zeta;
93  }
94 
95 
96  /**
97  * @brief Evaluate the Riemann zeta function @f$ \zeta(s) @f$
98  * by an alternate series for s > 0.
99  *
100  * The Riemann zeta function is defined by:
101  * \f[
102  * \zeta(s) = \sum_{k=1}^{\infty} \frac{1}{k^{s}} for s > 1
103  * \f]
104  * For s < 1 use the reflection formula:
105  * \f[
106  * \zeta(s) = 2^s \pi^{s-1} \Gamma(1-s) \zeta(1-s)
107  * \f]
108  */
109  template<typename _Tp>
110  _Tp
111  __riemann_zeta_alt(const _Tp __s)
112  {
113  _Tp __sgn = _Tp(1);
114  _Tp __zeta = _Tp(0);
115  for (unsigned int __i = 1; __i < 10000000; ++__i)
116  {
117  _Tp __term = __sgn / std::pow(__i, __s);
119  break;
120  __zeta += __term;
121  __sgn *= _Tp(-1);
122  }
123  __zeta /= _Tp(1) - std::pow(_Tp(2), _Tp(1) - __s);
124 
125  return __zeta;
126  }
127 
128 
129  /**
130  * @brief Evaluate the Riemann zeta function by series for all s != 1.
131  * Convergence is great until largish negative numbers.
132  * Then the convergence of the > 0 sum gets better.
133  *
134  * The series is:
135  * \f[
136  * \zeta(s) = \frac{1}{1-2^{1-s}}
137  * \sum_{n=0}^{\infty} \frac{1}{2^{n+1}}
138  * \sum_{k=0}^{n} (-1)^k \frac{n!}{(n-k)!k!} (k+1)^{-s}
139  * \f]
140  * Havil 2003, p. 206.
141  *
142  * The Riemann zeta function is defined by:
143  * \f[
144  * \zeta(s) = \sum_{k=1}^{\infty} \frac{1}{k^{s}} for s > 1
145  * \f]
146  * For s < 1 use the reflection formula:
147  * \f[
148  * \zeta(s) = 2^s \pi^{s-1} \Gamma(1-s) \zeta(1-s)
149  * \f]
150  */
151  template<typename _Tp>
152  _Tp
153  __riemann_zeta_glob(const _Tp __s)
154  {
155  _Tp __zeta = _Tp(0);
156 
157  const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
158  // Max e exponent before overflow.
159  const _Tp __max_bincoeff = std::numeric_limits<_Tp>::max_exponent10
160  * std::log(_Tp(10)) - _Tp(1);
161 
162  // This series works until the binomial coefficient blows up
163  // so use reflection.
164  if (__s < _Tp(0))
165  {
166 #if _GLIBCXX_USE_C99_MATH_TR1
167  if (std::tr1::fmod(__s,_Tp(2)) == _Tp(0))
168  return _Tp(0);
169  else
170 #endif
171  {
172  _Tp __zeta = __riemann_zeta_glob(_Tp(1) - __s);
173  __zeta *= std::pow(_Tp(2)
176 #if _GLIBCXX_USE_C99_MATH_TR1
177  * std::exp(std::tr1::lgamma(_Tp(1) - __s))
178 #else
179  * std::exp(__log_gamma(_Tp(1) - __s))
180 #endif
182  return __zeta;
183  }
184  }
185 
186  _Tp __num = _Tp(0.5L);
187  const unsigned int __maxit = 10000;
188  for (unsigned int __i = 0; __i < __maxit; ++__i)
189  {
190  bool __punt = false;
191  _Tp __sgn = _Tp(1);
192  _Tp __term = _Tp(0);
193  for (unsigned int __j = 0; __j <= __i; ++__j)
194  {
195 #if _GLIBCXX_USE_C99_MATH_TR1
196  _Tp __bincoeff = std::tr1::lgamma(_Tp(1 + __i))
197  - std::tr1::lgamma(_Tp(1 + __j))
198  - std::tr1::lgamma(_Tp(1 + __i - __j));
199 #else
200  _Tp __bincoeff = __log_gamma(_Tp(1 + __i))
201  - __log_gamma(_Tp(1 + __j))
202  - __log_gamma(_Tp(1 + __i - __j));
203 #endif
204  if (__bincoeff > __max_bincoeff)
205  {
206  // This only gets hit for x << 0.
207  __punt = true;
208  break;
209  }
210  __bincoeff = std::exp(__bincoeff);
211  __term += __sgn * __bincoeff * std::pow(_Tp(1 + __j), -__s);
212  __sgn *= _Tp(-1);
213  }
214  if (__punt)
215  break;
216  __term *= __num;
217  __zeta += __term;
218  if (std::abs(__term/__zeta) < __eps)
219  break;
220  __num *= _Tp(0.5L);
221  }
222 
223  __zeta /= _Tp(1) - std::pow(_Tp(2), _Tp(1) - __s);
224 
225  return __zeta;
226  }
227 
228 
229  /**
230  * @brief Compute the Riemann zeta function @f$ \zeta(s) @f$
231  * using the product over prime factors.
232  * \f[
233  * \zeta(s) = \Pi_{i=1}^\infty \frac{1}{1 - p_i^{-s}}
234  * \f]
235  * where @f$ {p_i} @f$ are the prime numbers.
236  *
237  * The Riemann zeta function is defined by:
238  * \f[
239  * \zeta(s) = \sum_{k=1}^{\infty} \frac{1}{k^{s}} for s > 1
240  * \f]
241  * For s < 1 use the reflection formula:
242  * \f[
243  * \zeta(s) = 2^s \pi^{s-1} \Gamma(1-s) \zeta(1-s)
244  * \f]
245  */
246  template<typename _Tp>
247  _Tp
248  __riemann_zeta_product(const _Tp __s)
249  {
250  static const _Tp __prime[] = {
251  _Tp(2), _Tp(3), _Tp(5), _Tp(7), _Tp(11), _Tp(13), _Tp(17), _Tp(19),
252  _Tp(23), _Tp(29), _Tp(31), _Tp(37), _Tp(41), _Tp(43), _Tp(47),
253  _Tp(53), _Tp(59), _Tp(61), _Tp(67), _Tp(71), _Tp(73), _Tp(79),
254  _Tp(83), _Tp(89), _Tp(97), _Tp(101), _Tp(103), _Tp(107), _Tp(109)
255  };
256  static const unsigned int __num_primes = sizeof(__prime) / sizeof(_Tp);
257 
258  _Tp __zeta = _Tp(1);
259  for (unsigned int __i = 0; __i < __num_primes; ++__i)
260  {
261  const _Tp __fact = _Tp(1) - std::pow(__prime[__i], -__s);
262  __zeta *= __fact;
263  if (_Tp(1) - __fact < std::numeric_limits<_Tp>::epsilon())
264  break;
265  }
266 
267  __zeta = _Tp(1) / __zeta;
268 
269  return __zeta;
270  }
271 
272 
273  /**
274  * @brief Return the Riemann zeta function @f$ \zeta(s) @f$.
275  *
276  * The Riemann zeta function is defined by:
277  * \f[
278  * \zeta(s) = \sum_{k=1}^{\infty} k^{-s} for s > 1
279  * \frac{(2\pi)^s}{pi} sin(\frac{\pi s}{2})
280  * \Gamma (1 - s) \zeta (1 - s) for s < 1
281  * \f]
282  * For s < 1 use the reflection formula:
283  * \f[
284  * \zeta(s) = 2^s \pi^{s-1} \Gamma(1-s) \zeta(1-s)
285  * \f]
286  */
287  template<typename _Tp>
288  _Tp
289  __riemann_zeta(const _Tp __s)
290  {
291  if (__isnan(__s))
293  else if (__s == _Tp(1))
295  else if (__s < -_Tp(19))
296  {
297  _Tp __zeta = __riemann_zeta_product(_Tp(1) - __s);
298  __zeta *= std::pow(_Tp(2) * __numeric_constants<_Tp>::__pi(), __s)
300 #if _GLIBCXX_USE_C99_MATH_TR1
301  * std::exp(std::tr1::lgamma(_Tp(1) - __s))
302 #else
303  * std::exp(__log_gamma(_Tp(1) - __s))
304 #endif
306  return __zeta;
307  }
308  else if (__s < _Tp(20))
309  {
310  // Global double sum or McLaurin?
311  bool __glob = true;
312  if (__glob)
313  return __riemann_zeta_glob(__s);
314  else
315  {
316  if (__s > _Tp(1))
317  return __riemann_zeta_sum(__s);
318  else
319  {
320  _Tp __zeta = std::pow(_Tp(2)
323 #if _GLIBCXX_USE_C99_MATH_TR1
324  * std::tr1::tgamma(_Tp(1) - __s)
325 #else
326  * std::exp(__log_gamma(_Tp(1) - __s))
327 #endif
328  * __riemann_zeta_sum(_Tp(1) - __s);
329  return __zeta;
330  }
331  }
332  }
333  else
334  return __riemann_zeta_product(__s);
335  }
336 
337 
338  /**
339  * @brief Return the Hurwitz zeta function @f$ \zeta(x,s) @f$
340  * for all s != 1 and x > -1.
341  *
342  * The Hurwitz zeta function is defined by:
343  * @f[
344  * \zeta(x,s) = \sum_{n=0}^{\infty} \frac{1}{(n + x)^s}
345  * @f]
346  * The Riemann zeta function is a special case:
347  * @f[
348  * \zeta(s) = \zeta(1,s)
349  * @f]
350  *
351  * This functions uses the double sum that converges for s != 1
352  * and x > -1:
353  * @f[
354  * \zeta(x,s) = \frac{1}{s-1}
355  * \sum_{n=0}^{\infty} \frac{1}{n + 1}
356  * \sum_{k=0}^{n} (-1)^k \frac{n!}{(n-k)!k!} (x+k)^{-s}
357  * @f]
358  */
359  template<typename _Tp>
360  _Tp
361  __hurwitz_zeta_glob(const _Tp __a, const _Tp __s)
362  {
363  _Tp __zeta = _Tp(0);
364 
365  const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
366  // Max e exponent before overflow.
367  const _Tp __max_bincoeff = std::numeric_limits<_Tp>::max_exponent10
368  * std::log(_Tp(10)) - _Tp(1);
369 
370  const unsigned int __maxit = 10000;
371  for (unsigned int __i = 0; __i < __maxit; ++__i)
372  {
373  bool __punt = false;
374  _Tp __sgn = _Tp(1);
375  _Tp __term = _Tp(0);
376  for (unsigned int __j = 0; __j <= __i; ++__j)
377  {
378 #if _GLIBCXX_USE_C99_MATH_TR1
379  _Tp __bincoeff = std::tr1::lgamma(_Tp(1 + __i))
380  - std::tr1::lgamma(_Tp(1 + __j))
381  - std::tr1::lgamma(_Tp(1 + __i - __j));
382 #else
383  _Tp __bincoeff = __log_gamma(_Tp(1 + __i))
384  - __log_gamma(_Tp(1 + __j))
385  - __log_gamma(_Tp(1 + __i - __j));
386 #endif
387  if (__bincoeff > __max_bincoeff)
388  {
389  // This only gets hit for x << 0.
390  __punt = true;
391  break;
392  }
393  __bincoeff = std::exp(__bincoeff);
394  __term += __sgn * __bincoeff * std::pow(_Tp(__a + __j), -__s);
395  __sgn *= _Tp(-1);
396  }
397  if (__punt)
398  break;
399  __term /= _Tp(__i + 1);
400  if (std::abs(__term / __zeta) < __eps)
401  break;
402  __zeta += __term;
403  }
404 
405  __zeta /= __s - _Tp(1);
406 
407  return __zeta;
408  }
409 
410 
411  /**
412  * @brief Return the Hurwitz zeta function @f$ \zeta(x,s) @f$
413  * for all s != 1 and x > -1.
414  *
415  * The Hurwitz zeta function is defined by:
416  * @f[
417  * \zeta(x,s) = \sum_{n=0}^{\infty} \frac{1}{(n + x)^s}
418  * @f]
419  * The Riemann zeta function is a special case:
420  * @f[
421  * \zeta(s) = \zeta(1,s)
422  * @f]
423  */
424  template<typename _Tp>
425  inline _Tp
426  __hurwitz_zeta(const _Tp __a, const _Tp __s)
427  {
428  return __hurwitz_zeta_glob(__a, __s);
429  }
430 
431  } // namespace std::tr1::__detail
432 }
433 }
434 
435 #endif // _GLIBCXX_TR1_RIEMANN_ZETA_TCC
complex< _Tp > log(const complex< _Tp > &)
Return complex natural logarithm of z.
Definition: complex:782
A structure for numeric constants.
_Tp __riemann_zeta_glob(const _Tp __s)
Evaluate the Riemann zeta function by series for all s != 1. Convergence is great until largish negat...
_Tp __hurwitz_zeta(const _Tp __a, const _Tp __s)
Return the Hurwitz zeta function for all s != 1 and x > -1.
_Tp __riemann_zeta_alt(const _Tp __s)
Evaluate the Riemann zeta function by an alternate series for s > 0.
_Tp __hurwitz_zeta_glob(const _Tp __a, const _Tp __s)
Return the Hurwitz zeta function for all s != 1 and x > -1.
_Tp __riemann_zeta(const _Tp __s)
Return the Riemann zeta function .
Properties of fundamental types.
Definition: limits:278
static _Tp epsilon()
Definition: limits:287
_Tp __riemann_zeta_sum(const _Tp __s)
Compute the Riemann zeta function by summation for s > 1.
_Tp __log_gamma(const _Tp __x)
Return . This will return values even for . To recover the sign of for any argument use __log_gamma_...
Definition: gamma.tcc:221
complex< _Tp > sin(const complex< _Tp > &)
Return complex sine of z.
Definition: complex:817
complex< _Tp > pow(const complex< _Tp > &, const _Tp &)
Return x to the y'th power.
Definition: complex:964
_Tp __riemann_zeta_product(const _Tp __s)
Compute the Riemann zeta function using the product over prime factors. where are the prime number...
_Tp abs(const complex< _Tp > &)
Return magnitude of z.
Definition: complex:594
static _Tp quiet_NaN()
Definition: limits:293
complex< _Tp > exp(const complex< _Tp > &)
Return complex base e exponential of z.
Definition: complex:755
static _Tp infinity()
Definition: limits:291