                                                         AN ANALYTIC METHOD FOR BOUNDING ψ(x).


                                                                                       JAN BÜTHE


                                                   Abstract. In this paper we present an analytic algorithm √    which calculates almost




arXiv:1511.02032v2 [math.NT] 22 Oct 2017
                                                   sharp bounds for the normalized remainder term (t − ψ(t))/ t for t ≤ x in expected
                                                   run time O(x1/2+ε ) for every ε > 0. The method has been implemented and used
                                                   to calculate such bounds for t ≤ 1019 . In particular, these imply that li(x) − π(x) is
                                                   positive for 2 ≤ x ≤ 1019 .




                                                                1. Introduction and statement of results
                                             This paper concerns the problem of calculating limited range approximations to the
                                           Chebyshov function                X
                                                                     ψ(x) =       log(p) = x + o(x).
                                                                                     pm ≤x
                                           More precisely, we are interested in calculating almost sharp bounds for the normalized
                                           error term
                                                                                         t − ψ(t)
                                           (1.1)                                Rψ (t) =    √
                                                                                              t
                                           in the prime number theorem for ψ(t) in finite intervals [x, Lx]. So far, such calculations
                                           seem to have been based on tabulating prime numbers (see e.g. [12, 13, 14]), whereas
                                           bounds for unlimited ranges are usually derived analytically (see e.g. [3, 4, 13, 15, 14, 16]).
                                           The elementary approach leads to a run time of Õ(x) for fixed L and x → ∞, where
                                           f (x) = Õ(g(x)) means there exists an A such that f (x) = O(g(x) log(x)A ). In this paper
                                           we present an analytic algorithm for this task, which satisfies the assertion of the following
                                           theorem.
                                           Theorem 1. For every triple (L, δ, θ) ∈ (1, ∞) × (0, ∞) × (0, 1/2], there exist effectively
                                           computable constants C1 = C1 (L, θ, δ) and C2 = C2 (L, θ, δ) and an algorithm √which takes
                                           x ≥ 2 and the zeros ρ of the Riemann zeta function with 0 < ℑ(ρ) ≤ C1 xθ log x with
                                                                                                    √
                                           an accuracy of x−C2 as input, performs Õ([1 + Ne (C1 xθ log x)]xθ ) arithmetic operations
                                           on Õ(xθ ) variables of size C2 log x, where Ne (T ) denotes the number of zeros with 0 <
                                           ℑ(ρ) ≤ T violating the Riemann Hypothesis, and outputs numbers ML+ (x) and ML− (x)
                                           satisfying
                                                                    t − ψ(t)                         t − ψ(t)
                                           (1.2)           sup ± √           ≤ ±ML± (x) ≤ sup ± √             + δx1/2−θ ,
                                                          x≤t≤Lx         t                 x≤t≤Lx         t
                                           where pluses and minuses are to be taken correspondingly.
                                             If L is sufficiently large the Riemann Hypothesis (RH) implies ([9, 15.2.1, Ex. 2b)])
                                                                                                 t − ψ(t)
                                           (1.3)                             lim inf sup ±          √     > 0.
                                                                              x→∞ x≤t≤Lx              t
                                           Furthermore, the expected run time for calculating zeros of the zeta function with imag-
                                           inary part up to T within an accuracy of T −c for any c > 0 is O(T 1+ε ) for every ε > 0,

                                             Date: October 24, 2017.
                                             2010 Mathematics Subject Classification. Primary 11N05, Secondary 11M26.
                                                                                             1
2                                       JAN BÜTHE


assuming RH and simplicity of the zeros (see [10]). So if we take θ = 1/2 and δ sufficiently
small, the algorithm calculates almost sharp bounds for Rψ (x) in [x, Lx] in expected run
time O(x1/2+ε ) for every ε > 0.
   The algorithm has been implemented and used to calculate analytic bounds for x ≤
1019 , using the zeros with imaginary part up to 1011 , whose calculation has been reported
in [5]. The calculated bounds also give rise to improved bounds for the functions
                      X              X                            X π(x1/k )
(1.4)         π(x) =      1, ϑ(x) =      log(p), and π ∗ (x) =                .
                                                                          k
                    p≤x              p≤x                         k≥1

The numerical results are summarized in the following theorem.
Theorem 2. The following estimates hold:
                             √
(1.5)       |x − ψ(x)| ≤ 0.94 x                                    for 11 < x ≤ 1019 ,
                             √
(1.6)         x − ϑ(x) ≤ 1.95 x                                  for 1423 ≤ x ≤ 1019 ,
                             √
(1.7)         x − ϑ(x) > 0.05 x                                      for 1 ≤ x ≤ 1019 ,
                          √
                 ∗          x
(1.8)  |li(x) − π (x)| <                                             for 2 ≤ x ≤ 1019 ,
                         log x
                           √ 
                             x          3.9    19.5 
(1.9)     li(x) − π(x) ≤        1.95 +      +                        for 2 ≤ x ≤ 1019 ,
                         log(x)        log x log(x)2
(1.10)      li(x) − π(x) > 0                                         for 2 ≤ x ≤ 1019 .
    In particular, this gives a new lower bound for the Skewes number, the number xs ∈
[2, ∞) where the first sign change of li(x) − π(x) occurs. The last published lower bound
appears to be xs ≥ 1.2 × 1017 in [12]. Furthermore, in the earlier paper [17] the second
author claims to have verified xs ≥ 1018 but no further explanation is given. In total, the
calculations took about 1, 200 hours on a 2.27 GHz Intel Xeon X7560 CPU.

                                   Acknowledgment
   The author wishes to thank Jens Franke, to whom he owes the idea of applying methods
for multiple evaluations of trigonometric sums to this problem and who suggested this
topic as part of the author’s PhD thesis.
   Furthermore, he wishes to thank the anonymous referees for their comments. He is
particularly grateful to the referee who spotted a mistake in the proof of Lemma 1, which
has been corrected in the final version.

                                      2. Notations
   In addition to the usual Landau O and Vinogradov ≪ notation, we frequently use
Turing’s big theta notation: g(t) = Θ(f (t)) for t ∈ U ⇔ g(t) ≤ |f (t)| for t ∈ U . Fur-
thermore, the notation f (t) ≍ g(t) is used for f (t) ≪ g(t) and g(t) ≪ f (t). Finally,
f± (t) := limhց0 f (t ± h) denotes the limit from the right, respectively left.

                            3. Description of the method
   The basic idea of the method presented in this paper is to use an explicit formula to
bound ψ(t) at sufficiently many well-distributed points in [x, Lx] which are then extended
to the whole interval by interpolation.
   To illustrate the first task recall the well-known approximate version of the von Man-
goldt explicit formula
                            X xρ                 1                 x          
(3.1)       ψ(x) = x −                − log 2π − log(1 − x−2 ) + O     log(x)2
                                   ρ             2                  T
                        |ℑ(ρ)|<T
                               BOUNDING ψ(x) ANALYTICALLY                                                3


