8

Possible Duplicate:
How to evaluate Riemann Zeta function

I have an amateur interest in the Zeta Function. I have read Edward's book on the topic, which is perhaps a little dated. I would like to know any modern methods for computing estimates of

$\zeta(s)$, specifically for the critical strip.

My goal is to convert the method into C++ code, and actually for my purposes the accuracy and convergence are less important than the ease of writing the code.

Thanks in advance.

Community
  • 1
JJG
  • 727
  • 7
  • 11

2 Answers2

8

This article of Gourdon and Sebah 'Numerical evaluation of the Riemann Zeta-function' should be a fine reading about different evaluation methods (the authors may have code and other stuff here).

You should not hope too much algorithmic evolution since Edwards' book : multiple-evaluation became faster but AFAIK the time required to evaluate a single value remains proportional to $\sqrt{|\Im(z)|}$ with the best method (Riemann-Siegel). A method easier to implement and more accurate for small values is Euler Maclaurin but in this case the time will be proportional to $|\Im(z)|$. These methods are all detailed in the first link so good reading !

For an implementation of Riemann-Siegel you may see this link that used some inspiration from Ralph Pugh's thesis that used Edward's book...

Raymond Manzoni
  • 43,021
  • 5
  • 86
  • 140
  • 1
    With respect to computations in the critical strip (and the numerical use of the Riemann-Siegel formula), one would do well to look into work by Andrew Odlyzko. – J. M. ain't a mathematician Aug 17 '12 at 16:05
  • 2
    @J.M.: references are in the first paper of course and Odlyzko has excellent and well known tables and papers here, – Raymond Manzoni Aug 17 '12 at 16:08
  • (Part 1 of 2) I've read a great deal of the paper provided and have a few comments to add. Firstly, please note there was a typo right after (9). It is written "z=2(t-m)-1", but it should say "z=2(-m)-1". Next, while this formula does seem to work, no sources seem to give a clear definition on how to efficiently evaluate the nth derivative of the function defined in the paper. You could use a recursive formula and evaluate it in terms of trig functions, but it would be terribly inefficient and prone to roundoffs as z goes to ±0.5. – Math Machine Aug 15 '20 at 17:18
  • (Part 2 of 3) You could also try using a Taylor's series expansion about t=1/2, but it would diverge if z is imaginary. Which isn't a problem if you're only evaluating in the critical strip, but you lose the benefit of being able to compute ζ(s) for all s of large imaginary part. The last thing I'd like to add about this paper is that some of the information provided is incomplete. More precisely, the author gives you a link so you can read more about some of it, but it requires paid access. I'd instead recommend visiting the wolfram alpha page. – Math Machine Aug 15 '20 at 17:24
  • (Part 3 of 3) My apologies, I didn't foresee comment 2 being this long. https://mathworld.wolfram.com/Riemann-SiegelFormula.html . Please note the function listed in Wolfram Alpha is slightly different than in the paper, namely the input is shifted by 1/2 and scaled down by 2. – Math Machine Aug 15 '20 at 17:27
  • @MathMachine: Thanks for the typo $t \to \tau$. Note that the link provided at the end of my answer has this part right and further details concerning your other problems (evaluation of derivatives of $\psi$ are often done as Taylor series using a CAS, only a fixed number $M$ of derivatives are evaluated since Riemann-Siegel is used for evaluation of zeros with reduced precision but with $t$ the imaginary part of $s=\frac 12+it$ large). – Raymond Manzoni Aug 15 '20 at 23:13
  • The Riemann-Siegel method is fairly easy to implement and understand. However, the Odlyzko–Schönhage algorithm (and variants of it developed by Bober and Hiary) are much faster than the Riemann-Siegel formula. I will make a new answer in the duplicated question. – Axion004 Nov 07 '21 at 23:38
5

Euler-Maclaurin Summation

This was one of the first techniques used to approximate the zeta function ans was in fact used by Euler to approximate $\zeta(2)$. However, this method is only used on the remainder after a certain number terms of the zeta series has been computed.

$$\zeta(s)=\sum_{n=1}^N \frac{1}{n^s}+\frac{N^{1-s}}{s-1}+\frac{N^{-s}}{2}+\sum_{r=1}^{q-1}\frac{B_{2r}}{(2r)!}s(s+1) \cdots(s+2r-2)N^{-s-2r+1}+\epsilon_{2q}(s)$$

where

$$|\epsilon_{2q}(s)| < \left|\frac{s(s+1) \cdots(s+2r-2)N^{-\operatorname{Re}[s]-2r+1}}{(2q)!(s+2q-1)}\right|$$

Alternating Series

The alternating zeta series is given by

$$\zeta_a (s)=\sum_{n=1}^\infty \frac{(-1)^{n-1}}{n^s}=(1-2^{1-s})\zeta(s)$$

Then, by accelerating this sum, we have

$$e_k = \sum_{j=k}^n {n \choose j}$$

$$\zeta(s)=\frac{1}{(1-2^{1-s})}\left(\sum_{k=1}^n \frac{(-1)^{k=1}}{k^s}+\frac{1}{2^n }\sum_{k=n+1}^{2n} \frac{(-1)^{k=1}e_k}{k^s}\right)+\epsilon_n(s)$$

where

$$\epsilon_n(s) < \frac{(1+|t/\operatorname{Re}[s]|)\exp(|t|\pi/2)}{8^n|1-2^{1-s}|}$$

Another, different acceleration gives a slightly faster, but more complicated result

$$d_k = n \sum_{j=k}^n \frac{(n + j − 1)!4^j}{(n − j)!(2j)!}$$

$$\zeta(s)=\frac{1}{d_0(1-2^{1-s})}\sum_{k=1}^n \frac{(-1)^{k-1}d_k}{k^s}+\epsilon_n(s)$$

where

$$|\epsilon_n(s)| \le \frac{2}{(3+\sqrt{8})^n |\Gamma(n)(1-2^{1-s})|}$$

This paper gives pseudocode (pg. 41) for the Zeta function using the Euler Maclaurin summation method and a Mathematica implementation (Appendix D, pg. 57).

Argon
  • 25,303
  • 2
    Notes: The first accelerated method under "Alternating Series" is essentially the Euler transformation applied to Dirichlet $\eta$, while the second method, which is based on Chebyshev polynomials, is due to Cohen, Rodriguez Villegas, and Zagier. – J. M. ain't a mathematician Aug 17 '12 at 16:00
  • When you wrote $(-1)^{k=1}$, did you mean to write $k-1$? And where did the $t$ in $\epsilon _n(s)$ come from? (My apologies if you wrote this answer too long ago to remember.) – a52 Dec 06 '17 at 09:35