An improved sieve of Eratosthenes
Harald Andres Helfgott
math.NT 11Y05, 11Y06, 11Y16


  We show how to carry out a sieve of Eratosthenes up to N in space
  O(N^1/3 (log N)^2/3) and time O(N log N).
  In comparison, the usual versions of the sieve
  take space about O(√(N)) and
  time at least linear on N. We can also apply our sieve to any subinterval
  of 1,N of length
  Ω(N^1/3) in time close to linear on the length
  of the interval. Before, such a thing was possible only for subintervals of
  1,N of length Ω(√(N)).
  
  Just as in (Galway, 2000), the approach here is
  related to Diophantine approximation, and also has close ties
  to Voronoï's work on the Dirichlet divisor problem. The advantage of the method here
  resides
  in the fact that, because the method we will give is based on the sieve of
  Eratosthenes, we will also be able to use it to factor integers,
  and not just to produce lists of consecutive primes.
  





A Random Block-Coordinate Douglas-Rachford Splitting Method with Low Computational Complexity for Binary Logistic Regression
This work was partly supported by the the CNRS MASTODONS project under grant 2016TABASCO.

    
Luis M. Briceño-Arias
Giovanni Chierchia
Emilie Chouzenoux
Jean-Christophe Pesquet

    Received  2009 month day; accepted  2009  month day
========================================================================================================================================================================================================================





§ INTRODUCTION


The sieve of Eratosthenes is a procedure for constructing all primes up
to N. More generally, such a sieve can be used
for factoring all integers up to N, or to compute the values f(n),
n≤ N, of arithmetical functions f that depend on the factorization
of integers. For instance, one can take f=μ, the Möbius function, or
f=λ, the Liouville function.

The point of the sieve of Eratosthenes is that it can be carried out in time
close to linear on N, even though determining whether an individual
integer is prime, let alone factoring it, takes much more than constant
time with current methods. Though a very naïve implementation of
the sieve would take space proportional to N, it is not hard to see how
to implement the sieve in space O(√(N)), simply by applying the
sieve to intervals of length O(√(N)) one at a time; the time taken
is still close to  linear[
  In its standard form,
  the sieve of Eratosthenes
  (segmented or not) takes time O(N loglog N); when used for sieving
  out primes, its time consumption can be reduced to O(N).

  We follow the convention that
arithmetic operations (e.g., adding or multiplying two numbers) take
O(1) time. This convention holds in the range in which using
a sieve is realistic. (Outside that range, such operations do take time
O((log N)^c), c close to 1.) We shall measure space in bits, again
reflecting how matters work in practice.]





  on N. (One could take shorter
intervals, but the algorithm would then become much less efficient.)
This is called the “segmented sieve”; see <cit.> for an early reference.

Of course, the output still takes space linear on N,
but that is of less importance: we can store the output in slower memory
(such as a hard drive), or give the output in batches to a program that
needs it and can handle it sequentially.
Then the space used is certainly O(√(N)).











There has been a long series of improvements to the basic segmented sieve.
Most of them improve the running time or space by a constant factor
or by a factor in the order of loglog N. Many work only when the
sieve is used to construct primes, as opposed to computing μ, say.
See <cit.>, which reviews the state of matters at the time
before improving on it; see also <cit.> for further references
and for a contemporary implementation combining some of the most useful
existing techniques.

In practice, saving space sometimes amounts to saving time, even when it seems, at
first, that there should be a trade-off.
As of the time of writing,
a good office computer can store about 10^12 integers in very slow memory
(a hard drive), about  10^9 integers in memory working at intermediate speed
(RAM), and about 10^6 integers in fast memory (cache); a program becomes
faster if it can run on cache, accessing RAM infrequently.
Having enough RAM is also an issue; sieves have been used, for instance, to
verify the binary Goldbach conjecture up to 4· 10^18 <cit.>,
and so we are close to the point at which O(√(N)) space
might not fit in RAM.
Space constraints
can become more severe if several processor cores work in parallel and share resources.

Moreover, finding all primes within a short interval can be useful in itself.
For instance, there are applications in which we verify a conjecture on
one interval at a time, and we neither need nor can store a very long interval
in memory. See the verification of Goldbach's (binary) conjecture up to 4· 10^18 in <cit.>, and, in particular <cit.>.