for T ≪ x, where the sum is taken over the non-trivial zeros of the Riemann zeta function.
Using this formula, approximating Rψ (t) within an accuracy of O(x1/2−θ ) can be done
by calculating the contribution of zeros with imaginary part up to T = Cxθ log(x)2 .
Extending these bounds to [x, Lx] with an error of size O(x1/2−θ ) can be achieved by
calculating approximations at O(xθ ) well-distributed points in [x, Lx], since if ξ ∈ [x, Lx]
and 0 < y ≪ x1−θ , then
                            ξ + y − (1 + O(x−θ ))ξ   ψ(ξ + y) − (1 + O(x−θ ))ψ(ξ)
    Rψ (ξ + y) − Rψ (ξ) =          √               −           √
                                      ξ+y                         ξ+y
                        = O(x1/2−θ ) + O(x−θ−1/2 )ψ(ξ) + O(x−1/2 )(ψ(ξ + y) − ψ(ξ))
                        = O(x1/2−θ )
by the prime number theorem and the Brun-Titchmarsh inequality.
   If the approximations are caltulated by directly evaluating the sum over zeros in (3.1)
this leads to a run time of Õ(x2θ ), which outperforms the naive method only for θ < 1/2.
This can be improved by using techniques for multiple evaluations of trigonometric sums,
such as the Odlyzko-Schönhage algorithm. These allow one to evaluate the contribution
of zeros on the critical line to the sum in (3.1) on geometric progressions of length T in
run time Õ(T ) = Õ(xθ ), reducing the run time of the algorithm to Õ(xθ ).
   In principle one could use an explicit version of (3.1), but we rather use a continuous
approximation to ψ(x) from [2] for which a similar explicit
                                                          √ formula exists. This decreases
the truncation bound in the sum over zeros to T = C ′ xθ log x, saving a factor log(x)3/2 .
Also, we rather use a simpler FFT method from [5] in place of the Odlyzko-Schönhage
algorithm for multiple evaluation of trigonometric sums.

3.1. Bounding ψ(x) analytically. Let
                                      1 X              1 X
                           ψ0 (x) =           log(p) +     log(p)
                                      2 pm <x          2 m
                                                           p ≤x

denote the normalized Chebyshov function. We intend to bound ψ(x) in terms of the
modified Chebyshov function
                                            X    1
(3.2)                ψc,ε (x) = ψ0 (x) +           Mx,c,ε (pm )
                                             m ε
                                                 m
                                         −ε  e   x<p <e x

introduced in [2], where
                                           Z log(t/x)
                    log t h ∗
(3.3) Mx,c,ε(t) =          χ         (t)            ηc,ε (τ )e−τ /2 dτ
                    λc,ε [x,exp(ε)x]        −ε
                                                                      Z ε                           i
                                                 − χ∗[exp(−ε)x,x] (t)            ηc,ε (τ )e−τ /2 dτ, ,
                                                                      log(t/x)

χ∗A denoting the normalized characteristic function which takes the value 1/2 on the
boundary of A,
                                              c         p
(3.4)                        ηc,ε (τ ) =           I0 (c 1 − (τ /ε)2 ),
                                         ε sinh(c)
         P∞         2n       2
I0 (y) =   n=0 (y/2) /(n!) denoting the 0-th modified Bessel function of the first kind,
            R ε         τ /2
and λc,ε = −ε ηc,ε (τ )e dτ .
    The function ψc,ε (x) is a continuous approximation to ψ(x) and we review some of its
properties. The first result provides bounds for ψ(x) in terms of ψc,ε .
4                                              JAN BÜTHE


Proposition 1 ([2, Proposition 4]). Let
                                   Rt
                                  − −∞ ηc,1 (τ ) dτ
                                                                       t<0
                         µc (t) = −µc (−t)                              t>0
                                  
                                    0                                   t=0
                                  

and let                                             Z t
                                         νc (t) =         µc (τ ) dτ.
                                                     −∞
Furthermore, let 0 ≤ α < 1, x > 100, and let 0 < ε < 10−2 , such that
                                              εxe−ε |νc (α)|
                                        B=                   >1
                                               2(µc )+ (α)
holds. We define
                                      h 2ε x |ν (α)|
                                               c            √    1            i
         A(x, c, ε, α) = e2ε log(eε x)               + 2.01ε x + log log(2x2 ) .
                                           log B                 2
Then we have
                           ψ(e−αε x) ≤ ψc,ε (x) + A(x, c, ε, α),
and
                            ψ(eαε x) ≥ ψc,ε (x) − A(x, c, ε, α).
   The modified Chebyshov function satisfies a similar explicit formula as ψ(x) but the
sum over zeros converges absolutely and is therefore more accessible to numerical calcu-
lations using a subset of the zeros of ζ(s). For the purpose of this paper, the following
approximate version will suffice.
Proposition 2. Let x ≥ 10, 0 < ε ≤ 10−4 and let
                                              p
                                       c sinh( c2 − (εt)2 )
(3.5)                    ℓc,ε (t) =         p
                                    sinh c    c2 − (εt)2
denote Logan’s function [7]. Then we have
                                              1    X ℓc,ε ( ρ − 1 )
(3.6)                  x − ψc,ε (x) =                       i   2i
                                                                    xρ + Θ(2).
                                         ℓc,ε (i/2) ρ        ρ

Proof. This is a corollary of [2, Proposition 2]: since ℓc,ε (−z) = ℓc,ε (z) we get
        X ℓc,ε ( ρ − 1 )                    ρ   1
                                 1 X ℓc,ε ( i − 2i )   ℓc,ε (i/2) X
                i    2i
                             =                       ≤              ℑ(ρ)−2 < 0.025 ℓc,ε(i/2)
         ρ
                 ρ               2 ρ   ρ(1 − ρ)             2     ρ

using the bijection ρ 7→ 1 − ρ of non-trivial zeros and [13, Lemma 17]. Furthermore, we
have γ/2 + 1 + log(π)/2 ≤ 1.87, − log(1 − x−2 )/2 ≤ 0.006 and 8ε|log ε| ≤ 0.008, so the
assertion follows.                                                                    
    We have the following tail bounds for truncating the sum over zeros.
Proposition 3 ([2, Proposition 3]). Let x > 1, 0 < ε ≤ 10−3 and c ≥ 3. Then we have
             X      ℓc,ε ( ρi − 2i
                                1
                                   ) xρ         x + 1 0.71√cε             c
(3.7)                                   ≤ 0.16         e      log(3c) log    .
                  c
                       ℓc,ε (i/2) ρ            sinh(c)                     ε
                |ℑ(ρ)|> ε

   Furthermore, if a ∈ (0, 1) such that a εc ≥ 103 holds, and if the Riemann Hypothesis
                                                     c
holds for all zeros with imaginary part in ( ac  ε , ε ], then we have
                X        ℓc,ε ( ρi − 2i
                                     1
                                        ) xρ   1 + 11cε         c  cosh(c√1 − a2 ) √
(3.8)                                        ≤             log                        x.
            ac        c
                            ℓc,ε (i/2) ρ         πca2            ε       sinh(c)
              ε <|ℑ(ρ)|≤ ε
                                 BOUNDING ψ(x) ANALYTICALLY                                5


