The sieve of Eratosthenes

A prime is a whole number above 1 whose only divisors are 1 and itself: 2, 3, 5, 7, 11, …. Every whole number above 1 is a product of primes in exactly one way, so primes are the building blocks of the whole numbers. This lesson tests single numbers, then finds all the primes up to n at once with a method more than 2,000 years old.

Testing one number

To test whether n is prime, try dividing it by 2, 3, 4, …. You can stop at √n, the square root of n. If n = a · b with a ≤ b, then a · a ≤ a · b = n, so a ≤ √n: every factor pair has one member at most √n.

from math import isqrt

def is_prime(n):
    if n < 2:
        return False
    for d in range(2, isqrt(n) + 1):
        if n % d == 0:
            return False
    return True

math.isqrt(n) is the square root rounded down, ⌊√n⌋, computed exactly for whole numbers of any size. int(n ** 0.5) goes through a float and can be off by one for large n. Trial division costs O(√n): for n = 10¹² that is up to a million divisions for one number.

Crossing out multiples

To find every prime up to n, Eratosthenes (about 240 BC) turned the question round. Instead of testing each number for divisors, take each prime and cross out its multiples. Whatever is never crossed out has no smaller prime divisor, so it is prime.

Two details make it fast.

• Start at p · p. A smaller multiple p · j, with j < p, has a prime factor smaller than p (a factor of j), so it was crossed out already.
• Stop at √n. Every crossed-out number up to n has a prime factor at most √n, by the factor-pair argument above. So after the primes up to √n, nothing new is crossed out.

from math import isqrt

def primes_upto(n):
    is_p = bytearray([1]) * (n + 1)   # one byte per number
    is_p[0] = is_p[1] = 0
    for p in range(2, isqrt(n) + 1):
        if is_p[p]:
            count = len(range(p * p, n + 1, p))
            is_p[p * p::p] = bytes(count)   # cross out
    return [i for i in range(n + 1) if is_p[i]]

The slice is_p[p * p::p] is every p-th entry from p² to the end, and bytes(count) is that many zero bytes. Assigning one to the other crosses out all the multiples in one step that runs in C, which is many times faster than a Python loop over them.

What the sieve costs

Prime p crosses out about n / p numbers, so the total is about n · (1/2 + 1/3 + 1/5 + 1/7 + …) over the primes up to √n. This sum of 1/p over primes grows like ln ln n, where ln is the natural logarithm, so the sieve is O(n log log n). That is nearly O(n): for n = 10⁹, ln ln n is only about 3.

Memory is the real limit: one byte per number is 100 MB for n = 10⁸. Storing only odd numbers halves it, and one bit per number (BITS) divides it by 8. How many primes are there? About n / ln n up to n (the prime number theorem): 78,498 up to a million, and 50,847,534 up to a billion.

Smallest factors

Instead of a yes-or-no mark, the sieve can store each number’s smallest prime factor: when prime p crosses out a number for the first time, p is its smallest prime factor. Then any number up to n can be broken into primes with one division per factor, as the second version of the figure shows for 45.

from math import isqrt

def smallest_factors(n):
    spf = list(range(n + 1))       # spf[x] = x until crossed out
    for p in range(2, isqrt(n) + 1):
        if spf[p] == p:            # p is prime
            for q in range(p * p, n + 1, p):
                if spf[q] == q:
                    spf[q] = p
    return spf

def factorise(x, spf):
    out = []
    while x > 1:
        out.append(spf[x])
        x //= spf[x]
    return out                     # factorise(45, spf) == [3, 3, 5]

Each division at least halves x, so factorising takes O(log x) steps after the sieve. A variant called the linear sieve crosses out each number exactly once, from its smallest prime factor, and runs in O(n).

Primes in a far window

What are the primes between 10¹² − 10⁶ and 10¹²? A sieve up to 10¹² needs a million million bytes. But only the window needs marks. Every composite number in it has a prime factor at most √10¹² = 10⁶. So:

1. sieve the small primes, up to √hi;
2. for each small prime p, cross out its multiples inside the window. The first multiple of p that is at least lo is (lo + p - 1) // p * p (lo rounded up to a multiple of p), but never start below p², or p itself would be crossed out.

This is a segmented sieve. Its memory is about √hi plus the width of the window, and its time is about the window’s sieve plus one division per small prime. MATH-P4 asks you to write it.

Now try MATH-P4, Primes in a window.