.PH ""
.nr Pt 1
.EQ
delim $$
gsize 12
define Cbar % size 9 bold up 14 | back 35 size 12 roman C %
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
.bp
.ls 1
.ce 99
.S 12
.B
New analytic algorithms in number theory
.R
.sp
.I
A. M. Odlyzko
.R
.sp .2
AT&T Bell Laboratories
Murray Hill, NJ 07974
USA
.ce 0
.sp
.H 1 "Introduction"
.P
.ls 2
Several new number theoretic algorithms are sketched here.
They all have the common feature that they rely on
bounded precision computations of analytic functions.
The main one of these algorithms is a new method, due
to A. Scho\*:nhage and the author, of calculating
values of the Riemann zeta function at multiple points.
This method enables one to verify the truth of the Riemann
Hypothesis (RH) for the first $n$ zeros in what is expected to
be $n sup {1 ^+^ o(1)}$ operations as $n ~->~ inf$, as
opposed to about $n sup 3/2$ operations for earlier methods.
This technique can also be extended to the evaluation of
other zeta and $L$-functions.
.pn 2
.PH "''- \\\\nP -''"
.P
Methods for the approximate computation of zeta and $L$-functions
have been shown by J. C. Lagarias and the author to lead to
fast algorithms for the exact computation of certain integer-valued number
theoretic functions.
It is shown below that these methods enable one to compute $pi (x)$, the
number of primes $<= ^x$, in $x sup {1/2 ^+^ o(1)}$ operations
as $x ~->~ inf$.
.P
Numerical computations of zeros of the zeta function
cannot ever prove the RH, but if there is
a counterexample to the RH, such computations might
find it.
This was undoubtedly the main motivation for the computations
that have been done.
The first few nontrivial zeros were computed by hand
by Riemann [Ed, Si], and further computations, by hand or with
the help of electromechanical technology, were carried
out in the early 1900's, culminating with the verification
of the RH for the first 1041 zeros (i.e., the 1041 zeros in
the upper half-plane $Im (s) ~>~ 0$ that are closest to the
real axis $Im (s) ~=~ 0$) by Titchmarsh and Comrie [25].
After World War II, these computations were extended using
electronic digital computers.
The latest result in this direction is the verification
of the RH for the first $1.5 ~cdot~ 10 sup 9$ zeros [16], a computation
that involved about two months on a modern supercomputer,
and used the same basic method for computing the zeta function
that was employed by Titchmarsh and Comrie.
.P
Computations of the zeros of the Riemann zeta function
and related functions (such as Dirichlet $L$-functions
and Epstein zeta functions) are also useful in providing
evidence for and against various other conjectures,
such as the Montgomery pair correlation conjecture
[17, 18],
and some related conjectures which predict
that the zeros of $zeta (s)$ ought to behave like eigenvalues
of random hermitian matrices, (whose distribution is known
relatively well, since they have been studied by physicists
who use them to model energy levels in many-particle systems).
For a discussion of these conjectures and the numerical
evidence about them, see [17, 18, 19, 20].
Another, similar, motivation for computing zeros of zeta
and $L$-functions comes from trying to understand better
the differences between those functions, like the Riemann
zeta function and Dirichlet $L$-function, that are expected
to satisfy the RH, and those, such as the Epstein zeta
functions, which are in most cases expected to violate
the RH.
See [2, 9] for some result on these problems.
Other applications for approximate values of the zeta
function arise in disproving certain conjectures, see [10, 21].
.P
As will be explained briefly in Section 2, the problem of
computing an approximate value of a zero of $zeta (s)$,
for example, and proving that this zero is on the critical
line, can be reduced to the problem of evaluating $zeta (s)$ to a certain
accuracy.
In Section 3, the previous methods of computing $zeta (s)$ will be
presented.
The fastest of them, the Riemann-Siegel formula, requires on
the order of $t sup 1/2$ operations to evaluate $zeta (1/2 ~+~ it )$ to
within $+- ~t sup -c$, say, for $t ~>~ 1, ~c ~>~ 0$.
In Section 4 new methods for evaluating $zeta (s)$ at multiple points
are presented.
The fastest of them, due to A. Scho\*:nhage and the author [22],
can compute any single value of $zeta ( 1/2 ~+~ it)$ with
$T ~<=~ t ~<=~ T ~+~ T sup 1/2$ to within $+- ~T sup -c$ in
$T sup o(1)$ operations (on numbers of $O( log ^T)$ bits), provided
a precomputation involving $T sup {1/2 ^+^ o(1)}$ operations is
carried out beforehand
and $T sup {1/2 ^+^ o(1)}$ storage locations are available.
This method leads to an algorithm for numerically verifying the RH
for the first $n$ zeros in what is expected to be $n sup {1 ^+^ o(1)}$
operations, as opposed to about $n sup 3/2$ operations using the
older methods.
These techniques extend to the computation of other zeta
and $L$-functions.
.P
The computation of zeros to a specified accuracy requires only the
computations of approximate values of the appropriate function.
It turns out that such approximate values can be
used to compute exactly various number theoretic functions,
such as $pi (x)$, the number of primes $<= ^x$.
The reason for this is related to the ``explicit formulas''
of Guinand [8] and Weil [29], which relate arithmetic functions
to zeros of the appropriate zeta and related functions.
An even older formula of Landau [15] (see also [7]) says that
for any fixed $y ~>~ 0$, as $T ~->~ inf$,
.DS 3
.EQ
sum from {0 ^<^ Im ( rho ) ^<^ T} ~e sup {rho y} ~=~ left {
matrix {
lcol {-^ T over {2 pi} ~ log ^p ~+~ O( log ^ T) above O( log ^T)}
lcol {roman if ~y ~=~ log ^ p sup m ^, above roman if ~y ~!=~ log ^p sup m ^,}
}
.EN
.DE
where $p$ denotes a prime and $m$ a positive integer, and $rho$
runs over the zeros of the zeta function.
The formula above can be used to test integers for
primality.
This idea has occurred to many people, but unfortunately it
does not seem to yield an efficient algorithm.
However, modifications of such ideas can be used to
compute relatively efficiently functions such as $pi (x)$, as was shown by
J. C. Lagarias and the author [13, 14].
Using the latest algorithm for evaluating the zeta function,
it is possible to compute $pi (x)$ in $x sup {1/2 ^+^ o(1)}$ steps
as $x ~->~ inf$, which is faster than the best known combinatorial
algorithm [12],
which requires around $x sup 2/3$ operations.
This new algorithm for $pi (x)$ and similar functions is sketched
briefly in Section 5.
.H 1 "Zeta function and the numerical verification of the RH"
The zeta function is defined for $Re (s) ~>~ 1$ by
.DS 3
.EQ (2.1)
zeta (s) ~=~ sum from n=1 to inf ~1 over {n sup s} ^,
.EN
.DE
but it can be extended by analytic continuation to an analytic
function throughout $Cbar ^\e^ "{" 1 "}"$, and it has a first
order pole at $s ~=~ 1$.
If we define, as usual,
.DS 3
.EQ (2.2)
xi (s) ~=~ pi sup {-^ s/2} ~GAMMA (s/2) ~zeta (s) ^,
.EN
.DE
then $xi (s)$ satisfies the functional equation
.DS 3
.EQ (2.3)
xi (s) ~=~ xi (1-s) ~~,~~ s ~member~ Cbar ^\e^ "{" 0,1"}" ~.
.EN
.DE
Since $xi (s)$ is real for real $s$, we have $xi ( s bar ) ~=~ {xi (s)} bar$,
and so for $t ~member~ roR$,
.DS 3
.EQ (2.4)
xi ( 1/2 ^+^ it) ~=~ xi (1 ^-^ ( 1/2 ^+^ it)) ~=~ xi ( 1/2 ^-^ it) ~=~
{xi (1/2 ^+^ it)} bar ^,
.EN
.DE
which implies that $xi ( 1/2 ^+^ it) ~member~ roR$.
Therefore if we find two values $t sub 1$ and $t sub 2$ such
that $xi ( t sub 1 ) ~<~ 0 ~<~ xi ( t sub 2 )$ (and the errors
in the computed values of $xi (t sub 1 )$ and $xi ( t sub 2 )$ are
smaller than those values, so we can be certain that the two
inequalities are valid), then we can conclude that there is a
value $t prime , ~t sub 1 ~<~ t prime ~<~ t sub 2$, such that
$zeta ( 1/2 ^+^ i t prime ) ~=~ 0$.
In this case we know that the zero is right on the critical line,
not just in some neighborhood of it.
.P
All the numerical verifications of the RH have relied on the use
of the intermediate value theorem as sketched above.
This theorem gives a lower bound for the number of zeros of
the zeta function that are on the critical line up to some height.
To assert that all of the zeros up to some height are on the
critical line, one can use the principle of the argument to
determine the total number of zeros up to that height and if
that agrees with the number of zeros that have been located on
the critical line, the verification is completed.
In practice it is not even necessary to use the principle of
the argument, since there is an elegant method of Turing [5, 20, 28],
based on a theorem of Littlewood about the zeta function,
which usually enables one to conclude that all the zeros up to some
height are on the critical line based only on computations of
$xi (s)$ on that line.
.P
It should be mentioned that the above strategy for numerically
verifying the RH is not guaranteed to work, even if the RH is
true, since there might be multiple zeros on the critical line.
However, such zeros have not been encountered yet and are not
expected to exist.
In practice, using the sophisticated strategies for locating
points $t$ at which to evaluate $xi ( 1/2 ^+^ it)$ that have
been developed [16], it takes on average about 1.2 evaluation
of $xi (s)$ to establish that a zero is on the critical line.
.P
The above discussion serves to reduce the problem of numerically
verifying the RH to that of computing $xi (s)$.
The next section discusses that problem.
.H 1 "Methods for calculating the zeta function"
For purposes of numerically verifying the RH or for computing
$pi (x)$ by the algorithm to presented in Section 5, it is
sufficient to be able to compute $zeta (s)$ to within $+- ~ | s | sup -c$
for various constants $c ~>~ 0$.
For some purposes (such as those of Section 5) it is necessary
to calculate $zeta (s)$ for $Re (s) ~>~ 1$, but to keep the
presentation simple we will deal only with $Re (s) ~=~ 1/2$ here.
(See [22] for more details.)
Since $GAMMA (s/2)$ and $pi sup {-^ s/2}$ are very easy to compute
accurately, computing
$xi (1/2 ^+^ it)$ is essentially the same as computing $zeta (1/2 ^+^ it)$.
.P
The simplest method for calculating $zeta (s)$ is to apply the
Euler-Maclaurin summation formula to (2.1).
This gives an analytic continuation of $zeta (s)$ to all of
$Cbar ^\e^ "{" 1 "}"$, and enables one to compute $zeta (s)$
anyplace to arbitrary accuracy by adjusting the parameters
in the formula appropriately.
However, this method is not very efficient, since to compute
$zeta ( 1/2 ^+^ it)$ to $+- ~t sup -1$, say, already takes
about $t$ steps.
Still, this was the method that was used by the first few people
to publish calculations of the zeros of the zeta function.
.P
A much more efficient method for computing $zeta ( 1/2 ^+^ it)$ was
known to Riemann, and was used by him to compute the first few
zeros.
It is one of what Hardy and Littlewood called approximate functional equation,
but it has the nice feature that the error term is given by
an asymptotic expansion.
Discovered in Riemann's unpublished papers by Siegel, it
became known as the Riemann-Siegel formula [5, 23].
It has been used for all large-scale computations of zeros
of the zeta function since the 1930's.
(The only other method to have been proposed so far, that of Turing [27],
was meant to be used at relatively small heights, where it turned out not
to be necessary.)
.P
Let
.DS 3
.EQ (3.1)
theta (t) ~mark =~ roman arg ^left [ pi sup {-^ it/2} ~GAMMA ( 1/4 ~+~ it /2) right ] ^,
.EN
.sp
.EQ (3.2)
Z(t) ~lineup =~ exp ^( i theta (t)) ~zeta ( 1/2 ~+~ it ) ^,
.EN
.DE
so that
.DS 3
.EQ
Z(t) ~=~ {xi (1/2 ^+^ it)} over {| pi sup {-^ 1/4 ^-^ it/2} ~GAMMA ( 1/4 ^+^ it/2)|} ~.
.EN
.DE
Also let
.DS 3
.EQ
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
Then the Riemann-Siegel formula [5, 6, 11, 20, 23, 25, 26] says that
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.3)
~~~~~
.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
R sub k ^( tau ) ~=~ O( tau sup {-^ (2k ^+^ 3)/4} ) ^,
.EN
.DE
and the $PHI sub j ^(z)$ are certain entire functions.
Explicit bounds for the $R sub k ^( tau )$ are known (the
best ones being due to Gabcke [6]), and bounds for the
error made in computing the $PHI sub j ^(x)$ are also known.
Thus the entire difficulty in the computation of $Z(t)$ resides
in the need to evaluate the sum of cosines in (3.3), which
takes on the order of $t sup 1/2$ steps.
.H 1 "New methods for evaluating the zeta function"
The Riemann-Siegel formula is still the fastest known method
for computing a single value of $zeta ( 1/2 ^+^ it)$ to medium
accuracy ($+- ~t sup -c$ for some $c ~>~ 0$) for large $t$.
Here we will show that if values of $zeta ( 1/2 ^+^ it)$ at a
large set of closely spaced $t$'s are desired, one can obtain
much faster algorithms.
.P
In the previous section we showed that the bottleneck in using
the Riemann-Siegel formula is the need to evaluate the sum of
about $t sup 1/2$ cosines.
(In the Euler-Maclaurin formula, one has to evaluate a similar
sum of about $t$ cosines.)
Note that
.DS 3
.EQ
2 ~ sum from n=1 to m ~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
so that it suffices to find an efficient method to compute the
complex-valued function
.DS 3
.EQ (4.1)
f(t) ~=~ sum from n=1 to m ~2 n sup {-^ 1/2} ~e sup {it ^ log ^n}
.EN
.DE
for the values of $t$ that are of interest.
We will count the number of elementary arithmetic operations on
numbers of $O( log ^t)$ digits.
(See [22] for more on the model of computation that's used.)
What makes the methods below work is that there is a lot
of structure in sums of the form (4.1), and if many of them
are to be computed, one can exploit this structure to lower
the average cost of an evaluation.
.P
The first observation to make is that if we are interested in computing
$f(t)$ (and thus $zeta ( 1/2 ^+^ it)$) for $T ~<=~ t ~<=~ T ~+~ T sup alpha$
for some $0 ~<~ alpha ~<=~ 1$, say, then it suffices to compute
$f(t)$ and a few similar sums at an evenly spaced grid of points,
$t ~=~ t sub 0 , ~t sub 1 ^,...,^ t sub H$, where
.DS 3
.EQ
t sub j ~=~ T ~+~ j delta ~~,~~ 0 ~<=~ j ~<=~ H ~=~ left floor T sup {alpha ^+^ eta} right floor ,
~delta ~=~ T sup {-^ eta} ^,
.EN
.DE
for some $eta ~epsilon~ (0,1 /10)$.
(Eventually $eta$ will go to zero.)
The reason for this is that
.DS 3
.EQ (4.2)
f sup (k) ^(t) ~=~ sum from n=1 to m ~2 ^i sup k ^n sup {-^ 1/2} ~
( log ^n) sup k ~e sup {it ^ log ^n} ^,
.EN
.DE
so that
.DS 3
.EQ
| f sup (k) ^(t) | ~<=~ 5 ( log ^ m ) sup k ~m sup 1/2 ~.
.EN
.DE
Therefore if we wish to compute $f(t)$ for some $t$,
$T ~<=~ t ~<=~ T ~+~ T sup alpha$, to within $+- ~T sup {-c} ,$ we can
use the Taylor series expansion around the nearest $t sub j$, provided
we have precomputed the values of $f sup (k) ^( t sub j )$ for
$k ~<=~ ( c ^+^ 100 ) eta sup -1 $, say.
Since the derivatives $f sup (k) ^(t)$ are of the same form as
$f(t)$, we have reduced the problem to that of computing series
of the form
.DS 3
.EQ (4.3)
g(t) ~=~ sum from n=1 to m ~a sub n ~e sup {it ^ log ^n}
.EN
.DE
at $t ~=~ t sub j ~=~ T ~+~ j delta$ for $0 ~<=~ j ~<=~ H$.
We should note that $m$ depends on $t$, $m ~=~ left floor (2 pi ) sup {-^ 1/2} ~t sup 1/2 right floor$.
.P
The trivial way to compute $g(t)$ at each of the $t sub j$ takes
about $m H$ operations.
One way to obtain a faster method is to use fast matrix multiplication
methods.
For simplicity, suppose that we are trying to compute $g( delta ^k)$
for $0 ~<=~ k ~<~ R sup 2$, where $delta ~<=~ 1$.
Consider the matrix $C ~=~ ( c sub jh )$, $0 ~<=~ j,h ~<=~ R-1$, where
.DS 3
.EQ (4.4)
c sub jh ~=~ left {
matrix {
lcol {a sub h ~exp ( ij (R-1) delta ~ log ^h) ^, above 0 ^,}
lcol {roman for ~h ~<=~ beta sub j ^, above roman for ~h ~>~ beta sub j ^,}
}
.EN
.DE
where the $beta sub j$ will be specified later.
Also, let $D ~=~ ( d sub hk )$, $0 ~<=~ h,k ~<=~ R-1$, be the
matrix with
.DS 3
.EQ (4.5)
d sub hk ~=~ exp ( i ^k^ delta ~ log ^h) ~.
.EN
.DE
Then the matrix $CD ~=~ ( b sub jk ) , ~0 ~<=~ j,k ~<=~ R-1$, has
.DS 3
.EQ (4.6)
b sub jk ~=~ sum from h=1 to {beta sub j} ~a sub h ~exp ( i [ j (R-1)
^+^ k] delta ~ log ^h) ~.
.EN
.DE
If we select $beta sub j ~=~ left floor (2 pi ) sup {-^ 1/2} ~j sup 1/2 ~( R-1 ) sup 1/2 right floor$,
then $b sub jk$ will differ from $g( delta ( j (R-1) ^+^ k))$ by
only a few terms, which can be easily computed.
The trivial method for multiplying $C$ and $D$ takes about $R sup 3$
operations, which is not much better than the ordinary $O( R sup 3 ^delta sup 1/2 )$
method of evaluating each sum of the form (4.6).
However, if we use fast matrix multiplication methods, of which
the latest one allows one to multiply two $n ~times~ n$ matrices
in time $O( n sup 2.479 )$ [24], we obtain a method with running
time about $R sup 2.479$.
Since in practice we would have $delta ~=~ R sup {-^ o(1)}, ~R ~=~ T sup {1/2 ^+^ o(1)}$, this
would let us compute all values of $zeta ( 1/2 ^+^ it)$ for $0 ~<=~ t ~<=~ T$
in time $T sup o(1)$ to within $+- ~T sup -c$ after an initial
precomputation that requires $<= ~T sup 1.24$ operations, and would thus be more
efficient than using the Riemann-Siegel formula.
(This method could also be modified to work on shorter intervals,
but not as efficiently.)
Unfortunately fast matrix multiplication methods are completely
impractical (except possibly for Strassen's original
$O ( n sup 2.807... )$ method for multiplying two $n ~times~ n$
matrices [1]), so that the above technique is of purely theoretical
significance.
.P
An algorithm that is both theoretically efficient and appears to be practical
can be obtained by a different approach.
Suppose that we wish to compute $g( t sub j ), ~0 ~<=~ j ~<=~ H$,
where $t sub j ~=~ T ~+~ j delta , ~0 ~<~ delta ~<~ 1/10$, and
$H ~<=~ T sup 1/2$.
Let $M ~=~ left floor (2 pi ) sup {-^ 1/2} ~T sup 1/2 right floor$, and
.DS 3
.EQ
g sup * ^(t) ~=~ sum from n=2 to M ~a sub n ^e sup {it ^ log ^n} ~.
.EN
.DE
Then $g( t sub j )$ and $g sup * ^( t sub j )$ differ at most
by two terms (for $T$ large), and so it suffices
to compute the $g sup * ^( t sub j )$ efficiently.
Choose $N ~=~ 2 sup r$ so that $H ~<~ N ~<=~ 2H$, and let
.DS 3
.EQ
omega ~=~ exp ( 2 ^ pi ^i /N ) ~.
.EN
.DE
Let
.DS 3
.EQ
h(k) ~=~ sum from j=0 to N-1 ~g sup * ^( t sub j ) ^omega sup jk ~~,~~
0 ~<=~ k ~<~ N ~.
.EN
.DE
The inverse Fourier transform gives
.DS 3
.EQ
g sup * ^( t sub j ) ~=~ 1 over N ~sum from k=0 to N-1 ~h(k) ^ omega sup {-^ jk}
~~,~~ 0 ~<=~ j ~<=~ H ^,
.EN
.DE
so if we can compute all the $h(k)$ accurately, we will be able
to recover the $g sup * ^( t sub j )$ in $O( N ^ log ^N)$ operations
using the Fast Fourier Transform [1].
(All the numbers that come up in this procedure are of moderate size,
so the required precision does not cause too much of a problem.)
.P
The $h(k)$ can be rewritten as
.DS 3
.EQ
h(k) ~mark =~ sum from n=2 to M ~a sub n ^e sup {T ^ log ^n}
~sum from j=0 to N-1 ~e sup {ij ^ delta ~ log ^n} ~omega sup jk
.EN
.sp
.EQ
~lineup =~ sum from n=2 to M ~{a sub n ^e sup {iT ^ log ^n}} over
{1 ~-~ omega sup k ~e sup {i delta ^ log ^n}} ~=~ sum from n=2 to M
~{-^ a sub n ^e sup {-^i ^( delta ^-^ T) ^ log ^n}} over {omega sup k ~-~ e sup {-^ i delta ^ log ^n}} ~.
.EN
.DE
(We are assuming here for simplicity that $omega sup k ~-~ exp (-^ i delta ~ log ^n)$
is not too small for all $k$ and $n$.
.P
Those $n$ for which the above difference is small or even
zero for some $k$ have to be treated separately, cf. [22].)
Thus evaluating $h(k)$ for $0 ~<=~ k ~<=~ N-1$ is no harder
than evaluating a rational function of the type
.DS 3
.EQ (4.7)
h sup * ^(z) ~=~ sum from n=1 to N ~{alpha sub n} over {z ^-^ beta sub n}
.EN
.DE
at the points $z ~=~ omega sup h ~,~ 0 ~<=~ k ~<=~ N-1$, where
the $alpha sub n$ and $beta sub n$ can complex numbers with
$| beta sub n | ~=~ 1$, the $alpha sub n$ are not too large,
and $| omega sup k ^-^ beta sub n |$ is not too small for all $k$ and $n$.
It is easy to show functions of the form (4.7) can be evaluated
in $O( N ^( log ^N ) sup 2 )$ operations over some finite
fields, by employing Fast Fourier Transforms methods.
What is more surprising is evaluations in $O( N ^( log ^N ) sup 2 )$
operations on numbers of $O( log ^N)$ bit, can also be performed.
This was shown by A. Scho\*:nhage, and is presented in [22].
The basic idea is to use Taylor series expansions of $h sup * ^(z)$
around a small set of points.
A set of functions $h sub p,q ^(z)$ is defined, each of which consists
of some of the terms in (4.7), such that the Taylor series of
each $h sub p,q ^(z)$ around a certain point converges in a relatively
large region, and such that for each $k$, $h sup * ^( omega sup k )$ can
be written as a sum of a subset of the $h sub p,q ^( omega sup k )$.
Full details can be found in [22].
.H 1 "Analytic algorithms for arithmetic functions"
In this section we will sketch the analytic algorithms for
certain arithmetic function that have been developed recently.
For simplicity of presentation, we will restrict ourselves to
the problem of computing $pi (x)$.
Full details and extensions to other functions can be found in [14].
.P
In spite of extensive interest in computing $pi (x)$, until recently
no algorithm with proved running time of $O( x sup {1 ^-^ epsilon} )$
for some $epsilon ~>~ 0$ was known.
(See [12, 13] for historical references.)
In [12] a modification of the Meissel-Lehmer algorithm, which uses
combinatorial sieving, was developed, which runs in time $x sup {2/3 ^+^ o(1)}$
and space $x sup {1/3 ^+^ o(1)}$.
This algorithm has been implemented and has been used to compute
values of $pi (x)$ for various $x$ up to $x ~=~ 4 ~cdot~ 10 sup 16$.
At about the same time, a new analytic algorithm for computing $pi (x)$
was developed [13, 14].
It relies on numerical integration of analytic transforms involving
the zeta function.
When one uses the Euler-Maclacurin formula to compute the zeta function,
the resulting algorithm requires time $x sup {2/3 ^+^ o(1)}$ and space
$x sup o(1)$.
With the Riemann-Siegel formula, the running time decreases to
$x sup {3/5 ^+^ o(1)}$ and the space requirement stays at $x sup o(1)$.
With the use of the new algorithm described in Section 4, time
drops to $x sup {1/2 ^+^ o(1)}$, but space required becomes
$x sup {1/4 ^+^ o(1)}$.
.P
The analytic algorithm for $pi (x)$ is based on an old formula of Riemann,
namely that
.DS 3
.EQ (5.1)
pi (x) ~+~ sum from j=2 to inf ~1 over j ^pi ^( x sup 1/j ) ~=~ 1 over
{2 ^ pi ^i} ~int from {2 ^-^ i inf} to {2 ^+^ i inf} ~{x sup s} over s
~ log ^zeta (s) ~ds
.EN
.DE
for $x ~!=~ p sup m$ for any prime $p$ and $m ~epsilon~ Z sup +$.
(If $x ~=~ p sup m , ~( 2m ) sup -1$ must be subtracted from the
left side of (5.1).)
One of the things which make the analytic algorithm work is that
$pi (x)$ is always an integer, and so it suffices to compute it
to within $+- ~1/3$, say.
The sum of the terms on the left side of (5.1) with $j ~>=~ 2$
can easily be computed to within $+- ~1/10$ in about $x sup 1/2$ steps.
(Even fewer are required if one uses the algorithm recursively.)
Thus if we could compute the integral on the right side of
(5.1) efficiently to within $+- ~1/10$, say, we would compute
$pi (x)$ rapidly.
Unfortunately the integral in (5.1) is not even absolutely convergent,
and no way has been found to compute it fast.
.P
To overcome the problem of slow convergence of the integral in (5.1),
we use a different formula, namely
.DS 3
.EQ (5.2)
sum from {p ^<=^ x} ~c(p) ~+~ sum from {{pile {p,m above p sup m ^<=^ x above m ^>=^ 2}}} ~1 over m ~c ( p sup m ) ~=~ 1 over {2^pi^i} ~int from {2 ^-^ i inf} to{2 ^+^ i inf}
~F(s) ~ log ~zeta (s) ~ds ^,
.EN
.DE
where
.DS 3
.EQ (5.3)
c(u) ~=~ 1 over {2^pi^i} ~int from {2 ^-^ i inf} to {2 ^+^ i inf} ~
F(s) ~u sup -s ~ds ^,
.EN
.DE
and where this formula holds whenever $F(s)$ is a sufficiently
nice function.
($F(s)$ analytic in $Re (s) ~>~ 1$, $| F(s) | ~=~ | s | sup -2$ as
$| s | ~->~ inf$, is sufficient but not necessary, for example.)
We would like to choose $F(s)$ so that (5.2) behaves like (5.1);
i.e., $c(p) ~=~ 1$ for $p ~<=~ x$ and $c (p) ~=~ 0$ for $p ~>~ x$.
Unfortunately, the uniqueness theorem for Mellin transforms
says that is impossible unless $F(s) ~=~ s sup -1 ~x sup s$.
We therefore select a somewhat different $F(s)$, in which the
sum on the left side of (5.2) is relatively close to that on
the left side of (5.1).
To be precise, we select $F(s)$ so that
.VL 8
.LI (i)
$F(2 ^+^ it) ~->~ 0$ rapidly as $| t | ~->~ inf$, so that the integral
in (5.2) is approximated well by a finite integral
.LE
.DS 3
.EQ (5.4)
1 over {2^ pi^ i} ~int from {2 ^-^ iT} to {2 ^+^ iT} ~F(s) ~ log ~zeta (s)
~ds ^;
.EN
.DE
.VL 8
.LI (ii)
Individual values of $F(s)$ can be computed rapidly, so that
the integral above can be evaluated by numerical integration;
.LI (iii)
$c(p) ~=~ 1$ for primes $p ~<=~ x-y$, $c(p) ~=~ 0$ for
primes $p ~>=~ x$, and $0 ~<=~ c(p) ~<=~ 1$ for primes $p ~epsilon~ (x-y ^,^ x)$;
.LI (iv)
Individual values of $c(k)$ can be computed efficiently.
.LE
Such $F(s)$ can be found for any $y$ and $T$ with the property
that $yT ~>=~ x sup {1 ^+^ o(1)}$, see [13, 14].
(The requirement that $y T ~>=~ x sup {1 ^+^ o(1)}$ is caused by
the uncertaintly principle of Fourier analysis, which says that
a function and its transform cannot both decrease too rapidly.)
The contribution of the sum
.DS 3
.EQ
sum from {{pile {p,m above p sup m ^<=^ x above m ^>=^ 2}}} ~1 over m
~c ( p sup m )
.EN
.DE
can be evaluated to within $+- ~1/100$ in $x sup {1/2 ^+^ o(1)}$ operations.
The difference
.DS 3
.EQ
pi (x) ~-~ sum from {p ^<=^ x} ~c(p) ~=~ sum from {x-y ^<^ p ^<=^ x} ~
(1 ~-~ c(p))
.EN
.DE
can be evaluated to within $+- ~1/100$ by finding all the primes
in $[x-y ^,^ x]$ and computing their contribution, and this can
be done in $y x sup o(1)$ operations.
Finally, the evaluation of the integral (5.4) can be reduced to evaluating
$zeta (2 ^+^ it)$ at a grid of points between $-T$ and $T$, spaced
$x sup {-^ o(1)}$ apart, and by the algorithm of Section 4 this can
be done in $T ^x sup o(1)$ operations.
Thus the running time of the algorithm is
.DS 3
.EQ
x sup {1/2 ^+^ o(1)} ~+~ y ^x sup o(1) ~+~ T ^x sup o(1) ^,
.EN
.DE
and selecting $y ~=~ T ~=~ x sup {1/2 ^+^ o(1)}$, we obtain
the desired running time bound.
Full details are presented in [14].
.bp
.ce
\f3REFERENCES\f1
.sp
.AL
.LI
A. V. Aho, J. E. Hopcroft, and J. D. Ullman, \f2The Design and Analysis of
Computer Algorithms\f1, Addison-Wesley, 1974.
.LI
E. Bombieri and D. A. Hejhal, manuscript in preparation.
.LI
D. Davies, An approximate functional equation for Dirichlet $L$-functions,
\f2Proc. Royal Soc. Ser.\f1 \f2A 284\f1 (1965), 224-236.
.LI
M. Deuring, Asymptotische Entwicklungen der
Dirichletschen L-Reihen, \f2Math. Annalen 168\f1
(1967), 1-30.
.LI
H. M. Edwards, \f2Riemann's Zeta Function\f1,
Academic Press, 1974.
.LI
W. Gabcke, Neue Herleitung und explizite Restabscha\*:tzung
der Riemann-Siegel-Formel, Ph.D. Dissertation,
Go\*:ttingen 1979.
.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
D. A. Hejhal, Zeros of Epstein zeta functions and supercomputers,
\f2Proc. 1986 Intern. Congress Math.\f1, to appear.
.LI
A. E. Ingham, On two conjectures in the theory of numbers,
\f2Amer. J. Math. 64\f1 (1942), 313-319.
.LI
A. Ivic\*v,
\f2The Riemann Zeta-function\f1, Wiley, 1985.
.LI
J. C. Lagarias, V. S. Miller, and A. M. Odlyzko,
Computing $pi (x)$: the Meissel-Lehmer
method, \f2Math. Comp\f1., \f244\f1 (1985), 537-560.
.LI
J. C. Lagarias and A. M. Odlyzko, New algorithms
for computing $pi (x)$, pp. 176-193 in \f2Number Theory: New York 1982\f1,
D. V. Chudnovsky, G. V. Chudnovsky, H. Cohn, and M. B. Nathanson, eds.,
Lecture Notes in Mathematics #1052, Springer-Verlag, 1984.
.LI
J. C. Lagarias and A. M. Odlyzko, Computing $pi (x)$:\0\0An analytic
method, \f2J. Algorithms\f1, to appear.
.LI
E. Landau, U\*:ber die Nullstellen der Zetafunktion,
\f2Math. Ann. 71\f1 (1911), 548-564.
.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
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
A. M. Odlyzko, On the distribution of spacings between zeros
of the zeta function, \f2Math. Comp.\f1, to appear.
.LI
A. M. Odlyzko, On the distribution of zeros of the Riemann zeta
function:\0\0Conjectures 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
C. L. Siegel, ${roman U} dotdot$ber Riemanns Nachlass zur
analytischen Zahlentheorie, \f2Quellen und
Studien zur Geschichte der Math. Astr. Phys.\f1 2
(1932), 45-80.
Reprinted in C.\ L. Siegel,
\f2Gesammelte Abhandlungen,\f1 Spring-Verlag, 1966, Vol. 1, pp. 275-310.
.LI
V. Strassen, Relative bilinear complexity and matrix multiplication,
to be published.
.LI
E. C. Titchmarsh, The zeros of the Riemann Zeta-function,
\f2Proc. Royal Soc. London 151\f1 (1935), 234-255 and
\f2157\f1 (1936), 261-263.
.LI
E. C. Titchmarsh, \f2The Theory of the Riemann
Zeta-Function\f1, Oxford Univ. Press, 1951.
.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. M. Turing, Some calculations of the Riemann Zeta-function,
\f2Proc. London Math. Soc.\f1 (3) \f23\f1 (1953), 99-117.
.LI
A. Weil, Sur les ``formules explicites'' de la theorie des
nombres premiers, \f2Comm. Sem. Math. Univ. Lund\f1,
tome supplementaire (1952), 252-265.
.LE