Modern formula for calculating Riemann Zeta Function [duplicate]

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.

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...

• 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. Aug 17, 2012 at 16:05
• @J.M.: references are in the first paper of course and Odlyzko has excellent and well known tables and papers here, Aug 17, 2012 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. Aug 15, 2020 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. Aug 15, 2020 at 17:24
• (Part 3 of 3) My apologies, I didn't foresee comment 2 being this long. 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. Aug 15, 2020 at 17:27

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).

• 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. Aug 17, 2012 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 6, 2017 at 9:35