Galway <cit.> found a way to sieve using space O(N^1/3) and
time O(N). Like the sieve in <cit.>, on which
it is based[Atkin and Bernstein's preprint was already available in 1999, as the bibliography in <cit.> states.], Galway's sieve is specific to finding prime numbers. There is also the algorithm in
<cit.>, specific, again, to finding primes:
under the assumption of the Generalized Riemann Hypothesis, it
finds all primes up to N in space O((log N)^3/loglog N) and
time O(N (log N)^2/loglog N); unconditionally, it runs in space
O(N^0.132) and time O(N^1.132). (It runs in time O(N log N) under the
assumption of a more specialized conjecture.)

We will show how to implement a sieve of Eratosthenes in space 
close to N^1/3 and still close to linear time. Our method is not limited
to finding prime numbers; it can be used to factor integers, or, of course,
to compute μ(n), λ(n) or other functions given by the
factorization of n.


  We can construct all primes p≤ N in 

    space  O(N^1/3 (log N)^2/3)        and      time   O(N log N).

We can also factor all integers n≤ N in 

    space  O(N^1/3 (log N)^5/3)     and      time   O(N log N).


  Moreover, for 
N^1/3 (log N)^2/3≤Δ≤ N, we can construct all primes in an interval
  N-Δ,N+Δ in

    space   O(Δ)     and      time   O(Δlog N)

and factor all integers in the same interval in

    space   O(Δlog N)      and      time   O(Δlog N).


Here we recall that space refers to the number of bits used, and
time to the number of operations of words used (on integers of size
O(N)). 

The main ideas come from elementary number theory.
In order for us to be
able to apply the sieve to an interval I of length O(N^1/3) without large
time inefficiencies, we need to be able to tell in advance which primes
(or integers) d
up to √(N) divide at least one integer in I, without testing each d
individually. We can do this by Diophantine approximation, followed by
a local linear
approximation to the function x↦ n/x for n fixed, and
then by solving what amounts to a linear equation 1.

The idea of using Diophantine approximation combined with a local linear approximation
is already present in <cit.>, where
it  was used to compute ∑_n≤ xτ(n) in time
O(x^1/3 (log x)^O(1)). (We write
τ(n) for the number of divisors of an integer n.)
The basic underlying idea in
<cit.> may be said to be the same as the one here:
we are speaking of a Diophantine
idea that stems ultimately from Voronoï's
work on the Dirichlet divisor problem <cit.> (in our work) and
Sierpinski's adaptation of the same method to the circle problem
<cit.> (in the case of <cit.>).[
To be precise, the immediate inspiration for <cit.> came
from Vinogradov's simplified version of Voronoï's method, as in
<cit.>. A 
bound of roughly the same quality as Voronoï's result was claimed long
before Voronoï in
<cit.>.
 While the proof there was apparently incomplete, Landau
later showed <cit.> that it could be made into an actual proof,
and sharpened to give a result matching Voronoï's.
Thanks are due to S. Patterson for pointing out Pfeiffer's result to the 
author.]
For that matter,  <cit.> already suggests that Voronoï's
work on the Dirichlet divisor problem could be used to make the
sieve of Eratosthenes in space about O(N^1/3) and close to linear time.
We should also make clear that Galway can sieve out efficiently
segments of length about x^1/3, just as we do.


One difference between this paper and Galway's
is that the relation to Voronoï's and Sierpinski's work in
Galway's paper may be said to be
more direct, in that Galway literally dissects a region
  between two circles, much as Sierpinski does. In the case of the
  present paper, we can say that Voronoï's 
 main idea originated in the context of giving
 elementary estimates (for ∑_n≤ xτ(n))
 and is now used to carry out
  an exact computation.

  Another precedent that must be mentioned is Oliveira e Silva's technique for
  efficient cache usage <cit.>, <cit.>.
    It seems that
  this technique can be combined with the algorithm here. Such a
  combination can be useful, for instance, when N is so large that N^1/3 bits fit in RAM but not in cache.

  Additional motivation. My initial interest
  in the problem stemmed from the fact that I had to compute values of μ(n)
  so as to check inequalities of the form[Mertens's conjecture
    states that the inequality ∑_n≤ xμ(n) ≤√(x) 
    holds for all x. That conjecture has been
    disproved <cit.>, but the inequality is known to hold for all
    x≤ 10^16
    <cit.>, and may hold in a far wider range.]
    ∑_n≤ xμ(n) ≤√(x),
    ∑_n≤ xμ(n)/n ≤ 2/√(x), etc., for all x less
  than some large finite N. The importance of such sums in number theory is
  clear. Explicit results on them for x bounded help complement analytic
  estimates, which are generally strong only when x is large.

  While we can obviously determine values of μ(n) by first factorizing
  n and then using the definition of μ, it is also possible and
  rather simple
  to save some space and time in practice by modifying the procedure
  we will give (Algorithm <ref>)
  so as to keep track of μ(n) instead of the list of factors.
  Time and space complexity remain the same.
  
  Conventions and notation.
  Integer operations are assumed to take constant time.
  As is usual,
  we write f(x) = O(g(x)) to mean that there exists a constant C>0
  such that |f(x)|≤ C g(x) for all large enough x. We write
  f(x) = O^*(g(x)) to mean that |f(x)|≤ g(x) for all x.
  We use either f(x) ≪ g(x) or g(x) ≫ f(x) to mean f(x) = O(g(x)).
    
  For n a non-zero integer and p a prime, v_p(n) denotes the largest
  k such that p^k|n. As is customary,
  we write (a,b) for the gcd (greatest common divisor)
  of a and b, provided that no confusion with the ordered pair (a,b)
  is possible.

  Given α∈ℝ, we denote by
 {α} the element of 0,1) congruent to
α modulo 1, i.e., modulo ℤ.

  Acknowledgements.
  The author is currently supported by funds from his Humboldt Professorship.
Thanks are due to Manuel Drehwald, for having written code implementing
an earlier version
of the algorithm, to Lola Thompson, for helpful suggestions, and to
an anonymous referee, for several useful remarks.




§ ANALYSIS OF THE PROBLEM


Let I =  n-Δ, n+Δ⊂ 0,x be
an interval. How can we tell which integers
m≤√(x) have at least one integer multiple in the interval I?

