.de TP
.rs
.sp |1i
..
.EQ
delim $$
gsize 11
.EN
.S 11
.ds HP 11 11 11
.ds HF 3 3
.am DE
.sp .25
..
.am DS
.ll 4.19i
..
.am FS
.ll 4.19i
..
.ll 7i
.nr W 7i
.po -.7i
.ce 99
SUPERCOMPUTERS AND THE RIEMANN ZETA FUNCTION
.sp 2
\f2A. M. Odlyzko\f1
.sp
AT&T Bell Laboratories
Murray Hill, New Jersey 07974
.ce 0
.sp 2
.po +.2i
.nr W 4.19i
.ll 4.19i
.de TP
.rs
.sp |.3i
..
.bp
.pl 9.8i
.ls 1
.P 0
\f2ABSTRACT\f1
.sp .25
.P
The Riemann Hypothesis, which specifies the location of zeros of the Riemann zeta function, and thus describes the behavior of primes, is one of the most famous
unsolved problems in mathematics, and extensive efforts have been made over more than a century
to check it numerically for large sets of cases.
Recently a new algorithm, invented by the speaker and
A. Scho\*:nhage, has been implemented,
and used to compute over 175 million zeros near zero number $10 sup 20$.
The new algorithm turned out to be over 5 orders of magnitude faster than older methods.
The crucial ingredients in it are a rational function evaluation method similar
to the Greengard-Rokhlin gravitational potential evaluation algorithm, the FFT, and band-limited function interpolation.
While the only present implementation is on a Cray, the algorithm can easily be parallelized.
.sp
.H 1 "Introduction and History"
The Riemann zeta function is defined for complex values of $s$ with $roman Re (s) ^>^1$ by
.DS 2
.EQ (1.1)
zeta (s) ~=~ sum from n=1 to inf ~ 1 over {n sup s} ^,
.EN
.DE
and it can be continued analytically to the entire complex plane with the exception of $s=1$, where it has a first order pole.
The zeta function was actually first defined by Euler in the first half of the eighteenth century\ [2,\|26].
Euler's work was motivated by the problem of evaluating
.DS 2
.EQ
sum from n=1 to inf ~ 1 over {n sup 2} ^,
.EN
.DE
(which is $zeta (2)$ in our notation), which was posed in the seventeenth century by Mengoli.
Euler eventually showed that $zeta (2)^=^ pi sup 2 /6$, but some of his initial efforts went into numerically evaluating $zeta (2)$ and involved development of what is now called the Euler-Maclaurin formula.
This was the first of many connections between numerical analysis and the zeta function.
.P
Euler found some
important relations
between the zeta function and primes, but it was only Riemann in the 1850's who showed the full extent of these connections.
In particular, Riemann showed that the distribution of primes is determined by the
\f2nontrivial zeros\f1 of
$zeta (s)$;
i.e., those zeros of $zeta (s)$ that lie in the strip $0^<^roman Re (s) ^<^1$ (called the critical strip).
Riemann further conjectured that all of the nontrivial zeros lie on the
\f2critical line\f1, $roman Re (s)^=^1/2$.
This conjecture, known as the
\f2Riemann Hypothesis\f1 (RH), is probably the most important unsolved problem
in mathematics.
(For general background and history, see [9].)
.P
Given the importance of the RH, it is not surprising that many attempts have been made to check it numerically for various sets of zeros.
Since $rho bar$ is a zero of $zeta (s)$ whenever $rho$ is, we consider only the nontrivial zeros $rho$ with $roman Im ( rho ) ^>^0$, and number them
$rho sub 1 ,^rho sub 2 ,^"..."$, so that
.DS 2
.EQ
0 ~<~ roman Im ( rho sub 1 ) ~<=~ roman Im ( rho sub 2 ) ~<=~... ^.
.EN
.DE
(In case of multiple zeros, they have to be counted according to their
multiplicity.)
The RH is the equivalent to the claim that each $rho sub n^=^1/2^+^i ^gamma sub n$,
where $gamma sub n$ is real.
A partial list of the most important numerical verifications of the RH is given in Table\ 1.
The number in parentheses
refers to the date of publication, and the $n$ entry denotes that the RH has been checked for the first $n$ zeros, so that $gamma sub 1 ,^gamma sub 2 ,^"...",^gamma sub n$ are real.
This means not just that the $gamma sub n$ are real to within some bound,
like $10 sup -20$, but that they are actually real.
It is possible to establish this rigorously, assuming of course that the programs and hardware are correct, due to some special properties of the zeta function.
For a fuller history of these computations, see [9,\|21].
.P
There were two main innovations that were introduced in the computations listed in Table\ 1.
One was a change in technology;
starting with the work of A. Turing, all computations were performed on
electronic digital computers, as opposed to paper and pencil or
the
electromechanical devices of earlier workers, with the latest result, the verification
of the RH for the first $1.5 times 10 sup 9$ zeros by
van\ de\ Lune,
te Riele, and Winter [17] taking about 1500 hours on a Cyber 205.
The other innovation was in algorithms, and was at least as important.
The computations of Gram, Backlund, and Hutchinson used a method for computing
the zeta function that is based on the Euler-Maclaurin formula.
It is effective, but not very efficient, since to compute a single value of $zeta (1/2 ^+^it)$ with this method takes on the order of $t$ steps.
(The \%$n$-th zero $1/2 ^+^i^gamma sub n$ has $gamma sub n ^wig ^2 pi n/ ( log ^n)$, with the $1.5 times 10 sup 9$-th zero having $gamma sub n ^approx^5 times 10 sup 8$.)
Starting with the work of Titchmarsh, a much more efficient algorithm has been employed.
As it turns out, this was not a new invention, but rather a discovery by C.\ L. Siegel [25] in Riemann's unpublished notes of a method
now known as the Riemann-Siegel formula.
This method was already used by Riemann to compute at least the first few zeros of the zeta function, and the results of these computations may have been a crucial factor in stimulating him to make the RH.
.P
The Euler-Maclaurin formula is a very general technique that applies to many summations.
In contrast, the Riemann-Siegel formula relies on special properties of the zeta function.
It is much more efficient in that only about $sqrt t$ operations are needed to evaluate $zeta ( 1/2 ^+^it )$, in contrast to the roughly $t$ operations needed by the Euler-Maclaurin method.
For computations near the \%$10 sup 9$-th zero, the gain in efficiency is about 4 orders of magnitude.
.P
Recently A. Scho\*:nhage and the author [19,\|24] have invented an even faster method for computing large sets of zeros of the zeta function.
In time roughly $sqrt T$, it enables one to compute up to $B$ zeros, where $B^<wig^sqrt T$ is approximately the storage capacity on moderately fast devices (such as magnetic disks) that is available.
This algorithm has been implemented on a Cray\ \%X-MP, and used to compute the zeros
listed in Table\ 2.
(The $N^=^10 sup 12$ entry, for example, means that 1,592,196 zeros were computed, starting with zero number $10 sup 12^-^6032$ and ending with zero number
$10 sup 12 ^+^1586163$.)
.pl 11.5i
An earlier implementation of the standard Riemann-Siegel formula method [20] calculated $10 sup 5$ zeros near zero number $10 sup 12$ in about 15 hours on the Cray\ \%X-MP.
Since $gamma sub n ^approx^2.6 times 10 sup 11$ for $n^=^10 sup 12$, and $gamma sub n ^approx^1.5 times 10 sup 19$ for $n^=^10 sup 20$, this indicates that computing $6.5 times 10 sup 7$ zeros near zero number $10 sup 20$ would have required about
.DS 2
.EQ
{6.5 times 10 sup 7} over {10 sup 5} ~ cdot ~ left ( {1.5 times 10 sup 19} over {2.6 times 10 sup 11} right ) sup 1/2 ~approx~ 5 times 10 sup 6
.EN
.DE
as long, or about $7.5 times 10 sup 7$ hours,
if one were to use the Riemann-Siegel method.
On the other hand, one run of the new algorithm determined those $6.5 times 10 sup 7$ zeros in about 250 hours (of otherwise idle time), which is about $3 times 10 sup 5$ times faster.
.P
A full description of the implementation of the new algorithm and the results that have been obtained with it are presented in [22].
The next section sketches
some of the main features of the new algorithm.
Now we will mention some of the main results from [22].
First of all, no counterexamples to the RH were found!
Further, extensive statistics about the zeros were collected.
The main purpose for this (and, in fact, for the whole project) was to
obtain more insight into the behavior of the zeros of the zeta function.
The fact that the RH holds for the first $1.5 times 10 sup 9$ zeros, as well as for over $3 times 10 sup 8$ zeros at heights as far up as the $2 times 10 sup 20$-th zero, is not entirely convincing evidence in favor of the RH, since there are many
examples of number theoretic conjectures that have been shown to be false,
but where the smallest known counterexamples are very high up (up to $exp ( exp (65))$ in one case [23]).
Thus it is important to examine more sophisticated heuristics as to whether
the RH is likely to be true,
and computational results are very valuable in this regard.
(See [21] for a full discussion.)
In particular, there is a set of conjectures, going back to Hilbert and Po\*'lya,
which imply the RH, and, when combined imaginatively with some observations of Dyson and H. Montgomery might lead one to expect that the zeros of the zeta function behave like eigenvalues of random Hermitian matrices.
More precisely, if we let
.DS 2
.EQ
delta sub n ~=~ ( gamma sub n+1 ^-^gamma sub n ) ~ {log ^( gamma sub n / ( 2 pi ))} over {2 pi} ^
.EN
.DE
be the normalized gap between consecutive zeros of the zeta function,
then the conjecture is that the $delta sub n$ behave like the gaps between
successive normalized eigenvalues in the Gaussian unitary ensemble (\f2GUE\f1),
which has been investigated very extensively by physicists [5,\|6,\|18].
In the case of the GUE, it is known that the $delta sub n$ have a particular
(complicated) distribution.
Fig.\ 1 compares the distribution of the $delta sub n$ for the first $10 sup 6$ zeros of the zeta function (scatterplot) to the GUE distribution (solid line).
As can be seen, the agreement is only moderate.
On the other hand, Fig.\ 2 shows a similar graph based on about $10 sup 6$ zeros near zero number $2 times 10 sup 20$, and here the agreement as almost perfect.
This might be taken as evidence that the highly speculative
Hilbert and Po\*'lya conjectures are indeed correct, and that the RH is true.
.P
The papers [21,\|22] contain many more statistical studies of the zeros and discussions of what they might mean for the truth of the RH.
At this point we will only mention that there are some interesting conjectural connections between
zeta function zeros and quantum chaos [4].
.H 1 "Algorithms and Implementations"
The new algorithm is based on the Riemann-Siegel formula.
In the standard implementation, almost all of the computing time is devoted to evaluating sums of the form
.DS 3
.EQ (2.1)
F(t) mark ~=~ 2 ~ sum from k=2 to {k sub 1} ~ k sup -1/2 ~ exp ( it ^log^k ) ^,
.EN
.sp .5
.nr Eq 1
.EQ where
~~
.EN
.sp .5
.nr Eq 0
.EQ (2.2)
k sub 1 lineup ~=~ left floor (t/( 2 pi )) sup 1/2 right floor ^.
.EN
.DE
Once $F(t)$ (or, actually, the $F(t) ^exp (-i theta (t))$, where $theta (t)$ is a certain easy to evaluate function) is computed,
$zeta ( 1/2 ^+^it )$ is very easy to obtain.
.P
For values of $t$ near zero number $10 sup 20$, $k sub 1 ^approx^1.5 times 10 sup 9$, so the sum in (2.1) is very long,
and it is desirable to avoid evaluating that sum term-by-term for each value of $t$ that one might wish to investigate.
The first crucial ingredient in the new algorithm is to observe that if one can
precompute $F(t)$ at an evenly spaced grid of points, $t^=^T$,
$T+ delta$, $T^+^2 delta ,^"...",^T^+^(R-1) delta$, which is dense
enough, then one can compute $F(t)$ at any point in that interval from the precomputed values.
The initial method that was proposed for doing this involved Taylor series
expansions [24].
However, a much more efficient method can be obtained by using band-limited
function interpolation, and it was used in the implementation [22].
It relies on the fact, well known to communication engineers and complex
analysts, that a band-limited function $G(t)$,
.DS 2
.EQ (2.3)
G(t) ~=~ int from {- tau} to tau ~ g(x) e sup ixt ^dx ^,
.EN
.DE
is determined by its samples at the points $n pi / tau$, $n$ running through the integers, provided only that it satisfies some mild growth conditions
[7,\|14,\|15].
In fact, one has the classical ``cardinal series''
.DS 2
.EQ (2.4)
G(t) ~=~ sum from {n= - inf} to inf ~ G left ( {n pi} over tau right ) ~
{sin ( tau t - n pi )} over {tau t - n pi} ^.
.EN
.DE
This series is not suitable for zeta function computations because of slow convergence, but it is possible to obtain a similar series with the $( sin ^u )/u$ kernel replaced by one that drops off much more rapidly,
provided one is willing to compute samples $G( n pi / beta )$ for some $beta ^>^tau$.
As is shown in [22], this yields a very efficient method for computing $G(t)$ from the $G( n pi / beta )$.
.P
We next consider the problem of evaluating $F(t)$, $F(T+ delta ) ,^"...",^F(T+ (R-1) delta )$.
We need to find a method that requires substantially fewer than $Rk sub 1$ operations to compute all these values.
The first step is to apply the discrete Fourier transform.
We let
.DS 2
.EQ (2.5)
h sub k ~=~ sum from j=0 to R-1 ~ F(T + j delta ) ~exp ( 2 pi ij / R) , ~~0^<=^k^<^R ^.
.EN
.DE
If $R$ is chosen appropriately (say $R^=^2 sup r$), then the $F(T+ j delta )$ can be computed efficiently from the $h sub k$ using the FFT.
(In the current implementation, the FFT routines that were used were those of Bailey [3].
They consume a negligible fraction of the total computing time.)
Thus we reduce to the problem of computing all the $h sub k$ efficiently.
.P
To compute the $h sub k$, which consumes most of the time, we note that when we put in the definition (2.1) of $F(t)$ into (2.5), we obtain
$h sub k ^=^ h ( exp ( 2 pi ik / R ))$, where
.DS 3
.EQ
h(z) mark ~=~ sum from k=2 to {k sub 1} ~ {a sub k} over {z- b sub k} ^,
.EN
.EQ (2.6)
~~
.EN
.EQ
b sub k lineup ~=~ exp ( i delta ^log^k ) ^,
.EN
.DE
and the $a sub k$ are certain complex constants.
The straightforward term-by-term evaluation of the series (2.6) takes $approx^k sub 1$ steps for each $k$, which takes $approx^Rk sub 1$ steps in all, and offers no saving over the standard evaluation of the Riemann-Siegel formula.
However, there is a way to compute all of $h sub k$ simultaneously in
$O(k sub 1 ( log^k sub 1 ) sup 2 )$ operations
(for $R^<wig^k sub 1$) [24].
This method relies on Taylor series expansions of partial sums of the sum in (2.1).
It is explained in detail in [22,\|24].
.P
The method of [24] for evaluating the rational function $h(z)$ of (2.6) simultaneously is very similar in spirit (although not in form) to the Greengard-Rokhlin algorithm for evaluating gravitational and Coulomb potentials [11] (see also [1,8,12,13]), and it can be adapted to handle such problems.
For gravitational potential problems, it appears that the algorithms of [11] and [24] are of comparable efficiences.
However, as is discussed in [22], it seems likely that for computing zeros of the zeta function,
the Greengard-Rokhlin algorithm could turn out to be faster, due to the special structure of the zeta function problem and space limitations.
.P
Although the current implementation of the new algorithm is over five orders of
magnitude faster than the standard Riemann-Siegel method, substantial further
improvements are possible.
One could come from using the Greengard-Rokhlin algorithm.
Others could come from more careful coding, including the use of assembly language.
(Right now all the programs are written in Fortran, and most of the important
loops are vectorized by the Cray Fortran compiler, with the most important loop operating at 100 Mflops/sec, but given the mix of arithmetic operations, this could
probably be increased by using assembly language and rearranging the computation.)
Perhaps the largest improvement that can be obtained would come from use of larger disk storage.
Main memory is not a major constraint (the basic program uses under 7 Mbytes), but lack of temporary
disk storage is.
By using larger devices (which, if one were willing to do extra work, could
even be magnetic tapes or optical disks) one could obtain substantially higher speeds.
.P
While the current program has only been implemented on a Cray \%X-MP (and its has also been tested on a Cray-2), it can easily be parallelized.
There are already some parallel implementations of the Greengard-Rokhlin algorithm
[10,\|16,\|27,\|28], and
so similar ones could be carried out for the zeta function algorithm.
.P
One interesting feature of the new computations is the problem of accuracy.
The current computations are not completely rigorous due to incomplete control
over roundoff errors.
This problem is due to the extremely large heights at which the computations are being carried out and to the 64-bit word size on the Cray.
Basically, double precision arithmetic is being used, but even that is not sufficient.
Larger word sizes (72 or 80 bits) would be very helpful here!
.sp
.HU "REFERENCES"
\0\0
.sp -1
.AL
.LI (1)
J. Ambrosiano, L. Greengard, and V. Rokhlin, The fast multipole method
for gridless particle simulations, \f2Computer Physics Comm. 48\f1 (1988),
117-125.
.LI (2)
R. Ayoub, Euler and the zeta function, \f2Amer. Math. Monthly 81\f1 (1974),
1067-1086.
.LI (3)
D. H. Bailey, A high-performance FFT algorithm
for vector supercomputers,
\f2Intern. J. Supercomputer Appl. 2, no. 1,\f1 (1988), 82-87.
.LI (4)
M. V. Berry, Quantum chaology, \f2Proc. Royal Soc. London A 413\f1 
(1987), 183-198.
.LI (5)
O. Bohigas and M.-J. Giannoni, Chaotic motion and random
matrix theories, pp. 1-99 in \f2Mathematical and Computational
Methods in Nuclear Physics\f1, J. S. Dehesa,
J. M. G. Gomez, and A. Polls, eds., Lecture Notes in Physics
#209, Springer 1984.
.LI (6)
T. A. Brody, J. Flores, J. P. French, P. A. Mello, A. Pandey,
and S. S. M. Wong, Random-matrix physics: spectrum and strength
fluctuations, \f2Rev. Modern Physics 53\f1 (1981), 385-479.
.LI (7)
P. L. Butzer, W. Splettsto\*:sser, and R. L. Stens, The sampling
theorem and linear prediction in signal analysis, 
\f2Jber. d. Deutschen Math.-Verein. 90\f1 (1988), 1-70.
.LI (8)
J. Carrier, L. Greengard, and V. Rokhlin, A fast adaptive algorithm
for particle simulations, \f2SIAM J. Stat. Sci. Comp. 9\f1 (1988), 669-686.
.LI (9)
H. M. Edwards, \f2Riemann's Zeta Function\f1, Academic Press, 1974.
.LI (10)
L. Greengard and W. D. Gropp, A parallel version of the fast
multipole algorithm, to be published.
.LI (11)
L. Greengard and V. Rokhlin, A fast algorithm for particle simulations,
\f2J. Computational Phys. 73\f1 (1987), 325-348.
.LI (12)
L. Greengard and V. Rokhlin, Rapid evaluation of potential fields in
three dimensions, to be published.
.LI (13)
L. Greengard and V. Rokhlin, On the efficient implementation of the
fast multipole algorithm, to be published.
.LI (14)
J. R. Higgins, Five short stories about the cardinal series,
\f2Bull. Amer. Math. Soc. 12\f1 (1985), 45-89.
.LI (15)
A. J. Jerri, The Shannon sampling theorem - its various
extensions and applications: a tutorial review, 
\f2Proc. IEEE 65\f1 (1977), 1565-1596.
.LI (16)
J. Katzenelson, Computational structure of the N-body problem,
\f2SIAM J. Stat. Sci. Comp.\f1, to appear.
.LI (17)
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., \f2Math. Comp. 46\f1 (1986), 667-681.
.LI (18)
M. L. Mehta, \f2Random Matrices\f1, Academic Press, 1967.
.LI (19)
A. M. Odlyzko, New analytic algorithms in number theory,
pp. 466-475 
in \f2Proc. Intern. Congress Math. 1986\f1, Amer. Math. Soc. 1987.
.LI (20)
A. M. Odlyzko, On the distribution of spacings between zeros of
the zeta function, \f2Math. Comp. 48\f1 (1987), 273-308.
.LI (21)
A. M. Odlyzko, Zeros of the Riemann
zeta function: Conjectures and computations, manuscript in preparation.
.LI (22)
A. M. Odlyzko, The 
$ 10 sup 20 $-th zero of the Riemann zeta function
and 175 million of its neighbors,
manuscript in preparation.
.LI (23)
A. M. Odlyzko and H. J. J. te Riele, Disproof of the
Mertens conjecture, \f2J. reine angew. Math. 357\f1 (1985), 138-160.
.LI (24)
A. M. Odlyzko and A. Scho\*:nhage, Fast algorithms for
multiple evaluations of the Riemann zeta function, \f2Trans.
Amer. Math. Soc. 309\f1 (1988), 797-809.
.LI (25)
C. L. Siegel, U\*:ber Riemanns Nachlass zur analytischen
Zahlentheorie, \f2Quellen und Studien zur Geschichte der Math. Astr. Phys. 2\f1
(1932), 45-80.
Reprinted in C. L. Siegel, \f2Gesammelte Abhandlungen\f1, Springer,
1966, Vol. 1, pp. 275-310.
.LI (26)
A. Weil, Prehistory of the zeta-function, 
in \f2Number Theory, Trace Formulas, and Discrete Groups,\f1
K. E. Aubert, E. Bombieri, and D. M. Goldfeld, eds.,
Proc. 1987 Selberg Symposium,
Academic Press, 1988, pp. 1-9.
.LI (27)
Feng Zhao, An $O(N)$ algorithm for three-dimensional $N$-body simulations,
MIT AI Lab. report # 995, October 1987.
.LI (28)
Feng Zhao and L. Johnsson, The parallel multipole method on the
Connection Machine, paper in preparation.
.LE
.sp
.DS
.TS
center;
c s
c s
c | c
l | n.
Table\ 1.\0 Numerical verifications of the Riemann Hypothesis
for the first $n$ zeros.
.sp
Investigator	$n$
.sp .3
_
.sp .3
Gram (1903)	10
Backlund (1914)	79
Hutchinson (1925)	138
Titchmarsh et al. (1936)	1,041
Turing (1953)	1,104
Lehmer (1956)	25,000
Meller (1958)	35,337
Lehman (1966)	250,000
Rosser et al. (1969)	3,500,000
Brent (1979)	81,000,001
van de Lune et al. (1986)	1,500,000,000
.TE
.DE
.sp 2
.DS
.TS
center;
c s s
c s s
c | c | c
c | c | c
c | n | l.
Table\ 2.\0 Large computed sets of zeros of the
Riemann zeta function.
.sp
		index of first
$N$	number of zeros	zero in set
.sp .3
_
.sp .4
$10 sup 12$	1,592,196	$N^-^ 6,032$
.sp .5
$10 sup 14$	1,685,452	$N^-^ 736$
.sp .5
$10 sup 16$	16,480,973	$N^-^5,946$
.sp .5
$10 sup 18$	16,671,047	$N^-^8,839$
.sp .5
$10 sup 19$	16,749,725	$N^-^13,607$
.sp .5
$10 sup 20$	175,587,726	$N^-^30,769,710$
.sp .5
$2 times 10 sup 20$	101,305,325	$N^-^633,984$
.TE
.DE
.sp
.HU "FIGURE CAPTIONS"
\0\0
.sp -1
.VL 8
.LI "Fig.\ 1."
Probability density of the normalized spacings $delta sub n$.
Solid line:
GUE prediction.
Scatterplot:\0 empirical data based on initial 1,000,000 zeros.
.sp
.LI "Fig.\ 2."
Probability density of the normalized spacings $delta sub n$.
Solid line:
GUE prediction.
Scatterplot:\0 empirical data based on
1,041,600 zeros near zero number $2 times 10 sup 20$.
.LE