How to evaluate Riemann Zeta function

How do I evaluate this function for given $s$?

$$\zeta(s) = \sum_{n=1}^\infty \frac1{n^s} = \frac{1}{1^s} + \frac{1}{2^s} + \frac{1}{3^s} + \cdots$$

Solution 1:

Assuming you're talking about numerical evaluation, along the critical strip $0<\Re s<1$ and "large" $\Im s$, (which is the region of interest for many) the Riemann-Siegel formula is standard; off the strip, what you can manage is a polyalgorithm.

For $\Re s\leq0$, one can use the reflection formula for $\zeta$,

$$\zeta(1-s)=\frac{2}{(2\pi)^s}\cos\left(\frac{s\pi}{2}\right)\Gamma(s)\zeta(s)$$

so we can consider evaluation for $\Re s>0$ in what follows.

Note that if $|s|$ is "large enough" (how large is "large" depends on the computing environment you're in), one can simply use the defining series $\zeta(s)=\sum\limits_{j=1}^\infty j^{-s}$, since its terms quickly diminish in magnitude. That leaves the problem of what to do for small to medium-sized $s$.

The key is to use the Dirichlet $\eta$ function:

$$\eta(s):=-\sum_{j=1}^\infty \frac{(-1)^j}{j^s}$$

which is related to $\zeta$ by the following identity:

$$\eta(s)=\left(1-2^{1-s}\right)\zeta(s)$$

The reason for our interest in $\eta$ is that although this sum is slowly convergent, it is an alternating series, and a number of algorithms exist for quickly finding the sum of an alternating series numerically.

One method is the Levin transformation; another one, which is one of the simplest methods for numerically summing alternating series (and my personal favorite) is the Cohen-Rodriguez Villegas-Zagier algorithm. The algorithm is a bit too long to describe here, so I will just have to point you to the original paper.

This is in fact identical to the approach taken by Borwein in this paper.

Solution 2:

You should look at the work of Wadim Zudilin. In particular you should look at "One of the numbers ζ(5), ζ(7), ζ(9), ζ(11) is irrational" (Turpion link, pdf 91k, gzip ps 80k) in Russian Math. Surveys 56:4 (2001), 774--776;

Solution 3:

Early implementations are based on variations of the Riemann-Siegel formula. A key figure in the development of numerical algorithms was Andrew Odlyzko. He wrote two fundamental papers on calculating zeros on the critical line - Fast algorithms for multiple evaluations of the Riemann zeta function and Supercomputers and the Riemann zeta function. A bunch of other people (including Alan Turing) have contributed to the verification of the Riemann hypothesis. A good summary is included in the report by Xavier Gourdon.

More recently, practitioners have been working on improvements to the Odlyzko–Schönhage algorithm (also see here). Two key figures are Jonathan Bober and Ghaith Hiary. Hiary was a PhD student of Odlyzko who works on fast computations of Dirichlet functions (his thesis was on a variation of the Odlyzko–Schönhage algorithm). There is an implementation of the Odlyzko–Schönhage algorithm on Bober's webpage. However, this is fairly difficult to read.

A bunch of students have written there own implementations of the Riemann-Siegel formula. Perhaps the most well known is the C implementation written by Glen Pugh (this is sometimes incorrectly cited as his PhD thesis. Glen wrote a wonderful PhD thesis about Lanzcos Gamma Approximations). Another computer implementation in C++ was written by Ken Takusagawa.

I wrote my own Riemann-Siegel code in Java (which wasn't a good idea. However, I worked as a Java programmer at the time and wasn't familiar with C/C++). I summarized my findings in a report which was later expanded in a Beamer presentation.

If you are interested in writing your own implementation then I would suggest reading Riemann's Zeta Function by Harold M. Edwards. An English translation of Bernhard Riemann's eight-page paper entitled "On the Number of Primes Less Than a Given Magnitude" is included at the end of the book.