Tuesday, January 19, 2010

Zeta evaluation with the Riemann-Siegel expansion

I'm very grateful to Juan Arias de Reyna who has contributed a module implementing zeta function evaluation using the Riemann-Siegel formula in mpmath. This finally permits rapid evaluation of the Riemann zeta function for arguments with large imaginary part.

To follow tradition on this blog, pictorial proof shall be given. Here is a plot of a segment of the critical line, ζ(1/2+ti) with t between 1010+50 and 1010+55:A complex plot of showing the critical strip, t ranging between 108+40 and 108+45 (note: the y scale is reversed):

Juan Arias de Reyna, who is a professor of mathematics at the University of Seville, has done a thorough job with this code. He has even proved rigorous error bounds for his algorithm (subject to assumptions that the underlying functions in mpmath being well-implemented). The news is that the bounds -- documented in an as-yet unpublished paper -- are valid off the critical line. The code also computes derivatives (up to 4th derivatives), although not as rigorously but still very robustly.

I integrated the code (and added a few optimizations) during the last few days. The zeta function in mpmath now automatically switches between Borwein's algorithm (close to the real line), Euler-Maclaurin summation, and the Riemann-Siegel expansion.

Some example values and timings on my laptop, computing ζ(1/2+10ni) for n up to 12:

The implementation also supports use of the fp context. Double precision unavoidably becomes insufficient as the imaginary part approaches 1015, but it has the advantage of speed in the range where it works:

We have done some comparison with Mathematica, and the mpmath version appears to be about as fast (a bit faster or a bit slower, sometimes substantially faster, depending on circumstance). The most expensive part of the computation occurs in a simple internal function that adds lots of n−s terms. I think for Sage, it will be very easy to switch to a Cython version of this function which should improve speed by a large factor.

But most importantly, Mathematica's Zeta[] is notoriously buggy for large imaginary arguments. As a first example, here Mathematica 7.0 computes two entirely different values at different precisions:

That's all for now. I'm back in school again, so maybe I won't have as much time for programming in the near future. But my schedule is flexible, so we'll see. A new release of mpmath shouldn't be far away.

Update: I should point out that the bugs in Mathematica's Zeta[] only seem to occur in recent versions. Mathematica 5 and older versions do not seem to have these problems.