Sieve of Eratosthenes in Python: Every Prime Under a Million, Fast

There are 78,498 prime numbers below one million, and a laptop can produce every one of them in a few milliseconds. It can only do that if you stop asking the obvious question.

The obvious question is "is this number prime?", asked once per number. Answer it by trial division and listing the primes below a million costs about 68 million modulo operations. The sieve of Eratosthenes asks a different question — "which numbers can I rule out?" — and answers it with 2.1 million writes and not a single division. Same list of primes, thirty-two times less work, and the gap widens with every extra digit.

Eratosthenes ran the library at Alexandria in the third century BC and is better known for measuring the circumference of the Earth. His method for primes is 2,200 years old and is still what production code uses when it needs every prime up to a bound. What follows proves the two optimisations that most implementations write without justifying, derives the O(n log log n) running time from the sum of the reciprocals of the primes rather than asserting it, and covers the two variants you will actually reach for: a segmented sieve for ranges too large to hold in memory, and a smallest-prime-factor table that turns factorisation into a lookup.

The idea

Write out every whole number from 1 to n. Cross out 1 straight away — a prime has exactly two distinct divisors and 1 has one, so it does not qualify.

Now repeat two steps until you run out of work:

1. Find the smallest number still standing. It is prime.
2. Cross out every multiple of it.

Step 1 is the part that deserves an argument, and the argument is three lines. Suppose k is the smallest number above 1 that is still standing. Every prime below k was picked in an earlier round, and each of those rounds crossed out all of its own multiples. So k is not a multiple of any prime smaller than itself. A composite number always has a prime factor smaller than itself. Therefore k is not composite.

Notice what step 2 does not involve. Crossing out the multiples of 7 means visiting 14, 21, 28, 35 — repeated addition with a stride of 7. There is no division anywhere in the algorithm and no test of any individual number. Trial division asks n separate questions, each costing up to the square root of n in divisions; the sieve replaces all of them with a handful of straight walks through an array.

Watching it work

Take the numbers 1 to 30. Cross out 1, then start.

Round 1 — the smallest survivor is 2, so 2 is prime. Cross out 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30. That is 14 cells written.

Round 2 — the smallest survivor is now 3. Cross out 9, 12, 15, 18, 21, 24, 27, 30 — 8 cells, and half of them (12, 18, 24, 30) were already crossed out by 2. The sieve does not check first. Reading a cell to find out whether it is already marked costs at least as much as writing it, and writing "already crossed out" again is harmless.

Round 3 — the smallest survivor is 5. Cross out 25 and 30. Only 25 is new; 10, 15 and 20 went in earlier rounds and 30 went twice already.

Round 4 — the smallest survivor is 7, and this is where you stop. 7 × 7 = 49 is already past 30, and every multiple of 7 that is 30 or less is 7 times something smaller than 7, so it carries a factor below 7 and went long ago. The same holds for 11, 13 and every later survivor: nothing is left to mark, so everything still standing is prime.

Ten primes — 2, 3, 5, 7, 11, 13, 17, 19, 23, 29 — for 24 writes and zero divisions.

The six-column layout is not decorative: every prime past 3 lands in the first or last column, because the other four columns are exactly the numbers divisible by 2 or by 3. Past 5, only 8 of every 30 consecutive numbers can be prime at all — the observation that wheel factorisation uses to shrink a sieve's memory.

The code

Start with the baseline, so the improvement is measurable rather than asserted.

from math import isqrt, log


def is_prime_by_trial_division(number: int) -> bool:
    """True if number is prime, decided by dividing by every candidate up to its square root."""
    if number < 2:
        return False
    divisor = 2
    # If number = a * b with a <= b, then a * a <= number, so any divisor above
    # the square root can only appear paired with one below it.
    while divisor * divisor <= number:
        if number % divisor == 0:
            return False
        divisor += 1
    return True


def primes_by_trial_division(limit: int) -> list[int]:
    """Every prime below limit, decided one number at a time."""
    return [number for number in range(2, limit) if is_prime_by_trial_division(number)]


print(primes_by_trial_division(30))
print(len(primes_by_trial_division(1000)))

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
168

