Proofs: Meissel-Lehmer algorithm

The Meissel-Lehmer algorithm returns pi(x) full

This proof shows: correctness.

Write $p_1 < p_2 < \cdots$ for the primes, $P_b = p_1 \cdots p_b$, and $\phi(t, b)$ for the number of integers $1 \le n \le t$ with no prime factor among $p_1, \ldots, p_b$. For $m \ge 2$, $\mathrm{lpf}(m)$ is its least prime factor. Two facts are used throughout. First, for real $t \ge 0$ and integers $u, v \ge 1$, $\lfloor \lfloor t/u \rfloor / v \rfloor = \lfloor t/(uv) \rfloor$. Second, the inclusion-exclusion form of $\phi$: $\phi(t, b) = \sum_{d \mid P_b} \mu(d) \lfloor t/d \rfloor$, because $\sum_{d \mid \gcd(n, P_b)} \mu(d)$ is 1 if $\gcd(n, P_b) = 1$ and 0 otherwise, and $\lfloor t/d \rfloor$ counts the $n \le t$ with $d \mid n$.

If $x < 2$, there is no prime up to $x$ and line L1 returns 0. Assume $x \ge 2$ from now on. Line L2 sets $y = \lfloor x^{1/3} \rfloor \ge 1$, so $y^3 \le x < (y+1)^3$, and $z = \lfloor x/y \rfloor$. Line L3 gives the primes up to $z$, the table $\pi(t)$ for $t \le z$, $\mathrm{lpf}$ on $[2, z]$ and $\mu$ on $[1, y]$; since $y \le z$, $a = \pi(y)$ is well defined, and $p_a \le y < p_{a+1}$. Since $y \le x^{1/3} \le x^{1/2}$, also $\lfloor \sqrt x \rfloor \le x / y$, hence $\lfloor \sqrt x \rfloor \le z$.

Step 1: $\pi(x) = \phi(x, a) + a - 1 - P_2$ with $P_2$ as in line L10. The numbers counted by $\phi(x, a)$ are the $1 \le n \le x$ all of whose prime factors are at least $p_{a+1} \ge y + 1 > x^{1/3}$. If such an $n$ had three prime factors counted with multiplicity, then $n > (x^{1/3})^3 = x$, which is impossible. So $n$ is 1, a prime $p$ with $y < p \le x$, or a product $pq$ of primes with $y < p \le q$ and $pq \le x$. The primes number $\pi(x) - a$. In a product $pq$ with $p \le q$, $p^2 \le pq \le x$, so $y < p \le \lfloor \sqrt x \rfloor$, and for such a $p$ the possible $q$ are the primes with $p \le q \le \lfloor x/p \rfloor$; as $p \le \lfloor x/p \rfloor$, they number $\pi(\lfloor x/p \rfloor) - \pi(p - 1) = \pi(\lfloor x/p \rfloor) - \pi(p) + 1$. Each product is counted once, for its smaller factor $p$. So $\phi(x, a) = 1 + (\pi(x) - a) + P_2$, which rearranges to the claim. All arguments of $\pi$ in line L10 are at most $z$: $\lfloor x/p \rfloor \le \lfloor x/y \rfloor = z$ for $p > y$, and $p \le \lfloor \sqrt x \rfloor \le z$. So line L10 computes $P_2$ from the table.

Step 2: $\phi(x, a) = S_0 + S$. By inclusion-exclusion, $\phi(x, a) = \sum_{d \mid P_a} \mu(d) \lfloor x/d \rfloor$. Split the sum by $d \le y$ and $d > y$.

The terms with $d \le y$. A squarefree $m \le y$ has all its prime factors at most $y$, hence at most $p_a$, so $m \mid P_a$; for a non-squarefree $m$, $\mu(m) = 0$. So these terms are $\sum_{m \le y} \mu(m) \lfloor x/m \rfloor = S_0$, which is the value set by line L4.

The terms with $d > y$. Let $\mathcal{D}$ be the set of $d \mid P_a$ with $d > y$, and $\mathcal{T}$ the set of triples $(m, b, r)$ with $0 \le b < a$, $m$ squarefree, $m \le y < m p_{b+1}$, every prime factor of $m$ larger than $p_{b+1}$, and $r \mid P_b$. We show that $(m, b, r) \mapsto m p_{b+1} r$ is a bijection from $\mathcal{T}$ onto $\mathcal{D}$.

