Linear Sieve, in depth: the exactly-once proof, lowest-prime-factor tables, multiplicative functions and measured speed against Eratosthenes

The one rule: mark by the smallest prime factor

Every composite x has a unique representation x = p * i in which p is the smallest prime factor of x. Then i = x / p, and every prime factor of i is at least p, so p <= lp[i]. The linear sieve generates composites from exactly these pairs. For each i from 2 to n, it walks the list of primes found so far in increasing order and marks i * p with lp = p. It stops as soon as p exceeds lp[i], or as soon as i * p exceeds n.

Each composite is marked at least once. Take a composite x with smallest prime factor p, and let i = x / p. Because i is smaller than x, the outer loop reaches i first. By then lp[i] is already known, because i is either prime (found when the loop reached it) or composite (marked earlier, which is the same argument applied to a smaller number). The prime p is no larger than i, so p is already in the prime list. Since p <= lp[i], the inner loop has not stopped before p, so it marks x.

Each composite is marked at most once. The pair (i, p) that marks x always satisfies p <= lp[i], which means p is the smallest prime factor of i * p. A number has only one smallest prime factor, and so it has only one such pair. No other (i, p) can produce x. The number of marks is therefore exactly the number of composites up to n. The inner loop also does one comparison that fails before it stops. That adds one step per i, so the total work is O(n), against about n ln ln n marks for Eratosthenes.

The code and a hand trace to 30

The code is short. The two break conditions are the whole algorithm, so read them first:

def linear_sieve(n):
    """lp[x] = smallest prime factor of x (lp[x] == x for primes); primes ascending."""
    lp = [0] * (n + 1)
    primes = []
    for i in range(2, n + 1):
        if lp[i] == 0:              # nobody marked i, so i is prime
            lp[i] = i
            primes.append(i)
        for p in primes:
            if p > lp[i] or i * p > n:
                break
            lp[i * p] = p           # p is the smallest prime factor of i*p
    return lp, primes

For n = 30 the trace is short enough to check by hand. The figure lists every i that marks something. Look at i = 6. Its lowest prime factor is 2, so it marks 12 and stops. It does not mark 18, because 18 = 6 * 3 and the smallest prime factor of 18 is 2, not 3. That number is marked later by i = 9 with p = 2. Likewise 30 is marked only by i = 15 with p = 2. Eratosthenes would also have crossed it at p = 3 and p = 5.

The output table for 2..30 is 2 3 2 5 2 7 2 3 2 11 2 13 2 3 2 17 2 19 2 3 2 23 2 5 2 3 2 29 2. Primes are the entries where lp[x] equals x. There are 10 primes and 19 composites up to 30, and the trace marks exactly those 19 composites.

Multiplicative functions in the same pass

A function f is multiplicative if f(ab) = f(a) f(b) whenever a and b are coprime. Examples are Euler's totient φ, the Möbius function μ, the divisor count d and the divisor sum σ. The linear sieve can fill all of them because each composite is created once, as m = i * p, and there are only two cases to handle.

• p < lp[i]: p does not divide i, so p and i are coprime and f(m) = f(i) f(p). This is the easy case.
• p == lp[i]: p already divides i, so f(m) is not f(i) f(p). You need the exponent of p in i, because the effect of raising that exponent depends on the function. Keep an array e[x] holding the exponent of lp[x] in x.

With e[] available, each function has a one-line update. For φ, raising the exponent of a prime that is already present multiplies the value by p. For μ, any repeated prime gives zero. For d, the factor (e + 1) becomes (e + 2):

def sieve_multiplicative(n):
    lp = [0] * (n + 1); primes = []
    phi = [0] * (n + 1); mu = [0] * (n + 1); d = [0] * (n + 1)
    e = [0] * (n + 1)                       # exponent of lp[x] in x
    phi[1] = mu[1] = d[1] = 1
    for i in range(2, n + 1):
        if lp[i] == 0:
            lp[i] = i; primes.append(i)
            phi[i] = i - 1; mu[i] = -1; d[i] = 2; e[i] = 1
        for p in primes:
            m = i * p
            if p > lp[i] or m > n:
                break
            lp[m] = p
            if p == lp[i]:                  # exponent of p grows by one
                e[m] = e[i] + 1
                phi[m] = phi[i] * p
                mu[m] = 0
                d[m] = d[i] // (e[i] + 1) * (e[i] + 2)
            else:                           # p is a new, smaller prime
                e[m] = 1
                phi[m] = phi[i] * (p - 1)
                mu[m] = -mu[i]
                d[m] = d[i] * 2
    return lp, primes, phi, mu, d

