(pdf version : pdf;
(postscript version : ps)
We expose techniques that permit to approximate the Riemann Zeta
function in different context. The motivation does not restrict only
to the numerical verification of the Riemann hypothesis
(see Numerical computations about the zeros of
the zeta function) :
in [5] for example, Lagarias and Odlyzko obtained the best
asymptotic known method to compute p(x), the number of primes less
than x, thanks to massive evaluation of the zetafunction; Odlyzko
and Te Riele in [7] used approximations of the first 2000
critical zeros of z(s) to disprove to Mertens conjecture;
finally, computations of z(n) at integer values of n also
consistute a motivation, together with series expressed in terms of
those values. (More details can be found on the motivation of
approximating z(s) in [2].)
The first section is dedicated to
numerical evaluations at a given point, the next section deals with
multievaluation.
1 Numerical evaluation of the Zeta function at a given point
In this section, we present different methods that permit to evaluate
numerically the Riemann Zeta function z(s) for complex numbers
s=s+it. The section is dedicated to general evaluation, and
does not specialize to evaluation of Zeta at integer argument for
example (see [2] where a more deep study is done).
As we shall see, the best method to approximate Zeta depends on the
purpose of evaluation : for evaluation at middle or high precision,
method of subsection 1.2 is suited, as for low
precision evaluation in the critical strip (for computations
concerning zeros) RiemannSiegel formula is more competitive.
This method is the historical first technique proposed to
numerically approximate z(s), and was developped by Euler
himself in his initial efforts to evaluate z(2).
The approach consists in applying
EulerMaclaurin summation not to the global z(s) series itself,
(see the section on analytic
continuation of the z(s) function where this technique was used), but
only to remainders å_{n > N} n^{s} of the series. The result is
z(s) =
N å
n=1
1
n^{s}
+
N^{1s}
s1
+
N^{s}
2
+
q1 å
r=1
B_{2r}
(2r)!
s(s+1)¼(s+2r2) N^{s2r+1} +E_{2q}(s)
(see [4] for details) for any complex number s=s+it
(with s = Â(s) > 2q+1). The terms B_{2r} are the Bernoulli numbers.
Backlung gave an estimate of the error term E_{2q}(s) in the form
E_{2q}(s) <
ê ê
s(s+1) ¼(s+2q1) B_{2q} N^{s2q+1}
(2q)!(s+2q1)
ê ê
.
However, this method is not very easy to use in the practice
since it requires evaluations of numerous
Bernoulli numbers, and the error term is difficult to control in this
nonconvergent expansion (see [2] where this difficulty is
presented).
Moreover, even with choosing carefully the parameters N and q, the
rate of convergence is not better than the one obtained with other
simpler techniques based on acceleration of alternating series
presented below.
A better method, more efficient and easier to implement, is based on the
alternating series convergence (see [1]).
It applies to complex numbers s=s+it with
s = Â(s) > 0. The starting point is the Zeta alternating series
z_{a}(s) =
¥ å
n=1
(1)^{n1}
n^{s}
= (12^{1s})z(s).
(1)
We apply the general technique of convergence of alternating series
of [3] described in Acceleration of the convergence of series. The technique applies
because the general term of the series 1/(n+1)^{s} is of the form
1
(n+1)^{s}
=
ó õ
1
0
f(x) x^{n} dx, with f(x) =
log(1/x)^{s1}
G(s)
Let P_{n}(x)=å_{k=0}^{n} p_{k} (x)^{k} be a polynomial such that
P_{n}(1) ¹ 0, and let define S_{n} as
S_{n} =
1
P_{n}(1)
n1 å
k=0
c_{k}
(1)^{k}
(k+1)^{s}
, c_{k}=
n å
j=k+1
p_{j}.
Then we have the identity
z_{a}(s)  S_{n} = x_{n}, x_{n} =
1
P_{n}(1)G(s)
ó õ
1
0
P_{n}(x)(log(1/x))^{s1}
1+x
dx.
(2)
As always, polynomials (P_{n}) are choosen such that
sup_{x Î [0,1]} P_{n}(x)/P_{n}(1) is small, making the value of
x_{n} small.
In [1], this identity is used with two families of
polynomials (P_{n}) : the Chebyshev polynomials shifted to [0,1],
and the polynomials x^{n}(1x)^{n}, giving the two following results :
Proposition 1
Let
d_{k} = n
n å
j=k
(n+j1)! 4^{j}
(nj)! (2j)!
then
z(s) =
1
d_{0} (12^{1s})
n å
k=1
(1)^{k1} d_{k}
k^{s}
+ g_{n}(s)
where for s=s+it with s ³ [ 1/2]
g_{n}(s) £
2
(3+Ö8)^{n}
1
G(s)
1
12^{1s}
£
3
(3+Ö8)^{n}
(1+2t) e^{tp/2}
12^{1s}
This approach, to compute z(s+it) with d decimal digits
of acurracy, requires a number of term n roughly equal to
1.3d+0.9t.
The next approach is based on the choice of the
polynomial x^{n}(1x)^{n}. Even if a little less efficient, the
technique is simpler. This time, the integral
in (2) is convergent whenever s = Â(s) > (n1), thus the technique can be also applied in this range.
Proposition 2
Let
e_{k} =
n å
j=k
æ è
n
j
ö ø
then
z(s) =
1
(12^{1s})
æ è
n å
k=1
(1)^{k1}
k^{s}
+
1
2^{n}
2n å
k=n+1
(1)^{k1} e_{kn}
k^{s}
ö ø
+ g_{n}(s)
where for s=s+it with s > 0
g_{n}(s) £
1
8^{n}
(1+t/s) e^{tp/2}
12^{1s}
.
If (n1) < s < 0 then
g_{n}(s) £
1
8^{n}
4^{s}
12^{1s} G(s)
.
Computational aspects
This technique is easy and efficient to compute z(s) with middle or
high accuracy, for values of s with small imaginary part
(in [1], computational experiments on the computation of
z(5) have been made proving
the efficiency of this approach compared to Zeta evaluation in classic
computer algebra systems like Maple, Mathematica or Pari).
In the practice, one should compute the terms sequentially, doing a
fixed number of additions and multiplications per term. One should
also notice that for integer values of s, the binarysplitting
technique can be used making these formula suitable for huge
computations in quasilinear time.
To compute values of z(1/2+it) for very large values of t,
(even with a low precision of 10 digits for example)
a number of needed terms proportional to t is required. As we
will see below, RiemannSiegel formula is more suited in this context.
The RiemannSiegel formula permits to approximate
z(s+it) in the critical strip 0 £ s £ 1 for
large values of t with a number of terms proportional to
Ö{t} (which is much better than previous methods that requires
a number of terms proportional to t).
It states in the general form as follows.
Theorem 1[RiemannSiegel formula]
Let x and y be positive real numbers such that 2pxy=t.
Then for s=s+it with 0 £ s £ 1, we have
z(s) =
å
n £ x
1
n^{s}
+ c(s)
å
n £ y
1
n^{1s}
+ O(x^{s}) +O(t^{1/2s} y^{s1}),
(3)
where the O holds uniformly for x > h and y > h for any h > 0.
We remind that c(s), often used to express functional equation of
Zeta, is defined by
c(s) = 2^{s} p^{s1} sin
æ è
ps
2
ö ø
G(1s).
Formula (3) looks
like the functional equation, for that reason it is also called the
approximate functional equation. A proof can be found
in [10] for example.
It is natural to use this formula when x=y=Ö{t/(2p)}, which
writes in the form
z(s) =
m å
n=1
1
n^{s}
+ c(s)
m å
n=1
1
n^{1s}
+ E_{m}(s), m =
ê ë
æ è
t
2p
ö ø
1/2
ú û
(4)
where the error term E_{m}(s) satisfies
E_{m}(s) = O(t^{s/2}).
An asymptotic expansion of the error term E_{m}(s) can be given
explicitly (see [10] or [8]), but
it is complicated in its general form and we just give its first
terms in the particular (but important) case s = 1/2 below.
RiemannSiegel formula on the critical line
It is interesting to specialize the RiemannSiegel formula on the
critical line s = 1/2.
We first define the RiemannSiegel Zfunction which permits a
convenient study of z(s) in this zone.
The RiemannSiegel Zfunction
The c(s) function satisfies c(s)c(1s)=1, thus for
s=[ 1/2]+it we have
c(
1
2
+it) c(
1
2
it) = 1
so that c([ 1/2]+it)=1. Therefore, we can define
implicitly the q(t) real function as
e^{iq(t)} = c(
1
2
+it)^{1/2},
(5)
with q(0)=0. Equivalently, it easily shown that
q(t) also satisfy
q(t) = arg( p^{it/2} G( 1/4 + i t/2 ) )
(6)
where the argument is defined by continuous variation of t starting
with the value 0 at t=0.
The RiemannSiegel Zfunction (also called Hardy Zfunction
in [2]) is defined for real values of t by
Z(t) = exp(i q(t)) z(
1
2
+it).
(7)
The functional equation of the Zeta function entails
z([ 1/2]it) = c([ 1/2]it) z([ 1/2]+it)
and we easily deduce that [`Z(t)] = Z(t).
In other words, the Zfunction is purely real and the equality
Z(t) = z([ 1/2]+it) holds. This wonderful property,
illustrated in figure 1,
makes this function convenient for the study of z(s) on the
critical line.
Figure 1: The Riemann Siegel function Z(t) for values of t near
0. In dashed, the value of z([ 1/2]+it).
RiemannSiegel formula applied on the Zfunction
Multiplying (4) by exp(iq(t)) for
s=[ 1/2]+it, we obtain
Z(t)
=
e^{iq(t)}
m å
n=1
n^{1/2it} + e^{iq(t)}
m å
n=1
n^{1/2+it} + O(t^{1/4})
=
2
m å
n=1
cos(q(t)tlogn)
Ön
+O(t^{1/4}).
(8)
This concise formula is very suitable for computation of Z(t) for
large values of t, as needed in computations relative to zeros of
the zeta function. For precise computations,
asymptotic formula of the error term is also available, and (using
formulation of [2]) one has for positive t
Z(t)
=
2
m å
n=1
cos(q(t)tlogn)
Ön
+
(1)^{m+1} t^{1/2}
M å
j=0
(1)^{j}t^{j} F_{j}(z) +R_{M}(t),
(9)
with R_{M}(t)=O(t^{(2M+3)/4}), where we used the notations
t =
æ Ö
t
2p
, m=ëtû, z = 2 (tm)1.
The first functions F_{j}(z) are defined by
F_{0}(z)
=
cos([ 1/2]pz^{2} + [ 3/8]p)
cos(pz)
F_{1}(z)
=
1
12p^{2}
F_{0}^{(3)}(z)
F_{2}(z)
=
1
16p^{2}
F_{0}^{(2)}(z) +
1
288p^{4}
F_{0}^{(6)}(z)
The general expression of F_{j}(z) for j > 2 is quite complicated
and we refer to [6] or [10] for it.
As exposed in [2], explicit bounds have been rigorously
obtained on the error term R_{M}(t), and for t ³ 200, one has
To be able to completly approximate Z(t) thanks to
formula (9), it remains to give an
approximation of
the q(t) function which is obtained from
expression (6) using Stirling series, giving
q(t) =
t
2
log
t
2p

t
2

p
8
+
1
48t
+
7
5760t^{3}
+ ¼
(10)
Practical approximation considerations
For practical purposes in computations relative to zeros of
z(s) it is not necessary to compute precisely the zeros but just
to locate them, and using M=1 or M=2 in
formula (9)
is usually sufficient. For t around 10^{10} for example,
the choice M=1 permits to obtain an absolute precision of Z(t)
smaller than 2×10^{14}, and with M=2 the precision is
smaller than 5×10^{20}. As for the number of terms involved
in the sum of (9), it is proportional
to Öt which is much better than previous approaches without
RiemannSiegel formula which required a number of terms of order
t. For t @ 10^{10} for example, RiemannSiegel formula only
requires m @ 40,000 terms, whereas approach of
proposition 1 requires at least
@ 9×10^{9} terms.
Quite recently, Odlyzko and Schönhage in [8]
gave an algorithm that is able to compute efficiently
values of z(s+it) at multiples values of of t.
More precisely and roughly speaking,
the algorithm is able to compute T^{1/2}
approximations of z(s+it) in a range of the form
T £ t £ T+T^{1/2} in time O(T^{1/2+e}) and using
O(T^{1/2+e}) bits of storage, whereas single evaluations at
the same values with techniques presented above would require a time
proportional to T.
We do not present the algorithm in its general form and we concentrate
on evaluations on the critical line s = 1/2 (and in fact on
evaluations of the RiemannSiegel function Z(t)) that are used for
numerically checking the Riemann hypothesis. In this context, the
general and original algorithm of [8] has been
refined for practical implementation by Odlyzko
(see [9]).
We sketch here the main ideas of this approach.
2.1 The OdlyzkoSchönhage algorithm for multievaluation of
Z(t) : generalities
The algorithm spends most of the time in the computation of
F(t) =
k_{1} å
n=1
1
Ön
exp(i tlogk), k_{1} =
ê ë
æ è
t
2p
ö ø
1/2
ú û
.
(11)
From RiemannSiegel formula (9), the value of
Z(t) is easily obtained from the value of F(t) as the real part of
2exp(iq(t)) F(t) and adding remainder terms.
We thus concentrate on multievaluations of F(t).
The algorithm is divided in two steps :
first, values of F(t) are computed on a regular grid of abscissa for
t (see 2.2).
From these values, an interpolation formula permits to obtain
efficiently any value of F(t) at a certain accurracy
(see 2.3).
We describe here how OdlyzkoSchönhage algorithm permits to evaluate
efficiently approximations of the values F(t) at evenly spaced
values t=T, T+d, ¼, T+(R1)d. We will later
precise what values of d and R we should choose.
Instead of computing directly the values F(T), F(T+d),
¼, F(T+(R1)d), the key idea
of [8] is to compute their discrete Fourier
transform, defined by
u_{k} =
R1 å
j=0
F(T+jd) w^{jk}, w = exp(2ip/R),
for 0 £ k < R. The property of the inverse Fourier transform gives
F(T+jd) =
1
R
R1 å
k=0
u_{k} w^{jk}.
In the algorithm, the value of R is choosen to be a
power of two, so the values F(T+jd) are efficiently obtained
from the (u_{k}) with an FFT transform. As described by
Odlyzko in [9], this FFT takes a negligeable
portion of the time. So the problem is now reduced to computing
efficiently the (u_{k}). Using the
definition (11) of F(t) and exchanging the
order of summations, we obtain that
u_{k} = w^{k} f(w^{k}), where f(z) is defined by
f(z) =
k_{1} å
k=1
a_{k}
zb_{k}
, b_{k}=e^{idlog k}, a_{k} =
e^{iTlogk}
k^{1/2}
(1e^{iRdlogk}),
and we are now lead to compute efficiently the complex values
f(w^{k}) for 0 £ k < R. The idea is to use Taylor series
expansions for different subset of the functions in the sum.
Given a complex number z_{0} on the unit circle and given L, 0 < L < 1,
consider for example the subset I of indices in {1,2,¼,k_{1}}
such that b_{k}z_{0} > L for all k Î I. We have the Taylor series expansion
f_{I}(z) º
å
k Î I
a_{k}
zb_{k}
=
å
n ³ 0
A_{n} (z_{0}z)^{n}, A_{n} =
å
k Î I
(z_{0}b_{k})^{n1}
valid for all z on the unit circle with zz_{0} < L. When
restricting to z with zz_{0} < L/2 for example, the convergence of
the Taylor series expansion is better than (1/2)^{n} and a fixed
number of terms (say V) in this expansion is sufficient to obtain a given
absolute precision. Approximations of values f_{I}(w^{j}) for
indices j such that w^{j}z_{0} < L/2 (let us call J this set
of indices) are then obtained from this
Taylor series expansion. When the number of elements in I is much
bigger than the number of terms V used in the series expansion, this
gives a much more efficient method to compute the values
(f(w^{j}))_{j Î J}. This
idea is used in [8] as a basic ingredient of a
sophisticated method that writes f(z) in the form åf_{Il},
the (I_{l}) being a partition of {1,2,¼,k_{1}} choosen such
that Taylor series convergence of f_{Il} is well controlled.
In the previous paragraph, we have seen how to compute efficiently an
approximation of the set of values F(T), F(T+d), ¼,
F(T+(R1)d). Our goal
here is to show how to compute efficiently F(t) for any t with T £ t £ T+(R1)d using an interpolation of these precomputed
values. The approach presented in [8] was to
compute several of the derivative F^{(h)}(T+jd) and to compute
desired value of F(t) by expanding in a Taylor series around the
nearest grid point. This approach has been significantly improved
in [9] where a bandlimited function interpolation is
used that permits to rely only on evaluations F(T+jd) and to
use a less dense grid.
The bandlimited function interpolation is classical in signal theory
for example. Under mild conditions, it states that a function with
a spectrum in [t,t]
G(t) =
ó õ
t
t
g(x) e^{ixt} dx
is determined by its samples at np/t with n an integer, and
more precisely we have the "cardinal series"
G(t) =
+¥ å
n=¥
G
æ è
np
t
ö ø
sin(ttnp)
ttnp
.
However, this series has a too slow convergence, and the idea of
Odlyzko was to use a denser sampling to obtain a faster convergence.
He applied his approach to the function
G(t) = F(t) e^{iat}, a =
logk_{1}
2
which spectrum is in [t,t] with t = [(logk_{1})/2],
and obtained the representation
F(t) =
l
b
+¥ å
n=¥
F
æ è
np
b
ö ø
e^{ia(np/bt)}
sinl(np/bt)
l(np/bt)
h(np/bt)
(12)
for any b > t, with l = (b+t)/2, and where h(u)
is a kernel function choosen as
h(u) =
c
sinh(c)
sinh(c^{2}e^{2}u^{2})^{1/2}
(c^{2}e^{2}u^{2})^{1/2}
, e = (bt)/2
where c is any fixed real positive constant.
In his computations, Odlyzko chose b @ 1.44 a (thus
d = p/b) and c=30, and observed that this choice permits
to approximate values of
F(t) with a precision of more than (roughly) 10 digits using
roughly 200 terms of the sum in (12).
The complete algorithm is presented in [9], where in
the context of zeros of the zetafunction on the critical line
s = 1/2, a lot of optimizations and practical considerations are
presented. Among the important improvements, Odlyzko applies the
multievaluation algorithm not directly to the function F(t) we have
presented, but to the function
F(k_{0},k_{1};t) =
k_{1} å
k=k_{0}
1
Ön
exp(i tlogk)
where k_{0} is a small integer that permits, in the bandlimited
function interpolation of Odlyzko, to decrease the size t of the
support spectrum of G(t) (and thus to optimize the interpolation),
with k_{0} small enough so that
the rest of the contribution å_{1 £ k < k0} takes a small
amount of time to compute. In [9] computations are done
near the zero number 10^{20}, for which a special care is
taken to the numerical imprecisions in all steps of the algorithm.
As for the addedvalue of this method compared to the classic
RiemannSiegel single evaluation method, Odlyzko estimates
in [9] that for computing zeros of z(s) on the critical
line near zero number 10^{20}, his implementation of this algorithm
is roughly 2×10^{5} faster, which proves the practical
efficiency of the method.
D. Borwein, P. Borwein, A. Jakinovski,
"An efficient algorithm for the Riemann Zeta function", (1995)
available from the CECM preprint server, URL=
http://www.cecm.sfu.ca/preprints/1995pp.html, CECM95043.
J. Borwein, D. Bradley, R. Crandall,
"Computational strategies for the Riemann Zeta Function", (1999)
available from the CECM preprint server, URL=
http://www.cecm.sfu.ca/preprints/1999pp.html, CECM98118.
A. Odlyzko, "The 10^{20}th zero of the Riemann
zeta function and 175 million of its neighbors", unpublished
manuscript, 1992. Available at URL=
http://www.research.att.com/ amo/unpublished/index.html.
J. van de Lune, H. J. J. te Riele and D. T. Winter,
On the zeros of the Riemann zeta function in the critical strip,
IV. Math. Comp. 46 (1986), 667681.