That is correct and it is the right tool when you have exactly one number to test. As a way to produce a whole list it throws away everything it learns: the moment it discovers that 4 is even, it forgets, and asks the same question again at 6.

Here is the sieve. limit is exclusive, matching range, so sieve_of_eratosthenes(30) returns the primes below 30.

def sieve_of_eratosthenes(limit: int) -> list[int]:
    """Return every prime strictly below limit, in ascending order.

    Rather than testing numbers, this marks composites: each prime it finds
    strikes off that prime's own multiples in a single sweep of additions.
    """
    if limit <= 2:
        return []

    is_prime = [True] * limit
    is_prime[0] = is_prime[1] = False

    candidate = 2
    # Every composite below limit has a prime factor no larger than its own
    # square root, so once candidate squared reaches limit nothing is left to mark.
    while candidate * candidate < limit:
        if is_prime[candidate]:
            # Multiples of candidate below candidate squared already carry a
            # smaller prime factor and were struck off on an earlier sweep.
            for multiple in range(candidate * candidate, limit, candidate):
                is_prime[multiple] = False
        candidate += 1

    return [number for number, prime in enumerate(is_prime) if prime]


print(sieve_of_eratosthenes(30))
print(sieve_of_eratosthenes(2), sieve_of_eratosthenes(3), sieve_of_eratosthenes(0))
print(sieve_of_eratosthenes(10_000) == primes_by_trial_division(10_000))

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
[] [2] []
True

Twenty-two lines, and the last one checks it against the baseline on every number below ten thousand.

In CPython the version above spends nearly all its time in the interpreter, one bytecode dispatch per write. Because the marking pattern is an arithmetic progression, you can hand the entire inner loop to C with a single slice assignment.

def sieve_flags(limit: int) -> bytearray:
    """The same sieve with the marking loop handed to C by slice assignment."""
    if limit <= 2:
        return bytearray(max(limit, 0))

    flags = bytearray(b"\x01") * limit
    flags[0:2] = b"\x00\x00"

    for candidate in range(2, isqrt(limit - 1) + 1):
        if flags[candidate]:
            start = candidate * candidate
            # bytearray(k) is k zero bytes, so one assignment clears the whole
            # arithmetic progression without a Python-level loop.
            flags[start::candidate] = bytearray(len(range(start, limit, candidate)))

    return flags


flags = sieve_flags(1_000_000)
print(sum(flags), flags[999_983], flags[999_984])

78498 1 0

Identical algorithm, identical number of writes, roughly forty times faster on a million-element sieve, and one byte per number instead of the eight a list of pointers costs. Converting the flags back into a list of primes then costs more than the sieve itself — if the question is "is 786,433 prime?", keep the flags and index into them.

How the code maps to the idea

The array is the sheet of paper, and the index is the number. is_prime[17] answers a question about 17: no hashing, no search, no comparison, just an offset into a contiguous block of memory. That is the decision the whole algorithm rests on, and it is also where the O(n) space cost comes from.

is_prime[0] = is_prime[1] = False handles the two numbers the loop never visits. Neither is prime and nothing crosses them out.

The outer loop stops at the square root. Any composite m can be written as a * b with a <= b, which forces a * a <= m, so m has a prime factor no larger than the square root of m — and therefore smaller than the square root of limit. Some earlier round has already crossed m out. Past that point every remaining round would sweep an empty range.

if is_prime[candidate] skips composites, and loses nothing by doing so. A composite candidate has a prime factor q smaller than itself, and every multiple of candidate is also a multiple of q, so q already dealt with all of them.

The inner loop starts at candidate * candidate, not 2 * candidate. Take any earlier multiple k * candidate with 2 <= k < candidate. Then k has a prime factor q with q <= k < candidate, so q ran in an earlier round, marking multiples of q from q * q upwards. And k * candidate is at least q * candidate, which exceeds q * q because candidate exceeds q. So it was already struck off: skipping it is correct, not merely convenient.

The inner loop step is candidate, which keeps the walk sequential in memory and free of arithmetic beyond one addition per cell.