Remark 1. It should be demonstrated that it is indeed more efficient to approximate
ψ(x) this way. Calculating ψ(x) within an accuracy of O(xδ ) via the modified Chebyshov
function can be done by choosing
              ε = xδ−1 log(x)1/2         and         c = (1 − δ) log(x) + 2 log log(x).
Since
                                                1
                                          |νc (0)| ∼ √
                                                2πc
for c → ∞ (see [2, Proposition 5]), Proposition 1 gives
                                        ψ(x) − ψc,ε (x) ≪ xδ
and from Propositions 2 and 3 we get
                                         1         X        ℓc,ε ( ρi − 2i
                                                                        1
                                                                           ) ρ
                  ψc,ε (x) − x =                                            x + O(xδ )
                                    ℓc,ε (i/2)                       ρ
                                                 |ℑ(ρ)|<T
                          √
with T = c/ε ∼ (1 − δ)x1−δ log x. If the same zeros are used in the von Mangoldt explicit
formula, the standard estimate (3.1) gives an error term which is larger by a factor of
size ≫ log(x)3/2 .
3.2. Interpolating bounds for ψ(x). Next, we give an estimate for the interpolation
error. For simplicity, we assume that t−ψ(t) changes sign in [x, Lx], or to be more precise:
we assume the upper, respectively, lower bound for Rψ (t) to be positive, respectively
negative. This is implied by RH if L is sufficiently large and has been the case in all
practical applications.
Proposition 4. Let 109 ≤ a < b and let
                                    a = x0 < x1 < · · · < xn = b
be a dissection of [a, b], whose maximal step size
                             ∆ = max{xk − xk−1 | k = 1, . . . , n}
                       −5
satisfies 10 ≤ ∆ ≤ 10 a. Then the following assertions hold:
    (1) Let M > 0 satisfy
                                          xk − ψ(xk )
                                             √        ≤M
                                               xk
         for k = 0, 1, . . . , n. Then
                                                                      
                    y − ψ(y)                    log a   ∆             2
                        √         ≤ 1.001 M + √            + log log(a )
                           y                       a log ∆
        holds for all y ∈ [a, b].
    (2) Let m < 0 satisfy
                                          xk − ψ(xk )
                                             √        ≥m
                                               xk
        for k = 0, 1, . . . , n. Then
                                                                  
                    y − ψ(y)              log a   ∆              2
                       √     ≥ 1.001  m −  √         + log log(a   )
                         y                   a log ∆
        holds for all y ∈ [a, b].
Proof. We start by proving
                                                                    
                                                    ∆              2
(3.9)             ψ(x) − ψ(x − y) ≤ log(x) 1.0001       + log log(x )
                                                  log ∆
6                                                  JAN BÜTHE

                                           ∆          10
for x ≥ a ≥ 109 and 0 ≤ y ≤ ∆  2 . Since log ∆ ≥ log 10 > 4 we may assume y ≥ 3. The
Brun-Titchmarsh inequality, as stated in [8], and the trivial estimate
                                                       Y
                       #{p | pm ∈ [X − Y, X]} ≤ 2 X 1/m + 1,
                                                       m
which holds for 0 < 2Y < X, yield
                                        X        1
           ψ(x) − ψ(x − y) ≤ log(x)
                                                 m
                                     x−y≤pm ≤x
                                                                           
                                                ⌊2 log x⌋
                                        2y        X          y          1 
                                                                          
                           ≤ log(x)         +             2 2 x1/m−1 +       .
                                       log y      m=2
                                                            m           m

Since
                                 ⌊2 log x⌋         Z 2 log x
                                   X         1                 dt
                                               ≤                  ≤ log log(x2 )
                                  m=2
                                             m      1           t
and
             ⌊2 log x⌋                            Z ∞
               X             y 1/m−1    y     2y      dt
                         2      x    ≤  √  +
              m=2
                             m2        2 x x2/3 2 t2
                                        y 1               y           y
                                     ≤ √     + x−1/6 < 0.6 √ ≤ 0.0002       ,
                                         x 2                 x        log y
this implies (3.9) since y 7→ logy y increases monotonically for y > e.
   Now let x ∈ {xk }nk=1 and let 0 ≤ y ≤ ∆/2. Then we have
             x − y − ψ(x − y)   x − ψ(x) ψ(x) − ψ(x − y)           y
(3.10)            √           = √         +      √           −√       .
                   x−y             x−y             x−y            x−y
Now if m < 0 and M > 0 satisfy the conditions in the theorem, then
                          √                        √
                            x        x − ψ(x)         x
              1.001m ≤ √        m≤ √           ≤ √       M ≤ 1.001M.
                          x−y          x−y         x−y
Furthermore (3.9) gives
                                                                           
                    ψ(x) − ψ(x − y)     log x             ∆
              0≤        √           ≤ √          1.0001       + log log(x2 )
                          x−y            x−y            log ∆
                                                                       
                                            log a    ∆
                                    ≤ 1.001 √             + log log(a2 ) .
                                               a log ∆
Since
                             y              ∆         log(a)∆
                       0≤ √        ≤ 1.0001 √ ≤ 1.001 √
                            x−y              a          a log ∆
the bound (3.10) yields
                                       x − y − ψ(x − y)
             log a  ∆
 1.001 m − √             + log log(a2 )   ≤      √
                a log ∆                           x−y
                                                 
                                                       log a  ∆                
                                                                              2
                                          ≤ 1.001 M + √            + log log(a ) .
                                                           a log ∆
The estimates
                                        x + y − ψ(x + y)
              log a  ∆
  1.001 m − √            + log log(a2 )   ≤      √
                 a log ∆                           x+y
                                                                                
                                                        log a  ∆
                                          ≤ 1.001 M + √            + log log(a2 )
                                                           a log ∆
             n−1
for x ∈ {xk }k=0 and 0 ≤ y ≤ ∆/2 are proven in an analogous way.                   
                                    BOUNDING ψ(x) ANALYTICALLY                                  7


3.3. Fundamental Theorem. We can now state the fundamental theorem for the an-
alytic method, which reduces the problem of bounding Rψ (t) on [x, Lx] to efficiently
approximating ψc,ε (t) at finitely many points. This is then dealt with in the next section.
Theorem 3. Let 0 < ε < 10−4 and let 0 ≤ α ≤ 1 satisfy
                                             ε x νc (α)
                                                        > 10.
                                            2(µc )+ (α)
Furthermore, let eαε 109 ≤ a < b and let
                                     a = x0 < x1 < · · · < xn = b
be a dissection of [a, b] whose maximal step size
                               ∆ = max{|xk − xk−1 | | k = 1, . . . , n}
                         −5
satisfies 10 ≤ ∆ ≤ 10 a.
   We define the error terms
                         √
           E1 = 1.001 α ε b,
                             √                                     log log(2a2 )
                                             ε b |ν (α)| −1                   
                                                     c
           E2 = 2.02 log(b) ε b |νc (α)| log                   +ε+       √
                                              2(µc )+ (α)              4 a