Our motivation for asking this question is that we will be taking intervals I
and sieving them by all integers, or all primes, m≤√(x).
Going over all integers ≤√(x) would take time at least √(x),
which could be much larger than Δ.








We will be able to sieve our intervals by m≤√(x) by producing
a list of those m which might have multiples in I,
without testing all m≤√(x) in succession. A quick probabilistic heuristic
hints that the number of such m should be proportional to Δlog x,
which, for Δ∼ x^1/3, is much smaller than √(x).

Let K≥ 2 be a parameter to be set later.
If m≤ K Δ,  

we simply sieve by m. Doing so takes time O(1+Δ/m). We could,
of course, test first whether m is prime; we will discuss this option later.

Assume henceforth that m>K Δ.
The interval I contains a multiple of m if and only if

    {n/m}∈ -Δ/m,Δ/m 1,

that is, {n/m}∈ -Δ/m,Δ/m + ℤ.

  Say we have already dealt with all m≤ M, where M≥ K Δ,
  and that we want to examine m close to m_0, where m_0>M.
Write m = m_0 + r.
The truncated Taylor expansion
    
    f(m) = f(m_0) + f'(m_0) r + f”(m_0 + θ r)/2 r^2

    (where θ is some element of 0,1) gives us, once
    we set f(x) = n/x,
    
    n/m = n/m_0 - n/m_0^2 r + O^*(n/
            m_-^3 r^2),

    where m_- = min(m,m_0).
    We will make sure that r is small enough that
    n r^2/m_-^3 ≲κΔ/m,
    where κ>0.
    (That is, we allow our error term to be not much larger than the interval
    we are trying to hit.) 
    We ensure r satisfies this condition by letting
    
    R = ⌊√(κΔ/n) M⌋,      
        m_0 = M + R.

    Then n r^2/m_-^3 ≤ n r^2/M^3 ≤κΔ/M for all m=m_0+r,
    r∈ -R,R.
    We let κ=1/4, since it is a value that is neither too large nor too
    small. We will also assume Δ≥ n^1/3, and so
    R≥⌊ K/2⌋≥ 1, since K≥ 2.

    
 
 
 
    
    We have thus reduced our problem to that
    of quickly finding all r in an interval -R,R such
    that P(r) ∈-η,η 1, where
    P(r) is the linear polynomial - (n/m_0^2) r + (n/m_0) and
    η = (1+κ) Δ/M = 5 Δ/4 M.
    In other words, we are being asked to find all approximate solutions
    to a linear equation in ℝ/ℤ efficiently. 
    Let us
    assume that K≥ 5/2, so that η≤ 1/2; otherwise there is not much
    to do.
    
    Let α_1 = - (n/m_0^2), α_0 = n/m_0.
  Given a rational number
  α, we can find – by means of continued fractions –
  an approximation a/q to α with q≤ Q and
  
    |a/q - α|≤1/q Q

  in time O(log Q).
 
 
  The procedure DiophAppr(α,Q) to obtain a/q is given in Algorithm <ref>; it runs in time O(log Q) and constant space.
  The fact that its output satisfies (<ref>)  follows from <cit.>
  or <cit.>.
  (In the notation of both sources:
      since any two consecutive approximants p_n/q_n, p_n+1/q_n+1
      to α satisfy |p_n/q_n - α|≤ 1/q_n q_n+1,
  the last p_n/q_n with q_n≤ Q satisfies |p_n/q_n-α|<1/q_n Q.)

  We invoke DiophAppr(α_1,2R), and obtain a rational
  a/q with (a,q)=1 and q≤ Q = 2 R satisfying (<ref>) for α = α_1.
  Hence, for r∈ - R,R,
  
    P(r) = α_1 r + α_0 = α_0 + a r/q + O^*(1/2q) ≡c + a r/q + O^*(1/q)  1

  for c = ⌊α_0 q + 1/2⌋.
  We have thus reduced our problem to a problem in ℤ/qℤ:
  if P(r) ∈-η,η, then
  
    c + a r ∈{-k-1,-k,…,k+1} q,

  where k = ⌊η q⌋. We must find the values of r
  satisfying (<ref>).

  The solutions to (<ref>) are clearly given by
  
    r ≡ - a^-1 (c + j)  q

  for -k-1≤ j≤ k+1. The multiplicative inverse a^-1 q
  of a can be computed in time O(log q) by the Euclidean algorithm.
In fact, we compute it
 at the same time as the continued fraction that gives us a/q: by <cit.>, two consecutive approximants
    p_n/q_n, p_n+1/q_n+1 satisfy p_n^-1 = (-1)^n+1 q_n+1 q_n;
moreover, p_n^-1 = (-1)^n+1 q_n-1 q_n. 

Thus, we simply need to go over all m of the form m_0+r, where
r goes over integers in -R,R satisfying (<ref>).
Since m_0-R = M > K Δ≥ 2 Δ, each such m has at most
one multiple in I= n-Δ,n+Δ. We do not miss any
m∈ m_0-R,m_0+R =  M,M+2 R that has a
multiple in I.

Due to the errors involved in the truncation of the Taylor expansion 
and in Diophantine approximation, we may obtain some m that have no
multiples within I. We simply ignore them, after checking that that is the
case.