Edge cases fall out of the single guard at the top. Limits of 0, 1 and 2 return an empty list before any array is allocated; a limit of 3 clears the first two flags, never enters the outer loop because 2 × 2 is not below 3, and answers [2].

Now count what the two approaches actually do, rather than trusting the description.

def count_trial_division(limit: int) -> int:
    """Modulo tests trial division performs to list every prime below limit."""
    tests = 0
    for number in range(2, limit):
        divisor = 2
        while divisor * divisor <= number:
            tests += 1
            if number % divisor == 0:
                break
            divisor += 1
    return tests


def count_sieve(limit: int, square_start: bool = True, stop_at_root: bool = True) -> int:
    """Cells the sieve writes to, with either optimisation switchable off."""
    writes = 0
    is_prime = [True] * limit
    is_prime[0] = is_prime[1] = False
    outer_bound = isqrt(limit - 1) + 1 if stop_at_root else limit
    for candidate in range(2, outer_bound):
        if is_prime[candidate]:
            start = candidate * candidate if square_start else 2 * candidate
            for multiple in range(start, limit, candidate):
                is_prime[multiple] = False
                writes += 1
    return writes


print(f"{'limit':>9}  {'modulo tests':>14}  {'sieve writes':>14}  {'ratio':>6}")
for limit in (1_000, 10_000, 100_000):
    tests = count_trial_division(limit)
    writes = count_sieve(limit)
    print(f"{limit:>9,}  {tests:>14,}  {writes:>14,}  {tests / writes:>6.1f}")

    limit    modulo tests    sieve writes   ratio
    1,000           5,287           1,409     3.8
   10,000         117,526          16,979     6.9
  100,000       2,745,693         193,076    14.2

The ratio roughly doubles for every tenfold increase in n, which is what you expect when one side grows like the square root of n and the other barely grows at all. Extend the same measurement to a million and trial division needs 67,740,403 modulo tests against the sieve's 2,122,046 writes — and a modulo is a much more expensive machine instruction than a byte store.

Complexity

Time: O(n log log n). Here is where that comes from.

For each prime p with p * p < n, the inner loop writes to the multiples of p from p * p up to n. That is (n - p * p) / p + 1 writes, which is a little under n / p. Add it up over every prime the outer loop reaches:

total writes  ≈  n × (1/2 + 1/3 + 1/5 + 1/7 + ... + 1/p)   for primes p ≤ √n

So the entire running time hangs on one quantity: the sum of the reciprocals of the primes up to a bound. That sum diverges, but astonishingly slowly. Mertens proved in 1874 that

1/2 + 1/3 + 1/5 + ... + 1/p  =  ln ln x + M + o(1)     for primes p ≤ x

where M ≈ 0.2615 is the Meissel-Mertens constant. The double logarithm is not a coincidence, and you can see where it comes from without the full proof. The prime number theorem says there are about x / ln x primes below x, which makes the k-th prime roughly k ln k in size. Summing reciprocals of the primes is then summing 1 / (k ln k) over k, and that sum behaves like the integral of 1 / (k ln k), which is ln ln k. One logarithm comes from the thinning-out of the primes; the other comes from integrating a reciprocal.

Substituting x = √n gives ln ln √n = ln(ln n / 2) = ln ln n - ln 2, so the count of writes is about n × (ln ln n - ln 2 + M), and dropping the constants leaves O(n log log n).

That is a claim you can check.

MEISSEL_MERTENS = 0.2614972128


def predicted_writes_per_number(limit: int) -> tuple[float, float]:
    """(plain n log log n estimate, same estimate corrected for the square start)."""
    base_primes = sieve_of_eratosthenes(isqrt(limit - 1) + 1)
    # Mertens: the sum of 1/p over primes p <= x is ln ln x + M.
    plain = log(log(limit)) - log(2) + MEISSEL_MERTENS
    # The square start skips the first p - 1 multiples of every base prime p.
    corrected = plain - (sum(base_primes) - len(base_primes)) / limit
    return plain, corrected