and
                                                      
                       log a         ∆
             E3 = 1.001 √                + log log(a2 ) .
                          a        log ∆
  Then the following assertions hold:
   (1) Let M > 0 satisfy
                                     xk − ψc,ε (xk )
                                         √           ≤M
                                           xk
       for k = 0, 1, . . . , n. Then
                               y − ψ(y)
                                  √     ≤ 1.01 (M + E1 + E2 + E3 )
                                   y
          holds for all y ∈ [eαε a, b].
      (2) Let m < 0 satisfy
                                           xk − ψc,ε (xk )
                                               √           ≥m
                                                 xk
          for k = 0, 1, . . . , n. Then
                                y − ψ(y)
                                   √     ≥ 1.01 (m − E1 − E2 − E3 )
                                    y
          holds for all y ∈ [a, e−αε b].
Proof. We start with the proof of the first assertion concerning the upper bound. Let
x̃ = eαε x. Then Proposition 1 yields
        x̃k − ψ(x̃k )   x̃k − xk   xk − ψc,ε (xk ) A(xk , c, ε, α)
             √        ≤ √        +     √          +     √
              x̃             x̃k         x̃k              x̃k
                                    √             x k − ψc,ε (xk )          A(xk , c, ε, α)
                      = 2 sinh(αε/2) xk + e−αε/2       √           + e−αε/2    √            .
                                                         xk                        xk
Under the suppositions of the proposition we have
                                            √
                                2 sinh(αε/2) xk ≤ E1 ,
8                                           JAN BÜTHE


and
                                      A(xk , c, ε, α)
                                         √            ≤ E2 .
                                             xk
Therefore,
                                   x̃k − ψ(x̃k )
                                       √         ≤ M + E1 + E2
                                         x̃k
for k = 0, . . . , n and Proposition 4 yields the desired estimate
                                                                                  !
               y − ψ(y)                                       ˜
                                                    log(ã)  ∆                 
                                                                              2
                    √     ≤ 1.001 M + E1 + E2 + √                 + log log(ã )
                     y                                 ã log ∆ ˜
                        ≤ 1.01 (M + E1 + E2 + E3 )
for all y ∈ [ã, b̃].
   The lower bound estimate follows in an analogous way by using
        x̂k − ψ(x̂k )                √           xk − ψc,ε (xk )         A(xk , c, ε, α)
            √         ≥ −2 sinh(αε/2) xk + eαε/2      √          − eαε/2    √            ,
              x̂k                                      x                        xk
where x̂k = e−αε xk .                                                                        

                                   4. Evaluation of ψc,ε
   We intend to evaluate the sum over zeros,
                                X∗ ℓc,ε (ρ/i − 1/2i)
(4.1)                                                xρ ,
                                             ρ
                                 |ℑ(ρ)|<T

in the explicit formula for ψc,ε for many values of x ∈ [x0 , Lx0 ]. If we take y = log(x),
denote non-trivial zeros of ζ(s) by ρ = β + iγ with β, γ ∈ R, normalize with the factor
e−y/2 , and remove possible violations of the RH, we encounter a trigonometric sum
                                              X
(4.2)                             Fψ,T (y) =      aρ eiyγ ,
                                                |γ|<T
                                                β=1/2

where
                                        ℓc,ε (ρ/i − 1/2i)
                                    aρ =                  .
                                                 ρ
   Such trigonometric sums can be evaluated efficiently on equidistant grids using the
Fast Fourier Transform (FFT). In this case we can calculate O(T ) values of Fψ,T (y)
using Õ(T ) arithmetic operations on variables of size O(log(T )) (see [10]). Furthermore,
the Fourier transform of F is supported on [−T, T ] so that F (y) can be recovered from
samples F (nπ/β) for some β > T by bandlimited function interpolation, where a single
evaluation can be done in Õ(1) (see [11]).
4.1. Multiple evaluations of trigonometric sums. Let
                                                N
                                                X
(4.3)                                 F (y) =         aj eiγj y
                                                j=1

with γj ∈ R and aj ∈ C. The first author of [5] proposed a simple method, based on the
FFT to evaluate F (y) simultaneously at integer values y ∈ [−Y, Y ] ∩ Z. The method is
similar to the Odlyzko-Schönhage algorithm [10].
   We briefly restate the algorithm and analyze the run time for the application in mind.
The algorithm is based on rounding eiγj onto the next Rth root of unity, where R = 2r
is a power of 2. Let nj ∈ Z such that
                                           2πnj     π
                                δj := γj −      =Θ       .
                                            R         R
                                    BOUNDING ψ(x) ANALYTICALLY                                                    9


Furthermore, let
                                          P (t) = b0 + · · · + bn tn
be a polynomial approximating f (t) = exp(it πY
                                              R ) in [−1, 1]. Then we have
                     N
                     X
           F (y) =         aj e2πinj y/R eiyδj
                     j=1
                     N
                     X                                                                         N
                                                                                                X             
                 =         aj e2πinj y/R P (yδj πY
                                                 R
                                                   )+Θ             f − P ; C 0 ([−1, 1])              |aj |
                     j=1                                                                        j=1

for y ∈ [−Y, Y ], where
                                      g; C 0 ([a, b]) := sup |g(t)|
                                                              t∈[a,b]

denotes the supremum norm on [a, b]. Now let
                               N               Rδ ℓ                            R
                                                 j
                                                        and fˆℓ (y) =
                               X                                                 X
              fℓ (k) =                   aj                                            fℓ (k)e2πiky/R .
                               j=1
                                               πY
                                                                                 k=1
                           nj ≡k mod R

Then we have
                     N
                     X                                      n
                                                            X              N
                                                                           X
                           aj e2πinj y/R P (xδj πY
                                                 R
                                                   )=             bℓ y ℓ         aj δjℓ e2πinj y/R
                     j=1                                    ℓ=1            j=1
                                                            n
                                                                  bℓ fˆℓ (y)y ℓ
                                                            X
                                                        =
                                                            l=1

and all values of fˆℓ on Z/RZ may be calculated appealing to the FFT. For the polyno-
mial P we choose the polynomial Pn of degree ≤ n which interpolates f (t) at the zeros
     2k−1
cos( 2n+2 π), k = 1, 2, . . . , n + 1, of the n + 1th Chebyshov polynomial. The standard error
estimate for polynomial interpolation then gives the bound
                                                        πY n+1 √8
                                          0
(4.4)                        f − Pn ; C ([−1, 1]) ≤                      .
                                                         2R     (n + 1)!
We then get the following result.
Proposition 5. Assume in (4.3) that there exist constants B, C ≥ 0 such that |aj | ≤ Bj C
for all j, and that γj ∈ [0, 2π), and let D, α > 0. Then there exists an N0 (α, B, C, D) such
that for all N, Y ∈ N satisfying N > N0 and log(N )−D ≤ N/Y ≤ log(N )D the algorithm
above takes each aj , j = 1, . . . , N, with an accuracy of N −2α−3 and each γj , j = 1, . . . , N,
with an accuracy of N −2α−C−4 as input and calculates F (y) for all y ∈ [−Y, Y ] ∩ Z within
an accuracy of N −α performing Õ(N ) arithmetic operations on Õ(N ) variables of size
O(log N ), where the implied constants depend on α, B, C and D only.
Proof. Let R = 2r denote the power of two which is closest to N . Then, in view of (4.4),
we have
                                         bℓ fˆℓ (y)y ℓ + O(N −2α )
                                  X