This version was checked against brute-force φ (count of k with gcd(k, x) = 1), brute-force d (count of divisors) and μ derived from factorisation, for every x up to 3,000, and all three agreed. The integer division in the d update is exact because d[i] contains the factor (e[i] + 1). For functions where the prime-power case has no simple rule, such as σ, keep a second array holding lp[x] ** e[x]. Then f(m) = f(m / pk) * f(pk) where pk is that power. The Möbius function article shows what μ is then used for, and the totient article covers φ beyond the sieve.

Factorising with the lp table

Once lp is built, factorising any x up to n is a loop that divides by lp[x] repeatedly. Every step at least halves x, so a query costs O(log x) and needs no trial division:

def factor(x, lp):
    out = []
    while x > 1:
        p, k = lp[x], 0
        while x % p == 0:
            x //= p; k += 1
        out.append((p, k))
    return out

# factor(360, lp)  -> [(2, 3), (3, 2), (5, 1)]
# factor(2310, lp) -> [(2, 1), (3, 1), (5, 1), (7, 1), (11, 1)]

This is the usual reason to build the table: millions of values below a known bound, factorised in a few array reads each after one O(n) build. For a single large number the table is useless, because you cannot build it up to 10^18. Use Pollard's rho for that.

Measured: marks versus time

Counting marks agrees with the proof: the linear sieve marks each composite once, and Eratosthenes, starting at p squared, does slightly more than twice as many:

The second measures time, and it is where the extra marks stop mattering. In Java 23, taking the best of three runs at n = 10^7, the linear sieve took about 68 ms. Eratosthenes on a boolean[] took about 55 ms. An odd-only Eratosthenes on a long[] bitset took about 28 ms. All three found 664,579 primes. At n = 10^8 the runs were noisy, with GC and paging swinging repeated runs by 2x or more. The first run showed about 0.9 s for the linear sieve, 0.8 s for the boolean array and 0.35 s for the bitset. The bitset won every repetition. The boolean array and the linear sieve swapped places once memory pressure set in, so treat their gap at that size as noise.

The reasons are mechanical. Eratosthenes writes with a fixed stride into a structure that is 4 to 64 times smaller, and its inner loop has no data-dependent branch. In the linear sieve, the inner loop length depends on lp[i], and every mark is a 4-byte store at an address that jumps around. Count the work in cache lines touched, not in marks, and measure on your own hardware.

Operational guidance: memory and integer width

Memory is the real limit. A 32-bit lp array costs 4(n + 1) bytes, which is about 400 MB at n = 10^8. Python lists cost far more, roughly 8 bytes per slot for the pointer plus the int objects. There are three standard ways to shrink it:

• Store lp only for composites. The smallest prime factor of a composite up to n is at most √n. Store 0 for primes and use a 16-bit array, which is valid while √n < 65,536, so up to about 4.29 * 10^9; factor() must then read 0 as prime. At n = 10^8 that halves the table to about 200 MB.
• Skip even numbers. Every even number has lp = 2, so store the table only for odd x. This halves memory again, at the cost of index arithmetic.
• Do you need the table at all? If the job is to list or count primes, a segmented odd-only bitset Eratosthenes runs in about √n memory plus one L1-sized segment. See the Eratosthenes article for that layout. The linear sieve does not segment naturally, because marking from i needs lp[i] from anywhere earlier in the range.

Integer width. i * p can overflow 32 bits before the comparison with n catches it when n is near 2^31. Compute it in 64 bits, as the Java benchmark does with (long) i * p, or compare p > n / i instead. Prime array size. In C or Java, size the prime array from an upper bound such as 1.26 n / ln n, valid for every n ≥ 2, or simply n / 2 + 10.

Failure modes

• Breaking at the wrong time. If the loop breaks before marking when p == lp[i], instead of after, prime squares such as 4, 9 and 25 are never marked and show up as primes. Test the first 30 values against a known list.
• Treating p == lp[i] as coprime. If the multiplicative update uses f(i) f(p) in both branches, you get φ(4) = 1 instead of 2 and μ(4) = 1 instead of 0. The brute-force cross-check catches this immediately.
• Choosing it for speed alone. If you only use the prime list, you pay for a table you throw away.

Trade-offs

For the last row, the prime-counting article covers the Meissel-Lehmer family, which counts primes without listing them.

What to do next

1. Type the 12-line sieve from memory and check lp[2..30] against the table above.
2. Add the multiplicative version and assert φ, μ and d against brute force up to a few thousand before you trust it anywhere.
3. Benchmark it against a bitset Eratosthenes in your production language at your real n. Keep the linear sieve only if you use lp or the functions.
4. Size the memory: 4 bytes per entry, or 2 bytes with the composite-only 16-bit layout. Then decide whether to compute the table at build time or at startup.
5. Read the Möbius article to see the sieve's μ output turned into coprime counting.