Prime Counting Function π(x), in depth: Legendre, Meissel-Lehmer and Lucy's O(x^3/4) algorithm, with tested code

How big is π(x)?

Before counting exactly, it helps to know roughly what answer to expect. The prime number theorem says π(x) behaves like x / ln x. The logarithmic integral li(x), the integral of 1 / ln t, is a much better estimate. Both are useful as sanity checks on an exact implementation, because a wrong answer is usually far off, not off by one.

li(x) is above π(x) in every row, and in every range anyone has computed directly. Littlewood proved that the sign does eventually flip infinitely often. The first crossing is far beyond any value you will meet, so a test asserting π(x) < li(x) is valid in practice. It is still not a theorem, so do not state it as one.

Legendre&#x27;s identity

Every composite number n ≤ x has a prime factor no larger than √x. So if we remove all multiples of the primes p₁, …, p_a up to √x, what remains is 1 and the primes between √x and x. Define φ(x, a) as the count of integers in [1, x] not divisible by any of the first a primes. Then:

π(x) = φ(x, a) + a − 1, where a = π(√x).

The + a adds back the a small primes that we removed. The − 1 removes the number 1, which survives the sieve but is not prime. φ has a simple recurrence. The survivors after a primes are the survivors after a − 1 primes, minus those that are multiples of p_a. Those multiples, divided by p_a, are exactly the survivors up to x / p_a:

φ(x, a) = φ(x, a − 1) − φ(⌊x / p_a⌋, a − 1), with φ(x, 0) = ⌊x⌋.

Worked example: π(100) by hand

Take x = 100. Then √100 = 10, the primes up to 10 are 2, 3, 5 and 7, so a = 4. Expanding the recurrence fully is inclusion and exclusion over every product of those primes:

So φ(100, 4) = 100 − 117 + 45 − 6 = 22. These 22 survivors are 1 and the 21 primes from 11 to 97. Then π(100) = 22 + 4 − 1 = 25, which is correct. Notice how many terms are zero. Whenever the product of primes exceeds x the term vanishes, and the recursion can stop. Even so, the plain recursion still explores far too many branches for large x. That is what the next two methods fix.

from functools import lru_cache
from math import isqrt

def primes_upto(n):
    is_p = bytearray([1]) * (n + 1)
    is_p[0:2] = b"\x00\x00"
    for i in range(2, isqrt(n) + 1):
        if is_p[i]:
            is_p[i * i::i] = bytearray(len(range(i * i, n + 1, i)))
    return [i for i in range(n + 1) if is_p[i]]

