Sieve of Eratosthenes: Start at p Squared, Stop at the Square Root

The sieve of Eratosthenes is more than two thousand years old and still the right answer whenever you need all the primes up to some limit. Instead of asking "is 91 prime?" it lets the primes it has already found cross out everything they can reach, and whatever is left standing is prime by construction. Two details make it fast, and they are exactly what interviewers ask about: why the crossing out starts at p * p, and why the outer loop can stop at the square root of n.

The problem it solves

Given n, list every prime up to n, or count them. The classic interview form is "count the primes less than n", and it also hides inside other questions: sum the primes in a range, check many numbers for primality, factorize many numbers quickly.

Testing each number separately with trial division checks divisors up to its square root, which costs about n√n in total for the whole range. For n in the millions that is billions of divisions. The sieve does the same job in close to linear time and uses no division at all.

The intuition

Write down every number from 2 to n, all unmarked. Then repeat:

• Take the smallest number nothing has crossed out. It is prime, because any smaller factor would be a smaller prime, and that prime would already have crossed it out. The proof is free.
• Cross out its multiples. Every one of them has this prime as a factor, so none of them is prime.

Two refinements turn that into the fast version:

• Start at p * p. A multiple of p smaller than p * p is p * m with m less than p, so it has a factor smaller than p, and a smaller prime already crossed it out. For 3, the multiple 6 was taken by 2, so 3 starts at 9. For 5, only 25 is left on a board that ends at 31.
• Stop when p * p passes n. Every composite number up to n has a factor no larger than its square root. Once p is past the square root of n, there is nothing left to cross out, and every unmarked number is prime.

The work per prime keeps shrinking: n / 2 strikes for 2, fewer than n / 3 for 3, and so on. Adding up n / p over the primes gives the n log log n total, which grows so slowly that it behaves almost like n.

Watch it run

The animation lays out thirty numbers, 2 to 31, six to a row. Testing each for divisors means many divisions, so it tests none of them. The first number nothing has crossed is 2, so 2 is prime, and every second cell goes: half the board is gone without a single division. The next untouched cell is 3, and untouched means no smaller number divides it. It starts at 9, because 6 was already taken by 2 and anything below 9 has a smaller factor that has already struck it. 5 survives, so 5 is prime, and only 25 is left for it to strike. The next candidate is 7, but 49 is off the board: past the square root of 31 there is nothing left to strike. Everything still standing is prime, eleven of them, found by crossing out rather than by testing.

The code

The lesson's sieve, with a counter for how many times a cell is crossed out:

import math

def primes_up_to(n):
    ok = [True] * (n + 1)
    ok[0] = ok[1] = False
    strikes = 0
    for p in range(2, int(n ** 0.5) + 1):
        if ok[p]:                            # p survived, so p is prime
            for k in range(p * p, n + 1, p): # below p * p, a smaller prime struck it
                ok[k] = False
                strikes += 1
    return [i for i, good in enumerate(ok) if good], strikes

primes, strikes = primes_up_to(31)
print(primes, strikes)       # [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31] 24
print(len(primes_up_to(100)[0]))                    # 25

for n in (10 ** 3, 10 ** 4, 10 ** 5):
    s = primes_up_to(n)[1]
    print(n, s, round(s / n, 2), round(math.log(math.log(n)), 2))
# 1000 1411 1.41 1.93
# 10000 16981 1.7 2.22
# 100000 193078 1.93 2.44

Strikes per number grow in step with log log n, a constant below it. The same algorithm with a byte array and slice assignment, which is how you would write it in Python for a million:

def sieve_fast(n):
    if n < 2:
        return []
    ok = bytearray([1]) * (n + 1)
    ok[0] = ok[1] = 0
    for p in range(2, math.isqrt(n) + 1):
        if ok[p]:
            ok[p * p::p] = bytes(len(range(p * p, n + 1, p)))
    return [i for i in range(n + 1) if ok[i]]

print(len(sieve_fast(10 ** 6)))                     # 78498

A common follow-up is to factorize many numbers. Store each number's smallest prime factor instead of a yes or no, and factorizing becomes repeated division by a table lookup:

def smallest_factors(n):
    spf = list(range(n + 1))
    for p in range(2, math.isqrt(n) + 1):
        if spf[p] == p:                              # p is prime
            for k in range(p * p, n + 1, p):
                if spf[k] == k:
                    spf[k] = p
    return spf

spf = smallest_factors(1000)

def factorize(x):
    out = []
    while x > 1:
        out.append(spf[x])
        x //= spf[x]
    return out

print(factorize(360), factorize(997))               # [2, 2, 2, 3, 3, 5] [997]

All three against trial division, on random limits up to 3,000 and every number up to 1,000:

import random

def is_prime(x):
    if x < 2:
        return False
    d = 2
    while d * d <= x:
        if x % d == 0:
            return False
        d += 1
    return True

random.seed(18)
ok = True
for _ in range(300):
    n = random.randint(0, 3000)
    want = [x for x in range(n + 1) if is_prime(x)]
    ok &= primes_up_to(n)[0] == want
    ok &= sieve_fast(n) == want
for x in range(2, 1001):
    f = factorize(x)
    ok &= math.prod(f) == x and all(is_prime(p) for p in f) and f == sorted(f)
print(ok)                                           # True

The complexity

• Time: O(n log log n). The strikes for each prime p number about n / p, and the sum of 1 / p over primes up to n grows like log log n.
• Space: O(n) for the table, one flag per number. A byte array keeps it to about a megabyte for a million.
• Trial division for every number: about O(n√n). Worth it only for a few isolated checks.
• Very large ranges that do not fit in memory use a segmented sieve: sieve the primes up to the square root once, then cross out one block of the range at a time.

Where it goes wrong

• Off-by-one on the limit. "Primes less than n" and "primes up to n" differ by whether n itself is included. Size the table for the one you were asked.
• Stopping the outer loop too early. range(2, int(n ** 0.5)) misses the square root itself; for n = 25 it never strikes 25. Add one, or use math.isqrt(n) + 1.
• Starting at 2 * p. Still correct, just slower: it re-strikes numbers smaller primes already handled.
• Forgetting 0 and 1. Neither is prime, and a table that starts all true reports them.
• Sieving to answer one question. For a single number, trial division up to its square root is simpler.

When it shows up in interviews

"Count primes less than n" is a well-known medium, and the sieve appears inside many number-theory questions. Interviewers ask why the inner loop starts at p * p, why the outer loop can stop at the square root, and what the complexity is. A good extension is the smallest-prime-factor table for fast factorization, or the segmented sieve when n is too large for memory.

How to say it in an interview

"I keep a boolean table from 0 to n, with 0 and 1 marked not prime. I walk p upward from 2. If p is still unmarked it is prime, because any smaller factor would already have crossed it out, so I cross out its multiples starting at p squared, since smaller multiples have a smaller factor and are already gone. I can stop when p squared passes n, because every composite up to n has a factor no bigger than its square root. What is left unmarked is prime. That is O(n log log n) time and O(n) space."

The divisibility facts behind it are in GCD and Euclid, and another way to answer many questions from one precomputed table is prefix sums.