§ DESCRIPTION OF THE ALGORITHM








We have already given a nearly full description of the algorithm. It only
remains to give the pseudocode (Algorithms <ref>–<ref>), with some commentary.
















All variables are integers, rationals, or tuples or sets with
integer or rational entries. Thus, all arithmetic operations can be
carried out exactly. We write (α) and (α) for
the numerator and denominator of α∈ℚ (written
in minimal terms: (α), (α) coprime, (α)>0).

In Algorithm <ref>, 
we are defining α_1 {- n/m_0^2}, α_0 {n/m_0},
rather than α_1  - n/m_0^2, α_0  n/m_0, as in the
exposition above.


The motivation is simply to keep variable sizes small. It is easy to see that
the values of α_0, α_1 matter only ℤ.



 §.§ Sieving for primes


The main function is NewSegSiev
(Algorithm <ref>).
It takes as inputs n and Δ, as well as the parameter K, which
affects time and space consumption.
The basic procedure is easy to summarize.
NewSegSiev uses SubSegSiev (Algorithm <ref>),
a segmented sieve of traditional type, to sieve for primes p≤ K Δ. 
It then proceeds as detailed in <ref> to find the integers m
in M, M+2 R that may divide an integer
in the interval n-Δ,n+Δ. Then it sieves by
them. The key step, in finding such m,
is to call DiophAppr, which
uses continued fractions in a standard way to find (i)
a rational approximation a/q to the input α_1, (ii)
the inverse a^-1 q.

We could try to improve on NewSegSieve 
by throwing out even values of m, say; as they are certainly not prime,
we need not sieve by them. More details are given in “Side notes on wheels”,
below.

Preexistent sieves as subroutines. Going back to SubSegSiev: we could avoid using it at all,
just by sieving by all integers m≤ K Δ, rather
than by all primes p≤ K Δ. That would give a running time
of O(K Δ + Δlog M) rather than O(K Δloglog M).
We take the slightly more complicated route in Algorithm <ref>,
not just because it is better for K = o((log M)/loglog M),
but also for the sake of exposition, in that we get to see several existing
forms of the sieve of Eratosthenes. Note, however, that none of this will
decrease the order of magnitude of the time taken by our entire algorithm,
since the total time will be at least in the order of Δlog M.

Function SimpleSiev is a relatively simple kind of sieve of
Eratosthenes. It sieves the integers up to N by the primes up to
√(N), where these primes are found by this very same process.
It is clear that we need to sieve only by the primes up to √(N),
since any composite number n≤ N has at least one prime factor
p≤√(N).
Sieving only by the primes, rather than by all integers,
is enough to take down the running time to O(N loglog N).
We also use a very primitive “wheel” in that we sieve using only
odd multiples of the primes. (Again, see the comments on wheels below.)

The basic segmented sieve is implemented as
SimpleSegSiev(n,Δ,M).
It uses SimpleSiev so as to determine the primes up to M.
Then it sieves the interval n,n+Δ by them.

We could use SimpleSegSiev instead of SubSegSiev.
The point of SubSegSiev is simply to reduce space consumption
by one more iteration: SubSegSiev determines the primes up to M by
SimpleSegSiev, taking only space O(√(M)) at a time; it sieves
n,n+Δ
by these primes as it goes along. The total time
taken by calls on SimpleSegSiev is O(M loglog M); to this
we add the time O(Δloglog M)
taken by sieving the interval n,n+Δ by the primes up to M.

Incidentally: it should be clear that all factors of loglog N, loglog M and the like are coming from Mertens's classical asymptotic statement

    ∑_p≤ N1/p = loglog N + O(1).

For instance, the number of times the
instructions S_n'-n 0, n'  n' + m
are executed in SimpleSegSiev or SubSegSiev
is at most

    ∑_p≤ M⌈Δ+1/p⌉   ≤Δ∑_p≤ M1/p + ∑_p≤ M 1 
       = O(Δloglog M + M).
 