(4.5)                   F (y) =
                                         1≤n≤log N

for every α > 0. It is easily seen from the discrete orthogonality of the Chebyshov
polynomials Tk that bℓ ≪ 3n ≪ N 2 , since the coefficients of Tk are bounded by 3k in
absolute value. Furthermore, we have fˆℓ (y) ≪ N C+1 /Y ℓ and trivially y ℓ ≪ Y ℓ . It
therefore suffices to calculate bℓ within an accuracy of O(N −2α−C−1 ), fˆℓ (y) within an
accuracy of O(Y −ℓ N −2α−2 ), and y ℓ within an accuracy of O(Y ℓ N −2α−C−3 ) in order to
calculate F (y) within an accuracy of O(N −2α ), which can all be carried out on variables of
10                                               JAN BÜTHE


size O(log N ). The calculation of bℓ takes Õ(1) arithmetic operations. For the calculation
of fˆℓ it suffices if the input variables aj and γj are given within an accuracy of N −2α−3 ,
respectively N −2α−C−4 and all values fˆℓ (y) are calculated via FFT performing Õ(N )
arithmetic operations on Õ(N ) variables. Calculating F (y) for a single value of y then
takes Õ(1) arithmetic operations, so the assertion follows.                               

4.2. Bandlimited function interpolation. The method outlined in the preceding sec-
tion is sufficient to obtain an algorithm satisfying Theorem 1. But for practical ap-
plications it can be necessary to reduce the memory requirement of the algorithm by
sub-dividing the sum over zeros. Then the number of evaluations is much larger than the
number of summands in the trigonometric sum and it is favorable to calculate sufficiently
many samples of the trigonometric sum to obtain intermediate values by bandlimited
function interpolation instead of repeatedly applying the method from the previous sec-
tion.
   We recall the interpolation formula from [11] which is a modification of the well-
known Shannon-Nyquist-Whittacker interpolation formula and give an explicit estimate
for truncating the infinite sum.

Proposition 6. Let
                                                      N
                                                      X
                                           F (y) =          aj eiγj y ,
                                                      j=1

where γj ∈ R and let τ = maxj {|γj |}. If β, λ and ε satisfy the inequalities

                                         τ ≤ λ − ε < λ + ε ≤ β,

then we have
                                                                   πn
                                    λX   πn  sin(λ(y − β ))          πn
(4.6)                     F (y) =      F               πn    ℓc,ε y −    ).
                                    β    β       λ(y − β )            β
                                     n∈Z

                           PN
Furthermore, if A =           j=1 |aj |, then we have


                                                πn
        λ      X              πn  sin(λ(y − β ))          πn 
(4.7)                     F                 πn    ℓc,ε y −
        β                      β      λ(y − β )            β
                 β |> ε
            |y− πn    c

                                                                                                   !
                                                                         2A      log e(c + 1)    2ε
                                                                     ≤                         +      .
                                                                       sinh(c)         π         β

Proof. The proof of (4.6) is outlined in [11], so we only prove the bound (4.7). We start
by estimating the contribution of summands with
                                                      nπ   c
                                                 y−      >
                                                      β    ε

to (4.6). Using the bounds

            sin(x)  1                              c       n      1 o
                   ≤ ,           |ℓc,1 (y)| ≤           min 1,         and |F (t)| ≤ A,
               x    x                           sinh(c)        |y| − c
                                   BOUNDING ψ(x) ANALYTICALLY                                  11


which hold for x 6= 0, y > c and t ∈ R, we get
                                                 πn
             λ     X          πn  sin(λ(y − β ))          πn
                          F                 πn    ℓc,ε y −     )
             β           c    β       λ(y − β )             β
                 y− nπ
                     β >ε
                                                                             
                        Ac    ε β                           1        1
                                                X
                  ≤                     +1 +
                                                                              
                                                               nπ     nπ
                     β sinh(c) c πε                       y − β ε(y − β ) − c
                                                                             
                                                  nπ  c+1
                                                  y− β >   ε
                               "             Z πβ (y− c+1
                                                       ε )
                                                                                  #
                      Ac       2ε   β                       dt
                 ≤                +    +              π
                   β sinh(c) c      πc   −∞      (y − β t)(ε(y − πβ t) − c)
                                              
                      A      2ε   1           
                 =              + log π(c + 1) .
                   sinh(c) β      π
   An analogous calculation gives the same estimate for the contribution of summands
                c
with y − nπ
          β < −ε.                                                                 

                                     5. Run time analysis
5.1. Proof of Theorem 1. Let L > 1, δ > 0 and θ ∈ [1/2, 1). We may assume x0 to be
sufficiently large, since the task can always be carried out trivially using the Eratosthenes
sieve in finite ranges. For simplicity we focus on proving the assertion concerning the
upper bound ML+ (x0 ). The considerations for ML− (x0 ) are almost the same.
   We first address the problem of bounding ψ(t) in I = [x0 , Lx0 ]. Let
                                                 X∗ ℓc,ε ( ρ − 1 ) xρ
                                                               i   2i
(5.1)                          ψ̃c,ε (x) = x −                          ,
                                                        ℓc,ε (i/2) ρ
                                             |γ|≤c/ε

and let η1 < θ. If we assume x0η1 −1 < ε < x−η   0
                                                   1
                                                     and take c = θ log x0 + log log x0 +
log log log x0 − log(δ/40), then Propositions 2 and 3 give the bound
                                                            δ 1−θ    δ 1−θ
(5.2)                   ψc,ε (t) − ψ̃c,ε (t) ≤ 2 + eo(1)     x    <   x
                                                           40 0     20 0
                                                                          √
for t ∈ I and x0 sufficiently large. Consequently, we may take ε = η2 x−θ
                                                                       0    log x0 for every
                                                       −1/2
η2 > 0, and since (µc )+ (0) = 1/2 and |νc (0)| ∼ (2πc)     , we may achieve
                                                                     δ 1−θ
(5.3)                        |ψ(t) − ψc,ε (t)| ≤ C(θ)η2 x1−θ
                                                         0   <        x
                                                                    20 0
for t ∈ I by use of Proposition 1. Now assume we may calculate t − ψ̃c,ε (t) for t ∈ I
within an accuracy < δx1−θ
                       0   /20 and denote this approximation by R(t). Then we get
                       R(t) 3δ 1/2−θ  t − ψ(t)   R(t) 3δ 1/2−θ
(5.4)                  √ − x0        ≤ √       ≤ √ + x0
                         t  20             t       t  20
for t ∈ I which we intend to interpolate. We cannot use Proposition 4 directly since
we assumed the bounds to have opposite sign and since this would also give a slightly
weaker result where applicable. Instead we estimate trivially, which increases the number
of grid points by a factor log x0 . Let S ⊂ I be a finite subset satisfying dist({s}, S \{s}) ≤
η3 x1−θ
    0   / log(x0 ) for all s ∈ S∪{x0 , 2x0 }. Now let s ∈ S, s±t ∈ I and |t| ≤ η3 x1−θ
                                                                                   0   / log(x0 ),
where η3 is sufficiently small. Then estimating as in (3.10) gives
                   (s ± t) − ψ(s ± t)   s − ψ(s)                 δ 1/2−θ