print(f"{'n':>9}  {'measured':>8}  {'ln ln n - ln 2 + M':>19}  {'corrected':>9}")
for limit in (10_000, 100_000, 1_000_000):
    measured = count_sieve(limit) / limit
    plain, corrected = predicted_writes_per_number(limit)
    print(f"{limit:>9,}  {measured:>8.3f}  {plain:>19.3f}  {corrected:>9.3f}")

        n  measured   ln ln n - ln 2 + M  corrected
   10,000     1.698                1.789      1.685
  100,000     1.931                2.012      1.920
1,000,000     2.122                2.194      2.118

The last column matches the measurement to within a fifth of a percent at a million. The formula is not an analogy; it is the count.

It is also worth seeing how small log log n is. At n = 1,000,000 the sieve does 2.12 writes per number. At n = 1,000,000,000,000 it does about 3.4. For any n you will ever sieve, log log n sits between 1 and 4 — which is why people describe the sieve as "basically linear" and are not being sloppy when they do.

The baseline it beats: O(n √n), or better in practice. Trial division tests each of n numbers with up to √n divisions, giving the plain O(n √n) bound. The measured cost is lower, because composites exit at their first divisor and only the primes pay the full square root — about n / ln n of them, for a true cost near n^1.5 / ln n. Either way it is a power of n, and the sieve is a logarithm of a logarithm.

Neither of the two famous optimisations changes the complexity class. They are constant-factor wins, and it is worth knowing how large.

print("writes at n = 1,000,000, by which optimisations are switched on")
print(f"  both:                 {count_sieve(1_000_000):>9,}")
print(f"  no square start:      {count_sieve(1_000_000, square_start=False):>9,}")
print(f"  neither:              {count_sieve(1_000_000, False, False):>9,}")

writes at n = 1,000,000, by which optimisations are switched on
  both:                 2,122,046
  no square start:      2,197,837
  neither:              2,775,208

Starting at p * p instead of 2 * p saves 3.6% of the writes — real, but far less than its fame suggests, and the share shrinks as n grows. Running the outer loop all the way to n while starting at 2 * p costs 31% more, because you are then summing 1/p over every prime below n rather than below √n, and the difference between those two sums is exactly ln 2. If you keep the square start, stopping at the square root saves no writes at all — the ranges past that point are empty — it saves the loop overhead of visiting the other n cells.

Space: O(n). One flag per number, whether or not you ever look at it. With a bytearray that is one byte per number: 1 MB for a million, 1 GB for a billion. A bit array cuts it eightfold and a wheel that skips multiples of 2, 3 and 5 cuts it by another factor of 30/8, but the shape of the cost does not change. This is the sieve's real limitation, and it is the reason the next section exists.

Two variants worth the extra code

Sieving a range that will not fit in memory

Suppose you want the primes just above one trillion. A flat sieve would need a trillion cells. But to sieve any block of numbers up to high, all you need are the primes up to √high — because that is the largest prime factor a composite in that range can be forced to have. For a trillion that is the 78,498 primes below a million, which fit comfortably in memory.

So: sieve the small primes once, then sweep the target range one block at a time, marking each block with those same primes.

