Sieve of Eratosthenes, in depth: a hand trace, the n log log n argument, odd-only and segmented sieves, memory layouts and cache-aware production use

The idea from first principles

Every composite number m has a prime factor no larger than the square root of m: if both factors of m = a x b were larger than the root, their product would exceed m. So if we remove every multiple of every prime up to the square root of n, whatever remains between 2 and n is prime.

The sieve does exactly that, in increasing order. Start with all numbers from 2 to n marked as candidates. Take the smallest unmarked number; it has no smaller prime factor, so it is prime. Cross off its multiples. Repeat with the next unmarked number. Stop when the next prime squared exceeds n.

Two optimisations follow from the same argument. Crossing off for prime p can start at p squared, because every smaller multiple k x p with k < p has a prime factor smaller than p and was already crossed off. And the outer loop can stop at the square root of n, because any composite left above that point would need a prime factor larger than its root.

Worked example: every prime up to 30

The survivors are 2, 3, 5, 7, 11, 13, 17, 19, 23 and 29: ten primes from 24 write operations. Notice that 12, 18, 24 and 30 were crossed twice. That repeated work is what the linear sieve removes, at the cost of storing the smallest prime factor of every number; for listing primes the plain sieve is usually faster in practice because its inner loop is a simple strided store.

Why the cost is n log log n

Prime p crosses off about n / p numbers. The total work is n times the sum of 1/p over primes up to the square root of n, and that sum grows like log log n (a result due to Mertens). Log log n is tiny: about 2.6 for a million and 3.0 for a billion. The sieve is linear for every practical purpose.

Measured, not estimated: for n = 1,000,000 the reference code below performs 2,122,048 crossing-off writes, about 2.1 writes per number. The bound n ln ln n gives 2.6 million, a little above the measured count, partly because starting at p squared skips work the bound counts. For comparison, trial division of every number up to a million, dividing only by primes up to its square root and stopping at the first factor, performs 13,927,402 divisions: about 6.6 times as many operations, each far more expensive than a byte store, and the gap widens as n grows.

If cost models are unfamiliar, Big-O notation in practice covers why constant factors still matter once the asymptotics are this good.

A reference implementation

In Python the inner loop should never be a Python loop. Slice assignment with a step crosses off all multiples of p in one call that runs in C. A bytearray stores one byte per number.

import math

def sieve(n):
    """Return a bytearray where is_prime[k] == 1 iff k is prime, for 0 <= k <= n."""
    is_prime = bytearray([1]) * (n + 1)
    is_prime[0:2] = b"\x00\x00"[: min(2, n + 1)]
    for p in range(2, math.isqrt(n) + 1):
        if is_prime[p]:
            is_prime[p * p :: p] = bytes(len(range(p * p, n + 1, p)))
    return is_prime

The second line handles n = 0 and n = 1 without special cases. On one laptop under CPython, sieve(10**8) took about 1.5 seconds and 100 MB. That is the limit of this version: memory is one byte per number, and every pass over the array streams the whole 100 MB through the cache.

Memory layouts for a billion numbers

Skipping even numbers is the cheapest big win, and the same bit-packing ideas appear in Bloom filters. The odd-only form maps index i to the number 2i + 1. A prime p is odd, so stepping p indices forward moves 2p in value, landing on the next odd multiple and skipping the even ones automatically.