(5.5)                    √            =    √     (1 + O(x−θ
                                                         0 )) +    x
                           s±t               s                  20
for x0 sufficiently large.
12                                      JAN BÜTHE


     Now let
                                            t − ψ(t)
(5.6)                              M0 = max    √     .
                                        t∈I      t
                                              √
Then in view of (5.4) the approximation R(s)/ s yields an upper bound M1 satisfying
                            s − ψ(s)            3δ 1/2−θ
(5.7)                          √     ≤ M1 ≤ M0 + x0
                                s               20
for s ∈ S. By (5.5) this extends to the bound
                    t − ψ(t)                        δ 1/2−θ
(5.8)                  √     ≤ M1 (1 + O(x−θ0 )) +    x       ≤ M0 + δx1−θ
                                                                        0
                         t                         20
                           √
for t ∈ I, since M0 = o( x0 ) for x0 → ∞.
    It remains to analyze the run time for evaluating
                                                    √ R(s) on such a set S. We may take
S = {exp(y0 + kh) | √  k ∈ Z} ∩ I, where y0 = log( Lx0 ) and h = η4 x−θ0 / log(x0 ). We take
T = c/ε ∼ C(δ, θ)xθ log x in (4.2) and consider the trigonometric sum F (y) = Fψ,T (y0 +
yh) which we intend to evaluate within an accuracy of δx−θ       0 /40 for y ∈ [−Y, Y ] ∩ Z,
where Y = max{|k| | exp(y0 + kh) ∈ I}. If x0 is sufficiently large, then F (y) satisfies
the suppositions of Proposition 5 (after reducing γh modulo 2π and evaluating aρ eiγy0
in (4.1), which is done in Õ(xθ )) with B = 1, C = 0 and D = 1. We have N ≍
(xθ0 log(x0 )3/2 ) and Y ≍ (xθ0 log(x0 )), so we may evaluate F (y) within an accuracy of
N −2θ using Õ(N ) = Õ(xθ ) arithmetic operations on Õ(xθ ) variables of size O(log x),
where the implied constants only depend on L, θ and δ. Furthermore, we may evaluate
the contribution of a single zero violating the Riemann hypothesis to the explicit formula
within sufficient accuracy performing Õ(xθ ) arithmetic operations on variables of size
O(log x). For x0 sufficiently large, this yields the desired accuracy and we can recover the
               √
values R(s)/ s with an error < δx1/2−θ /20.                                               
5.2. Reducing the memory requirement. One may reduce the space requirement of
the algorithm by splitting the sum over zeros, applying the method from section 4.1 to the
partial sums and using bandlimited function interpolation to calculate intermediate val-
ues. One then does not evaluate the full trigonometric sum anymore but rather calculates
upper and lower bounds for the partial sums which are subsequently used to calculate
bounds for the full trigonometric sum.
   More precisely, let (L, δ, θ) be an admissible triple in Theorem 1. Then we proceed
as in the proof of Theorem 1, but bound Fψ,T in the following way. Let N = ⌊xη ⌋, let
ρn = 1/2 + iγn be an enumeration of the zeros in the upper half plane satisfying RH
ordered by increasing absolute value and define
                                                X
(5.9)                      Fk (y) = e−iyτk              aρn eiyγ ,
                                          kN <n≤(k+1)N

where τk = (γ(k+1)N − γkN +1 )/2. Since γn ≍ n/ log n the functions Fk have bandwidth
≪ xγ and can thus be recovered from samples Fk (hℓ) where h ≫ x−η . In view of
Proposition 6 it thus suffices to calculate O(xη ) samples which by Proposition 5 can be
done performing Õ(xη ) arithmetic operations on Õ(xη ) variables of size O(log x). Now
                                             √
for each k with γ(k+1)N ≤ T = C1 (L, δ, θ)xθ log x the required Õ(xθ ) evaluations can be
done in Õ(xθ ) using the interpolation formula from Proposition 6. For each k only the
maximal and minimal values of ℜeiyτk Fk (y) are stored, from which one recovers upper
and lower bounds for ℜFψ,T (y). There are Õ(xθ−η ) values k to be considered, so in
total the algorithm performs Õ(x2θ−η ) arithmetic operations on Õ(xη ) variables of size
O(log x).
   It should be noted that the additional error from splitting the sum over zeros could be
avoided by adapting the method from [6] to this problem. This way one would split both
                              BOUNDING ψ(x) ANALYTICALLY                                  13


the trigonometric sum and the interval in question and use direct evaluation combined
with bandlimited function interpolation on every sub interval. For the calculations re-
ported in this paper this additional error was rather small (less than 1% of the calculated
bounds) and the author did not try out this method.

                                 6. Numerical results
   The algorithm has been implemented for L = 2, θ = 1/2 and variable δ and used to
calculate analytic bounds in the range between 1010 and 1019 .
   Function evaluations have been done using the multi-precision library MPFR and the
crucial calculations have been carried out using a 64-bit fixed point arithmetic.
   The calculations used the zeros with imaginary part up to 1011 whose calculation has
been reported in [5] and which were given within an accuracy of 2−64 . The amount of
memory was limited to 340 GB which required a sub-division of the sum over zeros for
x ≥ 4 × 1014 , the maximal amount of summands being 1.25 × 1010 . For the largest
calculation, concerning the interval [5.12 × 1018, 1.024 × 1019] the sum was divided into 13
pieces. This calculation took 290 hours on a 2.27 GHz Intel Xeon X7560 CPU. The run
time could have been reduced further by parallelizing the interpolation routine, which
accounted for half of the computing time. In total, the calculations took less than 1, 200
CPU hours.
   The largest value of a partial sum
                                  1       X ℓc,ε ( ρ − 1 )
(6.1)                                               i   2i
                                                           xρ−1/2
                             ℓc,ε (i/2)              ρ
                                        |ℑ(ρ)|<T
                                                                           √
in the explicit formula for the normalized remainder term (t − ψc,ε (t))/ t that occurred
in the calculations was 0.83545670 . . . at x = 36219716654216.6 . . . with c = 26, ε =
1.7 × 10−8 and T = 917, 647, 060 and the smallest value was −0.783738372378 at x =
1325006525152927089. . . . with c = 31, ε = 2.5 × 1010 and T = 3, 221, 225, 472. The
program aims to calculate the sum over zeros within an accuracy of 10−10 . This does not
include round-off errors, which could be larger but can still be shown to be bounded by
0.016 in these calculations [1]. In addition, the extremal values have been counter-checked
by direct evaluation of the sums in question and the largest deviation was < 6 × 10−12 .
A complete list of parameters and calculated values is given in the appendix to [1].
   The calculated bounds are listed in Table 1. In addition the bounds
(6.2)                              − 0.8 ≤ Rψ (t) ≤ 0.81
                    10
for 100 ≤ t ≤ 5 × 10 have been calculated using the Eratosthenes sieve. Together these
imply the bound (1.5), where the validity for 11 < t < 100 is easily checked by direct
evaluation.
6.1. Bounds for π(x), π ∗ (x), and ϑ(x). We provide several elementary lemmas for
deriving the bounds in Theorem 2 from the calculated bounds for ψ(x).
Lemma 1. Let 1 < a < b and suppose
                                         x − ψ(x)
