.PH ""
.nr Pt 1
.EQ
delim $$
gsize 12
define roR % { size 7 bold up 13 | back 13 size 12 roman R } %
.EN
.am DE
.sp
..
.ds HP 12 12 12 12 12 12 12
.ds HF 3 3 3 3 3 3 3
.ls 1
.ce 99
.B
.S 12 14
ON THE DISTRIBUTION OF SPACINGS BETWEEN
ZEROS OF THE ZETA FUNCTION
.R
.sp
.I
A. M. Odlyzko
.R
.sp .5
AT&T Bell Laboratories
Murray Hill, New Jersey
.sp 2
.B
ABSTRACT
.R
.ce 0
.sp 1.5
.P
.ls 2
A numerical study of the distribution of spacings between
zeros of the Riemann zeta function is presented.
It is based on values for the first $10 sup 5$ zeros and
for zeros number $10 sup 12 ~+~ 1$ to $10 sup 12 ~+~ 10 sup 5$
that are accurate to within $+-~ 10 sup -8$,
and which were calculated
on the Cray-1 and Cray X-MP computers.
This study tests the Montgomery
pair correlation conjecture as well as some further conjectures
that predict that the zeros of the zeta function behave
similarly to eigenvalues of random hermitian matrices.
Matrices of this type are used in modeling energy levels in
physics, and many statistical properties of their
eigenvalues are known.
The agreement between actual statistics for zeros of the zeta function and conjectured results is
generally good, and improves at larger heights.
Several initially unexpected phenomena
were found in the data and some were explained by relating
them to the primes.
.bp
.ls 1
.ce 99
.B
ON THE DISTRIBUTION OF SPACINGS BETWEEN
ZEROS OF THE ZETA FUNCTION
.R
.sp
.I
A. M. Odlyzko
.R
.sp .5
AT&T Bell Laboratories
Murray Hill, New Jersey
.ce 0
.sp
.H 1 "Introduction"
.ls 2
The Riemann Hypothesis (RH) has been of central interest
to number theorists for a long time, and many unsuccessful
attempts have been made to either prove or disprove it by analytic methods.
A series of large computations
have also been made, culminating in the recent verification
[52] that the RH holds for the first $1.5 ~times~ 10 sup 9$
zeros, an effort that involved over a thousand hours on
a modern supercomputer.
.pn 2
.PH "''- \\\\nP -''"
.P
A proof of the RH  would lead to tremendous improvements
in estimates for various arithmetic functions, such as the
difference between $pi (x)$, the number of primes $<= ^x$, and
$roman Li ^(x)$, the logarithmic integral of $x$.
However, many other questions would remain open, such as
the precise magnitude of the largest gap between consecutive
primes below a given bound.
Answers to such questions
depend on a much more detailed knowledge of the distribution of
zeros of the zeta function than is given by the RH.
Relatively little work has been devoted to the precise
distribution of the zeros.
The main reason for the lack of research in this area was undoubtedly the
feeling that there was little to be gained from studying problems
harder than the RH if the RH itself could not be proved.
.P
A major step towards a detailed study of the distribution of zeros
of the zeta function was made by H. L. Montgomery [56, 57].
Under the assumption of the RH, he showed that if we define
.DS 3
.EQ (1.1)
F( alpha ^,^ T) ~=~ 2 pi ( T ^ log ^ T ) sup -1 ~sum from
{{pile {0 ^<^ gamma ^<=^ T above 0 ^<^ gamma prime ^<=^ T}}} ~
T sup {i alpha ^( gamma - gamma prime )} ~4 over {4 ^+^ ( gamma - gamma prime
) sup 2}
.EN
.DE
for $alpha$ and $T$ real, $T ~>=~ 2$, where $1/2 ~+~ i gamma$ and
$1/2 ~+~ i gamma prime$ denote nontrivial zeros of the zeta
function, then
.DS 3
.EQ (1.2)
F( alpha ^,^ T) ~=~ (1 ^+^ o(1)) T sup {-^ 2 alpha} ~ log ^T ~+~ alpha ~+~
o(1) ~~roman as~~ T ~->~ inf ^,
.EN
.DE
uniformly for $0 ~<=~ alpha ~<=~ 1$.
Montgomery also presented heuristic arguments which suggested that
.DS 3
.EQ (1.3)
F( alpha ^,^ T) ~=~ 1 ^+^ o(1) ~~roman as~~ T ~->~ inf
.EN
.DE
uniformly for $alpha ~epsilon~ [a ^,^ b]$, where
$1 ~<=~ a ~<~ b ~<~ inf$ are any constants.
If the conjecture (1.3) were true, then the following estimate,
known as the Montgomery pair correlation conjecture [56] would follow:
.sp
.I
Conjecture.
For fixed $0 ~<~ alpha ~<~ beta ~<~ inf$,
.DS 3
.EQ (1.4)
{| "{" ( gamma ^,^ gamma prime ): ~0 ~<~ gamma ^,^ gamma prime ~<=~ T,
~2 pi alpha ( log ^T ) sup -1 ~<=~ gamma ^-^ gamma prime ~<=~ 2 pi beta
( log ^ T ) sup -1 "}" |} over {T over {2 pi} ~ log ^T}
.EN
.sp
.EQ
wig~
int sub alpha sup beta ~left ( 1 ~-~ left ( {sin~ pi u} over {pi u}
right ) sup 2 right ) ~du
.EN
.DE
as $T ~->~ inf$.
.R
.P
The above conjecture is very striking, particularly in predicting
that small gaps between zeros of the zeta function occur very infrequently.
This behavior is very different from that of many other number
theoretic functions.
The number of primes in short intervals, for example, is observed
experimentally and is conjectured theoretically to be
distributed like a Poisson random variable [29].
.P
The conjecture (1.4)
suggests some further topics for investigation.
In the language of mathematical physics, this conjecture
says that $1 ~-~ (( sin~ pi u ) /( pi u ) ) sup 2$
is the pair correlation function of zeros of the zeta
function.
F. J. Dyson pointed out [56] that the gaussian unitary
ensemble (GUE) has the
same pair correlation function.
The GUE has been studied very extensively in mathematical
physics together with the gaussian orthogonal ensemble (GOE)
and the gaussian symplectic ensemble (GSE) as models for
distribution of energy levels in many-particle systems.
The GUE consists of $n ^times^ n$ complex hermitian
matrices of the form $A ~=~ ( a sub jk )$, where
$a sub jj ~=~ 2 sup 1/2 ~sigma sub j,j$,
$a sub jk ~=~ sigma sub j,k ~+~ i ^eta sub j,k$ for
$j ~<~ k$, and $a sub j,k ~=~ {a bar} sub k,j ~=~ sigma sub k,j ~-~ i^ eta sub k,j$
for $j ~>~ k$, where the $sigma sub j,k$ and $eta sub j,k$ are independent
standard normal variables.
When $n ~->~ inf$ and the eigenvalues of the matrices
of the GUE are suitably normalized, their pair correlations
becomes $1 ~-~ (( sin~ pi u)/( pi u ) ) sup 2$.
(The pair correlations of the GOE and GSE are different.)
For more information about the GUE and the other ensembles,
the reader is referred to [2, 4, 11, 41, 54, 59, 63].
.P
The possible connection between zeros of the zeta function and
eigenvalues of random matrices is of interest in number theory
because of the Hilbert and Po\*'lya conjectures [56, 59] which say
that the zeros of the zeta function correspond to eigenvalues
of a positive linear operator.
If true, the Hilbert and Po\*'lya conjectures would imply the RH,
and some people feel that the most promising way to prove
the RH is by finding the right operator and establishing
the necessary results about it.
One can argue heuristically that if such an operator exists,
its eigenvalues might behave like those of a random operator,
which in turn might behave like the limit of a sequence of
random matrices, and that therefore the zeros of the zeta
function might be distributed like the eigenvalues of a large
random matrix.
For more on these very vague conjectures, see [3, 59].
Since Montgomery's result (1.1) is consistent with the GUE
distribution of eigenvalues, one could then argue that the
operator associated to the zeta function ought to be complex
hermitian,
and that GUE properties of eigenvalues might apply at
least approximately to the zeros of the zeta function.
.P
Since Montgomery's proof of (1.2) (which depends on the RH) was
published, several other results have been obtained.
Ozluk [62] has considered a $q$-analogue of Montgomery's method
and showed, roughly speaking, that if one considers a function
similar to the $F( alpha ^,^ T)$ of (1.1), but where one sums over
zeros of many Dirichlet $L$-functions, then the analogue of
Montgomery's conjecture (1.3) is true for $1 ~<=~ alpha ~<=~ 2$.
For other results about Montgomery's pair correlation conjecture
and its consequences, see [30, 31, 34, 35, 36, 42, 44, 59].
Overall, though, there is very little solid theoretical evidence
in favor of even the Montgomery pair correlation conjecture (1.4),
and essentially none for the speculative conjectures discussed
above that link zeros of the zeta function to eigenvalues
of random hermitian matrices.
.P
One feature of the speculative arguments mentioned above
is that they propose a connection between zeros of the zeta
function, about which relatively little is known theoretically
and the GUE, which has been investigated very
extensively and abounds in specific predictions.
This makes it possible to compare data for the zeros
of the zeta function against the theoretical predictions of
the GUE.
The purpose of this paper is to report
the results of such empirical tests of the pair correlation
conjecture and other GUE predictions.
.P
The agreement between the predictions of the GUE theory and
actual computed behavior of the zeros of the zeta function is
quite good, especially when one takes account of the fact that
the zeta function on the critical line approaches its asymptotic
behavior very slowly, as is explained in Section 2.
However, there are also several features of the data for
the zeros that are not predicted by the GUE theory.
One is that there are long-range correlations between spacings
of consecutive zeros.
These can be explained in terms of the primes through the use
of the ``explicit formulas'' of prime number theory.
Another is that small spacings between zeros are somewhat
more common among high zeros than is expected.
The deviation between the predicted and observed numbers is
not large, but it is surprising, since it shows up quite early.
It would be very desirable to carry out further computations
at much greater heights to see if this phenomenon persists there.
.P
This is the first study that computed large numbers of high
zeros of the zeta function to substantial accuracy.
The large scale efforts to verify the RH numerically were
designed only to prove that the zeros under investigation lie
on the critical line, and did not produce accurate values
for their positions.
Some very accurate values of zeros have also been computed
in order to use them in disproving a variety of number theoretic
conjectures, but they were limited to a small number of initial
zeros (such as the roughly 100 decimal digit values for
the first 2000 zeros of [60]).
For a numerical investigation of the Montgomery pair correlation
conjecture and the other conjectures alluded to above, it was
desirable to obtain a larger sample of zeros but only to medium
accuracy.
Furthermore, since there are many properties of the zeta
function that are true only asymptotically, and the asymptotic
behavior is often approached very slowly, it seemed desirable
to compute a sample of zeros very high up.
The computations on which this paper are based determined the
values of the first $10 sup 5$ zeros and also of zeros number
$n$, $10 sup 12 ~+~ 1 ~<=~ n ~<=~ 10 sup 12 ~+~ 10 sup 5$, each
to within $+-^ 10 sup -8$, using the \%Cray-1 and Cray \%X-MP computers, respectively.
In the second case it was also verified that
those zeros satisfy the RH.
Zero number $10 sup 12$ is at height approximately $2.7 ~times~ 10 sup 11$,
and given the available computational resources and the method
used it was about as high as one could go.
(A more efficient algorithm for evaluating the zeta function
was invented recently [61], but it is quite complicated and
has not yet been implemented.)
The computations used on the order 20\ hours on the
Cray \%X-MP, and the capabilities of this modern supercomputer
were essential for the project.
Several other samples of $10 sup 5$ zeros each at large heights
had been computed earlier on a \%Cray-1, but due to a defect in the manufacturer's
software, described at the end of Section 3, some of their
accuracy was lost, and so they are used for only some of the
statistical analyses.
.P
The comparison of the empirical data for the zeros of the zeta
function to the results that are proved for the GUE is made
in Section 2.
Section 3 is devoted to a description of the computation of
the zeros of the zeta function and the associated error analysis.
Section 4 gives some statistical information about the zeros
that were found, as well as of some other zeros.
This information is not directly relevant for the purposes
of Section\ 2, but are of interest because of the
comparison with the earlier computations of [10, 52].
Section 5 reports on several alternative techniques for computing
the zeros that might be of use for other computations.
Finally, Section 6 describes the computation of the GUE predictions.
.H 1 "Zeros and eigenvalues"
Before presenting the results of the empirical study of the zeros,
we first recall some standard notation, and then discuss some
known results about the zeta function and how they might
affect our interpretation of the data.
As usual, we consider only zeros $rho$ of the zeta function
with $Im ^( rho ) ~>~ 0$, and we number them $rho sub 1 , ^rho sub 2 ^,...$
(counting each according to its multiplicity) so that
$0 ~<~ Im ^( rho sub 1 ) ~<=~ Im ^( rho sub 2 ) ^<= ...^$.
All the zeros $rho sub n$ that have been computed so far are simple
and lie on the critical line, so we can write them as
$rho sub n ~=~ 1/2 ~+~ i ^gamma sub n$,
$gamma sub n ~epsilon~ roR sup +$.
We have $gamma sub 1 ~=~ 14.134 ^...,$
$gamma sub 2 ~=~ 21.022 ^...,$ etc.
.P
We let
.DS 3
.EQ (2.1)
theta (t) ~mark =~ roman arg ^ [ pi sup {-^ it/2} ~GAMMA ( 1/4 ~+~ it/2 ) ] ^,
.EN
.sp
.EQ (2.2)
S(t) ~lineup =~ pi sup -1 ~ roman arg ~zeta ( 1/2 ~+~ it ) ^,
.EN
.DE
where in both cases the argument is defined by continuous variation
of $s$ in $pi sup {-^ s/2} ~GAMMA ( s/2 )$ or $zeta (s)$, respectively,
starting at $s ~=~ 2$, going up vertically to $s ~=~ 2 ~+~ it$, and
then horizontally to $s ~=~ 1/2 ~+~ it$, and where we assume (in
the case of $S(t)$) that there are no zeros $rho$ with
$Im ( rho ) ~=~ t$.
If $N(t)$ denotes the number of zeros
$rho$ with $0 ~<~ Im ( rho ) ~<~ t$ (counted according to their
multiplicity), then [68; Ch. 9.3]
.DS 3
.EQ (2.3)
N(t) ~=~ 1 ~+~ pi sup -1 ~theta (t) ~+~ S(t) ~.
.EN
.DE
By Stirling's formula [40] we have
.DS 3
.EQ (2.4)
theta (t) ~=~ 1 over 2 ^t ~ log (t/( 2 pi e )) ~-~ pi /8 ~+~ O( t sup -1 )
~~roman as~~ t ~->~ inf ^,
.EN
.DE
so that
.DS 3
.EQ (2.5)
N(t) ~=~ t over {2 pi} ~ log ^ t over {2 pi e} ~+~ 7 over 8 ~+~ O( t sup -1 )
~+~ S(t) ~.
.EN
.DE
We will discuss the function $S(t)$ at greater length below.
For the moment, though, we will mention only that the best
unconditional bound that is known for $S(t)$ is [68; Thm. 9.4]
.DS 3
.EQ (2.6)
S(t) ~=~ O( log ^ t) ~~~roman as ~ t ~->~ inf ~.
.EN
.DE
.P
In view of the bound (2.6), we see by (2.5) that $N(t)$ is almost
exactly equal to $pi sup -1 ~theta (t)$.
In particular, zeros become denser and denser the higher up one
goes in the critical strip, and the average vertical spacing
between consecutive zeros at height $t$ is
$2 pi /( log^(t/(2 pi )))$.
Therefore, in order to study quantities that are largely invariant
with height, we define the normalized spacing between consecutive
zeros $1/2 ~+~ i ^gamma sub n$ and $1/2 ~+~ i ^gamma sub n+1$ (in
cases where both zeros satisfy the RH) to be
.DS 3
.EQ (2.7)
delta sub n ~=~  ( gamma sub n+1 ~-~ gamma sub n  ) ~{log ^ {gamma
sub n} over {2 pi} } over {2 pi} ~.
.EN
.DE
It then follows from (2.5) and (2.6) that the $delta sub n$ have mean
value 1 in the sense that for any positive integers $N$ and $M$,
.DS 3
.EQ (2.8)
sum from {n ^=^ N+1} to N+M ~delta sub n ~=~ M ~+~ O( log ^ NM ) ~.
.EN
.DE
In this notation, Montgomery's pair correlation conjecture can be
reformulated to say that for any fixed $alpha , ^beta$, $0 ~<~ alpha ~<~ beta ~<~ inf$,
.DS 3
.EQ (2.9)
N sup -1 ^|^ left { (n,k): ~1 ~<=~ n ~<=~ N, ~k ~>=~ 0, ~
delta sub n ^+^ delta sub n+1 ^+...+^ delta sub n+k ~epsilon ^ [
alpha , beta ] right } ^ |
.EN
.sp
.EQ
wig~ int sub alpha sup beta ^left ( 1 ^-^
left ( {sin~ pi u} over {pi u} right ) sup 2 right ) ~du
.EN
.DE
as $N ~->~ inf$.
If the pair correlation conjecture holds, we might even expect
something stronger to hold, namely that
.DS 3
.EQ
M sup -1 ^|^ left { (n,k): ~N ^+^ 1 ~<=~ n ~<=~ N ^+^ M, ~k ~>=~ 0, ~
delta sub n ^+...+^ delta sub n+k ~epsilon ^ [ alpha , beta ] right } ^ |
.EN
.sp
.EQ (2.10)
wig~ int sub alpha sup beta ~left ( 1 ^-^ left ( {sin~ pi u} over {pi u}
right ) sup 2 right ) ~du
.EN
.DE
as $N, ~M ~->~ inf$ with $M$ not too small compared to $N$,
say $M~>=~N sup eta$ for some $eta~>~0$.
.P
The definition of the $delta sub n$ makes it easy to compare
distribution of spacings between consecutive zeros of the zeta
function and the normalized spacings between consecutive eigenvalues
in the GUE, since the latter are also normalized to have mean value 1.
The \f2GUE hypothesis\f1 will be used to refer to the
conjecture that the distribution of the $delta sub n$ approaches
asymptotically the distribution known to hold for the GUE eigenvalues.
In particular, the GUE hypothesis leads to the prediction that
the $delta sub n$ should have a particular distribution function;
for any $alpha , beta$ with $0 ~<=~ alpha ~<~ beta ~<~ inf$,
.DS 3
.EQ (2.11)
M sup -1 ^|^ left { n: ~N ^+^ 1 ~<=~ n ~<=~ N ^+^ M, ~delta sub n ~epsilon ^
[ alpha , beta ] right } ^|~ wig~ int sub alpha sup beta ~p(0,u) ~du
.EN
.DE
as $M, ~N ~->~ inf$ with $M$ not too small compared to $N$, where
$p( 0,u)$ is a certain complicated probability density function
discussed in Section 6 and graphed (with the solid line curve) in
figures 3 and 4.
Similarly, one obtains the prediction that
.DS 3
.EQ (2.12)
M sup -1 ^|^ left { n: ~N ^+^ 1 ~<=~ n ~<=~ N ^+^ M, ~delta sub n ^+^ delta
sub n+1 ~epsilon ^ [ alpha , beta ] right } ^|~ wig~ int sub alpha sup beta
~p( 1,u )~du ^,
.EN
.DE
where $p(1,u)$ is another density function.
More generally, one obtains the prediction that the $delta sub n$
approach asymptotically the behavior of a stationary (but non-Markovian)
process, so that for any $k$, the empirical joint distribution function
of $delta sub n , ^delta sub n+1 ^,...,^ delta sub n+k$ for
$N ^+^ 1 ~<=~ n ~<=~ N^+^ M$ approaches that of the GUE process.
.P
Most of the studies to be reported below are based on the $delta sub n$
for $1 ~<=~ n ~<=~ 10 sup 5$ and
$10 sup 12 ^+^ 1 ~<=~ n ~<=~ 10 sup 12 ^+^ 10 sup 5$, although a few
other sets of $10 sup 5$ consecutive $delta sub n$ for $n ~<~ 10 sup 12$
will also be referred to occasionally.
The letter $N$ will be used to designate the starting point of a
data set, so that, for example, a graph or a table entry labelled
with $N ~=~ 10 sup 11$ will refer to data for $delta sub n$ or
$delta sub n ^+^ delta sub n+1$,
$10 sup 11 ^+^ 1 ~<=~ n ~<=~ 10 sup 11 ~+~ 10 sup 5$.
The reason for selecting sets of size $10 sup 5$ was to
obtain a sample so large that random sampling errors would
not be very significant.
The reason for not computing zeros higher than $rho sub n$
for $n ~=~ 10 sup 12 ~+~ 10 sup 5$ (approximately) is
that on the machine that was available to the author (a single-processor
Cray X-MP with 2 million words of memory) such computations
would have been prohibitively slow, since the author's program,
described in Section 3,
achieves high speed by using several auxiliary storage arrays, so that
the size of available memory was a major limit.
.P
Most of the studies of the available $delta sub n$ are concerned
with the pair correlation conjecture (2.10) and the conjectured
distributions (2.11) and (2.12) for $delta sub n$ and
$delta sub n ~+~ delta sub n+1$, respectively.
The reason for this is less the difficulty of computing the GUE
predictions for other quantities than the question of significance
of the comparison of the experimental data to the theoretical predictions.
By (2.5) and (2.6) we see that the interesting behavior of spacings
between zeros is due to $S(t)$.
The bound (2.6) is the best that is known unconditionally.
Even on the assumption of the RH, it is only known that
.DS 3
.EQ (2.13)
S(t) ~=~ O left ( {log ^t} over {log ^ log ^t} right ) ~~roman as~~ t ~->~ inf ~.
.EN
.DE
The true rate of growth is likely to be much smaller, though.
It has been known [68; Ch. 14.12] for a long time that on the RH,
.DS 3
.EQ
S(t) ~=~ OMEGA sub {+-} ~(( log ^ t ) sup {1/2 ^-^ epsilon} ) ~~roman as~~
t ~->~ inf
.EN
.DE
for any $epsilon ~>~ 0$,
and Montgomery\ [58] has shown more recently
that
.DS 3
.EQ (2.14)
S(t)~=~OMEGA sub {+-}~(( log^t) sup 1/2~( log^log^t) sup -1/2 )~~
roman as~~t~->~inf~.
.EN
.DE
It is thought likely that
.DS 3
.EQ (2.15)
S(t) ~=~ O (( log ^t ) sup {1/2 ^+^ epsilon} ) ~~roman as~~ t ~->~ inf
.EN
.DE
for every $epsilon ~>~ 0$ (cf. [49])
and Montgomery has even conjectured\ [58] that (2.14) represents
the right rate of growth of $S(t)$.
For statistical studies of the kind we are undertaking, though,
the typical size of $S(t)$ is more important.
Selberg [66, 68] has shown that for all $k ~epsilon~ Z sup +$,
.DS 3
.EQ (2.16)
int sub 0 sup T ~S(t) sup 2k ~dt ~wig~ {(2k)!} over {k! ^( 2 pi ) sup 2k} ~T
^( log ^ log ^ T ) sup k ~~roman as~~ T ~->~ inf ~.
.EN
.DE
This means that a typical value of $S(t)$ is on the order of
$( log ^ log ^ t ) sup 1/2$, an extremely small quantity.
Experimentally, it is known that $| S(t) | ~<~ 1$ for
$7 ~<~ t ~<=~ 280$, and $| S(t) | ~<~ 2$ for
$7 ~<~ t ~<=~ 6.82 ~cdot~ 10 sup 6$.
The largest value of $| S(t) |$ that has been observed so far appears
to be $S(t) ~=~ 2.3136 ^...$ for some $t$ near
$gamma sub n , ~n ~=~ 1,333, ~195, ~692$ [R].
(Large values of $S(t)$ are associated to violations of Gram's
and Rosser's rules, see [10, 23, 49, 59], which explains why some statistics about them are available.)
.P
Aside from (2.6) and Selberg's estimate (2.16), several other estimates
have been proved.
For example, if
.DS 3
.EQ
S sub 1 ^(t) ~=~ int sub 0 sup t ~S(u) ~du ^,
.EN
.DE
then $| S sub 1 ^(t) | ~=~ O( log ^t)$ unconditionally, and
$| S sub 1 ^(t) | ~=~ O(( log ^t) ~( log ^ log ^t) sup -2 )$ on
the RH [68].
The true maximal order of magnitude is probably again around
$( log ^t ) sup 1/2$.
Therefore there is a lot of cancellation in the local behavior
of $S(t)$.
Many results about the distribution of values of $S(t ^+^ h) ~-~ S(t)$
are also known [25, 26, 27, 32, 33, 44, 69].
For example, Theorem 13.2 of [69] says that if
$1/2 ~<~ eta ~<=~ 1$, then for $T sup eta ~<~ H ~<=~ T$,
$0 ~<~ h ~<~ 1$,
.DS 3
.EQ
int sub T sup T+H ~left { S(t^+^h) ~mark -~ S(t) right } sup 2k ^dt ~=~
H A sub k ~left { log ^( 2 ^+^ h ^ log ^T ) right } sup k
.EN
.sp
.EQ (2.17)
~~~~~
.EN
.sp
.EQ
~lineup +~ O left ( H ^c sup k ^k sup k ^ left { k sup k ~+~
( log ^(2 ^+^ h ^ log ^T ) ) sup k-1/2 right } right ) ^,
.EN
.DE
where $A sub k ~=~ (2k )! ^{size 16 /}^ ( 2 sup k ^pi sup 2k ^k!)$, and
$c ~>~ 0$ is a constant.
There are also theorems about distribution of
$S sub 1 ^(t^+^h) ~-~ S sub 1 ^(t)$, cf. [69].
.P
Since $S(t)$ grows very slowly, we can expect low zeros of the
zeta function to be very restricted in their behavior, and their
asymptotic behavior to be approached quite slowly.
Pseudo-random behavior of the zeros (like that predicted by the GUE
hypothesis) can only be expected over ranges of about
$( log ^T ) sup -1/2$ at height $T$; i.e., over a range of
about $( log ^n ) sup 1/2$ normalized spacings from zero number $n$.
Furthermore, we can expect gaps between next-to-nearest neighbors
(i.e., $delta sub n ^+^ delta sub n+1$) to be much more constrained
than those between nearest neighbors (i.e., $delta sub n$), and
to approach their asymptotic behavior even more slowly.
This is indeed what is observed, and it's the main reason for
restricting most of the present investigation to nearest or next-to-nearest
spacings.
We will see, in fact, that when we consider very distant spacings,
a totally different phenomenon occurs, whose explanation lies not
in the GUE, but rather in prime numbers.
.P
We now proceed to a discussion of the evidence.
Figures 1 and 2 show the data for the pair correlation conjecture
for the two sets of zeros $1/2 ~+~ i ^gamma sub n$,
$N ^+^ 1 ~<=~ n ~<=~ N ^+^ 10 sup 5$, and for $N ~=~ 0$ and
$N ~=~ 10 sup 12$, respectively.
For each interval $I ~=~ [ alpha , beta )$, with
$[ alpha , beta ) ~=~ [0, ~0.05)$, $[0.05, ~0.1) ^,...,^ [2.95, ~3)$,
a star is plotted at the point $x ~=~ ( alpha ^+^ beta )/2, ~y ~=~ a sub {alpha , beta}$,
where
.DS 3
.EQ
a sub {alpha , beta} ~=~ 20 over {10 sup 5} ^ |^  left { (n,k): ^N ^+^ 1 ~<=~ n
~<=~ N ^+^ 10 sup 5 , ^k ~>=~ 0, ~delta sub n ^+...+^ delta sub n+k ~epsilon
^[ alpha , beta )  right } ^ | ~.
.EN
.DE
The solid lines in both figures are the GUE prediction
$y ~=~ 1 ~-~ (( sin ^ pi^ x) / ( pi x )) sup 2$.
Note that for a typical value of $alpha~>=~1$,
$a sub {alpha ,^beta}$ is derived from
about 5000 points $delta sub n~+~"..."~+~delta sub n+k$ so the
random sampling error is expected to be on the order of 0.015.
All of the plots in this paper, as well as most of the statistical
computations, were done using the $S$ system [1].
.P
Figures 3 and 4 show the data on the distribution of the normalized
spacings $delta sub n$.
The stars again correspond to a histogram of the $delta sub n$;
for each interval $[ alpha , beta )$, $alpha ~=~ k/20, ~beta ~=~ alpha ~+~ 1/20$,
a star is plotted at $x ~=~ ( alpha ^+^ beta )/2$,
$y ~=~ b sub {alpha , beta}$, where
.DS 3
.EQ
b sub {alpha , beta} ~=~ 20 over {10 sup 5} ^|^ left { n:~ N ^+^ 1 ~<=~ n
~<=~ N ^+^ 10 sup 5 , ~delta sub n ~epsilon ^[ alpha , beta ) right } ^| ~.
.EN
.DE
The solid lines are the GUE predictions, $y ~=~ p( 0,x)$.
(As in explained in Section\ 6, no rigorous error analyses is
available for the computed values of $p(0,^x)$, but they
are thought to be quite accurate.)
Similarly, figures 5 and 6 show the data on the distribution of
$delta sub n ^+^ delta sub n+1$.
.P
The graphs at first glance show a satisfying amount of agreement
between the experimental data and the GUE predictions.
Graphs have roughly the same shape, and agreement between data
and prediction improves dramatically as one goes from the first
$10 sup 5$ zeros to the $10 sup 5$ zeros with
$10 sup 12 ^+^ 1 ~<=~ n ~<=~ 10 sup 12 ^+^ 10 sup 5$.
Moreover, the disagreement between data and prediction is concentrated
precisely where our discussion of the effect of the small size
of $S(t)$ would lead
one to expect it; there are fewer very large or very small values
of $delta sub n$ and $delta sub n ^+^ delta sub n+1$ in the data
than predicted, and more values near the mean.
.P
Further comparison of the distribution of zeros to the GUE
prediction is provided by Table 1, which shows values of
.DS 3
.EQ (2.18)
10 sup -5 ~sum from {n ^=^ N+1} to {N ^+^ 10 sup 5} ~( delta sub n ^-^ 1 ) sup k
.EN
.DE
and
.DS 3
.EQ (2.19)
10 sup -5 ~sum from {n ^=^ N+1} to {N ^+^ 10 sup 5} ~( delta sub n ^+^
delta sub n+1 ~-~ 2 ) sup k
.EN
.DE
for the two sets $N ~=~ 0$ and $N ~=~ 10 sup 12$, as well as the
GUE values.
(Numbers in Table 1, as well as others in this paper that don't have
"$ "..." $" at the end are usually rounded to the number of digits shown,
whereas those with "$ "..." $" are truncated.)
Table 2 shows moments of $log ^delta sub n , ~delta sub n sup -1$,
and $delta sub n sup -2$.
.P
The histograms of figures 1-6 and the moments of tables\ 1-2 provide
a very rough measure of the degree to which zeros of the zeta
function match the GUE predictions.
One can obtain a better quantitative measure of the degree of fit
between the observed and the predicted distributions using the
histogram data (e.g., one can use the $chi sup 2$-test [47] or
the asymptotic results of [24]).
In our case, since we know the predicted distribution quite well
(at least numerically), we will use the quantile-quantile $(q-q)$
plot to detect deviations between the empirical and theoretical
distributions.
Given a sample $x sub 1 ^,...,^ x sub n$, and a continuous
cumulative distribution function $F(z)$ for some distribution,
the theoretical $q-q$ plot is obtained by plotting $x sub (j)$ against
$q sub j$, where $x sub (1) ~<=~ x sub (2) ^<= ... <=^ x sub (n)$
are the $x sub i$ sorted in ascending order, and
the $q sub j$ are the theoretical quantiles defined by
$F( q sub j ) ~=~ (j ^-^ 1/2 ) /n$ [13].
The $q-q$ plot is a very sensitive tool for detecting differences
among distributions, especially near the tails.
If the $x sub i$ are drawn from a distribution corresponding to
$F(z)$, and the sample size $n$ is large, the $q-q$ plot will
be very close to the straight line $y ~=~ x$.
Figures 7-10 show segments of the $q-q$ plots of the $delta sub n$ and
$delta sub n ^+^ delta sub n+1$ for $N ~=~ 0$ and $N ~=~ 10 sup 12$.
.P
To obtain a quantitative measure of the agreement between
the observed distributions of the $delta sub n$ and the GUE
predictions, we use the Kolmogorov test [47; Section 30.49].
Given a sample $x sub 1 ^,...,^ x sub n$ drawn from a distribution
with the continuous cumulative distribution function $F(z)$,
we let $F sub e ^(z)$ be the sample distribution function:
.DS 3
.EQ
F sub e ^(z) ~=~ n sup -1 ^|^ left { k: ~1 ~<=~ k ~<=~ n, ~x sub k
~<=~ z right } ^| ~.
.EN
.DE
The Kolmogorov statistic is then
.DS 3
.EQ (2.20)
D ~=~ {roman {su p}} from z ~ | ~ F sub e ^(z) ~-~ F(z) ^| ~.
.EN
.DE
The asymptotic distribution of $D$ is known [47; Eq. 30.132];
if the $x sub i$ are drawn from the distribution defined by $F(z)$, then
.DS 3
.EQ (2.21)
lim from {n ^->^ inf} ~roman Prob ( D ~>~ u n sup {-^ 1/2} ) ~=~ g(u) ^,
.EN
.DE
where
.DS 3
.EQ (2.22)
g(u) ~=~ 2 ~sum from r=1 to inf ~( -1 ) sup r-1 ~ exp ( -^ 2 ^r sup 2 ^ u
sup 2 ) ~.
.EN
.DE
Table 3 gives the Kolmogorov statistic $D$ for $delta sub n$ and
$delta sub n ^+^ delta sub n+1$ for several blocks of $10 sup 5$
consecutive zeros.
Also given is an estimate of the probability that this statistic
would arise if the $delta sub n$ ($delta sub n ^+^ delta sub n+1$, resp.)
were drawn from the GUE distribution.
This estimate is obtained by evaluating
$g( D^ 10 sup 5/2 )$.
.P
The data presented so far are fairly consistent with the GUE
predictions.
The differences between computed values and the predicted ones
diminish as one studies higher zeros, and they are more
pronounced for $delta sub n ^+^ delta sub n+1$ than for $delta sub n$.
Moreover, the computed spacings $delta sub n$ and $delta sub n ^+^ delta sub n+1$
are more concentrated near their means than the GUE predictions,
which is again to be expected on the basis of the small size of
$S(t)$.
However, there are some indications, based on the tails of the
distribution of the $delta sub n$ and $delta sub n ^+^ delta sub n+1$,
that asymptotically the zeros of the zeta function might not obey
the GUE predictions.
Already in Table\ 2 we see that the mean value of $delta sub n sup -2$
for $N~=~10 sup 12$ exceeds the GUE prediction,
which indicates an excess of small $delta sub n$.
The $q-q$ plot of Fig. 8 (for $N ~=~ 10 sup 12$) is much
closer to a straight line than that of Fig. 7 (for $N ~=~ 0$), just as
that of Fig. 10 is closer to a straight line than that of Fig. 9.
However, the graphs in figures 9 and 10 are both below the
straight line $y ~=~ x$ (which corresponds to perfect agreement
between theory and prediction), which indicates that in both
cases there are fewer large spacings than the GUE hypothesis
predicts, which was to be expected.
The graph of Fig. 7 is above the straight line
$y ~=~ x$, which indicates fewer small spacings than predicted,
which again was to be expected,
while on the other hand the graph of Fig.\ 8 is below the line $y ~=~ x$, which
indicates an excess of small spacings,
which is quite unexpected.
The small size of $S(t)$ leads one to expect a substantial
deficiency in the number of small $delta sub n$ and
$delta sub n ^+^ delta sub n+1$, since the graph of $S(t)$ is
an uneven sawtooth curve, with $S(t)$ jumping up by 1 at $t ~=~ gamma sub n$,
then decreasing essentially linearly up to $t ~=~ gamma sub n+1$, jumping
up again by 1 at $t ~=~ gamma sub n+1$, and so on,
so that small gaps between zeros also correspond to
large values of $S(t)$.
Fig. 11 shows an initial part of the $q-q$ plot for the $delta sub n$
for $N ~=~ 2 ^cdot^ 10 sup 11$ which indicates behavior very similar to
that of Fig. 8.
The $q-q$ plot for $N ~=~ 10 sup 11$ is very
similar to those for $N ~=~ 2 ^cdot^ 10 sup 11$ and
$N ~=~ 10 sup 12$, while that of $N ~=~ 10 sup 9$ lies on
the other side
of the $y ~=~ x$ line.
(We should note again that the data for $N ~=~ 10 sup 9$,
$10 sup 11$, and $2 ^cdot^ 10 sup 11$ are not very accurate,
and so have to be treated with caution.
In fact, Fig. 11 provides a way to gauge the inaccuracy introduced
into those data sets by the bug described
in Section 3.
The $q-q$ plot of Fig. 11 has in places the appearance of a
staircase, with flat horizontal segments.
These flat segments arise from identical computed values of
the $delta sub n$.
The correct values of the $delta sub n$ are expected to be
all distinct and much more evenly distributed, as is the case for $N~=~0$ and $N~=~10 sup 12$.)
Table\ 4 shows the number of $delta sub n$ and $delta sub n ^+^ delta sub n+1$
that are very small or very large, as well as the number that
would be expected under the GUE for a sample of size $10 sup 5$.
Large $delta sub n$ and $delta sub n ^+^ delta sub n+1$ are generally
less numerous than the GUE predicts, which is to be expected.
The data for the high blocks of zeros 
show that
the number of $delta sub n ~<=~ 0.05$
and of $delta sub n~<=~0.1$
is in excess
of that predicted.
(We would have obtained the same results if we considered the pair
correlation conjecture (2.10) with $alpha ~=~ 0.01, ~beta ~=~ 0.05$, say.)
The excess is not very large, given that the expected number is
fairly small, but it may be significant that while the first few
sets (for $N ~=~ 0$, $5 ~cdot~ 10 sup 7$, and $10 sup 9$) all
show deficiencies, all the high sets of zeros ($N ~=~ 10 sup 11$,
$2 ~cdot~ 10 sup 11$, and $10 sup 12$) show an excess of
$delta sub n ~<=~ 0.05$ and of $delta sub n~<=~0.1$.
.P
Another way to investigate this question is to look at the
minimal values of $delta sub n$, shown in Table 5.
These values include some very small $delta sub n$ that were found
by Brent [10] and van de Lune et al. [52].
These investigators were not computing the $gamma sub n$, but
occasionally, when they had difficulty separating a close pair
of zeros, they noted the approximate locations of them.
In particular, Brent [10] found that $delta sub n ~=~ 0.001235 ^...$
for $n ~=~ 41, ~820, ~581$, and van de Lune et al. [52] found that
$delta sub n ~=~ 0.000310 ^...$ for $n ~=~ 1, ~048, ~449, ~114$.
Thus these values provide upper bounds for the minimal $delta sub n$
in the intervals investigated by those authors, and it is likely
that they are the minimal values.
.P
Since it is known [54, 55] that $p(0,t)$ has the Taylor series
expansion near $t ~=~ 0$ given by
.DS 3
.EQ (2.23)
p(0,t) ~=~ {pi sup 2} over 3 ~t sup 2 ~-~ {2 pi sup 4} over 45 ~t sup 4
~+~ {pi sup 6} over 315 ~t sup 6 ^+...^,
.EN
.DE
(for comparison, we note that $p(1,^t)~=~pi sup 6 t sup 7 /4050~+...~
[55])$,
the GUE prediction is that among $M$ consecutive $delta sub n$'s, the
probability that the minimal $delta sub n$ is $<=~ alpha ~M sup -1/3$
is approximately
.DS 3
.EQ (2.24)
1 ~-~ left ( 1 ^-^ {pi sup 2} over 9 ~{alpha sup 3} over M right ) sup M
~wig~ 1 ^-^ exp left ( -^ pi sup 2 ^ alpha sup 3 /9 right ) ~.
.EN
.DE
The function on the right side of (2.24) was used to compute
the last column in Table 5.
The probabilities of obtaining the minimal $delta sub n$ that are
observed from the GUE distribution are rather low, in some cases very low.
(We would have obtained comparable inconsistencies between data and predictions had we worked with the
pair correlation function instead of the distribution of the $delta sub n$.)
.P
It is hard to draw any firm conclusions from the data that is
available.
It would be very desirable to obtain data from much higher sets
of zeros to see whether they also have an excess of smaller
spacings compared to the GUE hypothesis.
If they do, this would raise doubt not only about the GUE
hypothesis but even about the Montgomery pair correlation conjecture.
This is an important possibility because the GUE hypothesis is very speculative,
as it relies on the hope that the eigenvalues of the hypothetical operator
associated to the zeta function would behave like those of a random
operator.
Thus it is easy to conceive of the GUE hypothesis being
false even when there is an appropriate operator.
On the other hand, the pair correlation conjecture depends only
on the assumption of randomness in the behavior of primes
(basically on error terms in the number of primes in
arithmetic progressions cancelling each other out).
Therefore falsity of the pair correlation conjecture would imply
some unusual behavior among the primes.
.P
It would also be desirable to obtain numerical data for zeros
of Dirichlet $L$-functions, using the generalization of the Riemann-Siegel
formula [20, 21] or the new method of [61].
Some data are available for zeros of Epstein zeta functions [8],
and they don't obey the GUE predictions, as they have many small
spacings between consecutive zeros.
This is not surprising, since Epstein zeta functions do not satisfy the RH
in general.
.P
We briefly discuss large $delta sub n$.
By [15],
.DS 3
.EQ (2.25)
log ~p(0,t) ~wig~ -^ {pi sup 2} over 8 ~t sup 2 ~~roman as~~ t ~->~ inf
.EN
.DE
(a result that had been derived earlier but less rigorously
by F. Dyson [22]).
Hence the GUE prediction is that
.DS 3
.EQ
max from {N ^+^1 ^<=^ n ^<=^ N ^+^M} ~delta sub n ~wig~ pi sup -1 ~( 8
^ log ^M) sup 1/2
.EN
.DE
as $N,M ~->~ inf$ with $M$ not too small compared with $N$, and this implies
.DS 3
.EQ (2.26)
max from {gamma sub n ^<^ T} ~( gamma sub n+1 ~-~ gamma sub n ) ~wig~
8 ^( 2 ^ log ^ T ) sup {-^ 1/2} ~.
.EN
.DE
Montgomery's conjecture [58; Eq. (11)] that
.DS 3
.EQ
| S(t) | ~=~ O( ( log ^ t ) sup 1/2 ~( log ^ log ^ t ) sup {-^ 1/2} )
.EN
.DE
as well as another argument of Joyner [45] suggest that
.DS 3
.EQ (2.27)
max from {gamma sub n ^<^ T} ~( gamma sub n+1 ~-~ gamma sub n ) ~=~
O( ( log ^ T ) sup {-^ 1/2} ~( log ^ log ^T ) sup {-^ 1/2} ) ~,
.EN
.DE
which contradicts (2.26).
Because of the slow growth of $( log ^ T ) sup 1/2$ and
of $( log ^ log ^ T ) sup 1/2$, our data do not allow us to conclude
anything about the truth of these conjectures.
(Some additional information about known large values of $delta sub n$
is presented in Section\ 4.)
.P
There are some measures of higher order correlation among zeros
that were computed and compared to the GUE prediction [7].
The GUE predictions are not known very accurately, but
the agreement seemed quite good
for short-range correlations.
Substantial deviations were found in the measure $sum sup 2 ^(r)$ of
[BHP] for large $r$, which equals the variance of the number
of (normalized) zeros in an interval of length $r$.
In the GUE case, $sum sup 2 ( r )$ is monotone increasing,
whereas for the zeros of the zeta function for $N~=~10 sup 12$,
this measure was increasing up to about $r~=~13$, and then started
to decrease slightly.
Given the bounds on $S(t)$, this is not very surprising.
.P
One feature of the behavior of the zero spacings which initially
seemed quite surprising involves very long-range correlations
among the $delta sub n$.
As is usual, we define the autocovariances of the $delta sub n$ by
.DS 3
.EQ (2.28)
c sub k ~=~ c sub k ^(N,M) ~=~ 1 over M ~sum from {n ^=^ N+1} to N+M
~( delta sub n ^-^ 1) ~( delta sub n+k ^-^ 1) ~.
.EN
.DE
It is not known exactly what the $c sub k$ are in the GUE, but
it has been conjectured by F. J. Dyson (unpublished) that for $k ~>~ 0$,
.DS 3
.EQ (2.29)
c sub k ~approx~ -1 over {2^ pi sup 2 ^k sup 2} ~,
.EN
.DE
where $approx$ in (2.29) indicates approximate rather
than asymptotic value.
This conjecture is intuitively appealing, since a large spacing
is expected to lead to smaller spacings nearby to preserve the
average distance, but this effect should apply less and less
the further apart the spacings one considers.
What one observes in practice, though, is quite different.
Table\ 6 shows the values of the $c sub k ^( N,M)$ for various
values of $k$, $M ~=~ 10 sup 5$, and 
for the two sets corresponding to
$N ~=~ 0$ and $N ~=~ 10 sup 12$.
Fig.\ 12 shows a graph of $c sub k ^( 10 sup 12 ^,^ 10 sup 5 )$
for $950~<=~k~<=~1000$.
Since the variance of the $delta sub n$ is around 1/6, random independent
$delta sub n$ would give values of $c sub k$ around
$(6 ~{sqrt {10 sup 5}} ) sup -1 ~approx~ 0.00053$.
The values in Table 5 are typically much larger than that, and
therefore statistically significant.
These values say, for example, that a large $delta sub n$ tends
to be associated with a small $delta sub n+963$ (for $N ~=~ 10 sup 12$).
Even more striking than the sizes of the $c sub k$ are the patterns
of signs that are visible in Table 6 and Fig. 12.
(These effects are much more striking for $N ~=~ 10 sup 12$ than
for $N ~=~ 0$ due to the ``averaging out'' of the spectrum of the
$delta sub n$ for $N ~=~ 0$ that is introduced by the normalization
(2.7), but we won't discuss this here.)
.P
The long-range correlations between the $delta sub n$, which will
be discussed below, should not obscure the fact that our data do
support the Dyson conjecture (2.29), at least in the weak form that
says that for every fixed $k$, as $N,M ~->~ inf$ with $M$ not too
small compared to $N$, $c sub k$ is approximated by the quantity
on the right side of (2.29).
Some support for this conjecture can be derived from Table 6.
The fact that (2.29) does not hold for $k$ large is to be expected
given the known results about $S(t)$.
The only thing that is surprising is that the $c sub k$ are
large for large $k$ and there are patterns to their behavior.
.P
An explanation for the long-range correlations between the $delta sub n$
can be obtained by looking at the spectrum of the $delta sub n$.
Fig. 13 shows a portion of the spectrum, i.e., a graph of
.DS 3
.EQ
f(x) ~=~ c ^ | ^ sum from {n ^=^ N+1} to N+M ~( delta sub n ^-^ 1)
~e sup {pi ^in^x} | sup 2
.EN
.DE
for a certain constant $c$, where $N ~=~ 10 sup 12$,
$M ~=~ 98303$.
No smoothing or tapering was done to the data.
The spectral lines are very sharp.
The line graph of Fig. 13 is based on approximately 12000 evenly
spaced values of $x$ between 0 and 1/4, so that the heights
of the peaks may not be very accurately portrayed, but their
locations are accurate.
(The function $f(x)$ is periodic with period 1.
For $0.25 ~<=~ x ~<=~ 1$, it has more sharp peaks, but
they get closer and closer together and eventually begin to
coalesce into a chaotic regime.)
To show more clearly the behavior of the spectrum of the $delta sub n$,
Fig. 14 shows a graph of $g(x) ~=~ max ( log ^f(x) , ~-~ 10)$.
(This definition is meant to deemphasize regions where $f(x)$ is small.)
.P
The peaks of $f(x)$ occur near $x ~=~ 2 ~( log ^N ) sup -1 ~ log ^q$,
where $q$ is a prime power, and the spikes in Fig. 13 are due to
2,3,4,5,7,8,9,11,13,16,17, and 19, in that order, with proper prime
powers having smaller peaks than primes.
The explanation for this behavior of the spectrum of the $delta sub n$
lies in the formulas that connect zeros of the zeta function
and prime numbers.
Some of them, such as the exact formula for $pi (x)$ in terms
of the zeros, or the more general ``explicit formulas'' of
Guinand [38] and Weil [72] are well known.
A related formula was proved by Landau [48].
(See [37] for a somewhat stronger version.)
It says, under the assumption of the RH (which we assume for
convenience, although Landau proved an unconditional result),
that for any $y ~>~ 0$ as $N~->~inf$,
.DS 3
.EQ
sum from n=1 to N ~e sup {i ^gamma sub n ^y} ~=~ left {
matrix {
lcol {-^ {gamma sub N} over {2 pi} ~exp (-y/2) ^ log ^ p ~+~ O( exp^(-y/2) log ^ N) ~
~~~roman if~ y ~=~ log ^ p sup m ^, above O( exp ( -y/2 ) ^ log ^N ) ~~~roman
if~ y ~!=~ log ^ p sup m ^,}
}
.EN
.DE
where $p$ denotes a prime and $m$ a positive integer.
Therefore we also expect that if
.DS 3
.EQ
h(y) ~=~ sum from {n ^=^ N+1} to N+M ~e sup {i ^ gamma sub n ^y} ^,
.EN
.DE
then $h(y)$ will be large precisely for $y$ in the vicinity of $log ^ p sup m$
and small elsewhere.
Fig. 15 shows a graph of $log ^ | h(y) |$ for $0.004 ~<=~ y ~<=~ 3$,
$N ~=~ 10 sup 12 , ~M ~=~ 4 ^cdot^ 10 sup 4$, which exhibits precisely
this behavior.
(This function $log^|h(y)|$ is not graphed for
$0~<=~y~<~0.004$, since it is very large there.)
Let us write
.DS 3
.EQ
gamma sub N+k ~=~ gamma sub N ~+~ k alpha ~+~ beta sub k ^,
.EN
.DE
where $alpha$ is the average spacing,
.DS 3
.EQ
alpha ~=~ 2 pi ~left ( log ^ {gamma sub N} over {2 pi} right ) sup -1 ~.
.EN
.DE
The $beta sub k$ are usually small, on the order of $1/ alpha$.
Then, for $y$ small, we can expect that
.DS 3
.EQ
h(y) ~mark =~ e sup {i ^gamma sub N ^y} ~sum from k=1 to M ~e sup
{i k alpha y ^+^ i beta sub k y}
.EN
.sp
.EQ (2.30)
~lineup approx~ e sup {i^ gamma sub N ^y} ~sum from k=1 to M ~e sup
{i k alpha y} ~+~ i ^y^ e sup {i^ gamma sub N ^y} ~sum from k=1 to M
~beta sub k ~e sup {ik alpha y}
~.
.EN
.DE
The first sum on the right side of (2.30) is small for $y$
away from integer multiples of $2 pi / alpha$, and so it is
the second sum that contributes the oscillations at
$y ~=~ log ^ p sup m$.
(For $y$ away from $log^p sup m$, and especially so for
$y$ small, the first sum on the right side of (2.30)
dominates, and since it equals the rapidly oscillating function
$|^sin^(M~alpha^y/2)|$
$| sin^( alpha y/2)| sup -1$
is absolute magnitude, it gives rise to the very
regular patterns of points visible in Fig.\ 15
because of the comparatively infrequent but regular spacing of sample
points.)
On the other hand,
.DS 3
.EQ
sum from {n^=^N+1} to N+M ~( delta sub n ^-^ 1) ~e sup {pi i nx} ~mark approx~
alpha sup -1 ~sum from {n^=^N+1} to N+M ~( gamma sub n+1 ^-^ gamma sub n
^-^ alpha ) ~e sup {pi i n x}
.EN
.sp
.EQ
~lineup =~ alpha sup -1 ~e sup {pi i N x} ~sum from k=1 to M ~( beta
sub k+1 ^-^ beta sub k ) ~e sup {pi i kx}
.EN
.sp
.EQ
~lineup =~ alpha sup -1 ~e sup {pi i N x} ~left { beta sub M+1 ~e sup
{pi i Mx} ^-^ 2 beta sub 1 ~e sup {pi i x} right }
.EN
.sp
.EQ
~lineup +~ alpha sup -1 ~e sup {pi i Nx} ~sum from k=1 to M ~beta sub k
~e sup {pi i kx} ^,
.EN
.DE
and so the spectrum of the $delta sub n$ 
can be expected to behave like the second sum on the right side
of (2.30) and to
have peaks
at frequencies $x ~=~ pi sup -1 ^alpha ~ log ~ p sup m$ for
primes $p$ and $m ~>=~ 1$.
.H 1 "Computation of the zeros"
This section describes the computations of the zeros that
were used in the comparisons with eigenvalues of random
hermitian matrices described in Section 2 and documents
the claim that the values that were found are accurate
to within $+-^10 sup -8$.
For simplicity we will restrict the discussion here to the
computations of $gamma sub n$ for $n ~=~ 10 sup 12 ~+~1$
to $n ~=~ 10 sup 12 ~+~ 101000$.
Exactly the same method was used to compute $gamma sub n$ for
$2001 ~<=~ n ~<=~ 101000$, but the error estimates in that
case are much easier.
(Values of $gamma sub n$ for $1 ~<=~ n ~<=~ 2000$ were taken
from the 105 decimal place values computed in [60].)
.P
The computations described here relied on the Riemann-Siegel
formula [23, 43, 68], which requires on the order of $t sup 1/2$
operation to compute $zeta ( 1/2 ~+~ it)$.
A new method has recently been invented [61] which enables
one to compute on the order of $T sup 1/2$ values of
$zeta ( 1/2 ~+~ it)$ for $T ~<=~ t ~<=~ T ~+~ T sup 1/2$ in a total
of $O( T sup {1/2 ^+^ epsilon} )$ operations for all $epsilon ~>~ 0$
(i.e., $O( T sup epsilon )$ per value), but this method is quite
involved and it has yet to be implemented.
.P
We first recall the Riemann-Siegel formula [10, 23, 28, 43, 59, 68].
Let
.DS 3
.EQ
theta (t) ~=~ roman arg [ pi sup {-^ it/2} ~GAMMA ( 1/4 ~+~ it/2 ) ] ^,
.EN
.DE
and
.DS 3
.EQ
Z(t) ~=~ exp ( i ^ theta ^(t)) ~zeta ( 1/2 ~+~ it ) ~.
.EN
.DE
Then $Z(t)$ is real for real $t$, and any sign change of
$Z(t)$ corresponds to a zero of $zeta (s)$ on the critical
line.
We define
.DS 3
.EQ (3.1)
tau ~=~ ( 2 pi ) sup -1 ^t, ~~m ~=~ left floor tau sup 1/2 right floor , ~~
z ~=~ 2( tau sup 1/2 ~-~ m) -1 ~.
.EN
.DE
(All the computations of zeros $gamma sub n$ for $n$ near
$10 sup 12$ that were carried out had $m ~=~ 206393$).
Then, for any $k ~>=~ 0$,
.DS 3
.EQ
Z(t) ~mark =~ 2 ^sum from n=1 to m ~n sup {-^ 1/2} ~ cos ( t ^ log ^(n)
~-~ theta (t))
.EN
.sp
.EQ (3.2)
~~~~~
.EN
.sp
.EQ
~lineup +~ ( -1 ) sup m+1 ~tau sup {-^ 1/4} ~ sum from j=0 to k ~
PHI sub j ^(z) ^( -1 ) sup j ~tau sup {-^ j/2} ~+~ R sub k ^( tau ) ^,
.EN
.DE
where
.DS 3
.EQ (3.3)
R sub k ^( tau ) ~=~ O( tau sup {-^( 2k+3)/4} ) ^,
.EN
.DE
and the $PHI sub j ^(z)$ are certain entire functions.
Gabcke [28] has derived very sharp explicit error estimates
for the remainder terms $R sub k ^( tau )$, and the one used
in the computations was
.DS 3
.EQ (3.4)
| R sub 4 ~( tau ) | ~<=~ 0.017 t sup -11/4 ~~roman for~~ t ~>=~ 200 ~.
.EN
.DE
The $PHI sub j ^(z)$ were computed using their Taylor series
expansion [17, 28].
.P
Before discussing the specific implementation of the Riemann-Siegel
formula that was used, we present the basic assumptions used
in the error analysis.
We will assume that the computer that was used (a Cray X-MP) performed as
usual, i.e., that there were no intermittent faults, and the same
program run on it at another time would produce the same results.
We also have to make some assumptions about the correct functioning
of the basic arithmetic units and the standard logarithm and cosine
routines supplied by the manufacturer.
It is often unsafe to trust the manufacturer's specifications
and several examples of this will be presented below.
.P
The model of computation we will use in the error analysis of
our algorithm is that developed by W. S. Brown [12].
It assumes that floating point numbers can be represented in
the signed, base $b ^, ~t$ digit form
.DS 3
.EQ (3.5)
+- ~b sup e ~sum from i=1 to t ~a sub i ~b sup -i ^,
.EN
.DE
where either $a sub 1 ^=...=^ a sub t ~=~ 0$ or $1 ~<=~ a sub 1 ~<~ b$,
$0 ~<=~ a sub i ~<~ b$ for $2 ~<=~ i ~<=~ t$, and $e sub min ~<=~ e ~<=~ e sub max$.
Brown [12] defines ``correct'' arithmetic for floating point
numbers in this model.
If a model number is one that is representable in the form
(3.5), then Brown's model can be briefly summarized by the following
rules:
.sp
.in +5n
If the result of applying an elementary operation $(+, ~-, ~times~ roman or~ /)$
to two model numbers is exactly a model number, then the implementation
must return that number.
.sp
Otherwise, the implementation may return anything between the
model numbers adjacent to the exact result.
.sp
Comparison of model numbers must be exact.
.sp
Comparison of non-model numbers $x$ and $y$ may report the result
of exact comparison between any pair of elements in the smallest
model intervals containing $x$ and $y$.
.in -5n
.sp
(For more details, see [12].)
.P
One advantage of Brown's model of computation is that it is
not necessary to worry about the internal details of the floating
point arithmetic unit, but instead it is sufficient to
test whether a particular machine obeys the specifications
of the model.
It is of course impractical to test all possible cases.
N. L. Schryer [64, 65] has implemented a very extensive
set of tests designed to reveal deviations from the rules of the
model and has used it to test many computers.
Schryer's tests are based on the fact that there tend to be
patterns in the way hardware designers make mistakes, and these
patterns can be detected by cleverly designed test arguments.
Schryer's tests have detected most of the previously known defects as
well as some new ones on a variety of machines.
The superficial description of the Cray-1 and the Cray X-MP (both of
which have 64-bit words) might lead one to expect that they
would satisfy Brown's model with $b ~=~ 2, ~t ~=~ 48$ for single
precision and $b ~=~ 2, ~t ~=~ 96$ for double precision.
However, Schryer's tests have revealed some faults.
For example, if $fl(x)$ denotes
the floating-point version of $x$, then $2 ~times~ fl ( 2 sup -1 ~+~ 2 sup -2 ^+...+^ 2 sup -48 ) ~=~ 2$
on both machines.
In fact, D.\ Winter has pointed out that it already follows from
the Cray hardware manual [19] that these machines do not
satisfy Brown's model for $b~=~2$, $t~=~48$.
On the other hand, very extensive tests by Schryer's programs have shown that these
machines appear to satisfy Brown's model for $b ~=~ 2, ~t ~=~ 47$
in single precision and $t ~=~ 94$ in double precision.
We will assume in our analysis that this is correct.
(Some defects in Cray software will be described at the end of
this section.)
.P
Using Brown's model, it is possible to bound the errors that are made
in various elementary arithmetic operations.
For example, if $x$ and $y$ are both model numbers, $2 sup k ~<~ x ~<~ 2 sup k+1$,
$0 ~<~ y ~<~ 2 sup k$, and $x +y ~<~ 2 sup k+1$, where
$1 ~<=~ k ~<~ 40$, say, then the single precision value $z$ that is computed
for $x+y$ will satisfy
.DS 3
.EQ
| z ~-~ (x+y) | ~<~ 2 sup k-46 ~.
.EN
.DE
.P
We now describe the computation of the sum of cosines
in (3.2), which is what consumes almost all of the computing time.
The argument $t$ is always maintained as a double precision
$(d p )$ variable.
Another $d p$ variable, $t sub 0$, is also maintained, which has
the property that $| t ~-~ t sub 0 | ~<=~ 3$.
Three arrays $d sub n , ~q sub n , ~u sub n$, $1 ~<=~ n ~<=~ m$,
are also used; $d sub n$ is the value of $log ^n$, computed
and stored in $d p$, $q sub n$ is the value of $2 ^n sup {-^ 1/2}$,
computed in $d p$ but stored in single precision $( s p )$
and $u sub n$ is the value of $t sub 0 ^ log ^n$ reduced modulo
$2 pi$, computed in $d p$ but again stored in $s p$.
(On the Cray-1 and Cray X-MP, conversion from $d p$ to $s p$ is
effected by truncation, and from $s p$ to $d p$ by padding with
zeros.)
To compute $Z(t)$ for a new value of $t$, a comparison of $t$
with $t sub 0$ is made.
If $| t ~-~ t sub 0 | ~>~ 3$, $t sub 0$ is set equal to $t$,
and the $u sub n$ are recomputed.
At that point we are in the remaining case of $| t ~-~ t sub 0 | ~<=~ 3$.
Here $delta$ is defined to be the $s p$ value of $t ~-~ t sub 0$,
and $t sub 1$ the $d p$ value of $t sub 0 ~+~ delta$ (so
that $| t sub 1 ~-~ t | ~<~ t sub 0 ~cdot~ 2 sup -90$, say).
Next, $theta ( t sub 1 )$ is computed in $d p$, reduced
modulo $2 pi$, and converted to a $s p$ variable $v$.
Then
the quantities
.DS 3
.EQ (3.6)
w sub n ~=~ q sub n ~ cos ( delta ^e sub n ~+~ ( u sub n ~-~ v)) ~~,~~
1 ~<=~ n ~<=~ m
.EN
.DE
where $e sub n$ is the $s p$ value of $d sub n$, are computed
in $s p$, using the manufacturer-supplied $s p$ cosine routine, and these values
are added up in a special way to be described below.
.P
The reason for the involved procedure described above is that for
$t ~>~ 10 sup 10$, say, $s p$ arithmetic is not sufficiently
accurate.
On the other hand, $d p$ arithmetic is much too slow
to be used all the time,
since on the Cray \%X-MP it is done in software and does not vectorize.
Since the costly step of recomputing the $u sub n$ is done very
infrequently (roughly once for every 100 evaluations of $Z(t)$,
on average, in the main, zero-locating run), the scheme above allows for relatively high accuracy at
reasonably high speed.
.P
The computation that is described above approximates $Z(t sub 1 )$
rather than $Z(t)$.
However, since $| t sub 1 ~-~ t | ~<~ t sub 0 ~cdot~ 2 sup -90$,
this difference is insignificant in practice,
since zeros were computed only to within $+-^ 5 ~cdot~ 10 sup -9$.
.P
We next estimate how much $w sub n$ differs from $2^ n sup {-^1/2} ~ cos ( t sub 1 ^ log ^n ~-~ theta ( t sub 1 ))$.
The function $theta ( t sub 1 )$ is computed entirely in $d p$ using
Stirling's formula [40].
The main term in that formula is $t sub 1 ^( log ^t sub 1 ~-~ 1 ~-~ log (2 pi ) )/2$.
No performance figures are supplied by the manufacturer for the
$d p$ logarithm routine.
However, tests performed with 1000 random inputs $epsilon ^( 0, ~10 sup 13 )$
(comparing the values obtained on the Cray X-MP with those
obtained with the Macsyma symbolic manipulation system [53],
which allows for variable precision) found a maximal error of
$<~ 2 ^cdot^ 10 sup -27$.
Therefore it was assumed that all $d p$ logarithm values
(for $2 ^cdot^ 10 sup 11 ~<~ t sub 1 ~<~ 3 ^cdot^ 10 sup 11$)
are correct to within $+- ~5 ^cdot^ 10 sup -27$.
Under this assumption, the $d p$ values of
$t sub 1 ^( log ^t sub 1 ~-~ 1 ~-~ log (2 pi ) )/2$
are accurate to within $+-~ 4 ^cdot^ 10 sup -27 ^t sub 1$.
The computation of the other terms in the expansion of $theta ( t sub 1 )$
is easily shown to introduce an error of at most $10 sup -17$.
Hence
after reduction modulo a $d p$ value of $2 pi$ (correct to
within $+-~ 10 sup -27$), the resulting $d p$ value of $theta ( t sub 1 )$ is in the
interval $(-0.01,~ 6.3)$, say, and (except for a possible additive factor
of $+- 2 pi$) is correct to within $+-~ 7 ^cdot^ 10 sup -27 ^t sub 1$.
The conversion to a $s p$ value introduces an error in
at most the last bit, and this error is therefore $<=~ 2 sup -44$.
Hence the $s p$ value $v$ of $theta ( t sub 1 )$ differs
from $theta ( t sub 1 )$ by at most
$2 sup -44 ~+~ 7 ^cdot^ 10 sup -27 ^t sub 1$ (and,
possibly, by an additional factor of exactly $2 pi$).
Similarly, the $s p$ value $u sub n$ differs from $t sub 0 ^ log ^ n$
by an integer multiple of $2 pi$ and at most
$2 sup -44 ~+~ 7 ^cdot^ 10 sup -27 ^t sub 1$
and again is in the range (-0.01, 6.3).
.P
Since $e sub n$ is obtained from $d sub n$ by truncation of the
bit pattern
it differs from $log ^n$ by at most
$( 2 sup -47 ~+~ 10 sup -24 )~ log ^n$, say.
Furthermore, the $s p$ value computed by multiplying $delta$ and
$e sub n$ in $s p$ can differ from $delta ^e sub n$ by at most
$( 2 sup -47 ~+~ 2 sup -93 ) ~delta ^e sub n$.
Hence the $s p$ value that is computed for $delta ^e sub n$ differs
from $delta ~ log ~n$ by at most $( 2 sup -45 ~+~ 10 sup -22 ) ~ log ~n$,
as $| delta | ~<=~ 3$.
Since $m ~<=~ 210000$, we have $| delta ^e sub n | ~<=~ 37$.
Subtraction of $v$ from $u sub n$ gives a $s p$ value that
is in the range (-7,7) but off by at most $2 sup -44$ from the
correct value, and adding it to the $s p$ value of $delta ^e sub n$
yields a value that is in the interval (-45,45), but is up to
$2 sup -41 ~+~ 2 sup -44$ away from the correct sum of the computed
values of $delta ^e sub n$, $u sub n$, and $-v$.
Therefore the $s p$ value of $r sub n ~=~ delta e sub n ~+~ ( u sub n ^-^ v )$ differs from
$t sub 1 ^ log ^ n ~-~ theta ( t sub 1 )$ by an integer multiple
of $2 pi$ and an additive factor that is
.DS 3
.EQ (3.7)
<=~ 2 sup -41 ~+~ 2 sup -44 ~+~ ( 2 sup -45 ^+^ 10 sup -22 ) ^ log ^n ~+~
2 sup -43 ~+~ 1.4 ^cdot^ 10 sup -26 ^t sub 1 ~<=~ 10 sup -12
.EN
.DE
in absolute magnitude, for $t sub 1 ~<=~ 2.7 ^cdot^ 10 sup 11$
and $n ~<=~ 210000$.
Hence
.DS 3
.EQ (3.8)
| cos ( r sub n ) ~-~ cos ( t sub 1 ^ log ^ n ~-~ theta ( t sub 1 )) |
~<=~ 10 sup -12
~.
.EN
.DE
.P
No rigorous bounds are known for the error made by the
Cray X-MP $s p$ cosine routine.
Some estimates are supplied [18; Appendix B] based on a random sample
of $10 sup 4$ arguments; for example, for $x ~member~ ( -^ pi ^,^ pi )$, the
maximal difference observed between $cos (x)$ and the computed
value is given as $1.01 ^cdot^ 10 sup -14$, and the standard
deviation as $2.31 ~cdot~ ~10 sup -15$.
However, these figures appear to come from an old edition of the
manual, whereas substantial changes have been made to the software
in the meantime.
A sample of $10 sup 4$ points $x$, drawn uniformly at random from
the interval $(-^ pi ^,^ pi )$ showed that the maximal difference
between values of $cos (x)$ computed in $d p$ and $s p$
was only $5.34 ^cdot^ 10 sup -15$, and the standard deviation
$1.4 ^cdot^ 10 sup -15$, which is considerably better than
claimed in [18].
A similar sample of $10 sup 6$ points $x$, drawn uniformly
at random from (-50,50) showed that the maximal difference
was $<=~ 1.35 ^cdot^ 10 sup -14$ and the standard deviation
$<=~ 1.65 ^cdot^ 10 sup -15$.
A comparison of several hundred values computed with the
$s p$ cosine routine with the same values computed
in the Macsyma [53] system showed similar results.
Therefore it seemed safe to assume that for $x ~member~ (-50, 50)$,
$s p$ values of $cos (x)$ are accurate to within
$+-~ 5 ^cdot^ 10 sup -14$.
With this assumption, the computed value of $cos ( r sub n )$
differs from the correct value $cos ( t sub 1 ^ log ^ n ~-~ theta ( t sub 1 ))$
by at most $1.05 ^cdot^ 10 sup -12$.
.P
Since $q sub n$ is obtained by converting the $d p$ value of
$2^ n sup {-^ 1/2}$ to $s p$, it is correct to within
$+-~ ( 2 sup -47 ~+~ 2 sup -90 ) ~2 ^ n sup {-^ 1/2}$, say.
Therefore the $s p$ value of $q sub n ~ cos ( r sub n )$ that
is computed differs from the desired value of
$2^ n sup {-^ 1/2} ~ cos ( t sub 1 ^ log ^ n ~-~ theta ( t sub 1 ))$
by at most
.DS 3
.EQ (3.9)
2^ n sup {-^ 1/2} ~( 1.05 ^cdot^ 10 sup -12 ~+~ 2 sup -46 ~+~ 2 sup -80 )
~<=~ 2.14 ^cdot^ 10 sup -12 ~n sup {-^1/2} ~.
.EN
.DE
The sum of the term on the right sides of (3.9) for
$1 ~<=~ n ~<=~ m ~<~ 207000$ is $<~ 1.95 ^cdot^ 10 sup -9$.
.P
To add up the computed $s p$ values of $q sub n ~ cos ( r sub n )$,
the values for $1 ~<=~ n ~<=~ 1280$ were converted to $d p$
and added in $d p$.
The $n$ with $1280 ~<~ n ~<=~ m$ were partitioned into sets $B$
such that for each $B$, $B$ had $<=~ 120$ elements and
.DS 3
.EQ
sum from {n ^epsilon^ B} ~2 ^ n sup {-^1/2} ~<~ 1/2 ~.
.EN
.DE
For each $B$, the sum $S sub B$ of the computed values of
$q sub n ~ cos ( r sub n )$ for $n ~epsilon~ B$ was
evaluated in $s p$, making
an error of at most $2 sup -41$, converted to $d p$, and
added in $d p$ to the sums of the other sets $B$ and of
the initial 1280 values.
There were $<~ 3600$ sets $B$, so the total error made in this
addition is
.DS 3
.EQ
<~ 3600 ~cdot~ 2 sup -41 ~+~ 2 sup -60 ~<~ 1.7 ^cdot^ 10 sup -9
~.
.EN
.DE
.P
Gabcke's bound (3.4) for the remainder term in the Riemann-Siegel
formula and the error terms given in [17, 28] for the Taylor
series expansions of the $PHI sub j ^(z)$ guarantee that these
sources contribute an additional error of at most $5 ~cdot~ 10 sup -11$.
We conclude therefore that the computed values of $Z( t sub 1 )$
are correct to within $+-^ 3.7 ~cdot~ 10 sup -9$.
.P
The actual computation of zeros was done in two stages.
In the first stage, the cumbersome procedure described
above for computing the sum of the $w sub n$ was replaced
by a simpler one which gave approximately a 10% speedup
of the entire routine for evaluating $Z(t)$.
This modified routine was used to locate the zeros to
a nominal accuracy (i.e., disregarding any errors in computation
or from neglecting remainder terms in the Riemann-Siegel formula)
of $+-~ 5 ^cdot^ 10 sup -9$.
The procedure was the standard one [10, 52] of locating Gram
blocks and searching for the expected number of sign changes
of $Z(t)$ in them.
Once these sign changes were found, zeros were computed more
accurately using the Brent algorithm [9] for locating zeros
of functions by a combination of linear and quadratic interpolation.
On average, eight evaluations of $Z(t)$ were used for each zero.
.P
In the second stage the more accurate algorithm described above
was used.
For each value $gamma$ of a zero that was computed in the first
stage, $Z(t)$ was evaluated at the two points $t ~=~ gamma ~+-~ 8 ^cdot^ 10 sup -9$,
the signs of the computed values were verified to be opposite,
and their absolute magnitudes were checked to be greater than
the upper bound $3.7 ^cdot^ 10 sup -9$ for the error that we
obtained above.
In a few cases this check failed.
For these zeros a more accurate version of the program, operating
almost entirely in double precision, was used to verify the
correctness of their locations.
.P
The time for the first stage of the computation of the zeros
was approximately 16 hours on a single-processor Cray X-MP.
The same program runs only about half as fast on the Cray-1.
Since the cycle time on the X-MP is only about 30% faster
than on the Cray-1, most of the gain seems to be due to the
larger number of memory ports on the X-MP.
.P
The second stage, the careful verification that the computed
values were as accurate as claimed, was carried out twice,
once using the CFT 1.14 Fortran compiler, and once using an early
unofficial version of the CFT 1.15 computer.
Each run took about 5 hours.
(These runs were slower than the initial zero-locating
ones not only because of the less efficient method for
adding up the $w sub n$ that was used, but also because
the double precision array recomputations were relatively
more frequent.)
The reason for undertaking two runs was that both compilers
have been reported by other users to have
bugs which have occasionally led to incorrect answers when operating
with long vectors.
The exact circumstances that lead to such wrong computations are
for the most part not known very well (most seemed to be
associated with long complex vectors or various optimization features
that were not used in our program), and the bugs affecting the
two compilers seemed different.
The results of the two computations were compared, and they seemed
identical to the full $s p$ accuracy.
.P
One very annoying bug in the Cray-1 and Cray X-MP software was
discovered during the computations.
When $D$ is declared to be a $d p$ variable, and one uses
the free-format statement
.DS 3
.EQ
READ ~*, ~D
.EN
.DE
to read a number, this number is first converted into a $s p$
value, and only then converted to $d p$, so that $D$ agrees
with the value being read in only about 48 bits.
The values that had been computed earlier for zeros $gamma sub n$,
$10 sup 9 ~<=~ n ~<=~ 10 sup 9 ~+~ 10 sup 5$,
$10 sup 11 ~<=~ n ~<=~ 10 sup 11 ~+~ 10 sup 5$, and
$2 ^cdot^ 10 sup 11 ~<=~ n ~<=~ 2 ^cdot^ 10 sup 11 ~+~ 10 sup 5$
had been converted for archival storage using a program that used
the above free-format input.
As a result, the values that are available for those zeros are
not very accurate and were not used for most of the analyses in Section 2.
.H 1 "Statistics of zeros"
The main computation discovered 101111 zeros between Gram points
$g sub n$ with $n ~=~ 10 sup 12 ~-~ 14$ and
$n ~=~ 10 sup 12 ~+~ 101097$.
Since the interval $[g sub n , ~g sub n+12 )$ for $n ~=~ 10 sup 12 ~-~ 14$
is the union of 6 Gram blocks of length 1 and 3 Gram blocks of
length 2,
and the interval $[g sub n , ^g sub n+12 )$ for $n ~=~ 10 sup 12 ^+^ 10 sup 5 ^+^ 1085$
is the union of 8 Gram blocks of length 1 and 2 Gram blocks of length 2,
we can conclude by the now standard Turing method
(Theorem 3.2 of [10], for example) that
the 101087 zeros that were found in the interval
$[g sub p ~,~ g sub q ), ~p ~=~ 10 sup 12 ^-^ 2$,
$q ~=~ 10 sup 12 ~+~ 10 sup 5 ~+~ 1085$ are all the
zeros of the zeta function in that range, and are indeed
the zeros $gamma sub n$ with $10 sup 12  ~<=~ n ~<=~ 10 sup 12 ~+~ 10 sup 5 ~+~ 1086$.
.P
The interval $[g sub p ^,^ g sub q )$ for $p ~=~ 10 sup 12 ^-^ 2$,
$q ~=~ 10 sup 12 ^+^ 10 sup 5 ^+^ 1$ appears to consist of 65265
Gram blocks of length 1, 10686 of length 2, 2850 of length 3,
847 of length 4, 217 of length 5, 49 of length 6, and 7 of length 7.
(We say ``appears'' because the values of $Z(t)$ at Gram points
were not checked to verify that their signs were unambiguous, except
near violations of Rosser's rule.)
.P
Four exceptions to Rosser's rule were found during the computation
of the zeros $gamma sub n$, $10 sup 12 ^+^ 1 ~<=~ n ~<=~ 10 sup 12 ^+^ 10 sup 5$.
In the now-standard terminology of [52], all were
of length 2, and two were of type 1
($n ~=~ 10 sup 12 ^+^ 15934$ and $n ~=~ 10 sup 12 ^+^ 93436$) and
two were of type 2 ($n ~=~ 10 sup 12 ^+^ 13452$ and
$n ~=~ 10 sup 12 ^+^ 17086$).
Earlier computations had found 2 exceptions to Rosser's rule
among zeros $gamma sub n$ with
$2 ^cdot^ 10 sup 11 ~+~ 1 ~<=~ n ~<=~ 2 ^cdot^ 10 sup 11 ~+~ 10 sup 5$
(both of length 2, type 2, for $n ~=~ 2 ^cdot^ 10 sup 11 ~+~ 5469$ and
$n ~=~ 2 ^cdot^ 10 sup 11 ~+~ 52084$), one exception among the
$gamma sub n$ with $10 sup 11 ~+~ 1 ~<=~ n ~<=~ 10 sup 11 ~+~ 10 sup 5$
(for $n ~=~ 10 sup 11 ~+~ 10477$ and of length 2 and type 2), and no exceptions
for $10 sup 9 ~+~ 1 ~<=~ n ~<=~ 10 sup 9 ~+~ 10 sup 5$.
Thus the frequency with which Rosser's rule is violated appears
to increase very rapidly with increasing height (cf. [10, 52]).
.P
The largest $delta sub n$ that this author is aware of is
$delta sub n ~=~ 4.2626^...$ for $n ~=~ 184, ~155, ~671, ~040$.
(We have $gamma sub n ~wig~ 5.3 ^cdot^ 10 sup 10$ here.)
$Z(t)$ atains a value of approximately \-214 near
$gamma sub n$, and there is a violation of Rosser's rule of length 2, type 2 in the
vincinity of $gamma sub n$.
This zero was discovered by an application of a method for
constructing values of $t$ for which $| Z(t) |$ is large
that was based on the Lovasz lattice basis reduction algorithm
in a manner somewhat similar to that of [60].
Several other values of $t~<wig~10 sup 11$ with $|Z(t)|$ large
were constructed that way, and a substantial fraction of them
corresponded to violations of Rosser's rule,
although all the violations discovered that way were of known type.
It is clear that this method can be used to construct $t$ with
extremely large values of $|Z(t)|$, and the main barrier to its
use is the inability to compute the corresponding values of
$Z(t)$ efficiently.
(It is hoped that the algorithm of [61] will be implemented in the
near future, which would enable one to explore the behavior
of $Z(t)$ near such special points where it is large.)
.P
Other methods of constructing values of $t$ with interesting
behavior of $Z(t)$ are presented in [46] and [50].
Van de Lune [50] has found several values of $t$ for which $| Z(t) |$
is very large.
Thanks to 
D.\ Hejhal, the National Science Foundation, and the Minnesota Supercomputer Institute,
this author was able to use a \%Cray-2 to study some of these points.
The huge memory of the \%Cray-2 made it possible to utilize
the algorithm described in Section\ 3 to compute $Z(t)$ for the most
interesting of the values of $t$ found by van\ de\ Lune.
(The highest point mentioned below required about 45\ million words of memory.)
Rigorous error analysis was not done, and the
computations spanned blocks of 20\ Gram intervals centered approximately
at the points found by van\ de\ Lune, so it is impossible to assert
that all of the zeros in those intervals have been found and that they
satisfy the RH, but it seems very likely that they do.
Table\ 7 presents the results of these computations.
For each point $t$ that was checked, the approximate value of $Z(t)$
is given (taken from [50], since the computations reported have used a different
grid of points at which to evaluate $Z(t)$), the type of violation of
Rosser's rule that was found (if any), and the maximal $delta sub n$
that was found in that neighborhood.
The value of $Z(t)~approx~-453.9$ at the last point in the list gives the
largest values of $|Z(t)|$ that has been found so far.
An entry of ``no'' in the ``Rosser rule'' column means that no violation
of Rosser's rule was found.
The entry (2,2) means that a violation of length 2 and type 2
[52] was found.
The entry ``a'' means a previously undetected type of violation
of length\ 3, with the two Gram intervals preceding the block with zero
pattern (0,\ 1,\ 0) having zero pattern (2,2).
The entry ``b'' denotes another new type of violation,
this time of length 2, with Gram block with zero patterns (0,0) being
followed by a block with zero pattern (2,\ 1,\ 1,\ 2).
.H 1 "Other methods for computing the zeta function"
The program described in the preceding section was designed for
relatively accurate and rapid computation of the zeros of the
zeta function at fairly large weight.
On the single-processor Cray X-MP with 2 million words of memory
on which this program was run, it could only be used for zeros
$gamma sub n$ with $n$ not much larger than $3 ~cdot~ 10 sup 12$ due to
the limitations on memory storage.
In terms of speed, on the Cray X-MP
it ran at about two thirds of the speed of the latest version [73] of the Cyber-205 program (written
partially in assembly language) that was used for the verification
of the RH for the first $1.5 ^times^ 10 sup 9$ zeros [52],
since on average it required about $3.5 ^cdot^ 10 sup -3$ seconds
to evaluate a sum of the form
.DS 3
.EQ (5.1)
sum from {n ~=~ K+1} to {K ~+~ 10 sup 4} ~2 ^n sup {-^1/2} ~ cos
( t ^ log ^(n) ~-~ theta (t)) ^,
.EN
.DE
whereas the data in [73] indicate that a comparable computation
with the program outlined there requires about $2.5 ~cdot~ 10 sup -3$ seconds.
The program described in Section 4 has the advantages of being
more portable (since it is written exclusively in Fortran) and
more accurate.
The error estimates in it are
not completely rigorous, as the error bounds assumed for the
manufacturer's cosine routine are based on statistical sampling.
The main program described in [51, 52, 73] used the manufacturer's
$d p$ cosine routines only to compute a few values to medium accuracy.
These values were then used inside the program to obtain $s p$ values
of the cosine by linear interpolation.
Thus the [51, 52, 73] program can be regarded as somewhat more
trustworthy than the one described here.
(In some cases, though, where $Z(t)$ was small, a program in which
all operations were carried out using the manufacturer's $d p$ routines
was used in [51, 52], and its trustworthiness is comparable to ours.)
.P
It is possible to write programs that are much faster than that described
in Section 3, are almost as accurate, are written in Fortran, and whose
errors can be analyzed rigorously.
Note that most of the effort in computing $Z(t)$ by the Riemann-Siegel
formula (3.2) is spent in evaluating
.DS 3
.EQ
sum from n=1 to m ~2 ^ n sup {-^ 1/2} ~ cos ( t ^ log ^n ~-~ theta (t))
~=~ Re ~e sup {-^ i theta (t)} ~sum from n=1 to m ~2^ n sup {-^1/2}
~e sup {it ^ log ^n} ~.
.EN
.DE
Since $theta (t)$ is easy to compute, it would suffice to find an efficient method to compute
.DS 3
.EQ (5.2)
f(t) ~=~ sum from n=1 to m ~2^ n sup {-^1/2} ~e sup {it ^ log ^n} ~.
.EN
.DE
Now $m ~=~ left floor (2 pi ) sup {-^ 1/2} ~t sup 1/2 right floor$ is
constant for large stretches of values of $t$.
Suppose that $m$ is constant over the interval $t sub 0 ~<=~ t ~<=~ t sub 1$,
and suppose we wish to find sign changes of
$Z(t)$ in this interval.
Without much loss of generality we can assume that $t sub 0$
is a Gram point, $t sub 0 ~=~ g sub k$.
We compute the complex numbers
.DS 3
.EQ
a sub n ~=~ 2^ n sup {-^1/2} ~e sup {i t sub 0 ^ log ^n} ~~,~~ 1 ~<=~
n ~<=~ m
^,
.EN
.DE
using some very accurate and rigorously analyzable method and
store them.
If $delta ~=~ g sub k+1 ~-~ g sub k$, then we also compute
.DS 3
.EQ
b sub n,j ~=~ e sup {i ^ delta 2 sup -j ^ log ^n} ~~,~~ 1 ~<=~ n ~<=~ m,
~~~0 ~<=~ j ~<=~ J
^,
.EN
.DE
($J$ depending on accuracy to which zeros are to be computed),
and again store them.
To compute $Z( t sub 0 ) ~=~ Z( g sub k )$, we evaluate
$a sub 1 ~+~ a sub 2 ^+...+^ a sub m$.
To compute $Z( g sub k+1 )$, we compute
.DS 3
.EQ (5.3)
a sub n sup (1) ~=~ a sub n ~cdot~ b sub n,0 ~~,~~ 1 ~<=~ n ~<=~ m ^,
.EN
.DE
store them, and then evaluate $a sub 1 sup (1) ^+...+^ a sub m sup (1)$.
To compute $Z(t)$ for some $t$ very close to $g sub k+2$, we compute
.DS 3
.EQ
a sub n sup (2) ~=~ a sub n sup (1) ~cdot~ b sub n,0 ~~,~~ 1 ~<=~ n ~<=~
m ^,
.EN
.DE
store them, and then evaluate $a sub 1 sup (2) ^+...+^ a sub m sup (2)$.
(The differences $g sub l+1 ~-~ g sub l$ change very slowly with $l$.)
To compute $Z(t)$ for some $t$ close to $( g sub k+3 ~+~ g sub k+2 )/2$,
we compute
.DS 3
.EQ
a sub n sup (3) ~=~ a sub n sup (3) ~cdot~ b sub n,1 ~~,~~ 1 ~<=~ n ~<=~
m ^,
.EN
.DE
store them, and evaluate $a sub 1 sup (3) ^+...+^ a sub m sup (3)$,
and so on.
When we reach the point when the accumulated roundoff errors in
the $a sub n sup (h)$ get too big, we reset $t sub 0$ to
the value of a Gram point near $t sub 0$ and recompute the $a sub n$.
Since the recomputation of the $a sub n$ would be very infrequent
in computations such as that aimed at numerically verifying the RH
for a large set of consecutive zeros (where $Z(t)$ is evaluated at
all Gram points), the running time of the algorithm would be dominated
by the time to compute the component by component product
$c ~=~ ( c sub 1 ^,...,^ c sub m ) ~=~ (a sub 1 ~b sub 1 ^,...,^ a sub m ~b sub m )$
of two complex vectors $a ~=~ ( a sub 1 ^,...,^ a sub m )$ and
$b ~=~ ( b sub 1 ^,...,^ b sub m )$ and then evaluate the complex
sum $c sub 1 ^+...+^ c sub m$.
These operations are vectorized automatically by the Cray X-MP
Fortran compiler, and for $m ~=~ 10 sup 4$ require only about
$0.8 ^times^ 10 sup -3$ seconds, which would yield a program
running about 4 times faster than that described in Section 3.
Such a program has not been implemented due to the large storage
requirements, which would not have allowed computations with
zeros around the $10 sup 12$-th zero.
However, for verifying the RH around the $10 sup 10$-th zero such
a Fortran program ought to be practical and, since most of the
computation there involves evaluating $Z(t)$ at Gram points, it
would probably run about 3 times faster on the Cray X-MP than
the program described in [73] run on the Cyber-205.
Some test runs performed on the Cyber-205 by H. te Riele
have shown that the ideas described above
appear to lead to
algorithms that are faster than those of [73] on that machine as well.
.H 1 "Distribution of eigenvalues computations"
In the standard terminology [15, 54, 55] we let $p(k,u)$ denote
the probability density of $k$-th neighbor spacing in the GUE;
i.e., the probability that the (normalized) difference between
an eigenvalue of the GUE and the $(k+1)$-st smallest eigenvalue
of those that are larger than it lies in the interval $(a,b)$ is
.DS 3
.EQ (6.1)
int sub a sup b ~p( k,u ) ~du ~.
.EN
.DE
(This is referred to as $p sub 2 ^(k;u)$ in [15, 54, 55],
where the subscript 2 refers to the GUE.)
.P
Define $lambda sub n$ and $f sub n ^(x)$ to be the eigenvalues
and eigenfunctions of the integral equation (finite Fourier transform)
.DS 3
.EQ (6.2)
lambda ^f(x) ~=~ int sub -1 sup 1 ~{ sin "{" pi t(x-y)/2 "}" } over
{pi (x-y)} ~f(y) ~dy ^,
.EN
.DE
where $t ~member~ roR sup +$ is a parameter, and where the
$f sub n ^(x)$ satisfy the orthonormality condition
.DS 3
.EQ (6.3)
int sub -1 sup 1 ~f sub k ^(x) ~f sub m ^(x) ~dx ~=~ delta sub km ~.
.EN
.DE
.P
The $lambda sub n$ and $f sub n ^(x)$ can be renumbered so that
.DS 3
.EQ
1 ~>=~ lambda sub 0 ~>=~ lambda sub 1 ^>=...>=^ 0 ~.
.EN
.DE
Mehta and des Cloizeaux [55] have shown that
.DS 3
.EQ (6.4)
p(k,t) ~=~ 4 t sup -2 ~prod from n ~( 1- lambda sub n ) ~sum from
{0 ^<=^ j sub 1 ^<...<^ j sub k+2} ~V( j sub 1 ^,...,^ j sub k+2 ) ^,
.EN
.DE
where
.DS 3
.EQ (6.5)
V( j sub 1 ^,...,^ j sub k+2 ) ~=~ 4 ~prod from i=1 to k+2 ~{lambda sub {j sub i}}
over {1- lambda sub {j sub i}} ~sum from { {pile {1 ^<=^ scrL ^<^ m ^<=^ k+2 above
j sub scrL ^+^ j sub m ^==^ 1 ^( mod 2) above scrL ^<^ m}}} ~
f sub {j sub scrL} ^( 1 ) sup 2 ~f sub {j sub m} ^ ( 1 ) sup 2 ~.
.EN
.DE
(Note that notations differ from reference to reference by various
scaling factors.)
The $p(k,t)$ were computed using formulas (6.4) and (6.5), which are
much better than some of the older methods (cf. [54, 55]) which
relied on the formula
.DS 3
.EQ (6.6)
p(0,t) ~=~ {d sup 2} over {dt sup 2} ~prod from n ~( 1- lambda sub n ) ~.
.EN
.DE
For example, the values for $p(0,t)$ given in [54; Appendix A.12,
Table A.1], which were apparently derived using (6.6), seem to be
rather inaccurate, since for $t ~=~ 0.891$, the value given there is
0.943, whereas (6.4) and (6.5) gave the value of 0.933.
The discrepancy was most pronounced near the peak of the graph
of $p(0,t)$, which is where one would expect numerical differentiation
routines to be least accurate.
.P
The $f sub n ^(x)$ are essentially the linear prolate functions.
They were computed using the program of Van Buren [71], with slight
modifications introduced by this author and S. P. Lloyd.
In Van Buren's notation, the $lambda sub n$ defined above become
$lambda sub n ^( pi t/2)$, and the $f sub n ^(x)$ become
($lambda sub n ^( pi t/2 ) ) sup -1/2 ~psi sub n ^( pi t / 2,x)$.
.P
Van Buren's program uses a complicated combination of a variational
procedure and expansion in terms of Legendre and spherical
Bessel functions.
There is no rigorous error analysis for it.
However, it appears to be fairly accurate,
as is shown by comparing its output to that of other calculations.
Furthermore, tests such as numerical integration of
.DS 3
.EQ
int from 0 to inf~p(k,^t)~dt
.EN
.DE
to verify that we obtain values very close to 1 were all
positive.
.HU "Acknowledgments"
The author thanks
P.\ Barrucand, R.\ A. Becker,
O. Bohigas, R. P. Brent, W. S. Cleveland,
F. J. Dyson,
P. X. Gallagher,
E.\ Grosswald,
D. Hejhal, B. Kleiner, H. J. Landau, S. P. Lloyd, C. L. Mallows,
H. L. Montgomery, A. Pandey, H. J. J. te Riele, L. Schoenfeld, N. L. Schryer, D. S. Slepian, and D.\ T. Winter
for providing helpful information.
.bp
.ce
.B
References
.R
.sp
.AL
.LI
R. A. Becker and J. M. Chambers,
\f2S: An Interactive Environment for Data Analysis and Graphics\f1,
Wadsworth, 1984.
.LI
M. V. Berry, Semiclassical theory of spectral rigidity,
\f2Proc. Royal Soc. London A 400\f1 (1985), 229-251.
.LI
M. V. Berry, Riemann's zeta function:\0\0A model for
quantum chaos?, in \f2Proceedings of Second International Conference
on Quantum Chaos\f1, T. Seligman, ed., Springer, 1986.
(To appear.)
.LI
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
O. Bohigas, M. J. Giannoni, and C. Schmit, Characterization
of chaotic quantum spectra and universality of level fluctuation
laws, \f2Physical Rev. Letters 52\f1 (1984), 1-4.
.LI
O. Bohigas, M. J. Giannoni, and C. Schmit, Spectral properties
of the Laplacian and random matrix theories, \f2J. Physique-Lettres\f1,
to appear.
.LI
O. Bohigas, R. U. Haq, and A. Pandey, Higher-order correlations
in spectra of complex systems, \f2Phys. Rev. Letters 54\f1
(1985), 1645-1648.
.LI
E. Bombieri and D. Hejhal, manuscript in preparation.
.LI
R. P. Brent, \f2Algorithms for Minimization without Derivatives\f1,
Prentice-Hall, 1973.
.LI
R. P. Brent, On the zeros of the Riemann zeta function in
the critical strip, \f2Math. Comp., 33\f1 (1979), 1361-1372.
.LI
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
W. S. Brown, A simple but realistic model of floating-point
computations, \f2ACM Trans. Math. Software\f1 \f27\f1 (1981), 445-480.
.LI
J. M. Chambers, W. S. Cleveland, B. Kleiner, and P. A. Tukey,
\f2Graphical Methods for Data Analysis\f1, Wadsworth, 1983.
.LI
J. des Cloizeaux and M. L. Mehta, Some asymptotic expressions
for prolate spheroidal functions and for the eigenvalues of
differential and integral equations of which they are solutions,
\f2J. Math. Phys. 13\f1 (1972), 1745-1754.
.LI
J. des Cloizeaux and M. L. Mehta, Asymptotic behavior of
spacing distributions for the eigenvalues of random matrices,
\f2J. Math. Phys. 14\f1 (1973), 1648-1650.
.LI
J. B. Conrey, A. Ghosh, D. A. Goldston, S. M. Gonek and
D. R. Heath-Brown, Distribution of gaps between zeros of the zeta
function, \f2Quarterly J. of Math. Oxford\f1 (2), \f236\f1 (1985), 43-51.
.LI
F. D. Crary and J. B. Rosser, High precision coefficients
related to the zeta function, MRC Technical Summary Report #1344,
Univ. of Wisconsin, Madison, May 1975, 171 pp.; reviewed by
R. P. Brent in \f2Math. Comp. 31\f1 (1977), 803-804.
.LI
Cray Research, Inc., \f2CRAY X-MP and CRAY-1 Computer Systems;
Library Reference Manual SR-0014\f1, Revision I, Dec. 1984.
.LI
Cray Research, Inc., \f2CRAY-1 Computer Systems, S Series
Mainframe Reference Manual HR-0029\f1, Nov. 1982.
.LI
D.\ Davies, An approximate functional equation for
Dirichlet L-functions,
\f2Proc. Royal Soc. Ser. A 284\f1 (1965), 224-236.
.LI
M. Deuring, Asymptotische Entwicklungen der Dirichletschen
L-Reihen, \f2Math. Ann. 168\f1 (1967), 1-30.
.LI
F. J. Dyson, Statistical theory of the energy levels of
complex systems. II, \f2J. Math. Phys. 3\f1 (1962), 157-165.
.LI
H. M. Edwards, \f2Riemann's Zeta Function\f1, Academic Press, 1974.
.LI
D. Freedman and P. Diaconis, On the histogram as a density
estimator: $L sub 2$ theory,
\f2Z. Wahrscheinlichkeitstheorie v. Geb. 57\f1 (1981), 453-476.
.LI
A. Fujii, On the zeros of Dirichlet $L$-functions.
I. \f2Trans. Amer. Math. Soc.\f1 \f2196\f1 (1974), 225-235.
.LI
A. Fujii, On the uniformity of the distribution of zeros
of the Riemann zeta function, \f2J. reine angew. Math. 302\f1
(1978), 167-205.
.LI
A. Fujii, On the zeros of Dirichlet $L$-functions.
II (With corrections to ``On the zeros of Dirichlet $L$-function. I''
and the subsequent papers), \f2Trans. Amer. Math. Soc.\f1 \f2267\f1 (1981),
33-40.
.LI
W. Gabcke, Neue Herleitung und explicite Restabscha\*:tzung
der Riemann-Siegel-Formel, Ph.D. Dissertation, Go\*:ttingen, 1979.
.LI
P. X. Gallagher, On the distribution of primes in short
intervals, \f2Mathematika 23\f1 (1976), 4-9.
.LI
P. X. Gallagher, Pair correlation of zeros of the zeta
function, \f2J. reine\f1 \f2angew. Math.\f1 \f2362\f1 (1985), 72-86.
.LI
P. X. Gallagher and J. H. Mueller, Primes and zeros in short
intervals, \f2J. reine angew. Math.\f1 303/304 (1978), 205-220.
.LI
A. Ghosh, On Riemann's zeta-function--sign changes of $S(T)$,
pp. 25-46 in \f2Recent Progress in Analytic Number Theory\f1,
vol. 1, H. Halberstam and C. Hooley, eds., Academic Press, 1981.
.LI
A. Ghosh, On the Riemann zeta-function--mean value theorems
and the distribution of $| S(t) |$,
\f2J. Number Theory 17\f1 (1983), 93-102.
.LI
D. A. Goldston, Prime numbers and the pair correlation of
zeros of the zeta function, 
pp. 82-91 in \f2Topics in Analytic Number Theory\f1,
S.\ W. Graham and J.\ D. Vaaler, eds., Univ. Texas Press, 1985.
.LI
D. A. Goldston and D. R. Heath-Brown, A note on the differences
between consecutive primes, \f2Math. Ann.\f1 \f2266\f1 (1984), 317-320.
.LI
D. A. Goldston and H. L. Montgomery, Pair correlation of zeros
and primes in short intervals, Proc. 1984 Stillwater Conference on
Analytic Number Theory and Diophantine Problems, Birkhauser Verlag
(to appear).
.LI
S. M. Gonek, A formula of Landau and mean values of $zeta (s)$,
pp. 92-97 in \f2Topics in Analytic Number Theory\f1, S. W. Graham
and J. D. Vaaler, eds., Univ. Texas Press, 1985.
.LI
A. P. Guinand, A summation formula in the theory of prime
numbers, \f2Proc. London Math. Soc.\f1 (2) \f250\f1 (1948), 107-119.
.LI
M. C. Gutzwiller, Stochastic behavior in quantum scattering,
\f2Physica 7D\f1 (1983), 341-355.
.LI
\f2Handbook of Mathematical Functions\f1, M. Abramowitz and
I. A. Stegun, eds., National Bureau of Standards, 9th printing, 1970.
.LI
R. U. Haq, A. Pandey, and O. Bohigas, Fluctuation properties
of nuclear energy levels: Do theory and experiment agree?
\f2Physical Rev. Letters 48\f1 (1982), 1086-1089.
.LI
D. R. Heath-Brown, Gaps between primes and the pair correlation
of zeros of the zeta-function, \f2Acta Arith.\f1 \f241\f1 (1982), 85-99.
.LI
A. Ivic\*v,
\f2The Riemann Zeta-function\f1, Wiley, 1985.
.LI
D. Joyner, Distribution theorems for $L$-functions, to be published.
.LI
D. Joyner, On the Dyson-Montgomery hypothesis, to be published.
.LI
E. Karkoschka and P. Werner, Einige Ausnahmen zur
Rosserschen Regel in der Theorie der Riemannschen Zetafunktion,
\f2Computing 27\f1 (1981), 57-69.
.LI
M. G. Kendall and A. Stuart, \f2The Advanced Theory of Statistics\f1,
Griffin, 1981.
.LI
E. Landau, U\*:ber die Nullstellen der Zetafunktion,
\f2Math. Ann. 71\f1 (1911), 548-564.
.LI
R. S. Lehman, On the distribution of zeros of the Riemann zeta
function, \f2Proc. London Math. Soc. 20\f1 (1970), 303-320.
.LI
J. van de Lune, Some observations concerning the zero-curves
of the real and imaginary parts of Riemann's zeta function,
Report ZW 201/83, Mathematical Center Amsterdam, December 1983.
.LI
J. van de Lune, H. J. J. te Riele, and D. T. Winter,
Rigorous high speed separation of zeros of Riemann's zeta
function, Report NW 113/81, Mathematical Center, Amsterdam, 1981.
.LI
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
MATHLAB Group, \f2MACSYMA Reference Manual\f1, MIT Laboratory
for Computer Science, 1977.
.LI
M. L. Mehta, \f2Random Matrices\f1, Academic Press, 1967.
.LI
M. L. Mehta and J. des Cloizeaux, The probabilities for
several consecutive eigenvalues of a random matrix,
\f2Indian J. Pure Appl. Math. 3\f1 (1972), 329-351.
.LI
H. L. Montgomery, The pair correlation of zeros of the zeta
function, Proc. Symp. Pure Math. \f224\f1, A.M.S., Providence
1973, 181-193.
.LI
H. L. Montgomery, Distribution of zeros of the Riemann
zeta function, \f2Proc. Int. Congress Math. Vancouver\f1 (1974), 379-381.
.LI
H. L. Montgomery, Extreme values of the Riemann zeta function,
\f2Comm. Math. Helv.\f1 \f252\f1 (1977), 511-518.
.LI
A. M. Odlyzko, Distribution of zeros of the Riemann
zeta function: Conjectures and computations, manuscript in preparation.
.LI
A. M. Odlyzko and H. J. J. te Riele, Disproof of the
Mertens conjecture, \f2J. reine angew. Math. 357\f1 (1985), 138-160.
.LI
A. M. Odlyzko and A. Scho\*:nhage, Fast algorithms for
multiple evaluations of the Riemann zeta function, to be published.
.LI
A. E. Ozluk, Pair correlation of zeros of Dirichlet
L-functions, Ph.D. dissertation, Univ. of Michigan, Ann Arbor,
Mich., 1982.
.LI
C. E. Porter, ed., \f2Statistical Theories of Spectra: Fluctuations\f1,
Academic Press, 1965.
.LI
N. L. Schryer, A test of a computer's floating-point
arithmetic unit, AT&T Bell Laboratories Computing Science Technical
Report #89, 1981.
.LI
N. L. Schryer, manuscript in preparation.
.LI
A. Selberg, Contributions to the theory of the Riemann
zeta-function, \f2Arch. for Math. og Naturv. B, 48\f1 (1946), 89-155.
.LI
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
E. L. Titchmarsh, \f2The Theory of the Riemann Zeta-function\f1,
Oxford Univ. Press, 1951.
.LI
K.-M. Tsang, \f2The Distribution of the Values of the Riemann
Zeta-function\f1, Ph.D. Dissertation, Princeton, 1984.
.LI
A. M. Turing, A method for the calculation of the zeta-function,
\f2Proc. London Math. Soc. ser. 2\f1, \f248\f1 (1943), 180-197.
.LI
A. L. Van Buren, A Fortran computer program for calculating the
linear prolate functions, Report 7994, Naval Research Laboratory,
Washington, May 1976.
.LI
A. Weil, Sur les ``formules explicites'' de la theorie des
nombres premiers, \f2Comm. Sem. Math. Univ. Lund\f1, tome
supplementaire (1952), 252-265.
.LI
D. Winter and H. te Riele, Optimization of a program for the
verification of the Riemann hypothesis, \f2Supercomputer 5\f1 (1985),
29-32.
.LE
.bp
.DS
.TS
center;
c s s s s s s
c c s s c s s
c || c | c | c || c | c | c
n || n | n | n || n | n | n.
.ls 1
\f2Table 1.\f1\0\0Moments of $delta sub n ^-^ 1$ and of $delta sub n ^+^ delta sub n+1 ^-^ 2$, for blocks of the form $N^+^1~<=~n~<=~N~+~10 sup 5$.
.sp
	$delta sub n~-~1$	$delta sub n ^+^ delta sub n+1~-~2$
.sp 2
$k$	$N~=~0$	$N~=~10 sup 12$	GUE	$N~=~0$	$N~=~10 sup 12$	GUE
.sp .2
_
.sp
2	0.161	0.176	0.180	0.207	0.236	0.249
.sp .5
3	0.031	0.035	0.038	0.028	0.027	0.030
.sp .5
4	0.081	0.096	0.101	0.123	0.168	0.185
.sp .5
5	0.046	0.057	0.066	0.047	0.063	0.073
.sp .5
6	0.075	0.098	0.111	0.119	0.206	0.237
.sp .5
7	0.072	0.103	0.124	0.078	0.150	0.185
.sp .5
8	0.103	0.160	0.197	0.155	0.367	0.451
.sp .5
9	0.126	0.223	0.290	0.142	0.409	0.544
.sp .5
10	0.181	0.361	0.488	0.252	0.871	1.178
.TE
.DE
.bp
.DS
.ce 2
\f2Table 2\f1. Moments of log $delta sub n$, $delta sub n sup -1$, and $delta sub n sup -2$,
.sp .5
for blocks of the form $N~+~1~<=~n~<=~N~+~10 sup 5$.
.sp
.TS
center, tab(:);
c | r | r | r.
moments
of:$N~=~0$:$N~=~10 sup 12$:GUE
.sp .3
_
.sp .3
log $delta sub n$:\-0.0912:\-0.1021:\-0.1035
.sp .5
$delta sub n sup -1$:1.2363:1.2744:1.2758
.sp .5
$delta sub n sup -2$:2.2235:2.6149:2.5633
.TE
.DE
.bp
.DS
.TS
center;
c s s s s
c s s s s
c c s c s
c || c | c || c | c
c || n | c || n | c.
.ls 1
\f2Table 3\f1.\0\0Kolmogorov statistic for $delta sub n$ and $delta sub n ^+^ delta sub n+1$,
.sp .5
for blocks of zeros of the form $N^+^1 ~<=~ n ~<=~ N^+^ 10 sup 5$.
.sp
	$delta sub n$	$delta sub n ^+^ delta sub n+1$
.sp 2
$N$	$D$	prob.	$D$	prob.
.sp .5
_
.sp
0	0.0207	$10 sup -37$	0.0278	$10 sup -67$
.sp .5
$10 sup 9$	0.0081	$4 ^cdot^ 10 sup -6$	0.0138	$10 sup -16$
.sp .5
$10 sup 11$	0.0052	$9 ^cdot^ 10 sup -3$	0.0099	$6 ^cdot^ 10 sup -9$
.sp .5
$2 ^cdot^ 10 sup 11$	0.0058	$2 ^cdot^ 10 sup -3$	0.0084	$1.5 ^cdot^ 10 sup -6$
.sp .5
$10 sup 12$	0.0045	$3.5 ^cdot^ 10 sup -2$	0.0091	$1.3 ^cdot^ 10 sup -6$
.TE
.DE
.bp
.DS
.TS
center;
c s s s s s
c s s s s s
c s s s s s
c | c | c | c | c | c
c | c | c | c | c | c
c | n | n | n | n | n.
.ls 1
\f2Table 4\f1.\0\0Numbers of very small and very large
.sp .2
$delta sub n$ and $delta sub n ^+^ delta sub n+1$ in blocks of the form
.sp .2
$N^+^ 1 ~<=~ n ~<=~ N ^+^ 10 sup 5$, and the GUE predictions.
.sp
	number of	number of	number of	number of	number of
.sp .2
	$delta sub n ~<=~ 0.05$	$delta sub n~<=~0.1$	$delta sub n ~>=~ 2.8$	$delta sub n ^+^ delta sub n+1 ~<=~ 0.6$	$delta sub n ^+^ delta sub n+1 ~>=~ 4$
.sp .5
_
.sp
$N~=~ 0$	11	79	1	2	0
.sp .5
$N~=~ 5 ^cdot^ 10 sup 7$	9	103	9	24	3
.sp .5
$N~=~ 10 sup 9$	11	97	14	33	5
.sp .5
$N~=~ 10 sup 11$	16	107	11	41	8
.sp .5
$N~=~ 2^cdot^ 10 sup 11$	20	123	20	25	9
.sp .5
$N~=~ 10 sup 12$	18	124	14	34	6
.sp .5
GUE	13.68	108.80	19.68	38.63	13.57
.TE
.DE
.bp
.DS
.TS
center;
c s s s s s s
c s s s s s s
c s s s s s s
c s s s s s s
c | c | c | c | c | c | c
c | c | c | c | c | c | c
c | c | c | n | n | n | c.
.ls 1
\f2Table 5\f1.\0\0Extremal values of $delta sub n$ and $delta sub n ^+^ delta sub n+1$ among zeros
.sp .2
number $n$, $N^+^1 ~<=~ n ~<=~ N^+^M$, and the probability
.sp .2
that the minimum value of $delta sub n$ in the GUE
.sp .2
would not exceed the value in the third column.
.sp
						prob.
.sp .2
$N$	$M$	$min~delta sub n$	$max~delta sub n$	$min~delta sub n ^+^ delta sub n+1$	$max~delta sub n ^+^ delta sub n+1$	$min~delta sub n$
.sp .5
_
.sp
0	$10 sup 5$	0.02186	2.8052	0.57647	3.9022	0.682
.sp .2
$5 ^cdot^ 10 sup 7$	$10 sup 5$	0.00734	3.1551	0.41841	4.3058	0.042
.sp .2
$10 sup 9$	$10 sup 5$	0.01679	3.0680	0.37404	4.1650	0.405
.sp .2
$10 sup 11$	$10 sup 5$	0.00738	3.2811	0.30784	4.4753	0.043
.sp .2
$2 ^cdot^ 10 sup 11$	$10 sup 5$	0.01515	3.1109	0.39744	4.3121	0.317
.sp .2
$10 sup 12$	$10 sup 5$	0.01103	3.1893	0.44260	4.3371	0.137
.sp .5
_
.sp .5
0	$10 sup 6$	0.00970	3.2143	0.26344	4.1517	0.632
.sp .2
0	$7.5 ^cdot^ 10 sup 7$	$<=~0.00124$				$<=~0.145$
.sp .2
0	$1.5 ^cdot^ 10 sup 9$	$<=~0.000310$				$<=~0.048$
.TE
.DE
.bp
.DS
.TS
center;
c s s
c | c | c
n | n | n.
.ls 1
\f2Table 6\f1.\0\0Autocovariances of the $delta sub n$, $N^+^1 ~<=~ n ~<=~ N^+^ 10 sup 5$.
.sp
$k$	$N~=~0$	$N~=~ 10 sup 12$
.sp .5
_
.sp
0	.1607429	.1762933
1	-.0574023	-.0582498
2	-.0126083	-.0143371
3	-.0065874	-.0047603
4	-.0045317	-.0035270
5	-.0031454	-.0011576
6	-.0011362	-.0015305
7	-.0007084	-.0008964
8	-.0013904	-.0009071
9	.0013483	-.0011673
10	.0034456	-.0001692
11	.0018714	-.0009839
12	-.0002503	-.0002831
13	-.0005412	-.0003967
14	.0025227	-.0011130
15	.0046388	-.0002101
16	.0025451	-.0005048
17	.0010829	.0003850
.sp .5
_
.sp .5
55	-.0039299	-.0102228
56	-.0013716	-.0069591
57	.0032256	-.0011287
58	.0016816	.0012409
59	-.0014840	.0020236
60	-.0041143	.0017857
61	-.0017425	.0011242
62	.0019009	.0015995
63	.0000127	.0012965
.sp .5
_
.sp .5
960	.0002534	.0047463
961	-.0002862	.0013039
962	.0003855	-.0039218
963	-.0002286	-.0061988
964	.0001077	-.0015785
965	-.0000796	-.0011387
966	-.0002891	-.0050911
967	.0000783	-.0053232
968	-.0001676	.0008364
969	-.0000848	.0044897
.TE
.DE
.bp
.DS
.TS
center;
c s s s
c | c | c | c
c | c | c | c
r | n | c | n.
.ls 1
\f2Table 7.\f1\0\0Zeta function near selected van de Lune points.
.sp
		Rosser
.sp .2
$t$	$Z(t)$	rule	$max ~delta sub n$
.sp .5
_
.sp .5
18,132,299,244.660	-133.2	no	2.95
.sp .2
18,139,553,794.750	142.2	no	3.08
.sp .2
907,663,606,940.329	229.3	no	3.34
.sp .2
9,065,450,718,497.579	-253.5	a	3.02
.sp .2
45,323,986,866,893.743	-320.7	(2,2)	3.56
.sp .2
67,263,231,798,214.250	441.4	b	4.11
.sp .2
725,177,880,629,981.915	-453.9	no	3.43
.TE
.DE
.bp
.VL 11
.LI "\f2Fig. 1.\f1"
Pair correlation of zeros of the zeta function.
Solid line:\0\0GUE prediction.
Scatter plot:\0\0empirical data based on zeros $gamma sub n$,
$1 ~<=~ n ~<=~ 10 sup 5$.
.LI "\f2Fig. 2.\f1"
Pair correlation of zeros of the zeta function.
Solid line:\0\0GUE prediction.
Scatter plot:\0\0empirical data based on zeros
$gamma sub n$, $10 sup 12 ^+^ 1 ~<=~ n ~<=~ 10 sup 12 ^+^ 10 sup 5$.
.LI "\f2Fig. 3.\f1"
Probability density of the normalized spacings $delta sub n$.
Solid line:\0\0GUE prediction.
Scatter plot:\0\0empirical data based on zeros $gamma sub n$,
$1 ~<=~ n ~<=~ 10 sup 5$.
.LI "\f2Fig. 4.\f1"
Probability density of the normalized spacings $delta sub n$.
Solid line:\0\0GUE prediction.
Scatter plot:\0\0empirical data based on zeros $gamma sub n$,
$10 sup 12 ^+^1 ~<=~ n ~<=~ 10 sup 12 ^+^ 10 sup 5$.
.LI "\f2Fig. 5.\f1"
Probability density of the normalized spacings $delta sub n ^+^ delta sub n+1$.
Solid line:\0\0GUE prediction.
Scatter plot:\0\0empirical data based on zeros $gamma sub n$,
$1 ~<=~ n ~<=~ 10 sup 5$.
.LI "\f2Fig. 6.\f1"
Probability density of the normalized spacings $delta sub n ^+^ delta sub n+1$.
Solid line:\0\0GUE prediction.
Scatter plot:\0\0empirical data based on zeros $gamma sub n$,
$10 sup 12 ^+^1 ~<=~ n ~<=~ 10 sup 12 ^+^ 10 sup 5$.
.LI "\f2Fig. 7.\f1"
Initial segment of the quantile-quantile plot of the normalized
spacings $delta sub n$ against the GUE prediction.
Data based on zeros $gamma sub n$, $1 ~<=~ n ~<=~ 10 sup 5$.
Straight line $y ~=~ x$ drawn to facilitate comparisons.
.LI "\f2Fig. 8.\f1"
Initial segment of the quantile-quantile plot of the normalized
spacings $delta sub n$ against the GUE prediction.
Data based on zeros $gamma sub n$, $10 sup 12 ~+~ 1 ~<=~ n ~<=~ 10 sup 12 ~+~ 10 sup 5$.
Straight line $y ~=~ x$ drawn to facilitate comparisons.
.LI "\f2Fig. 9.\f1"
Final segment of the quantile-quantile plot of the normalized
spacings $delta sub n ~+~ delta sub n+1$ against the GUE prediction.
Data based on zeros $gamma sub n$, $1 ~<=~ n ~<=~ 10 sup 5$.
Straight line $y ~=~ x$ drawn to facilitate comparisons.
.LI "\f2Fig. 10.\f1"
Final segment of the quantile-quantile plot of the normalized
spacings $delta sub n ~+~ delta sub n+1$ against the GUE prediction.
Data based on zeros $gamma sub n$, $10 sup 12 ~+~ 1 ~<=~ n ~<=~ 10 sup 12 ~+~ 10 sup 5$.
Straight line $y ~=~ x$ drawn to facilitate comparisons.
.LI "\f2Fig. 11.\f1"
Initial segment of the quantile-quantile plot of the normalized
spacings $delta sub n$ against the GUE prediction.
Data based on zeros $gamma sub n$, $2 ^cdot^ 10 sup 11 ~+~ 1 ~<=~ n ~<=~ 2 ^cdot^ 10 sup 11 ~+~ 10 sup 5$.
Straight line $y ~=~ x$ drawn to facilitate comparisons.
.LI "\f2Fig. 12.\f1"
Autocovariance of the $delta sub n$.
The entry for $k$ is the mean value of
$( delta sub n ^-^ 1 ) ~( delta sub n+k ^-^ 1)$,
averaged over $10 sup 12 ~+~ 1 ~<=~ n ~<=~ 10 sup 12 ~+~ 10 sup 5$.
.LI "\f2Fig. 13.\f1"
Spectrum of the $delta sub n$.
Value plotted for $x$ is $c^ | sum from n ~( delta sub n ^-^ 1) ~exp ( pi i n x ) | sup 2$,
where $c$ is a constant, and $n$ runs over
$10 sup 12 ~+~ 1 ~<=~ n ~<=~ 10 sup 12 ~+~ 98303$.
(Heights of peaks are not represented accurately due to limited sampling.)
.LI "\f2Fig. 14.\f1"
Logarithm of the spectrum of the $delta sub n$.
Value plotted for $x$ is $2^log^ ^ | sum from n ~( delta sub n ^-^ 1) ~exp ( pi ^in^ x) | ~+~ c$ for
a constant $c$, where $n$ runs over $10 sup 12 ~+~ 1 ~<=~ n ~<=~ 10 sup 12 ~+~ 98303$,
and where values $< ~-~ 10$ are replaced by -10.
.LI "\f2Fig. 15.\f1"
Graph of $2^ log ^ | sum ~exp ( i ^gamma sub n ^y) |$ for $0.004 ~<=~ y ~<=~ 3$, where
$n$ runs over $10 sup 12 ~+~ 1 ~<=~ n ~<=~ 10 sup 12 ~+~ 40000$,
and values $< 0$ are replaced by 0.
.LE