Sieve of Eratosthenes — Prime Number Generation

The Sieve of Eratosthenes is 2300 years old and still the fastest way to generate all primes up to n in most practical scenarios. That longevity is not nostalgia — it is a testament to how completely it solves the problem. In competitive programming, you hit this algorithm in the first 10 minutes of almost any number theory problem. In production, you hit it when building RSA key generation utilities, cryptographic prime tables, or anything involving primality at scale.

What makes it elegant is not the algorithm itself — it's obvious once explained — but the complexity. O(n log log n) time is almost linear. For n=10^7, a correctly implemented sieve runs in under 100ms in Python and under 5ms in C. Compare that to trial division on each number: O(n√n) = O(10^10.5) for n=10^7 — three orders of magnitude slower. The choice is not subtle.

The Basic Sieve

The implementation is four lines of meaningful code. The subtlety is in where the inner loop starts.

Starting the inner loop at pp instead of 2p is the critical optimisation. Every composite number smaller than p² already has a prime factor smaller than p and has therefore already been marked. When we reach prime p=5, multiples 10, 15, 20 were already marked by 2 and 3. The first unmarked multiple of 5 is 5×5=25. This halves the total work.

The second optimisation is memory layout. A Python list of booleans uses 28 bytes per element (Python object overhead). Use a bytearray or numpy bool array instead — 28x less memory and dramatically better cache performance. For n=10^7, a Python list uses 268MB; a bytearray uses 10MB. Cache misses dominate the runtime at large n, so this is not premature optimisation.

def sieve(n: int) -> list[int]:
    """Return all primes <= n using Sieve of Eratosthenes."""
    is_prime = [True] * (n + 1)
    is_prime[0] = is_prime[1] = False
    p = 2
    while p * p <= n:
        if is_prime[p]:
            # Mark multiples of p starting from p*p
            for multiple in range(p * p, n + 1, p):
                is_prime[multiple] = False
        p += 1
    return [i for i, prime in enumerate(is_prime) if prime]

print(sieve(50))
# [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47]
print(f'Primes up to 100: {len(sieve(100))}')  # 25

Why Start Marking at p²?

The harmonic series of primes 1/2 + 1/3 + 1/5 + 1/7 + ... grows as log(log(n)) — this is Mertens' theorem from 1874. Each prime p contributes n/p marking operations. Summing over all primes up to √n gives the total work: n × Σ(1/p) ≈ n × log(log(n)).

The practical implication: for n=10^9, log(log(n)) ≈ 3.04. The sieve does roughly 3× the work of a simple O(n) pass. This is why the sieve feels "almost free" in practice — the constant is tiny.

Memory limits in production: At n=10^9, a bit-packed sieve uses 125MB. A bytearray uses 1GB. Bit packing (using Python's bitarray library or numpy with uint8) is the practical choice above n=10^8. The segmented sieve reduces this further — work on √n-sized chunks, keeping only the current segment in cache at any time.

Segmented Sieve — Large Ranges

The segmented sieve is what you reach for when n doesn't fit in RAM or when you need primes in a specific range [L, R] where R is large but R-L is manageable.

Real scenario: Finding all primes in [10^12, 10^12 + 10^6]. You obviously cannot sieve up to 10^12. But √(10^12) = 10^6 — you can sieve all primes up to 10^6 (trivial: ~78,498 primes), then use them to sieve the segment. The segment has only 10^6 numbers — fits easily in L1/L2 cache.

This is the pattern used in distributed prime-finding projects (like primegrid.com) and in production RSA implementations that need to test primality of large numbers using a quick sieve pre-filter before Miller-Rabin.

import math

def segmented_sieve(n: int) -> list[int]:
    """Memory-efficient sieve for large n using O(sqrt(n)) space."""
    limit = int(math.sqrt(n)) + 1
    base_primes = sieve(limit)
    primes = list(base_primes)
    low, high = limit + 1, min(2 * limit, n)
    while low <= n:
        segment = [True] * (high - low + 1)
        for p in base_primes:
            start = max(p * p, ((low + p - 1) // p) * p)
            for j in range(start, high + 1, p):
                segment[j - low] = False
        for i, is_p in enumerate(segment):
            if is_p:
                primes.append(low + i)
        low = high + 1
        high = min(high + limit, n)
    return primes

print(len(segmented_sieve(1_000_000)))  # 78498

Smallest Prime Factor Sieve

A variant that stores the smallest prime factor (SPF) for each number instead of just prime/composite. Enables O(log n) factorisation of any number after O(n log log n) preprocessing.

def spf_sieve(n: int) -> list[int]:
    """Smallest prime factor sieve."""
    spf = list(range(n + 1))  # spf[i] = i initially
    p = 2
    while p * p <= n:
        if spf[p] == p:  # p is prime
            for multiple in range(p * p, n + 1, p):
                if spf[multiple] == multiple:
                    spf[multiple] = p
        p += 1
    return spf

def factorise(n: int, spf: list) -> list[int]:
    factors = []
    while n > 1:
        factors.append(spf[n])
        n //= spf[n]
    return factors

spf = spf_sieve(100)
print(factorise(84, spf))   # [2, 2, 3, 7]
print(factorise(360, spf))  # [2, 2, 2, 3, 3, 5]

Frequently Asked Questions