.nr Pt 1
.ds HP 10 10 10 10 10
.ds HF 3 3 3 3 3 3
.S 11
.EQ
delim $$
gsize 11
define ref % "" sup %
define lcg % "\b'\(lc\(bv'" %
define rcg % "\b'\(rc\(bv'" %
define lfl % "\b'\(bv\(lf'" %
define rfl % "\b'\(bv\(rf'" %
.EN
.am DE
.sp
..
.ND ""
.TL
Computing \(*p(x): An Analytic Method
.AU "J. C. Lagarias" JL MH 11218 2966 2C-373
.AU "A. M. Odlyzko" AMO MH
.AS 1
.ls 2
.P
This paper presents a class of algorithms for
computing $pi (x)$, the number of primes $<=~ x$.
For each $b$ with $0 ~<=~ b ~<=~ 1/4$ and each
$epsilon ~>~ 0$ there is an algorithm $A(b, epsilon )$
which computes $pi (x)$ using $O( x sup {3/5 ^-^ 2b/5 ~+~ epsilon} )$
bit operations and uses $O( x sup {b ~+~ epsilon} )$ bits
of storage.
In particular one can compute $pi (x)$ using $O( x sup {3/5 ~+~ epsilon} )$
bit operations and $O( x sup epsilon )$ space, or using
$O( x sup {1/2 ~+~ epsilon} )$ bit operations and
$O( x sup {1/4 ~+~ epsilon} )$ space.
These algorithms are based on numerical
integration of certain integral transforms
of the Riemann $zeta$-function.
This technique can be generalized to evaluate
many other arithmetic functions, including
the functions
$pi (x ; ^ k , ^ scrL )$ counting the
number of primes $p~==~scrL$ (mod $k$) with
$p~<=~x$ and the function $M(x)$
which is the partial sum of the Mo\*:bius
function $mu (n)$ for all $n~<=~x$.
.AE
.ls 1
.MT 4
.ls 2
.H 1 "Introduction"
The problem of computing $pi (x)$, the number
of primes $p~<=~x$, is a very old one.
The sieve of Eratosthenes, in its usual form,
finds all primes $p~<=~x$ in
.DS 2
.EQ
O(x~log~log~x)
.EN
.DE
steps, where we consider
for our model of computation a random access
machine with about $x~lcg log sub 2 x rcg$-bit
storage locations, and where a ``step'' is defined as any elementary
arithmetic operation on integers $<=~x$.
Recently an improved version of this sieve
was found\ [18], which requires only
.DS 2
.EQ
O( x over {log~log~x} )
.EN
.DE
steps on such a machine.
On the other hand, by the Prime Number Theorem
.DS 2
.EQ (1.1)
pi (x)~wig~x over {log~x}~~roman as~~x~->~inf~,
.EN
.DE
so one cannot hope to find an algorithm
that finds all primes $p~<=~x$ in fewer than
about $x( log~x) sup -1$ operations.
.P
For a long time no method for computing
$pi (x)$ faster than finding all primes
$p~<=~x$ was known.
For example, Gauss was led to conjecture
the Prime Number Theorem\ (1.1) by inspection
of counts of primes in various intervals
that had been compiled by various authors,
but those counts were apparently obtained
by the sieve of Eratosthenes.
Legendre\ [5] used the principle of
inclusion-exclusion to obtain the formula
.DS 2
.EQ
pi (x) - pi ( sqrt x )+1~=~[x]-
sum from {p <= x sup 1/2} [ x over p ]+
sum from {p < q <= x sup 1/2} [ x over pq ]-...,
.EN
.DE
where $p$ and $q$ run over primes.
The sum on the right side above contains
approximately $6 pi sup -2 (1- log~2)x$
nonzero terms, and so is not suitable
for computation.
The first substantial progress towards an
efficient algorithm was made in the second
half of the 19-th century by the
astronomer Meissel, who found a combinatorial
method\ [15] that enabled him to compute
several large values of $pi (x)$.
Meissel's work culminated in his publication
of the value of $pi (10 sup 9 )$\ [16]
(which turned out to be too small by 56).
Many other mathematicians came up with
variants of Meissel's method (see\ [5]), but in
the era of hand computation none were willing
to try their methods on $pi (x)$ for
$x~>~10 sup 9$.
(It was pointed out by B. Berndt that there
is a curious table of values of $pi (x)$
for a few $x~<=~10 sup 8$ in a notebook
of Ramanujan\ [19].
They are all either correct or else very
close to the true value, but there is no
indication how they were derived.
They are more accurate than the values given
by Ramanujan's approximate formulas
for $pi (x)$ [9].)
.P
Further progress on computing $pi (x)$
was made by D.\ H. Lehmer\ [12], who
significantly simplified and generalized
Meissel's method and used it to compute
$pi (10 sup 10 )$ on an electronic
computer.
(His value of $pi ( 10 sup 10 )$ turned out to be too high
by 1\ [2].)
Further computations were carried out
by Mapes\ [14] and Bohman\ [2], who computed
several values up to and including
$pi ( 10 sup 13 )$.
(Bohman's value of $pi (10 sup 13 )$ turned
out to be too low by 941 [10].)
.P
The asymptotic running times of Meissel's and Lehmer's algorithms
for computing $pi (x)$ are rather difficult
to estimate, but as is pointed out in\ [10,11],
all the published versions seem to require
time $>=~c( epsilon ) x sup {1- epsilon}$
for any $epsilon~>~0$.
However, V. Miller and the authors of this paper
obtained a new version of that
method which requires only $O( x sup {2/3+ epsilon} )$
time and $O(x sup {1/3+ epsilon })$ space to compute $pi (x)$,
see [10].
This version of the Meissel-Lehmer
method has been used to compute values of
$pi (x)$ up to $x~=~4 times 10 sup 16$.
.P
In this paper we present new algorithms for
computing $pi (x)$, based on entirely
different ideas, which require $O(x sup {3/5+ epsilon })$
time and $O(x sup epsilon )$ space
for any $epsilon~>~0$
in one version and $O( x sup {1/2+ epsilon} )$ time and
$O( x sup {1/4+ epsilon} )$ space in another version.
We call this class of algorithms \f2analytic\f1
$pi (x)$-\f2algorithms\f1.
.R
An informal description of the algorithms
is presented in Section\ 2.
That section also discusses generalizations
of these algorithms to compute other
arithmetical functions.
Section\ 3 provides a more formal
definition of the algorithms together with
a rigorous estimate of their running time,
while Section\ 4 contains the proof of
the lemmas needed in Section\ 3.
.P
The announcement [11] of analytic algorithms for
computing $pi (x)$
presented only the $O( x sup {3/5+ epsilon} )$
time and $O( x sup epsilon )$ space algorithm.
A new method of computing multiple values of the Riemann zeta function [17] enables us
to describe here, for each $b$ with $0~<=~ b ~<=~ 1/4$, an algorithm
for computing $pi (x)$ that runs in time $O( x sup {(3 ^-^ 2b) / 5+ epsilon} )$
and space $O( x sup {b+ epsilon} )$ for every $epsilon > 0$.
.P
Our work on this subject was stimulated by
H.\ S. Wilf's note\ [23], which discusses
the question of what is a satisfactory
answer to an enumerative problem.
.H 1 "Informal Description and Discussion"
We shall carry out the complexity analysis of
analytic $pi (x)$-algorithms using for
a model of computation a random access
machine (RAM) with $O( log ^x)$ bit addresses.
In fact the analysis can be carried over to apply to
a multi-tape
Turing machine, since the $x sup epsilon$
factors are more than sufficient to take
care of the extra complexity introduced
by counting bit operations, and for the specific application
that is made in this paper, the algorithm of [17] can be
implemented on a multi-tape Turing machine.
This result contrasts with the Meissel-Lehmer-type
algorithm of \ [10] which relies on
sieving in an essential way, and sieving
cannot be implemented quickly on a
Turing machine.
For simplicity of presentation, we will only count
the number of elementary operations on $O( log ^x)$ bit
numbers that our algorithms perform.
.P
The analytic $pi (x)$-algorithms are motivated by the exact
formulas of prime number theory, which
express many arithmetical functions in terms
of zeros of the Riemann zeta and related
functions.
For example, Riemann's formula\ [6] states
that if
.DS 2
.EQ
pi sup star (x)~=~sum from p~
{sum prime} from {cpile {m=1 above p sup m <= x}} to inf~
1 over m~,
.EN
.DE
where $p$ runs over the primes, and the
prime in the summation indicates that
if $x~=~p sup m$ for some prime $p$,
then one should take $(2m) sup -1$ instead
of $m sup -1$ in the sum, then
.DS 2
.EQ (2.1)
pi sup star (x)~=~Li(x)- sum from rho~Li(x sup rho )~-~
log~2~+~int from x to inf~dt over {t(t sup 2 -1) roman log~t}~,
.EN
.DE
where $Li(u)$ is given by
.DS 2
.EQ
Li(u)~=~int from 0 to u~dt over {log~t}~,
.EN
.DE
the value of the integral being its Cauchy principal
value, and $rho$ runs over
zeros of $zeta (s)$ with $roman Re ( rho )~>~0$.
.P
The class of analytic $pi (x)$-algorithms is actually developed from integral
transform formulas from which ``explicit formulas''
like (2.1) are derived.
For example, (2.1) comes from the identity
.DS 2
.EQ (2.2)
pi sup star (x)~=~sum from p~
{sum prime} from {cpile {m=1 above p sup m <= x}} to inf~
1 over m~=~1 over {2 pi i}~int from (a)~
{x sup s} over s~log~zeta (s)~ds~,
.EN
.DE
where $(a)$ for $a~>~1$ denotes the straight
line $s~=~a+it$,
$- inf~<~t~<~inf$.
Eq.\ (2.2) follows immediately from
the Euler product for $zeta (s)$, which
yields for $roman Re (s)~>~1$
.DS 2
.EQ (2.3)
log~zeta (s)~=~log~prod from p (1-p sup -s ) sup -1~=~
sum from p ~log (1-p sup -s ) sup -1~=~
sum from p~sum from m=1 to inf~1 over m (p sup m ) sup -s~,
.EN
.DE
and the classical formula [3]
.DS 2
.EQ (2.4)
1 over {2 pi i}~int from (a)~{(x/n) sup s} over s~ds~=~
left "{"
matrix {
lcol {1 above 1/2 above 0}
ccol {"" above "" above ""}
ccol {"" above "" above ""}
lcol {n~<~x, above n~=~x, above n~>~x~.}
}
.EN
.DE
The sum on the left side of (2.2) is almost
exactly what we want to compute; it differs
from $pi (x)$ (for $x~nomem~Z$,
say) only by the contribution of primes
$p~<=~sqrt x$, and those terms can
be computed in $O(x sup {1/2+ epsilon })$
steps to within $+- 1/10$, say.
(These terms could be computed much faster
by the recursive use of an analytic $pi (x)$-algorithm,
but that is not important.)
Since $pi (x)$ is an integer, to know
it exactly it would therefore suffice
to compute the integral on the right
side of (2.2) to within $+- 1/5$.
Unfortunately that is not easy to do,
since on any line $s~=~a+it$,
$|x sup s |~=~x sup a$, and $log~zeta (s)$
oscillates boundedly, so the integral
is not absolutely convergent.
The rate of convergence can be estimated\ [3],
but it is not sufficient for our purposes.
This seems to be the reason that
Riesel and Go\*:hl\ [20], who investigated
formulas like (2.1) numerically,
did not obtain very good results.
.P
The main idea of the analytic $pi (x)$-algorithms is
to use, in place of (2.3), another of
the family of Mellin transform identities
to which (2.2) belongs, namely
.DS 2
.EQ (2.5)
sum from p~sum from m=1 to inf~
{c(p sup m )} over m~=~1 over {2 pi i}~
int from (a)~F(s)~log~zeta (s) ~ds~,
.EN
.DE
where $c(u)$ and $F(s)$ are
a Mellin transform pair:
.DS 3
.EQ (2.6)
c(u)~mark =~1 over {2 pi i}~int from (a)~
F(s)u sup -s ds~,
.EN
.sp
.EQ (2.7)
F(s)~lineup =~int from 0 to inf~c(u)u sup s-1 du~.
.EN
.DE
(Here $c(u)$ and $F(s)$ have to be sufficiently ``nice''
for the integrals in (2.6) and (2.7) to
converge, but this will easily hold for
the functions we will use.)
Ideally we would like to utilize the previously
used step function $c sub 0 (u)$
with $c sub 0 (u)~=~1$ for $u~<~x$ and
$c sub 0 (u)~=~0$ for $u~>~x$, but then,
by (2.7), we would again have $F(s)~=~s sup -1 x sup s$,
which makes the integral difficult to compute.
Instead, we use a function $c(u)$ that
is very close to $c sub 0 (u)$; for an
appropriately chosen parameter $y~member~(1,x)$,
we will have $c(u)~=~1$ for $0~<=~u~<=~x-y$,
$c(u)~=~0$ for $u~>=~x$, and $0~<=~c(u)~<=~1$ for
$x-y~<=~u~<=~x$.
In the interval $(x-y,x)$,
$c(u)$ will be chosen so as to make
$F(s)$ decrease very rapidly as
$| roman Im (s)|~->~inf$, which will make it
possible to evaluate the integral in
(2.4) numerically.
The integral, of course, will no
longer equal $pi sup star (x)$.
However, the difference will be
.DS 2
.EQ (2.8)
pi sup star (x) - sum from p~sum from m=1 to inf~
{c(p sup m )} over m~
~=~ sum from {cpile {p above x-y}} sum from
{cpile {m=1 above < p sup m < x}} to inf~
{1-c(p sup m )} over m~,
.EN
.DE
which depends only on the
$p sup m~member~(x-y,x)$.
All these prime powers can be determined
(using primality testing, or sieving) in
about $y$ steps $(O(yx sup epsilon )$,
to be precise).
Therefore, if $c(u)$ is easy to compute
accurately, we can obtain $pi sup star (x)$,
and hence $pi (x)$, in about $y$ steps
(for $y~>~x sup 1/2$, say) once we compute
the integral in (2.5).
It turns out that one can choose $c(u)$
so that in the integral in (2.5) the
contribution of the region with
$| roman Im (s)|~>=~T$ is negligible for
$T~>wig~xy sup -1$.
Now to compute the integral over the range
$| roman Im (s)|~<=~T$ numerically requires
$O(Tx sup epsilon )$ evaluations of the integrand,
which, provided $F(s)$ is easy to compute,
reduces to evaluating $zeta (s)$.
Given any $beta$, $0 ~<=~ beta ~<=~ 1/2$, the algorithm
of [17] computes all the necessary
values of $zeta (a+it)$ for $u ~<=~ t ~<=~ u ~+~ u sup beta$
in about $u sup 1/2$ operations using $u sup beta$ storage,
so the computation of the integral takes about $T sup {3/2 ^-^ beta}$
steps and space $T sup beta$.
(For technical reasons, in some ranges we
use the Euler-Maclaurin formula to compute
$zeta (s)$, but this does not affect the
running time analysis.)
Thus the total running time of the resulting algorithm
is on the order of
.DS 2
.EQ
y~+~T sup {3/2 ^-^ beta} ~,
.EN
.DE
where $yT~>wig~x$, and the space requirement is
on the order of $T sup beta$.
We select $beta ~=~ 5b/ (2(1+b))$,
$T ~=~ x sup {b/ beta}$, and $y$ on the order of $x/T$,
and obtain the desired result.
.P
In the next two sections we show that
the heuristic analysis presented above
can be made rigorous.
Three main ingredients are used
which make the new algorithm fast:
1) the ability to choose
$c(u)$ so that both $c(u)$ and $F(s)$
are easy to compute, and $F(s)$ decreases
rapidly, while $c(u)$ is nearly a step function,
2) the ability to find
all the primes and prime powers in the interval
$(x-y,x)$ in time not much larger
than $y$, and 3) the ability to compute
the integral numerically quickly,
using the
algorithm of [17]
to compute values of $zeta (s)$.
It is not possible to improve on the first two of
these ingredients.
Since a function and its Mellin transform
cannot both be rapidly decreasing
(this is related to the uncertainty
principle of quantum mechanics), we cannot
choose $c(u)$ with $c(u)~=~1$
for $0~<=~u~<=~x-y$ and $c(u)~=~0$
for $u~>=~x$ for which $F(s)$ will
decrease significantly faster than
the one we utilize (cf.\ [13]).
As far as the second ingredient is concerned, there are at least $y sup {1- epsilon}$
primes and prime powers in the interval
$(x-y,x)$, while we can find them in time
$O(y sup {1+ epsilon} )$,
for any $epsilon~>~0$, so little improvement is possible in this area.
Improvements in the algorithm could come from
improving the third ingredient.
The new algorithm of [17] enables one to compute
individual values of $zeta (s)$ to reasonable accuracy
(within $| s | sup -c$ for any $c > 0$) in average
time $O( | s | sup epsilon )$ (although at the cost
of requiring $O( | s | sup 1/2 )$ storage), which is
essentially optimal.
However,
the new algorithm only requires the
computation of the integral in (2.5),
and it is conceivable that
there are ways to compute the whole
integral directly using various
functional equations or integral
representations that would be faster than
using quadrature formulae
but we have
not found any way to do this.
.P
In terms of space requirements, numerical
integration requires only $O(x sup epsilon )$
locations.
Finding all the primes $p~member~(x-y,x)$
requires more than $O(x sup epsilon )$
locations if one uses sieving,
but the recent primality testing method
of Adleman et\ al.\ [1] makes it possible
to bypass that problem and implement the entire
algorithm in $O(x sup epsilon )$ storage
locations.
Thus the only place where nontrivial storage
is required is in the implementation of the algorithm
of [17] for computing the zeta function.
It should also be noted that analytic $pi (x)$-algorithms
lend themselves to parallel implementation, as
both numerical integration and
primality testing can be carried out by breaking
the ranges of operation into blocks and
performing the operation independently on each block.
Thus with $m$ processors one can
achieve essentially $m$-fold speedup.
.P
It is possible to replace the $O(x sup epsilon )$
terms in the bounds on both the running
times and the space requirements
in these algorithms by smaller and more
explicit expressions.
To simplify the exposition, we
have not done so.
In order to implement these algorithms,
it would be necessary to optimize many parameters
and to explicitly determine various constants,
and this would require
substantial additional work.
For example, one would need to determine
explicit estimates for the remainder
terms in the Riemann-Siegel formula in
the region $roman Re (s)~>~1/2$, similar to those
that have been obtained for $roman Re (s)~=~1/2$\ [7].
It would also be necessary to investigate
the various possible choices of the
function $c(u)$, balancing the need for
ease of computing $c(u)$ and $F(s)$
against the requirement that $F(s)$
decrease rapidly as $| Im (s)|~->~inf$.
(See\ [13] for some possible functions
to use.)
In practice, the integral in (2.5) would
probably be estimated by deforming
the contour of integration to go partially
through the critical strip, which would
create additional computational difficulties.
Instead of first estimating $pi sup star (x)$
and then computing $pi (x)$, it might
be advantageous to estimate $pi (x)$
itself, which involves using a formula
somewhat different from (2.2).
It would also be necessary to investigate
various quadrature rules for numerical
integration to determine which is best.
Finally, it would be necessary to implement the
algorithm of [17].
In summary, implementing an analytic $pi (x)$ algorithm that is competitive
with the Meissel-Lehmer method \f2in practice\f1
would require substantial additional analysis.
.P
The basic ideas underlying analytic $pi (x)$-algorithms
can be used
for computing other arithmetical functions,
such as
.DS 2
.EQ
M(x)~=~sum from {n <= x} mu (n)~,
.EN
.DE
where $mu (n)$ is the
Mo\*:bius
$mu$-function, or $pi (x;q,a)$, the number
of primes $p~<=~x$,
$p~==~a$ (mod $q$).
Generally speaking, the analytic $pi (x)$-algorithm can
be adapted to compute
.DS 2
.EQ
G(x)~=~sum from {n <= x}~g(n)~,
.EN
.DE
where $g(n)$ is an integer-valued function (so that
computing $G(x)$ exactly makes sense)
such that the Dirichlet series
.DS 2
.EQ
sum from n=1 to inf~g(n) over {n sup s}
.EN
.DE
is easy to compute for large $roman Re (s)$ (in
the case of $M(x)$, this series equals
$zeta (s) sup -1$, whereas for
$pi (x;q,a)$ it equals a combination of logarithms
of Dirichlet $L$-functions
$L(s, chi )$, which can be computed by a generalized
Riemann-Siegel type formula\ [4] and the method of [17]), and
where the values of $g(n)$ for
$x-y~<~n~<~x$ can be computed fast.
The exact time and storage requirements depend
on the method of computing $G(x)$
and the $g(n)$ for $x-y~<~n~<~x$.
.H 1 "The Algorithm"
Given any $b member [0, 1/4]$ and say positive $delta$,
we describe an algorithm for computing $pi (x)$,
which we call Algorithm $A( b, delta )$,
and show that it computes $pi (x)$ in
$O(x sup {(3 ^-^ 2b) /5 + delta} )$
elementary operations on numbers of $O( log ^x)$ bits using $O(x sup {b+ delta} )$ space as
$x~->~inf$.
(The constants implied by the O-notation
depend on $delta$).
The verification of the correctness of
the algorithm and of the bounds for
the storage and running time depend on a number of lemmas,
which are stated in this section
and proved in Section\ 4.
At the end of the section we indicate how the algorithms
$A( b, delta )$ may be combined to give
a single algorithm,
Algorithm $A(b)$, which computes $pi (x)$
in time $O( x sup {(3 ^-^ 2b) /5 + epsilon} )$
and space $O(x sup {b+ epsilon} )$ as
$x~->~inf$ for any $epsilon~>~0$.
.P
Algorithm $A( b, delta )$ is given run-time exponents $b$ and
$delta$ and a test value $x$ as input,
and produces $pi (x)$ as output.
The run-time exponent $delta$ is used to
compute a parameter $k~=~[200~delta sup -1 ]~+~1$,
which is held fixed in the rest of the algorithm.
In particular the kernel functions
$F sub x,y (s)$ and $g(w)$
and the auxiliary function
$h sub a,x,y (t)$ defined below depend on
the parameter $k$.
.P
We first describe the kernel function
$F(s)$ that we use.
We start with the polynomial
.DS 2
.EQ (3.1)
f(v)~=~f sub k (v)~=~c sub k ^ v sup k (1-v ) sup k~,
.EN
.DE
where $c sub k$ is a normalization constant chosen to
make
.DS 2
.EQ (3.2)
int sub 0 sup 1 f sub k (v)~dv~=~1~.
.EN
.DE
(Note that $c sub k$ is rational.)
We define the parametrized family
$F sub x,y (s)$ of kernel functions by
.DS 2
.EQ (3.3)
F sub x,y (s)~=~s sup -1~
int sub 0 sup 1 (x-vy) sup s
f sub k (v) dv~.
.EN
.DE
The integral defines $F sub x,y (s)$ as
an analytic function of $s$ for
$sigma~=~Re (s)~>~0$, provided $0~<=~y~<=~x$.
We will use $F sub x,y (s)$ and
its inverse Mellin transform
after making an appropriate choice of
the parameter $y$, which turns out to be
$y~=~x sup {(3 ^-^ 2b) / (5 ^-^ 2b)}$.
Our first lemma computes this Mellin
transform pair explicitly and bounds
the growth of $F sub x,y (s)$.
.sp
.I
Lemma 1.  (1) We have
.R
.DS 2
.EQ (3.4)
F sub x,y (s)~=~[y sup k s(s+1)...(s+k)] sup -1~
sum from r=1 to k+1~
{a sub r x sup s+k+r~-~b sub r (x-y) sup s+k+r} over
{y sup r (s+k+1)...(s+k+r)}~,
.EN
.DE
.I where
$~~a sub r~=~f sub k sup (k+r-1) (0)~,~~~~
b sub r~=~f sub k sup (k+r-1) (1)~=~
(-1) sup k+r-1 a sub r~.$
.sp
(2)  \f2For\f1  $sigma~=~Re (s)~>~0$
and $|s|~>=~1$,
.DS 2
.EQ (3.5)
|F sub x,y (s)|~<=~(2k) sup k
x sup {sigma +k} y sup -k |s| sup -k-1~.
.EN
.DE
.I
(3)  For $u~>~0$ and $a~>~0$
.R
.DS 2
.EQ (3.6)
1 over {2 pi i} int sub {a-i inf} sup {a+i inf}
F sub x,y (s) u sup -s ds~=~
g sub k left (
x-u over y right )~,
.EN
.DE
.I where
.DS 2
.EQ (3.7)
g sub k (w)~=~left {
lpile {
0 above {int sub 0 sup w f sub k (v)dv}
above 1}~~~~
lpile {
w~<=~0~, above
0~<~w~<~1~, above w~>=~1~~.
}
.EN
.DE
.P
Our object is to compute $pi (x)$ by evaluating
the contour integral
.DS 2
.EQ (3.8)
P(x,y)~=~1 over {2 pi i} ^ int
sub {a-i inf} sup {a+i inf}
F sub x,y (s)~log ^ zeta (s)~~ds
.EN
.DE
in two different ways,
to within an error of $+- 1 over 10$.
We choose the vertical line of integration
to be $a~=~3 over 2$.
.P
First we evaluate $P(x,y)$ in terms of
$pi (x)$.
Using the inverse Mellin transform, we have
by (2.3) and (3.6) that
.DS 2
.EQ (3.9)
P(x,y)~=~pi (x)~+~
sum from m=2 to inf 1 over m
sum from p from {\s+2p sup m~<=~x\s0}
1 ~-~sum from m=1 to inf
~1 over m ~ sum from p from
{\s+2x-y<p sup m <x\f0}
left ( 1~-~g
left ( {x-p sup m} over y right ) right )~.
.EN
.DE
The difference between $P(x,y)$ and $pi (x)$
can be computed quickly, as is shown in
the next lemma.
.sp
.I
Lemma 2.  A machine with $O(x sup epsilon )$
bits of storage can compute
.R
.DS 2
.EQ (3.10)
sum from m=2 to inf 1 over m
sum from p from {\s+2p sup m <= x\s0} 1
.EN
.DE
.I
to within $+- 10 sup -3$ in $O(x sup {1/2+ epsilon /10} )$
operations for $x~>=~x sub 0 ( epsilon )$.
Such a machine can compute
.R
.DS 2
.EQ (3.11)
sum from m=1 to inf ~1 over m
~sum from p
from {x-y < p sup m < x}~
left ( 1~-~g left ( {x-p sup m} over y
right ) right )
.EN
.DE
.I
to within $+- 10 sup -3$ in $O(yx sup {epsilon /10} )$
operations for $x~>=~x sub 0 ( epsilon )$.
.R
.P
Second we evaluate
$P(x,y)$ in (3.8) directly to within $+- 10 sup -1$.
By the bound (3.5) of Lemma 1, for
$T~>=~10$, say,
.DS 2
.EQ (3.12)
left | 1 over {2 pi i}~int from 3/2+iT to {3/2+i inf}
F sub x,y (s)~log~zeta (s)~ds right |~=~
O(x sup 3/2+k y sup -k T sup -k )~,
.EN
.DE
so if we choose
.DS 2
.EQ (3.13)
T~=~y sup -1 x sup {1+ delta /10}~,
.EN
.DE
then the quantity on the right side of
(3.12) will be $O(x sup -8 )$, since
$k~>=~200 ^ delta sup -1$.
Therefore, for the above choice of $T$,
we will have
.DS 2
.EQ (3.14)
|P(x,y)~-~1 over {2 pi i}~int from 3/2-iT to 3/2+iT~
F sub x,y (s)~log~zeta (s)~ds|~=~O(x sup -8 )~.
.EN
.DE
.P
We have thus reduced the problem of
computing $P(x,y)$ to that of estimating
a finite integral.
In our case we have ${F sub x,y} bar~=~F sub x,y ( s bar )$
and ${zeta (s)} bar~=~zeta ( s bar )$, so
the finite integral is
.DS 2
.EQ (3.15)
1 over {2 pi i}~int from 3/2-iT to 3/2+iT
F sub x,y (s)~
log~zeta (s)~ds~=~1 over pi~int from 0 to T~
h sub 3/2 (t)dt~,
.EN
.DE
where
.DS 2
.EQ (3.16)
h sub a (t)~=~
h sub a,x,y (t)~=~
Re ^ "{" F sub x,y (a+it)~log~zeta (a+it) "}"~.
.EN
.DE
We evaluate this last integral
numerically using the
Euler-Maclaurin
summation formula
([8], Eqn. 25.4.7).
.sp
.I
Lemma 3.
\0Let $g(t)$ be complex-valued,
with its first $2m$ derivatives continuous on
[0,T].
Then for any positive integer $n$ we have
.R
.DS 3
.EQ
int sub 0 sup T g(t)dt~= mark~
T over n left [
sum from j=0 to n g ( jT over n )~-~
{g(0)^+^g(T)} over 2
right ]
.EN
.sp
.EQ (3.17)
lineup~-~sum from i=1 to m-1~
{B sub 2i} over (2i)!~
( T over n ) sup 2i-1
[g sup (2i-1) (T)~-~g sup (2i-1) (0) ]
~+~R sub 2m~,
.EN
.DE
.I
where the $B sub 2i$ are the Bernoulli
numbers, and the remainder
$R sub 2m$ is bounded by
.R
.DS 2
.EQ (3.18)
|R sub 2m |~<=~n~{|B sub 2m |} over
(2m)!~( T over n ) sup 2m+1~
max from {0 <= t <= T}~|g sup (2m) (t)|~.
.EN
.DE
.P
We will apply Lemma 3 for
$g(t)~=~h sub a (t)$
with $a~=~3/2$.
To estimate the remainder term,
we observe that $h sub a,x,y (t)~epsilon~
C sup inf (- inf , inf )$ and
its derivatives are easily bounded,
as follows.
.sp
.I
Lemma 4.  For integers $m~>=~1$,
.R
.DS 2
.EQ (3.19)
left | {d sup m} over {dt sup m}~
h sub a,x,y (t) right |~=~
O left (
m! ( a-1 over 2 ) sup -m~
x sup {(a+1)/2^+^k}
right )~,
.EN
.DE
.I
and the constant implied by the
O-symbol depends on $k$ and $a$.
.R
.P
Now we apply Lemma\ 3, choosing
$a~=~3/2$ with number of steps
.DS 2
.EQ (3.20)
n~=~left [
T^x sup {delta /10}
right ]~=~left [
y sup -1 ^x sup {1^+^ delta /5}
right ]
.EN
.DE
and number of derivatives $2m~=~2k sup 2$.
Using Lemma\ 4 to estimate the remainder term
gives
.DS 2
.EQ (3.21)
|R sub {2k sup 2} |~=~
O(x sup -5 )~,
.EN
.DE
where the implied constant depends on
$k~=~[200^ delta sup -1 ]~+~1$.
.R
.P
It remains to evaluate the main term in
the Euler-Maclaurin summation formula (3.17).
Since $F(3/2+it)$ and all its derivatives
up to some fixed order can be computed
to within $O(x sup -10 )$ in
$O(x sup {delta /10 })$ operations,
and they are all
$O(x sup 3/2 )$ in size,
it suffices to show that we can compute
$log~zeta (s)$ and its derivatives at
$s~=~3/2+it$ to within $O(x sup -10 )$
efficiently.
Since $b sup -1~<=~| zeta (3/2+it)|~<=~b$ for some fixed
$b~>~0$, it suffices
to show that we can compute $zeta (3/2+it)$ and
its derivatives up to some fixed order
efficiently.
The Euler-Maclaurin summation formula gives
the following result for evaluating
$zeta (s)$ and its derivatives.
.sp
.I
Lemma 5. For any positive integer $m$ and any
$epsilon~>~0$ and for $1~<~a~<~2$,
there is an algorithm to compute
.R
.DS 2
.EQ
{partial sup j} over {partial t sup j} ~zeta (a+it)
.EN
.DE
.R
.I
for any $t$ with
$0~<=~t~<=~x$
and any $j$ with $0~<=~j~<=~2m$
with an error
of $O(x sup -10 )$ in $O(tx sup {epsilon /10} )$
operations using $O(x sup epsilon )$ space,
for $x~>=~x sub 0 ( epsilon )$.
The constant in the O-symbol depends on
$m$, $epsilon$, and $a$.
.R
.P
The Riemann-Siegel formula can be used to compute values
of $zeta (s)$ even faster.
It is difficult to obtain good bounds for the error in that
formula for small $t$, but one easily obtains a result valid
for $t ~>=~ x sup 1/10$.
.sp
.I
Lemma 6.
For $x sup 1/10 ~<=~ t ~<=~ x$ and $1 ~<=~ a ~<=~ 2$, the Riemann-Siegel
formula can be used to compute $zeta (a+it)$ with an error of
$O( x sup -10 )$ in $O( t sup 1/2 x sup {epsilon /10} )$
operations using $O( x sup epsilon )$ space, for $x ~>=~ x sub 0 ^( epsilon )$.
.R
.P
This bound already suffices to give the desired Algorithm
$A( b, delta )$ in the case $b ~=~ 0$.
For the case of general $b$ in $0 ~<~ b ~<=~ 1 over 4$ we need
a yet faster method to compute many values of $zeta (s)$, and
to this end we use the method of [17].
.sp
.I
Lemma 7.
If $0 ~<=~ beta ~<=~ 1/2$, $1 ~<=~ a ~<=~ 2$, and $epsilon > 0$
are given, then there is a constant $c$ and an algorithm such
that for every $x ~>=~ 1$ and $T ~>=~ x sup 1/10$ it will perform
$<=~ c T sup {1/2 + epsilon}$ operations using
$<=~ c T sup {beta + epsilon}$ bits of storage and will then
be capable of computing any value $zeta (a+it)$,
$T ~<=~ t ~<=~ T ~+~ T sup beta$ to within $+-~ x sup -10$ in
$<=~ c T sup {epsilon}$ operations using
the precomputed values.
.R
.P
We use the algorithm of Lemma 5 to compute
$h sub 3/2 ( jT over n )$ for
$jT over n~<=~x sup 1/10$ and also to compute
$h sup (j) (0)$ and $h sup (j) (T)$ for
$1 <= j~<=~2k sup 2$.
We then use either the algorithm of Lemma\ 7 to or that of Lemma 6 in the case $b=0$ to compute
$h( jT over n )$ for
$x sup 1/10~<=~jT/n~<=~T$.
By Lemma\ 3 and (3.21), this gives us
$int from 0 to T~h sub 3/2 (t)dt$ to within an error
of $O(x sup -8 )$ using $O(T sup {beta + epsilon} )$ space in the computation
and using
.DS 2
.EQ
O(x sup {1/5 + delta /10} + T sup {3/2 ^-^ beta ^+^ delta / 10} )
.EN
.DE
operations.
Combining this with the previous estimates
shows that we can compute $pi (x)$ to
within $+- 10 sup -1$ in
.DS 2
.EQ
O( "{" x sup 1/2 + y + x sup 1/5 + T sup {3/2 ^-^ beta} "}" x sup {delta /10} )
.EN
.DE
operations and $O(T sup {beta + delta} )$ space.
By (3.13), this yields a running time for Algorithm
$A( b, delta )$ of $O( x sup {(3 ^-^ 2 beta ) /(5 ^-^ 2 beta ) ^+^ 2 delta /5)}$
and space requirement of $O( x sup {2 beta /(5 ^-^ 2 beta ) ^+^ delta})$
if we take
$y~=~x sup {(3 ^-^ 2 beta ) / (5 ^-^ 2 beta )}$,
where the constant implied by the O-symbol depends on
the given $delta$.
If we now choose
$beta ~=~ 5b/(2(1+b))$, we obtain the desired estimate.
.P
Finally, we observe that these algorithms may be
combined to obtain a single Algorithm $A(b)$
which computes $pi (x)$ in
$O(x sup {(3-2b)/5+ epsilon} )$ steps using
$O(x sup {b+ epsilon} )$ space for any $epsilon~>~0$.
When given $x$, Algorithm $A(b)$ computes a value
$delta~=~delta (x)$ depending on $x$ and then
runs Algorithm $A( b, delta )$.
If the value $delta (x)$ goes to zero sufficiently
slowly as $x~->~inf$, Algorithm $A(b)$ will have
the required running time property.
Note that Algorithm $A(b)$ must compute various
quantities such as $f sub k (r)$ which are treated
as fixed overhead in Algorithm $A( b, delta )$.
An analysis which we omit of the dependence on
$delta$ of the fixed overhead and the running times
of various steps in Algorithm $A( b, delta )$ indicates
that choosing $delta (x)~=~( log^log^x) sup -1$
is adequate to obtain the stated running time
bound for Algorithm $A(b)$.
.H 1 "Proofs of Lemmas"
.I
Proof of Lemma 1.
.R
Integration by parts shows that
.DS 2
.EQ (4.1)
F sub x,y (s)~=~1 over {y sup k s(s+1)...(s+k)}~int from 0 to 1~
(x-vy) sup s+k f sub k sup (k) (v)dv~,
.EN
.DE
since $f sub k sup (r) (1)~=~f sub k sup (r) (0)~=~0$ for
$0~<=~r~<=~k-1$.
Since
.DS 2
.EQ
|f sub k sup (k) (r)|~<=~(2k) sup k~,~~~~
0~<=~r~<=~1~.
.EN
.DE
this implies (3.5).
Another $k$ integrations by parts of (4.1)
yield (3.4).
.P
To prove (3), we use (3.3) and an easily justified
exchange of the order of integration to
obtain for $a~>~1$ that
.DS 2
.EQ
1 over {2 pi i} int from {a-i inf} to {a+i inf}
F sub x,y (s)^ u sup -s^ ds~=~
int sub 0 sup 1
left [
1 over {2 ^pi^i} ~int from {a-i inf} to {a+i inf}
~1 over s~ ( x-vy over u ) sup s ds
right ] f sub k (v)dv~.
.EN
.DE
We obtain the desired result by using (2.4)
on the inner integrand.
\(sq
.sp
\f2Proof of Lemma 2\f1.\0\0We make use of the
recent algorithm of Adleman, Pomerance,
and Rumely\ [1], which tests whether an
integer $n$ is prime in
.DS 2
.EQ
O(( log~n) sup {c~log~log~log~n} )
.EN
.DE
operations (and also storage locations)
on a random access machine for some
$c~>~0$.
(This algorithm can also be implemented
on a Turing machine.)
.P
Any prime power $p sup m$ that contributes
to the sum (3.10) must have $p~<=~sqrt x$.
Therefore we initially set $S~=~0$, and
proceed to test each $n~<=~sqrt x$ for
primality in $O(x sup {epsilon /100} )$
bit operations.
If $n$ is composite, we go on to test
$n~+~1$.
If $n$ is prime, we compute successively
$n,n sup 2 ,...,n sup M$, where $M$
is the smallest integer such that $n sup M~>~x$
(note that $M~<=~1+ log sub 2~ x)$,
and then increment $S$ by the floating
point approximations to 1/2,1/3,...,1/($M-$1),
each accurate to $+- x sup -1$.
Thus all arithmetic will be performed on
$O( log~x)$-bit binary numbers.
At the end, $S$ will differ from the sum (3.10)
by at most $O(x sup -1/2 )$, since $S$ is
obtained by adding up $O(x sup 1/2 )$ terms,
each of which differs from the desired
amount by a most $x sup -1$.
This proves the first part of the lemma.
.P
To prove the second part of the lemma, we
determine, for each $n~member~[x-y,x]$, the
largest $m$ such that $n~=~p sup m$ for
some $p~member~Z sup +$.
There are few values of $m$ to test, since
$m~<=~log sub 2~ x$, and for each m
we can check whether $n$ is an $m$-th power
or not by computing $t~=~n sup 1/m$ using
$O( log~x)$-bit floating point approximations,
and then testing whether $[t-1] sup m$,
$[t] sup m$, or $[t+1] sup m$ equals $n$.
Thus this part takes $O(x sup {epsilon /100} )$
operations.
Once $p$ is determined, we test it for
primality.
If it is not prime, we go on to $n+1$.
If $p$ is prime, we compute
.DS 2
.EQ
1 over m left (
1-g left ( {x-p sup m} over y right ) right )
.EN
.DE
to within $+- x sup -2$ using
$O( log~x)$-bit floating point approximations,
where we use the fact that $g(w)$ is a polynomial
of degree $2k+1$, and thus its values are
easy to compute, so that each computation
takes $O(x sup {epsilon /100} )$ steps.
Since there are $<=~y+1$ terms, this
completes the proof of
the lemma.
\(sq
.sp
\f2Proof of Lemma 4\f1.\0\0Let $a~>~1$, and recall that
.DS 2
.EQ
h sub a,x,y (t)~=~roman Re~"{" H (t) "}"~,
.EN
.DE
where
.DS 2
.EQ
H(t)~=~H sub a,x,y (t)~=~F sub x,y (a+it)~log~zeta (a+it)
.EN
.DE
is an analytic function of $t$ for
$| roman Im (t)|~<~a-1$.
Now by Cauchy's formula, for real $t$ we have
.DS 2
.EQ
{d sup m} over {dt sup m}~H(t)~
=~m! over {2 pi i}~
int from {|z-t|= a-1 over 2}~H(z)dz over {(z-t) sup m+1}~.
.EN
.DE
Hence
.DS 2
.EQ (4.2)
left ""
left |
{d sup m} over {dt sup m} right ""~H(t) right | ~<=~
m!^ left ( a-1 over 2 right ) sup -m
left ( max from {|s-(a+it)|= a-1 over 2}~ |F sub x,y (s)| right )
left (
max from {|s-(a+it)| = a-1 over 2}~
| log~zeta (s)|
right )~.
.EN
.DE
Now
.DS 2
.EQ
| log~zeta (s) | ~=~O(1)
.EN
.DE
for $Re (s)~>=~a-(a-1)/2~=~(a+1)/2$ and by Lemma\ 1 we have
.DS 2
.EQ
max from {|s-(a+it)| <= a-1 over 2}~
|F sub x,y (s)|~<=~(2k) sup k
x sup {(a+1)/2 ^+^ k}~
y sup -k left | a+1 over 2 right | sup -k-1~.
.EN
.DE
Combining these inequalities with (4.2) proves the lemma.
\(sq
.sp
\f2Proof of Lemma 5.\f1\0\0This follows
from the Euler-Maclaurin formula
as applied to
$zeta (a+it)$ and its derivatives.
We present in detail only the analysis
for $zeta (s)$.
It is well known [6,\ p.\ 114], [21] that
for any positive integers $q$ and $n$
and for $s~=~sigma~+~it$
with $sigma~>~1$
we have
.DS 3
.EQ (4.3)
zeta (s)~= mark~sum from n=1 to q-1~n sup -s~+~1 over 2~
q sup -s~+~{q sup 1-s} over s-1~
+~sum from r=1 to n~{B sub 2r} over (2r)!~
q sup 1-s-2r~prod from j=0 to 2r-2~(s+j)~
+~R sub 2n sup * (q,s),
.EN
.DE
where
.DS 2
.EQ (4.4)
R sub 2n sup * (q,s)~=~-~
{s(s+1)...(s+2n)} over (2n+1)!~
int sub q sup inf
B bar sub 2n+1 (t)t sup -s-2n dt~,
.EN
.DE
and $B bar sub 2n (t)$ is a periodic function with
period\ 1, which equals the
Bernoulli polynomial of order $2n$ for $0~<=~t~<=~1$.
We have the estimate [6,\ p.\ 115]
.DS 2
.EQ
|R sub 2n sup * (q,s)|~=~
O left ( {|B sub 2n+2 |} over (2n+2)!~
q sup {-1- sigma - 2n}~|s+2n+1| over {sigma +2n+1}~
prod from j=0 to 2n~|s+j| right )~,
.EN
.DE
where the constant implied by the O-notation
is independent of $n$.
Since
.DS 2
.EQ
|B sub 2r |~=~2(2r)! over {(2 pi ) sup 2r}~zeta (2r)~,
.EN
.DE
the $O$-term above is, for $1~<~sigma~<=~2$,
.DS 2
.EQ
O left ( {|s+2n+2| sup 2n+2} over {q sup 2+2n} right )~.
.EN
.DE
We choose $n~=~[200^ epsilon sup -1 ]~+~1$
and $q~wig~tx sup {epsilon /20}$,
and this term becomes
$O( 2 sup 2n+2 (2n+4) sup 2n+2 x sup -10 )$.
The time needed to compute the remaining
terms is proportional to $q$ multiplied
by the time needed to compute each term
to an accuracy of $O(x sup -11 )$, which
is $O(x sup {epsilon /20 }$).
The only term whose computation time depends
on $a$ is $q sup 1-s over s-1$,
and this only for $s$ very near 1.
A similar analysis can be given for computing
${partial sup j} over {partial s sup j}~
zeta (s)$ from the formulae (4.3) and (4.4)
differentiated $j$ times,
for $1~<=~j~<=~m$.
\(sq
.sp
\f2Proof of Lemma 6\f1.
We use the Riemann-Siegel formula [21].
(This formula is also derived in [6,22], but for a restricted
range of values.)
If $s ~=~ sigma +it$, $1 ~<~ sigma ~<~ 2$,
$m ~=~ left floor (2 pi ) sup {- 1/2} t sup 1/2 right floor$,
$N ~member~ Z sup +$, and $t ~>=~ c ~N$ for a certain fixed
$c ~>~ 0$, then
.DS 3
.EQ (4.5)
zeta (s) ~mark =~ sum from n=1 to m ~n sup -s ~+~ chi (x) ~sum from
n=1 to m ^n sup s-1
.EN
.sp
.EQ
~lineup +~ (-1 ) sup m-1 ^(2 pi t ) sup (s-1)/2 exp (-it/2-i pi /8-i pi
(s-1)/2) GAMMA (1-s) ^ "{" S sub N + O ( e sup {- delta prime ^t} ) ~+~
O(( c prime N t sup -1 ) sup N/6 ) "}" ~,
.EN
.DE
where $delta prime ^,^ c prime ~>~ 0$ are fixed constants,
.DS 3
.EQ
chi (s) ~=~ 2 sup s ^pi sup s-1 sin ( {pi s} over 2 ) GAMMA (1-s) ~,
.EN
.DE
and
.DS 3
.EQ
S sub N ~=~ S sub N (s) ~=~ sum from n=0 to N-1 ~sum from {r <= n/2} ~
{n ! i sup r-n} over {r ! (n-2r) ! 2 sup n} ~( 2 over pi ) sup n/2-r
^a sub n (s) PSI sup (n-2r) ^(2 "{" t sup 1/2 ( 2 pi ) sup {-1/2} "}" ) ~,
.EN
.DE
with
.DS 3
.EQ (4.6)
PSI (u) ~=~ {cos pi ( 1 over 2 u sup 2 -u- 1 over 8 )} over {cos ~pi u} ~,
.EN
.DE
and the coefficients $a sub n ^(s)$ are given by the recurrence
.DS 3
.EQ
a sub n ^(s) ~=~ 1, ~a sub 1 ^(s) ~=~ {sigma -1} over {sqrt t} ^, ~~
a sub 2 ^(s) ~=~ {( sigma -1) ^( sigma -2)} over 2t ~,
.EN
.sp
.EQ (4.7)
(n+1) ^sqrt t ^a sub n+1 ^(s) ~=~ ( sigma -n-1) a sub n ^(s) ~+~ i ~
a sub n-2 ^(s) ~~roman for ~n ~>=~ 2 ~,
.EN
.DE
and satisfy the bound
.DS 3
.EQ
| a sub n ^(s) | ~=~ O( t sup -n/6 ) ~.
.EN
.DE
(In the definition of $S sub N , ~"{" u "}"$ denotes the fractional
part of $u$.)
.P
We take $N ~=~ 600$, say.
Since [8, Eq. (6.1.45)]
.DS 3
.EQ
lim from {| y | ^->^ inf} ~(2 pi ) sup -1/2 | GAMMA (x+iy) | e sup {1/2 pi | y |}
^| y | sup {1/2 -x} ~=~ 1
.EN
.DE
we have
.DS 3
.EQ
| ( 2 pi t ) sup (s-1)/2 ^exp (-it/2-i pi /8-i pi (s-1)/2) GAMMA (1-s)
| ~mark =~ (2 pi t ) sup {( sigma -1)/2} e sup {pi t/2} | GAMMA (1-s) |
.EN
.sp
.EQ
~lineup wig (2 pi ) sup {sigma /2} t sup {- sigma /2} ~~roman as~~ t ~->~ inf
.EN
.DE
Hence when $| t | ~>=~ x sup 1/10$ and $1 ~<=~ sigma ~<~ 2$,
the formula (4.5) gives
.DS 3
.EQ (4.8)
zeta (s) ~mark =~ sum from n=1 to m ^n sup -s ~+~ chi (x) ~sum from n=1 to m
^n sup s-1
.EN
.sp
.EQ
~lineup +~ (-1 ) sup m-1 ^(2 pi t ) sup (s-1)/2 ^exp left ( -~ it over 2
~-~ {i pi} over 8 ~-~ {i pi (s-1)} over 2 right ) GAMMA (1-s) S sub N ^(s)
~+~ O( x sup -10 ) ~.
.EN
.DE
.P
We first observe that for any fixed $N$ and $x sup 1/10 ~<=~ t ~<=~ x$,
the term containing $S sub N ^(s)$ can be evaluated to accuracy
$O( x sup -10 )$ in $O( x sup epsilon )$ operations.
To see this, we first note that by Stirling's formula [8; Eq. 6.1.42],
.DS 3
.EQ
log ~GAMMA (1-s) ~=~ (1/2 -s) ~ mark log ~(1-s) ~-~ (1-s) ~+~ 1 over 2
log (2 pi )
.EN
.sp
.EQ (4.9)
~ lineup + ~ sum from h=1 to M ~{B sub 2h} over {2h ^(2h -1)} ~
(1-s ) sup 1-2h ~+~ O( t sup -1-2M ) ^,
.EN
.DE
where we choose $M ~=~ 20$ and the constant implied by the $O$-notation
is (as will be the case for all other such constants) dependent
only on $sigma$ and $delta$.
Hence we can compute $log ^GAMMA (1-s)$ using (4.9) to an accuracy
of $O( x sup -20 )$ in $O( x sup epsilon )$ operations.
Similarly, we can compute
.DS 3
.EQ
w ~=~ -i pi (s-1)/2 ~-~ it/2 ~-~ i pi /8
.EN
.DE
to within $O( x sup -20 )$ in $O( x sup epsilon )$ operations.
Since [8; Eq. 6.1.45]
.DS 3
.EQ
| GAMMA (1-s) | ~=~ O( t sup {1/2 - sigma} ~e sup {- pi t/2} )
.EN
.DE
(a result that follows from (4.8) also), we have
.DS 3
.EQ
| ( 2 pi t ) sup (s-1)/2 ~e sup w ~GAMMA (1-s) | ~=~ O( t sup {- sigma /2} ) ^,
.EN
.DE
and we can compute it to an accuracy of $O( x sup -20 )$ in
$O( x sup epsilon )$ operations.
.P
Next we consider $S sub N ^(s)$, where we selected $N ~=~ 600$.
Note that $PSI (u)$ is actually an entire function, and its
Taylor series expansion around $u ~=~ 1/2$ can be used to
evaluate $PSI (u)$ and its first $N$ derivatives in a neighborhood
of $u ~=~ 1/2$ (and similarly in a neighborhood of $u ~=~ -1/2$),
while away from $u ~=~ +- ^1/2$, we can differentiate the expansion
(4.7) and evaluate individual terms.
Each of the first $N$ derivatives takes $O( x sup epsilon )$
operations to compute to within $O( x sup -20 )$.
Finally the $a sub n ^(s)$ can be computed using the recurrence
(4.7) to within $O( x sup -20 )$ in $O( x sup epsilon )$
operations.
Combining all these estimates, we conclude that the entire
term in (4.8) containing $S sub N ^(s)$ can be computed to
within $O( x sup -20 )$ in $O( x sup epsilon )$
operations, using $O( x sup epsilon )$ storage.
.P
Stirling's formula (4.9) allows us to compute $chi (s)$ to
accuracy $O( x sup -12 )$ in a similar way.
For $x sup 1/10 ~<=~ t ~<=~ x$ the two sums in (4.8) are
evaluated to accuracy $O( x sup -10 )$ by evaluation term-by-term
in time $O( m x sup epsilon ) ~=~ O( t sup {1/2 + epsilon} )$.
Lemma 6 follows.
\(sq
.bp
.ce
\f3REFERENCES\f1
.sp 3
.RL
.LI
L. M. Adleman, C. Pomerance, and R. S. Rumely,
On distinguishing prime numbers from
composite numbers, \f2Annals Math.\f1 117 (1983), 173-206.
.LI
J. Bohman, On the number of primes less
than a given limit, \f2BIT 12\f1 (1972), 576-577.
.LI
H. Davenport, \f2Multiplicative Number Theory\f1,
2nd ed. (revised by H.\ L. Montgomery),
Springer-Verlag, 1980.
.LI
M. Deuring, Asymptotische Entwicklungen der
Dirichletschen L-Reihen, \f2Math. Annalen 168\f1
(1967), 1-30.
.LI
L. E. Dickson, \f2History of the Theory of Numbers\f1,
Chelsea reprint. Vol. 1, Chapter XVIII.
.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
\f2Handbook of Mathematical Functions\f1, M. Abramowitz
and I.\ A. Stegun, eds., National Bureau
of Standards, 9th printing, 1970.
.LI
G. H. Hardy, \f2Ramanujan\f1, Cambridge Univ. Press,
1940.
(Reprinted by Chelsea.)
.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
D. H. Lehmer, On the exact number of primes
less than a given limit, \f2Illinois J. Math. 3\f1
(1959), 381-388.
.LI
B. F. Logan, Note on ideal sharp-cutoff filters,
preprint.
.LI
D. C. Mapes, Fast method of computing the
number of primes less than a given
limit, \f2Math. Comp. 17\f1 (1963), 179-185.
.LI
E. D. F. Meissel, ${roman U} dotdot$ber die Bestimmung der
Primzahlmenge innerhalb gegebener Grenzen,
\f2Math. Annalen, 2\f1 (1870), 636-642.
.LI
E. D. F. Meissel, Berechrung der Menge der
Primzahlen, welche innerhalb der ersten
Milliarde na\*:turlichen Zahlen vorkommen,
\f2Math. Annalen, 25\f1 (1885), 251-257.
.LI
A. M. Odlyzko and A. Scho\*:nhage,
Fast algorithms for multiple evaluations of the Riemann
zeta function, to be published.
.LI
P. Pritchard, A sublinear additive sieve
for finding prime numbers, \f2Comm. ACM 24\f1
(1981), 18-23.
.LI
S. Ramanujan, \f2Notebooks\f1, Tata Institute,
1957, 3rd notebook, p. 371.
.LI
H. Riesel and G. Go\*:hl, Some calculations related
to Riemann's prime number formula,
\f2Math. Comp. 24\f1 (1970), 969-983.
.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
E. C. Titchmarsh, \f2The Theory of the Riemann
Zeta-Function\f1, Oxford Univ. Press, 1951.
.LI
H. S. Wilf, What is an answer?
\f2Am. Math. Monthly, 89\f1 (1982), 289-292.
.LE
.ls 1
.CS