def segmented_sieve(low: int, high: int) -> list[int]:
    """Every prime in the half-open range [low, high), one block at a time.

    Memory holds the primes up to the square root of high plus a single block,
    so high may be far beyond anything a flat sieve could store.
    """
    if high <= 2:
        return []
    low = max(low, 2)
    base_primes = sieve_of_eratosthenes(isqrt(high - 1) + 1)

    found: list[int] = []
    block_size = 1 << 16
    for block_start in range(low, high, block_size):
        block_stop = min(block_start + block_size, high)
        flags = bytearray(b"\x01") * (block_stop - block_start)

        for prime in base_primes:
            # First multiple of prime at or above block_start, but never below
            # prime squared, or a base prime would strike itself out.
            first = max(prime * prime, -(-block_start // prime) * prime)
            for multiple in range(first, block_stop, prime):
                flags[multiple - block_start] = 0

        found.extend(block_start + offset for offset, flag in enumerate(flags) if flag)
    return found


print(segmented_sieve(0, 30) == sieve_of_eratosthenes(30))
print(segmented_sieve(1_000_000, 1_000_100))
print(len(segmented_sieve(10**12, 10**12 + 10_000)))

True
[1000003, 1000033, 1000037, 1000039, 1000081, 1000099]
335

There are 335 primes in the ten thousand numbers starting at a trillion, found without ever allocating more than 65,536 flags plus the base primes. Memory is O(√high + block size). The max(prime * prime, ...) term is the same square-start rule as before, and it is what stops the block containing 2 from crossing 2 out.

Two costs to be aware of. Every block pays a pass over all 78,498 base primes even when most of them have no multiple inside it, so very small blocks are wasteful — 32 KB to 256 KB is the usual choice, sized to fit the CPU cache, which is also where the speed comes from. And you cannot ask for primes near 10^30 this way: √high grows too, and the base sieve becomes the bottleneck.

Sieving the smallest prime factor

Change one thing: instead of storing a flag, store the smallest prime that divides each index. The loop is the same shape, plus one condition so that the first prime to reach a cell keeps it.

def smallest_prime_factor_sieve(limit: int) -> list[int]:
    """spf[number] is the smallest prime dividing number, for every number below limit."""
    spf = list(range(limit))
    candidate = 2
    while candidate * candidate < limit:
        if spf[candidate] == candidate:  # nothing smaller divides it, so it is prime
            for multiple in range(candidate * candidate, limit, candidate):
                if spf[multiple] == multiple:  # not yet claimed by a smaller prime
                    spf[multiple] = candidate
        candidate += 1
    return spf


def factorise(number: int, spf: list[int]) -> list[int]:
    """Prime factors of number with multiplicity. Needs 1 <= number < len(spf)."""
    factors: list[int] = []
    while number > 1:
        prime = spf[number]
        factors.append(prime)
        number //= prime  # peel off one copy of that prime and continue
    return factors


spf = smallest_prime_factor_sieve(1_000_000)
print(factorise(999_999, spf))
print(factorise(999_983, spf))
print(factorise(360, spf), factorise(1, spf))

[3, 3, 3, 7, 11, 13, 37]
[999983]
[2, 2, 2, 3, 3, 5] []

A cell keeps the first prime that reaches it, and rounds run in increasing order of prime, so that first arrival is the smallest prime factor. spf[number] == number is then exactly the test for primality, which is why the table also tells you 999,983 is prime.

Factorising afterwards is O(log n): every step divides the number by a prime of at least 2, so it cannot take more than log2(n) steps — at most 20 for a number below a million. Trial-division factorisation of the same number takes up to √n divisions, so about a thousand. If you need to factorise many numbers below a fixed bound, this table is the answer, and the same sweep can build Euler's totient, the Möbius function or a divisor count at no extra asymptotic cost.

The price is memory: a Python list of a million distinct integers is roughly 40 MB, against 1 MB for the bytearray of flags. array("i", range(limit)) from the standard library brings it down to 4 MB.

When to use it, and when not to

Use it when you need every prime below a fixed bound, or a per-number table derived from factorisation, and that bound fits in memory. This is precomputation: pay O(n log log n) once at startup, then answer any number of "is it prime" or "factorise it" questions in constant or logarithmic time.

Do not use it to test a single number, especially a large one. Asking whether one 60-digit number is prime by sieving is not slow, it is impossible — the array would need more cells than there are atoms in the observable universe. The right tool is the Miller-Rabin test, which is a handful of modular exponentiations (see fast exponentiation for the log-time squaring trick it depends on) and is deterministic for every 64-bit input with a fixed set of seven bases. For one modest number, trial division by 2, 3 and then 6k ± 1 is fine and needs no memory at all.

Do not use a flat sieve for a narrow window far out. Primes between 10^12 and 10^12 + 10^4 want the segmented sieve above, which costs O(√high) memory instead of O(high).

Be careful with the "linear sieve". There is a well-known O(n) variant that marks each composite exactly once by pairing it with its smallest prime factor. The bound is genuinely better, but in practice it is often slower than the classic sieve, because it accesses memory out of order and loses the sequential cache behaviour that makes the plain version fast. Measure before you adopt it.

Python's standard library has no prime function at all. math.isqrt and math.gcd are there; primality and factorisation are not. If you would rather not maintain this yourself, SymPy provides sympy.sieve, sympy.primerange, sympy.isprime and sympy.factorint.

Where it shows up in the real world

RSA and Diffie-Hellman key generation. Producing a 2048-bit RSA key means finding two 1024-bit primes, and the only way to find one is to test random odd candidates until one passes. A Miller-Rabin round on a 1024-bit number is a full modular exponentiation — thousands of times more expensive than a small-division check. So implementations sieve first. OpenSSL ships a static table of the first 2,048 primes (the largest is 17,863), and its candidate generator computes the candidate's remainder against every one of them, then walks forward through candidates by updating those remainders — a segmented sieve over the candidate window. By Mertens' third theorem, that table alone eliminates about 88% of odd candidates before a single expensive test runs.

Competitive programming and Project Euler. "Find the sum of all the primes below two million" is Project Euler problem 10, and it is a two-line answer with a sieve and an hour of waiting without one. Contest solutions routinely sieve to 10^6 or 10^7 during setup and then answer thousands of queries by lookup; the smallest-prime-factor table is the standard way to factorise a hundred thousand inputs inside a time limit.

Hash table sizing. The GNU C++ standard library picks the bucket count for unordered_map from a hard-coded list of primes, growing to the next entry when the load factor is exceeded. A prime bucket count keeps a poorly distributed hash function from collapsing onto a few buckets. That list was generated offline by a sieve.

Prime counting and number theory software. Kim Walisch's primesieve — a heavily optimised segmented, wheel-factorised sieve — is the standard tool for generating or counting primes anywhere in the 64-bit range, and is what research code reaches for instead of writing its own.

Common mistakes

Starting the inner loop at 2 * candidate and calling it a bug. It is not a bug, it is 3.6% more writes. The actual bug is starting at 2 * candidate and running the outer loop to n, which costs 31% more and is what most first attempts do.

Testing before writing. if is_prime[multiple]: is_prime[multiple] = False looks like it avoids work. It adds a read to every write and avoids nothing, because the write is idempotent.

Using a set of composites instead of a flat array. A set turns each mark into a hash computation and a probe into scattered memory, and costs an order of magnitude more space. The sieve is fast precisely because the index is the number and the writes are sequential.

Off-by-one at the limit. Decide whether your limit is inclusive or exclusive and put it in the docstring. range-style exclusive is the convention above; sieve_of_eratosthenes(30) therefore does not consider 30 itself.

Forgetting 0 and 1. [True] * limit claims both are prime. Clear them explicitly; nothing in the loop ever visits them.

Sieving inside a loop. Building a fresh sieve every time you need a primality check turns an O(n log log n) precomputation into a per-query cost. Sieve once, keep the array, index into it.

Practice

1. Sum every prime below two million, and check that the total is what you would expect from a table of prime sums.
2. Count the twin primes below one million — pairs p and p + 2 that are both prime.
3. Extend the sieve so that it also records, for every number below the limit, how many distinct prime factors it has.
4. Write nth_prime(k) by sieving up to the bound k * (ln k + ln ln k), which is known to exceed the k-th prime for k of at least 6, then indexing the result.
5. Use the segmented sieve to find the first gap of at least 100 between consecutive primes, and report where it starts.

Summary

The sieve of Eratosthenes wins by inverting the question. Testing numbers costs a square root per number; ruling out multiples costs one write per composite discovery, and the total number of those writes is n times the sum of the reciprocals of the primes up to √n — which Mertens tells us is ln ln n - ln 2 + M, a number between 1 and 4 for every input you will ever use. That is the whole argument, and the measured writes agree with it to a fraction of a percent.

Write it once, keep the array, and stop asking numbers whether they are prime.

Keep reading

• The Euclidean Algorithm — the other 2,000-year-old algorithm still in daily use, and the one behind RSA key setup.
• Fast Exponentiation — modular exponentiation in log n steps, which is what Miller-Rabin runs after the sieve has done the cheap rejections.
• Big O Notation — how to build the counting arguments used above from scratch.
• Arrays and Dynamic Arrays — why indexing a contiguous block is the fastest lookup there is, which is the sieve's whole advantage.