It maps into $\mathcal{D}$: $m$, $p_{b+1}$ and $r$ are squarefree, with prime factors in the disjoint ranges $(p_{b+1}, p_a]$ (as $m \le y$), $\{p_{b+1}\}$ and $[p_1, p_b]$; so their product is squarefree, divides $P_a$, and is at least $m p_{b+1} > y$.

It is onto and one to one. Let $d \in \mathcal{D}$ have prime factors $q_1 > q_2 > \cdots > q_k$. The prefix products $1, q_1, q_1 q_2, \ldots, d$ increase, the first is $1 \le y$ and the last is $d > y$, so there is a unique $j \ge 1$ with $q_1 \cdots q_{j-1} \le y < q_1 \cdots q_j$. Put $m = q_1 \cdots q_{j-1}$, $p_{b+1} = q_j$ and $r = q_{j+1} \cdots q_k$. Then $(m, b, r) \in \mathcal{T}$: $b < a$ since $q_j \le p_a$, the prime factors of $m$ exceed $q_j$, and $r \mid P_b$ since its factors are below $q_j$. Conversely, for any $(m, b, r) \in \mathcal{T}$ with product $d$, the decreasing list of the prime factors of $d$ is that of $m$, then $p_{b+1}$, then that of $r$, so the construction recovers exactly $m$, $b$ and $r$: the prefix products through the factors of $m$ end at $m \le y$, and the next one is $m p_{b+1} > y$.

Therefore, using $\mu(m p_{b+1} r) = \mu(m) \cdot (-1) \cdot \mu(r)$ for these coprime squarefree factors and the floor identity,

where the sum over $m$ runs over the squarefree $m$ with $y / p_{b+1} < m \le y$ whose prime factors exceed $p_{b+1}$. For an integer $m$, $y / p_{b+1} < m$ holds if and only if $\lfloor y / p_{b+1} \rfloor < m$. The condition "every prime factor exceeds $p_{b+1}$" is $\mathrm{lpf}(m) > p_{b+1}$ for $m \ge 2$; and $m = 1$ never occurs, because $p_{b+1} \le p_a \le y$ gives $\lfloor y / p_{b+1} \rfloor \ge 1$. Squarefree is $\mu(m) \ne 0$. So the pairs $(b, m)$ visited by line L7 are exactly the terms of $S$.

Step 3: the Fenwick tree. Let $A_b$ be the set of $n \in [1, z]$ with no prime factor among $p_1, \ldots, p_b$. We claim that at the start of round $b$ of line L6, the array behind $F$ holds 1 at the positions in $A_b$ and 0 elsewhere. For $b = 0$, line L5 sets all of $[1, z]$, which is $A_0$. If the claim holds at the start of round $b$, line L9 sets to 0 every multiple of $p_{b+1}$ in $[1, z]$ and changes nothing else, so afterwards the 1 positions are $A_b$ minus the multiples of $p_{b+1}$, which is $A_{b+1}$. The test $F[k] = 1$ only skips multiples that are already 0. During round $b$, lines L7 and L8 run before line L9 and do not change $F$, so each query of round $b$ sees $A_b$, and $F.\mathrm{prefix}(n) = |A_b \cap [1, n]| = \phi(n, b)$ for $0 \le n \le z$. The query index is $n = \lfloor x / (m p_{b+1}) \rfloor$ with $m p_{b+1} \ge y + 1$, so $n \le \lfloor x/y \rfloor = z$, and the query is in range. Hence line L8 adds exactly the term $-\mu(m) \phi(\lfloor x / (m p_{b+1}) \rfloor, b)$ of $S$, once for each term.

Conclusion. After the loop of line L6, the variable $\phi$ equals $S_0 + S = \phi(x, a)$ by Step 2. Line L10 computes $P_2$, and line L11 returns $\phi(x, a) + a - 1 - P_2 = \pi(x)$ by Step 1. $\blacksquare$