def odd_sieve(n):
    """Odd numbers only: index i stands for 2*i + 1. Half the memory, half the work."""
    if n < 2:
        return []
    size = (n - 1) // 2 + 1                 # indices for 1, 3, 5, ..., up to n
    odd = bytearray([1]) * size
    odd[0] = 0                              # 1 is not prime
    for i in range(1, (math.isqrt(n) - 1) // 2 + 1):
        if odd[i]:
            p = 2 * i + 1
            start = (p * p - 1) // 2        # index of p*p
            odd[start::p] = bytes(len(range(start, size, p)))
    return [2] + [2 * i + 1 for i, f in enumerate(odd) if f]

Wheels generalise the trick: removing multiples of 2, 3 and 5 leaves 8 residues out of every 30, and a bit per residue packs 30 numbers into one byte. Beyond 2, 3, 5 and 7 the returns diminish while the index arithmetic gets harder.

The segmented sieve

The plain sieve touches the whole array for each prime. Once the array is larger than the CPU cache, nearly every crossing-off write for large p is a cache miss, and the algorithm becomes memory-bound. The segmented sieve fixes this, and makes memory independent of n.

First sieve the base primes up to the square root of the upper limit. Then walk the range in windows of S numbers. For each window, cross off multiples of each base prime, starting at the first multiple inside the window, and collect the survivors. Each window is small enough to stay in cache while every base prime passes over it.

def segmented_primes(lo, hi, segment=1 << 18):
    """Yield primes in [lo, hi) using O(sqrt(hi) + segment) memory."""
    base = [p for p, f in enumerate(sieve(math.isqrt(hi - 1) + 1)) if f]
    for seg_lo in range(max(lo, 2), hi, segment):
        seg_hi = min(seg_lo + segment, hi)
        mark = bytearray([1]) * (seg_hi - seg_lo)
        for p in base:
            if p * p >= seg_hi:
                break
            start = max(p * p, (seg_lo + p - 1) // p * p)
            mark[start - seg_lo :: p] = bytes(len(range(start, seg_hi, p)))
        for i, f in enumerate(mark):
            if f:
                yield seg_lo + i

The function sieves any range, not only one starting at zero: there are 36,249 primes between 10¹² and 10¹² + 10⁶, and finding them needs only the 78,498 base primes below a million plus one window. Listing primes near 10¹² with a plain sieve would need a terabyte.

In CPython this version is slower than the plain one, about 6.5 seconds to 10⁸ on the same laptop, because the survivor loop yields one Python object per number. That is an interpreter cost, not an algorithmic one. In C, C++ or Rust, where the survivor scan is a tight loop, a segmented sieve with a window sized to the L1 or L2 data cache is typically several times faster than the plain sieve at large n, and windows can be handed to separate threads with no shared writes. The primesieve library by Kim Walisch is the reference implementation of these ideas if you need production speed rather than your own code.

Testing a sieve

Sieve bugs are almost always off-by-one errors at the edges: the limit itself, the first window, the first multiple inside a window. Test against independent implementations and against known prime counts.

for n in (0, 1, 2, 3, 4, 25, 10**6):
    a = [i for i, f in enumerate(sieve(n)) if f]
    assert a == odd_sieve(n) == list(segmented_primes(0, n + 1, 1000))
assert sum(sieve(10**6)) == 78498
assert len(odd_sieve(10**7)) == 664579
assert sum(1 for _ in segmented_primes(10**12, 10**12 + 10**6)) == 36249

Known values of the prime counting function are the best oracles: there are 78,498 primes below 10⁶, 664,579 below 10⁷ and 5,761,455 below 10⁸. Include tiny limits (0, 1, 2, 3, 4) and perfect squares of primes such as 25, since those exercise the stopping condition.

Failure modes

• Excluding the limit. Allocating n cells instead of n + 1 silently drops n when it is prime.
• Overflow of p squared. In 32-bit integer languages, p * p overflows for p above 46,340, and the loop starts at a negative index or never terminates. Use 64-bit arithmetic or compare p with the integer square root.
• Wrong first multiple in a window. Starting at ceil(lo/p)*p without the p*p floor crosses off p itself when p lies inside the window.
• Using a sieve for one big number. To test whether a single 60-digit number is prime, a sieve is useless; use Miller-Rabin.
• Allocating the full range by habit. A service that sieves to 10¹⁰ per request will be killed by the memory limit before it is slow.
• Rebuilding on every call. If many requests ask for primes below the same limit, sieve once at start-up and share the read-only table.

Choosing the right tool

What to do next

1. Implement the plain sieve from memory and check it against the trace to 30.
2. Add the odd-only version and confirm both return 78,498 primes below one million.
3. Implement the segmented version and test it on a window that does not start at zero.
4. Run the test block, including limits 0 to 4 and 25.
5. Port the segmented sieve to a compiled language, vary the window size from 16 KB to 4 MB, and plot run time against window size to find your cache sweet spot.
6. Replace any per-request sieve in your code with a shared table or a Miller-Rabin check, whichever matches the access pattern.