def legendre(x):
    ps = primes_upto(isqrt(x))

    @lru_cache(maxsize=None)
    def phi(y, a):
        if a == 0:
            return y
        return phi(y, a - 1) - phi(y // ps[a - 1], a - 1)

    a = len(ps)
    return phi(x, a) + a - 1          # legendre(100) == 25

Meissel-Lehmer: sieve less, correct with P2

Meissel's improvement is to sieve with fewer primes. Choose a = π(x^(1/3)) instead of π(√x). Then φ(x, a) counts 1, the primes above x^(1/3), and also some composites that have no prime factor up to x^(1/3). Such a composite cannot have three prime factors, because their product would exceed x. So each one is exactly p × q with x^(1/3) < p ≤ q. Call the count of those products P2(x, a). Then:

π(x) = φ(x, a) + a − 1 − P2(x, a)

P2 is easy to count. For each prime p between x^(1/3) and √x, the valid partners q are the primes in [p, x / p], and there are π(x / p) − π(p) + 1 of them. Summing over the b − a primes p_a+1…p_b with b = π(√x) gives the formula used in the code. The look-ups need π(y) for y up to x^(2/3), which a single sieve provides. Sieving to x^(2/3) instead of x is the whole saving. At x = 10^12 that is a sieve of 10^8, which is small.

from bisect import bisect_right

def icbrt(x):                          # exact integer cube root; never trust x ** (1/3)
    c = round(x ** (1 / 3))
    while c ** 3 > x:
        c -= 1
    while (c + 1) ** 3 <= x:
        c += 1
    return c

def meissel_lehmer(x):
    if x < 2:
        return 0
    c, r = icbrt(x), isqrt(x)
    limit = max(r, x // (c + 1))      # the largest x // p we will look up
    primes = primes_upto(limit)
    pi = lambda y: bisect_right(primes, y)
    a, b = pi(c), pi(r)

    @lru_cache(maxsize=None)
    def phi(y, k):
        if k == 0 or y == 0:
            return y
        return phi(y, k - 1) - phi(y // primes[k - 1], k - 1)

    # i is the 0-based index of p, so pi(p) - 1 == i
    p2 = sum(pi(x // primes[i]) - i for i in range(a, b))
    return phi(x, a) + a - 1 - p2

The P2 term in the code uses π(x / p) − i, where i is p's zero-based index. That is π(x / p) − (π(p) − 1), the same count as above. Lehmer's version also subtracts a P3 term for products of three primes, which lets a be even smaller. The idea is the same, with more bookkeeping.

Lucy&#x27;s algorithm

The method most competitive programmers use was posted by the user Lucy_Hedgehog on the Project Euler forum, so it is called Lucy's algorithm. It rests on one observation. During a sieve, we only ever need counts at values of the form ⌊x / k⌋, and there are fewer than 2√x distinct such values. For x = 100 they are 100, 50, 33, 25, 20, 16, 14, 12, 11, 10 and then 9 down to 1.

Let S[v] be the count of numbers in [2, v] that survive sieving by the primes processed so far. Start with S[v] = v − 1. When we process prime p, the numbers removed from [2, v] are the composites whose smallest prime factor is p. Each is p × m, where m ≤ v / p survives sieving by the primes below p and m ≥ p. Their count is S[⌊v/p⌋] − S[p − 1]. That gives the update:

S[v] ← S[v] − (S[⌊v/p⌋] − S[p − 1]), for every v ≥ p².

Values below p² are not affected, because the smallest such composite is p² itself. That is why the inner loop can stop early. Since ⌊⌊x/a⌋/b⌋ = ⌊x/(ab)⌋, every look-up lands on a value we already store. After processing every prime up to √x, S[x] = π(x). The total work is O(x^(3/4)) and memory is O(√x).

def lucy(n):
    r = isqrt(n)
    V = [n // i for i in range(1, r + 1)]
    V += list(range(V[-1] - 1, 0, -1))     # descending: n, n//2, ..., 1
    S = {v: v - 1 for v in V}
    for p in range(2, r + 1):
        if S[p] > S[p - 1]:                 # p survived, so p is prime
            sp = S[p - 1]                   # primes below p
            p2 = p * p
            for v in V:
                if v < p2:
                    break                   # V is descending: nothing smaller changes
                S[v] -= S[v // p] - sp
    return S[n]

In pure Python this computes π(10^10) = 455,052,511 in about three seconds on an ordinary desktop, and π(10^11) = 4,118,054,813 in under a minute. Replace the dictionary with two arrays, one indexed by v for small v and one by n // v for large v, then port it to C++ with 64-bit integers, and 10^13 becomes a matter of seconds rather than minutes. The same skeleton also computes sums of primes, if each S[v] starts as the sum 2 + … + v and each p is a weight.

Beyond x^(3/4)

Lucy's algorithm is the right tool up to about 10^13. Beyond that, the combinatorial line continues. Lagarias, Miller and Odlyzko (LMO, 1985) restructured Meissel-Lehmer to run in O(x^(2/3) / log x) time with a segmented sieve, so memory stays small. Deléglise and Rivat (1996) reduced the time to O(x^(2/3) / log² x). Xavier Gourdon's 2001 variant improved the constants and parallelised it. There is also an analytic method, based on the Riemann zeta function, with better asymptotic complexity. It is much harder to implement correctly and has been used mainly to compute and cross-check a few record values.

For production, do not write any of these yourself. Kim Walisch's open-source primecount library implements Deléglise-Rivat and Gourdon with multi-threading, and its documentation lists its own verified results. Use it as both your engine and your oracle. Your own implementation is valuable for learning and for variants the library does not support, such as sums of primes or primes in a residue class.

Testing a prime counter

Prime counting fails quietly. An off-by-one in a floor or a bad cube root returns a plausible number. Test it like this:

• Brute force for small x. Compare against a sieve for every x from 0 to a few thousand. That catches boundary bugs at perfect squares and cubes, which random tests almost never hit.
• Two methods against each other. Run Lucy and Meissel-Lehmer on hundreds of random x up to 10^7. They share almost no code, so agreement is strong evidence.
• Known values. Assert π(10^k) for k = 1 to 11: 4, 25, 168, 1,229, 9,592, 78,498, 664,579, 5,761,455, 50,847,534, 455,052,511 and 4,118,054,813.
• Approximation bounds. Assert x / ln x < π(x) < li(x) for every tested x ≥ 17. This cheap check catches gross errors on inputs too large to verify exactly.

Failure modes

• Floating-point roots. int(x ** 0.5) and round(x ** (1/3)) are wrong for large perfect powers. Use math.isqrt and an integer-corrected cube root, as above.
• Overflow in C or C++. p², and products such as p × q near x, overflow 64-bit integers once x passes about 9.2 × 10^18. Use __int128 for intermediates, or cap the input range and assert it.
• Unbounded φ recursion. Plain Legendre explores about 2^a branches. Memoise it, cut branches where y = 0 or y < p_k, and precompute φ(y, k) for small k using the periodicity of the first few primes.
• Dictionary memory. The dictionary in the Lucy code above stores about 2√x entries, which is fine at 10^12 (2 million), but switch to arrays before 10^14.
• Counting 1 as prime, or forgetting 2. Every formula has a ± 1 term. Test x = 0, 1, 2 and 3 explicitly.

Trade-offs

What to do next

1. Trace φ(100, 4) by hand, as in the table above, until the + a − 1 term feels obvious.
2. Type in lucy() and check it against a sieve for every x below 5,000.
3. Add the known-values test for 10^1 to 10^11 and the x / ln x < π(x) < li(x) check.
4. Convert Lucy's dictionary into two flat arrays and time π(10^11) before and after.
5. Modify it to return the sum of primes up to x, and verify the sum below 10 is 17 and below 2 million is 142,913,828,922.
6. Install primecount, compare its result with yours at 10^12, then read the Deléglise-Rivat paper with a working baseline in hand.
7. Continue with advanced number theory and Euler's totient function, whose prefix sums use the same ⌊x / k⌋ trick.