(6.3)                               c≤      √     ≤C
                                             x
holds for x ∈ [a, b]. Then
            x − ϑ(x)                             x1/3 + x1/5 + 2 log(x) x1/13
(6.4)         √      ≤ C + 1 − c x−1/4 + 1.03883             √
                x                                              x
and
                               x − ϑ(x)
(6.5)                            √      ≥ c + 1 − C x−1/4
                                   x
14                                            JAN BÜTHE


              Table 1. Upper and lower bounds Mψ± (x) for t−ψ(t)
                                                            √
                                                             t
                                                                 in [x, 2x]



                           x    Mψ− (x) Mψ+ (x)                   x     Mψ− (x) Mψ+ (x)
                           10                                     12
                      10          −.77      .85            10             −.80      .81
                  2 × 1010        −.75      .64        2 × 1012           −.79      .76
                  4 × 1010        −.73      .80        4 × 1012           −.73      .73
                  8 × 1010        −.80      .86        8 × 1012           −.80      .76
                 16 × 1010        −.88      .68       16 × 1012           −.80      .68
                 32 × 1010        −.88      .78       32 × 1012           −.67      .93
                 64 × 1010        −.66      .74       64 × 1012           −.78      .77

                        x Mψ− (x) Mψ+ (x)                     x Mψ− (x) Mψ+ (x)
                     1014   −.79      .72                  1016   −.88      .74
                 2 × 1014   −.60      .76              2 × 1016   −.87      .70
                 4 × 1014   −.65      .73              4 × 1016   −.65      .73
                 8 × 1014   −.81      .88              8 × 1016   −.82      .77
                16 × 1014   −.66      .86             16 × 1016   −.71      .92
                32 × 1014   −.74      .86             32 × 1016   −.78      .71
                64 × 1014   −.73      .66             64 × 1016   −.94      .82
                                                     128 × 1016   −.94      .75
                                                     256 × 1016   −.82      .86
                                                     512 × 1016   −.83      .94


hold for x ∈ [a2 , b].
Proof. We need to bound ϑ(x) in terms of ψ(t). To this end we use
                                   ∞                      ⌊2 log x⌋
                                   X                        X
                                               1/k
(6.6)                     ϑ(x) =         µ(k)ψ(x     )=               µ(k)ψ(x1/k ),
                                   k=1                      k=1

and the bounds
(6.7)                    ψ(x) ≤ x log x                                      for x ≥ 1,
(6.8)                    ψ(x) < 1.03883 x                                    for x > 0,

and

(6.9)                    ψ(x) ≥ 0.82 x                                  for x ≥ 100.
The first bound is trivial, the second is proven
                                             Pn in [14, Theorem 12] and the third bound
follows from [14, Theorem 10]. Now, since k=4 µ(k) ≤ 0 for n < 39 and since ψ(x1/k )
decreases monotonously with increasing k, we get
                                                   ⌊2 log x⌋
                                     √        1/3
                                                     X
                     ϑ(x) ≤ ψ(x) − ψ( x) − ψ(x ) +           ψ(x1/n )
                                                                      n=39
                                     √                  2 1/39
                       ≤ ψ(x) − ψ( x) − 0.82 x1/3 +        x     log(x)2
                                                       39
from (6.6) for x ≥ 106 , where we used (6.7) and (6.9) on the second line. The term
              2 1/39
−0.82 x1/3 + 39 x    log(x)2 is easily seen to be negative for x ≥ 106 , so we get
                                                    √
(6.10)                            ϑ(x) ≤ ψ(x) − ψ( x)
                                  BOUNDING ψ(x) ANALYTICALLY                                        15


first for x ≥ 106 , and then by directly checking the remaining values even for x ≥ 0. For
the lower bound we proceed in a similar way, using nk=6 µ(k) ≥ 0 for n < 13, which
                                                        P
gives
                                                           ⌊2 log x⌋
                               √                             X
               ϑ(x) ≥ ψ(x) − ψ( x) − ψ(x1/3 ) − ψ(x1/5 ) −           ψ(x1/n )
                                                                       n=13
                                √
(6.11)               ≥ ψ(x) − ψ( x) − 1.03883 (x1/3 + x1/5 + 2 log(x)x1/13 )
                                                                               √
for x ≥ 1, where we used (6.8) on the second line. Putting ϑ(x) = ψ(x) − ψ( x) + r(x),
the inequalities (6.5) and (6.4) now easily follow by inserting (6.10), respectively (6.11)
and (6.3) into
                                                    √        √
                  x − ϑ(x)    x − ψ(x)                x − ψ( x) r(x)
                           =           + 1 − x−1/4                − √ .
                      x           x                     x1/4           x
                                                                                         
  In order to prove (1.6) and (1.7) we first apply Lemma 1 with a = 100, b = 5 × 1010
and −c = C = 0.81, which gives (1.6) and (1.7) for 107 ≤ x ≤ 5 × 1010 . Switching the
parameters to b = 32 × 1012 and −c = C = 0.88 extends them to 5 × 108 ≤ x ≤ 32 × 1012
and taking b = 1019 and −c = C = 0.94 gives them for 32 × 1012 ≤ x ≤ 1019 . For the
remaining values smaller than 107 the bounds have been verified by a direct computation.
Lemma 2. Let b > 107 , 12 < a < b, let c < 0 and C > 0 satisfy
                                      x − ψ(x)
(6.12)                          c≤       √     ≤C
                                          x
for all x ∈ [a, b], and let
                                                          a − ψ(a)
                                  A = π ∗ (a) − li(a) +            .
                                                            log a
Then we have
                  li(x) − π ∗ (x)   x − ψ(x)    2C       5    log x
                    √             ≤    √     +       1+       +A √ ,
                      x/ log x          x      log x    log x      x
and
                 li(x) − π ∗ (x)   x − ψ(x)     2c      5    log x
                   √             ≥    √     +       1+       +A √
                     x/ log x          x      log x    log x      x
                       7
for all x ∈ [max{a, 10 }, b].
Proof. Partial summation gives
                                                                              Z x
                                                   x − ψ(x) a − ψ(a)                t − ψ(t)
(6.13)       π ∗ (x) − π ∗ (a) = li(x) − li(a) −           +         −                        dt.
                                                     log x    log a            a    t log(t)2
It thus suffices to prove
                               Z x                    √
                                         dt             x        5 
(6.14)                    0≤         √          ≤ 2       2
                                                            1 +
                                a     t log(t)2     log(x)      log x
for a ≥ 12 and x ≥ 107 . Applying the substitution u = log t gives
                Z x               Z log x u/2
                         dt               e
                    √           =              du
                 a    t log(t)2
                                   log a   u2
                                  Z log a+i∞ u/2        Z log x+i∞ u/2
                                              e                   e
                                =                2
                                                   du −                du.
                                   log a       u         log x     u2
For any α > 0 we get
                         Z α+i∞                                  Z ∞
                                  eu/2        2   8                       eit/2
(6.15)           e−α/2               2
                                       du = − 2 − 3 + 24i                         dt
                          α        u         α   α                0     (α + it)4