Side note on wheels. In general, a “wheel” is just
(ℤ/P ℤ)^*, where P = ∏_p≤ C p for some constant
C. We would use it by sieving only by multiples m· p' of the
primes p',
where m reduces P to an element of the wheel. Obviously m  P
should be constantly updated by shifting as m increases, rather than be
determined by division each time; hence the “wheel”. It is possible (and
common, at least in theoretical analyses) to choose C of size
δlog N for a small constant δ
(say, δ=1/4), with the consequence that P∼ N^δ.
Then only a proportion ∏_p≤ C (1-1/p) ∼e^-γ/log C
= O(1/loglog N) of the elements of
ℤ/Pℤ (here γ is Euler's constant)
are in (ℤ/P ℤ)^*. With appropriate coding, this fact can
be used to reduce the running time of a simple sieve or a segmented
sieve for primes (such as SimpleSiev or SimpleSegSiev)
by a factor of loglog N <cit.>.
We will not bother including wheels in our pseudocode, since they would
reduce only a second-order term of the total time bound in this fashion; the 
time complexity of NewSegSiev would remain the same. They
can obviously be added in implementation.

It is tempting to introduce a wheel in a different place, namely,
to make sure that the variable m
in NewSegSiev has no prime factors
p<δlog N. The hope there would be to reduce the main term in the
total time by a factor of loglog N. However, as we will later see,
q will be usually close to 2 R; thus, most of the time, it will not
make sense to use a large wheel as m goes over integers in
(m_0+r_0 + qℤ)∩ M,M+2R, as there will be
few such integers. Using a wheel on q instead – something that would
take time ≫ q to set up – would make no sense:
we do not go over ℤ/qℤ several times, but rather
go over a few elements of it once.

What can make sense is introduce a very small wheel, of bounded size, to
attempt to gain a constant factor. For instance, a wheel of size 2
would work as follows: if q is even, then, in the loop on j
in NewSegSiev, we consider only values of j such that
c+j is odd; if q is odd, then we hop over every other value of m
in a given congruence class mod q within M,M+2R,
considering only odd values m. This sort of modification will reduce total
time consumption only by a constant factor, and so, for the sake of
simplicity, we do not include it in the pseudocode, or use it further.

Variable size. We should bound the size of our variables in case
we choose to implement our algorithm using fixed-size integers and rationals. It is easy
to see that our integers will be of size ≤ n+Δ. We should
also bound the size of our rational variables.
It is trivial to modify function NewSegSiev
so that m_0 is always ≤√(n+Δ).


It is then easy to show that
we work entirely with integers – in numerators, in denominators
or on their own – of size
≤min(n+Δ,(2 R)^2)≤ n+Δ.

For n very large and Δ not as large, it may be helpful to
store some variables (such as a, a^-1, q, k, j and r_0)
in smaller integer types
(64 bits, say); these variables are all bounded by
2 R ≤ M √(Δ/n)≤√(Δ (n+Δ)/n).

















 §.§ Sieving for factorization


We will now see how to modify our algorithm so that it factorizes all integers
in the interval n-Δ,n+Δ, rather than simply finding all
primes in that interval. Time and space usage will not be much greater
than when we are just finding primes.
It goes without saying that this makes it possible
to compute various arithmetic functions (the Möbius function μ(m),
the Liouville function Λ(m), etc.) for all m∈ n-Δ,n+Δ.

For the sake of clarity, we first give a well-known procedure for
factoring all integers in an interval n,n+Δ by
means of the sieve of Eratosthenes (Algorithm <ref>), just
as we went over a traditional segmented sieve (Algorithm <ref>)
before describing our sieve for primes.
We will later reuse most of the subroutines.




Our new sieve for factoring
(NewSegSievFac,
Algorithm <ref>), designed for intervals around x
of length Δ≫ x^1/3, is very similar
to our sieve for primes (Algorithm <ref>).
We use a classical segmented
sieve (SubSegSievFac, Algorithm <ref>)
to find all factors p≤ KΔ of n+j for -Δ≤ j≤Δ.
After the call to SubSegSievFac, the variable Π_j contains
∏_p≤ K Δ p^v_p(n+j).
Since K Δ > (2 n)^1/3 > (n+Δ)^1/3,
we see that (n+j)/Π_j is either 1 or the product of at most two
primes > K Δ (not necessarily distinct).
In the innermost loop of SubSegSievFac,
when we come across an m>KΔ that divides not just n+j, but
(n+j)/Π_j,
we multiply Π_j by m, or by its square, if m^2|(n+j)/Π_j,
and include m (or its square) in the factorization.
Note that m
has to be a prime, or else Π_j would have already been
multiplied by some factor p^v_p(n+j), p|m, p<m, either earlier in the loop
or in SubSegSievFac, thus making m|(n+j)/Π_j impossible.

We should also explain the purpose of the final loop in
NewSegSievFac. (The loop in the classical procedure
SegSievFac is identical, and plays the same role.)
Once we take care of all m≤√(n+Δ) (and possibly some beyond),
what we have, for each -Δ≤ j≤Δ, is either
that Π_j = 1, and F_j contains a full factorization of n+j,
or Π_j  1, and F_j is missing a single large prime factor p,
which has to be equal to n/Π_j. We include that prime factor in the
factorization and are done.







§ TIME ANALYSIS. PARAMETER CHOICE.


The space and time consumption of Algorithms
<ref>–<ref>
is clearly as stated.






Time consumption of main algorithm.
Let us analyze the time consumption of Algorithm <ref>.
(Its space consumption, namely, O(Δ + √(K Δ)), will be clear.)
Sieving by primes p≤ KΔ gets done
by the traditional segmented-procedure
SubSegSieve(n,2Δ,KΔ), which takes time
O(KΔloglog K Δ) and space O(√(KΔ) +Δ).
We must analyze now how much time it takes to sieve by integers
KΔ<m≤√(n+Δ).
(Our algorithm sieves by integers KΔ<m≤√(n+Δ),
not just by primes, simply because the Diophantine-approximation
algorithm cannot tell in advance which of the integers it outputs
will be prime. As we discussed before, we could apply a small wheel in the
hope of saving a constant factor in time.)



























































We will use the main ideas of Vinogradov's proof of the bound

    ∑_n≤ xτ(n) = x log x + O(x^1/3 (log x)^5/3),

as given in <cit.>[Actually,
  in the given reference, Vinogradov gives a bound of
  O(x^1/3 (log x)^2) on the error term in (<ref>).
  The reason is simply that he did not choose his parameters optimally:
  if the parameter τ in his proof is set to
  if the value of τ in the
(A log A)^1/3 rather
  than A, the resulting bound is indeed as in (<ref>).], though
we will not use the bound (<ref>) itself. Our treatment will
be self-contained.
  
Algorithm <ref> does not correspond
extremely closely to Vinogradov's approach – 
we use our approximations n/m_0, -n/m_0^2 on the relatively broad intervals
on which they are useful, whereas Vinogradov's procedure changes 
approximations constantly. Nevertheless, we will be able to use the basic
approach that he took to bound an error term, though we of course will use
it to bound time consumption.

Let us look at an interval M,M+2R.
Finding a, a^-1, q by Diophantine approximation (function
DiophAppr) takes, as we know, time O(logmax(Q,m_0^2)) =
O(log n). The number of iterations of the
outer loop in NewSegSieve is

    ≪√(4 n/Δ)·log√(n)/K Δ≤√(n/Δ)log n
    ,
 and thus the total time taken
by the instructions inside the outer loop but outside the inner loop is
O(√(n/Δ) (log n)^2).

The time it takes to go over and sieve by all m∈ M,M+2R
congruent to m_0 - a^-1 (c+j)  q for all
-k-1≤ j≤ k+1 (where k=⌊η q⌋) is at most

    ⌈2 R+1/q⌉· (⌊η q⌋+ 3)    ≤(2 R/q + 2) (η q + 3) ≪
    (R+q) η + R/q  + 1
       ≪Δ^3/2/√(n) + R/q + 1,

since R = ⌊√(Δ/4 n)· M⌋, 
η = 5 Δ/4 M and q≤ R.


Since the number of times the main loop is executed – that is, the number
of intervals M,M+2R we consider – is
given by (<ref>),
the total contribution of the first and last terms in the last line
of (<ref>) is

    ≪(Δ^3/2/√(n) + 1)
    √(n/Δ)log n
    ≪Δlog n,

since Δ≥ n^1/3. It remains to account for the contribution
of R/q.
We may assume M is large enough for R to be ≥ 3, as otherwise
the contribution of R/q is bounded by thrice the contribution of the last
term 1.

Now we proceed much as in <cit.>.
We will examine how α_1 = {- n/m_0^2} changes as m_0 increases;
we will then be able to tell how often Diophantine approximations a/q
to α_1 with given q can occur. To be precise:
we will see by how much m_0 has to increase for -n/m_0^2 to increase by
1 or more, and we also want to know for how long -n/m_0^2 can have
 a given, fixed Diophantine approximation a/q as m_0 increases.

Consider two intervals m_0-R,m_0+R,
m_0'-R',m_0'+R', where m_0<m_0' and R≤ R'. Then

    (-n/(m_0')^2) - (-n/m_0^2)
    = n (m_0')^2-m_0^2/m_0^2 (m_0')^2.


If α_1 = -n/m_0^2 is a/q + O^*(1/q R) and
α_1' = - n/(m_0')^2 is a/q + O^*(1/q R'), it follows that
n/m_0^2 - n/(m_0')^2 is O^*(1/q R + 1/q R') = O^*(2/q R).
Suppose that this is the case.

Since Δ<n, we know that
M' := m_0' - R' ≤√(n+Δ)< √(2 n),
R'≤ M' √(Δ/4 n) < M'/2
and m_0' = M'+R' < (√(2)+1/2) √(n)< 2√(n); in the same way,
m_0 = M+R < 3 M/2. 
Clearly, m_0'≤ 2 m_0, as otherwise
n/m_0^2 - n/(m_0')^2 ≥ 3 n/4 m_0^2≥ 3 n/(m_0')^2 > 3/4,
giving us a contradiction to 2/q R ≤ 2/R ≤ 2/3.
Hence 
    2/q R≥
    n (m_0')^2-m_0^2/m_0^2 (m_0')^2 = n ·m_0'-m_0/m_0^2·m_0'+m_0/(m_0')^2≥3 n/4·m_0'-m_0/m_0^3.

In other words, when the same approximant a/q is valid at m_0 and m_0',

    m_0'-m_0 ≤8/3m_0^3/q R n.








The number of intervals for which α_1 = -n/m_0^2 has a given
approximation a/q is thus at most

    8/3m_0^3/q R (2 R + 1) n + 1   ≪m_0^3/q (R + 1)^2 n + 1≤m_0^3/q M^2 Δ/4 + 1
       < (3/2)^2 4 m_0/q Δ + 1 ≪m_0/q Δ + 1.




We should also see when n/m_0^2-n/(m_0')^2 ≥ 1, or rather
when n/m_0^2-n/(m_0')^2 ≥ 1-2/R, as then n/m_0^2 and n/(m_0')^2
may have Diophantine approximations that, while distinct, are congruent
1, i.e., differ by an integer. 
By (<ref>), the inequality
n/m_0^2-n/(m_0')^2 ≥ 1-2/R
is fulfilled exactly when

    n ((m_0')^2-m_0^2) ≥ m_0^2 (m_0')^2 (1-2/R),

and, since R≥ 3, that implies

    6 n (m_0' - m_0) > m_0^3.

Hence, for given m_0, it makes sense to consider all following intervals with
m_0'≤ m_0 + m_0^3/6 n. In that range, n/m_0^2 and n/(m_0')^2 can
have the same Diophantine approximation 1 only if they in fact have the
same Diophantine approximation.

Since R increases and the intervals are of width 2 R,
there will be at most
m_0^3/12 R n + 1 ≪ m_0^2/√(Δ n) + 1≪ m_0^2/√(Δ n)
intervals with m_0'≤ m_0+m_0^3/6 n.
(Note that m_0^2/√(Δ n)≥ K^2 Δ^2/√(Δ n)> 1.)
Among those intervals,
≪ m_0'/q Δ + 1 ≪ m_0/q Δ + 1 will have a given approximation a/q  1.

As we said before, we have to account for the total contribution of R/q.
Since 1/q and 1/q^2 decrease as q increases, the worst-case scenario
is for there to be as many m_0's as possible for which the approximation
has q as small as possible. We would have
⌊ O(m_0/q Δ+1) ϕ(q) ⌋ = O(m_0/Δ+q) values of a/q with q≤ Q, and none for
q>Q, where Q≪ m_0/(Δ n)^1/4.
(To bound Q, we apply
∑_q≤ Qϕ(q) = |{1≤ q_1,q_2≤ Q: q_1, q_2 coprime}| ≫ Q^2.)
The contribution of R/q for
all intervals
with m_0'≤ m_0+m_0^3/6 n is thus 

    ≪ R ∑_q≤ Q(m_0/Δ q + 1)
    ≪ R ·(m_0/Δlog Q + Q).

Now split M_0, 2 M_0 into
O(n/M_0^2) chunks of the form m_0, m_0 (1+ m_0^2/6 n).
We obtain that the contribution for all intervals with m_0 inside
M_0,2 M_0 is

    ≪n/M_0^2· 
    2 M_0 √(Δ/4 n)·(M_0/Δlog n +
    M_0/(Δ n)^1/4)
       ≪√(n)/√(Δ)log n +(n Δ)^1/4,

and thus the total contribution will be

    ≪√(n)/√(Δ) (log n)^2 +(n Δ)^1/4log n.


We conclude that the total time consumption of Algorithm
<ref> is

    ≪Δloglog K Δ + Δlog n + 
    √(n)/√(Δ) (log n)^2 +(n Δ)^1/4log n
       ≪Δlog n + √(n/Δ) (log n)^2,

since Δ≥ n^1/3. Under the stronger assumption
Δ≥ n^1/3log^2/3 n, we obtain that the total time consumption
is

    O(Δlog n).








Time consumption

  of Algorithm <ref>. We must now analyze
Algorithm <ref>,
which factors integers
in the interval n-Δ,n+Δ. 
Everything is much the same as for Algorithm <ref>,
except in one respect. The total space taken
will be

    O(Δlog n)

rather than O(Δ), simply because storing the list of prime factors
of an integer in n-Δ,n+Δ
takes O(log n) bits. (To wit: the number of bits taken to store a factor
p^α is O(log p) + O(log(α+1)), which is O(αlog p).
Hence, storing all the factors p_i^α_i of an integer
n = p_1^α_1… p_k^α_k takes
O(∑_i α_i log p_i) = O(log n) bits.)






















The total time estimate is still

    O(Δlog n).

The claims in (<ref>) and
(<ref>) follow immediately once one splits the interval
1,N into intervals of length 2Δ.
The main theorem is thus proved.



§ FURTHER PERSPECTIVES

   
   It is tempting to try to improve on this algorithm by taking a longer truncated Taylor
   expansion in (<ref>). However, this would require us to
   find the solutions x∈{-R,-R+1,…,R} to 
   a P(x)∈ -η,ηℤ,
   where P is a polynomial with
   P ≥ 2 and η≪ 1/R or thereabouts.

Solving a quadratic modular equation
a_2 x^2 + a_1 x + a_0 ≡ 0  q is not the main difficulty in our
context. We can obviously reduce such a problem to taking square-roots
q. Now, the problem of finding square-roots modulo q
(for q arbitrary) is well-known to be equivalent to
   factoring q <cit.>.
(There is a classical algorithm (Tonelli-Shanks) for finding square-roots 
to prime modulus.) This is not a problem, since we can factor
all integers q≤ x^1/4 (say) in advance, before we start sieving.

The problem is that it is not at all clear how to reduce 
finding solutions to P(x)∈ -η,ηℤ,
P=2,
to finding solutions to quadratic equations q. We can try to find
rational approximations with the same denominator q (simultaneous
Diophantine approximation) to the non-constant coefficients of P,
but such approximations will be generally worse than when we approximate a single
real number, and so matters do not work out: we need a more, not less,
precise approximation than before to the leading coefficient, since it
is now the coefficient of a quadratic term.

Concretely, for P(x) = α_2 x^2 + α_1 x + α_0,
in order to reduce the problem of finding solutions to
P(r) ∈ -η,ηℤ with
r∈ R,R∩ℤ to the problem of solving
a quadratic modular equation,
 we would need
a_1, a_2, q with |α_2 - a_2/q|≤ 1/q R^2,
|α_1 - a_1/q|≤ 1/q R. In general, we can do such a thing only for 
q≫ R^3, and that is much too large. For one thing, for ϵ∼ 1/R,
we would have to solve ϵ q ∼ R^2 distinct equations, and that is
obviously too many.

   

 §.§ Geometric interpretation

As we have seen, sieving integers in the interval 
n-Δ,n+Δ reduces to finding integers m such
that {n/m} is close to 0 modulo 1.
This task is equivalent to finding integer points
close to a hyperbola x↦ n/x. Seen from this perspective, our
approach consists simply in approximating a hyperbola locally by linear
functions. Such a geometric perspective is already present (and dominant) in
 <cit.>.

What we did in <ref> then amounts to finding points close
to a line, starting with an approximation of the slope by a rational, and then
proceeding by modular arithmetic. 

Taking one more term in the Taylor expansion would be
 the same as approximating a hyperbola by segments of parabolas, rather than 
by segments of lines. We could then take longer segments, and thus
hope to capture points near the hyperbola efficiently even when 
Δ is considerably smaller than n^1/3. (We can hope to
capture them efficiently because the segments are long enough that
we expect at least one point in each segment.) The problem of how to find
the points remains. It seems difficult to do so without reducing a
non-linear problem 1 to a problem q, and we do not yet know how to carry out such a reduction
well.




§ A FEW WORDS ON THE IMPLEMENTATION


I have written and tested a simple 
implementation as proof-of-purpose, that is, mainly so as to check that the algorithms work as described.
On the higher end of what can be done in 64-bit
computer arithmetic (say: n=5· 10^18 and Δ = 2· 10^7 or
Δ = 4· 10^7),
my implementation of algorithm NewSegSiev
runs substantially faster than my own implementation of the traditional
algorithm SegSiev. Still, on those same inputs, my implementation of
NewSegSiev is clearly slower (by a factor between 2 and
2.5) than a publicly
available, highly optimized implementation <cit.> of the traditional
algorithm (basically SegSiev, but improved in all the ways
detailed below, and some other ones) on the same interval.
Diophantine approximation (Algorithm <ref>)
turns out to take less than 2 percent of
the running time of NewSegSiev, so
the fact that <cit.> runs faster is due to better
coding, and not to the overhead of Diophantine approximation. For the values of
n and Δ above, we are talking about
total running times of only a few seconds in all cases.







The advantage of NewSegSieve over existing methods
should become clearer as n grows larger –
meaning substantially larger than 2^64. 
However, on the range n>2^64,
it seems harder to find highly optimized implementations
of the traditional algorithm for comparison.

Here are a few hints for the reader who would like to write a more serious
program on his or her own. Most of these tricks are standard, but are scattered
here and there in the literature (and in code).



  * Obviously, 
  we can save on space by storing the sieve as a bit array. Saving on space
  leads to better cache usage, and hence often to time savings as well.
  Our very conventions for space usage (expressed in terms of bits, rather
  than words) reflect this fact. Of course, in some ranges, it can save time
  to be a little more wasteful and store some integer data as words
  (e.g., prime factors, in the case of the sieve applied to factorization).

  * We can first apply a simple sieve to the integers between
  1 and M = ∏_p≤ p_0 p (where p_0=17, say), taking only the
  effect of primes p≤ p_0 into account (whether we are sieving for
  primality, or computing the Möbius function
  μ, or factoring numbers:  in the nth entry, we would store whether n
  is coprime to P, or, instead, store μ((n,M^2)), if we
  are computing μ,
  or the set 
  {p≤ p_0: p|n}, if we are factoring numbers). We then initialize
  our sieve by repeating that block of length M, and so we do not need
  to sieve by the primes p≤ p_0 ever again.
  
  
  
  
  
  

































  * We can of course implement a sieve in parallel in a trivial sense,
  by letting different processors look at different intervals (of length
  at least n^3 (log n)^2/3, in our case). There seems to be a small
  literature on parallel implementations of sieves; see, e.g.,
  <cit.>.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  


  * At the end of <ref>, we went briefly over the issue of
  “false alarms”, that is, values of m that lie in a valid congruence
  class m≡ m_0+r_0  q but do not actually have
  multiples within the interval I we are sieving. Such false alarms do not
  change the outcome of the algorithm, but they do waste
  some time.  It is possible to avoid them; the details are given in the
  first version of the present paper, available on arXiv.org. There, in fact,
  the algorithm without “false alarms” is called NewSegSiev,
  and presented as the main version of
  the algorithm.
  
  
  
  Unfortunately, at least in my implementation, the process of eliminating
  false alarms takes more time than it saves.
  
  
  
  

  * As we said in the introduction, Oliveira e Silva has shown
  how to use cache efficiently when implementing a sieve of Eratosthenes
  <cit.>, <cit.>. In effect, when sieving an interval
  of length L, he needs not much more than √(L) units of memory in
  cache. There seems to be no reason why Oliveira e Silva's technique
  can't be combined with the algorithms put forward here. Thus we could
  hope to sieve intervals of the form n-Δ, n+Δ,
  Δ∼ n^1/3 (log n)^2/3, while using no more than
  O(n^1/6 (log n)^1/3) units of memory in cache at a time (and
  O(n^1/3 (log n)^2/3) units of memory in total). The range up to about
  n∼ 10^36 could then become accessible, at least for short or medium-sized intervals. Note that we can still use
  128-bit-integer arithmetic in that range.
  
alpha