16                                          JAN BÜTHE


       Table 2. Upper and lower bounds Mπ±∗ (x) for (li(t) − π ∗ (t)) log
                                                                       √ t in [x, 2x].
                                                                         t




                       x    Mπ−∗ (x) Mπ+∗ (x)                x    Mπ−∗ (x) Mπ+∗ (x)
                       10                                   12
                    10        −.87        .95          10           −.88        .89
                2 × 1010      −.84        .73      2 × 1012         −.87        .84
                4 × 1010      −.82        .89      4 × 1012         −.81        .81
                8 × 1010      −.89        .95      8 × 1012         −.87        .84
               16 × 1010      −.97        .76     16 × 1012         −.87        .76
               32 × 1010      −.96        .86     32 × 1012         −.74          1
               64 × 1010      −.74        .82     64 × 1012         −.85        .84

                       x    Mπ−∗ (x) Mπ+∗ (x)              x Mπ−∗ (x) Mπ+∗ (x)
                       14                                    16
                   10         −.86        .79           10     −.94        .80
               2 × 1014       −.67        .83       2 × 1016   −.93        .76
               4 × 1014       −.72        .80       4 × 1016   −.71        .79
               8 × 1014       −.87        .95       8 × 1016   −.88        .83
              16 × 1014       −.72        .93      16 × 1016   −.77        .98
              32 × 1014       −.80        .92      32 × 1016   −.84        .77
              64 × 1014       −.79        .72      64 × 1016      −1       .88
                                                  128 × 1016      −1       .80
                                                  256 × 1016   −.87        .91
                                                  512 × 1016   −.88        .99


by repeated integration by parts. Here, the last integral on the right hand side is bounded
in absolute value by             Z α        Z ∞
                                       dt        dt      4
                                        4
                                          +        4
                                                      =    3
                                                             ,
                                   0   α     α   t      3α
and we get
                      Z α+i∞ u/2
                               e             eα/2       4        16 
                                     du = −2         1 +    +  Θ        .
                        α       u2            α2         α         α2
Thus, the integral on the left hand side is negative for α ≥ log(12). Furthermore, we have
16                    7
 α ≤ 1 for α ≥ log(10 ) so we get (6.14).                                                 
  Choosing a = 100 in Lemma 2 and using the bounds from (6.2) and Table 1 gives the
bounds listed in Table 2. Similarly, one obtains the bound (1.8) for x ≥ 107 and the
remaining values can again be checked by a direct computation.
Lemma 3. Let b > 107 , 12 < a < b, c ≤ 0 and C ≥ 0 such that
                                     x − ϑ(x)
(6.16)                          c≤      √     ≤C
                                          x
holds for all x ∈ [a, b], and let
                                                         a − ϑ(a)
                                    A = π(a) − li(a) +            .
                                                           log a
Then we have
                   li(x) − π(x)   x − ϑ(x)    2C       5    log x
                     √          ≤   √      +       1+       +A √ ,
                       x/ log x       x      log x    log x      x
and
                   li(x) − π(x)   x − ϑ(x)     2c      5    log x
                     √          ≥   √      +       1+       +A √
                       x/ log x       x      log x    log x      x
                                  BOUNDING ψ(x) ANALYTICALLY                                          17


for all x ∈ [max{a, 107}, b]. Furthermore, the implication
               t − ϑ(t) > 0 for 2 ≤ t ≤ T        ⇒     li(t) − π(t) > 0 for 2 ≤ t ≤ T
holds.
Proof. The first assertion follows from
                                                                Z x
                                         x − ϑ(x) a − ϑ(a)           t − ϑ(t)
(6.17)          π(x) − π(a) = li(x) − li(a) −     +           −             2
                                                                              dt
                                           log x      log a      a t log(t)
in the same way as in the proof of Lemma 2. The second part is well-known and follows,
e.g., by taking a = 10 in (6.17) since
                                                     10 − ϑ(10)
                                π(10) − li(10) +                > 0.1.
                                                       log(10)
                                                                                                      
  Choosing a = 1, 500 in Lemma 3 and using (1.6) gives (1.9) for 107 ≤ x ≤ 1019 and the
remaining values have again been checked directly. The bound (1.10) follows from (1.7)
and [14, Theorem 19].

                                            References
 1. J. Büthe, Untersuchung der Primzahlzählfunktion und verwandter Funktionen, Ph.D. thesis, Bonn
    University, March 2015.
 2. Jan Büthe, Estimating π(x) and related functions under partial RH assumptions, Math. Comp. 85
    (2016), no. 301, 2483–2498.
 3. P. Dusart, Autour de la fonction qui compte le nombre des nombres primiers, Ph.D. thesis, Université
    de Limoges, 1998.
 4. Laura Faber and Habiba Kadiri, New bounds for ψ(x), Math. Comp. 84 (2015), no. 293, 1339–1357.
 5. J. Franke, Th. Kleinjung, J. Büthe, and A. Jost, A practical analytic method for calculating π(x),
    Math. Comp. (to appear).
 6. Ghaith A. Hiary, An amortized-complexity method to compute the Riemann zeta function, Math.
    Comp. 80 (2011), no. 275, 1785–1796.
 7. B. F. Logan, Bounds for the tails of sharp-cutoff filter kernels, SIAM J. Math. Anal. 19 (1988),
    no. 2, 372–376.
 8. H. L. Montgomery and R. C. Vaughan, The large sieve, Mathematika 20 (1973), 119–134.
 9.          , Multiplicative Number Theory I. Classical Theory, Cambridge University Press, 2006.
10. A. M. Odlyzko and A. Schönhage, Fast algorithms for multiple evaluations of the Riemann zeta
    function, Trans. Amer. Math. Soc. 309 (1988), no. 2, 797–809.
11. A.M. Odlyzko, The 1020 -th zero of the Riemann zeta function and 175 million of its neighbors,
    http://www.dtc.umn.edu/ odlyzko/unpublished/zeta.10to20.1992.ps, 1992.
12. D. J. Platt and T. S. Trudgian, On the first sign change of θ(x)− x, Math. Comp. 85 (2016), no. 299,
    1539–1547.
13. Barkley Rosser, Explicit bounds for some functions of prime numbers, Amer. J. Math. 63 (1941),
    211–232.
14. J. Barkley Rosser and Lowell Schoenfeld, Approximate formulas for some functions of prime numbers,
    Illinois J. Math. 6 (1962), 64–94.
15.          , Sharper bounds for the Chebyshev functions θ(x) and ψ(x), Math. Comp. 29 (1975), 243–
    269.
16. Lowell Schoenfeld, Sharper bounds for the Chebyshev functions θ(x) and ψ(x). II, Math. Comp. 30
    (1976), no. 134, 337–360.
                                                                                                     13
17. Douglas A. Stoll and Patrick Demichel, The impact of ζ(s) complex zeros on π(x) for x < 1010 ,
    Math. Comp. 80 (2011), no. 276, 2381–2394.

   Hausdorff Center for Mathematics, Endenicher Allee 62, 53115 Bonn
   E-mail address: jan.buethe@hcm.uni-bonn.de
