Meissel–Lehmer Algorithm

The “Meissel–Lehmer algorithm” is an algorithm that counts the primes in $1\sim n$ in sublinear time.

Notation

$\left[x\right]$ denotes the floor of $x$. \n$p_k$ denotes the $k$-th prime, with $p_1=2$. \n$\pi\left(x\right)$ denotes the number of primes in $1\sim x$. \n$\mu\left(x\right)$ denotes the Möbius function. \nFor a set $S$, $\# S$ denotes the size of $S$. \n$\delta\left(x\right)$ denotes the smallest prime factor of $x$. \n$P^+\left(x\right)$ denotes the largest prime factor of $x$.

The Meissel–Lehmer algorithm for computing π(x)

Define $\phi\left(x,a\right)$ as the number of positive integers less than $x$ whose every prime factor is greater than $p_a$, namely:

$$\phi\left(x,a\right)=\#\big\{n\le x\mid n\bmod p=0 \implies p>p_a\big\}\tag{1}$$

We further define $P_k\left(x,a\right)$ as the number of positive integers less than $x$ that have exactly $k$ prime factors counted with multiplicity and whose every prime factor is greater than $p_a$, namely:

$$P_k\left(x,a\right)=\#\big\{n\le x\mid n=q_1q_2\cdots q_k \implies \forall i,q_i>p_a\big\}\tag{2}$$

In particular, define $P_0\left(x,a\right)=1$, so that:

$$\phi\left(x,a\right)=P_0\left(x,a\right)+P_1\left(x,a\right)+\cdots+P_k\left(x,a\right)+\cdots$$

This infinite sum can in fact be written as a finite sum, because when $p_a^k>x$, we have $P_k\left(x,a\right)=0$.

Let $y$ be an integer satisfying $x^{1/3}\le y\le x^{1/2}$, and write $a=\pi\left(y\right)$.

For $k\ge 3$, we have $P_1\left(x,a\right)=\pi\left(x\right)-a$ and $P_k\left(x,a\right)=0$, from which we obtain:

$$\pi\left(x\right)=\phi\left(x,a\right)+a-1-P_2\left(x,a\right)\tag{3}$$

Thus, computing $\pi\left(x\right)$ can be reduced to computing $\phi\left(x,a\right)$ and $P_2\left(x,a\right)$.

Computing P₂(x,a)

By equation $\left(2\right)$, the number $P_2\left(x,a\right)$ equals the number of prime pairs satisfying $y<p\le q$ and $pq\le x$, where the pair is $\left(p,q\right)$.

First note that $p\in \left[y+1,\sqrt{x}\right]$. In addition, for each $p$, we have $q\in\left[p,x/p\right]$. Therefore:

$$P_2\left(x,a\right)=\sum_{y<p\le \sqrt{x}}{\left(\pi\left(\dfrac{x}{p}\right)-\pi\left(p\right)+1\right)}\tag{4}$$

When $p\in \left[y+1,\sqrt{x}\right]$, we have $\dfrac{x}{p}\in \left[1,\dfrac{x}{y}\right]$. Thus, we can sieve the interval $\left[1,\dfrac{x}{y}\right]$, and then, for every prime $p\in \left[y+1,\sqrt{x}\right]$, compute $\pi\left(\dfrac{x}{p}\right)-\pi\left(p\right)+1$. To reduce the space complexity of this algorithm, we can use blocking with block length $L$. If $L=y$, then $O\left(\dfrac{x}{y}\log{\log{x}}\right)$ can be computed in $O\left(y\right)$ time and $P_2\left(x,a\right)$ space.

Computing ϕ(x,a)

For $b\le a$, consider all positive integers not exceeding $x$ whose every prime factor is greater than $p_{b-1}$. These numbers fall into two classes:

1. Those divisible by $p_b$;
2. Those not divisible by $p_b$.

Class $1$ contains $\phi\left(\dfrac{x}{p_b},b-1\right)$ numbers, while class two contains $\phi\left(x,b\right)$ numbers.

Therefore, we conclude:

Theorem $5.1$: The function $\phi$ satisfies the following properties

$$\phi\left(u,0\right)=\left[u\right]\tag{5}$$ $$\phi\left(x,b\right)=\phi\left(x,b-1\right)-\phi\left(\dfrac{x}{p_b},b-1\right)\tag{6}$$

A simple method for computing $\phi\left(x,a\right)$ can be derived from this theorem: repeatedly apply equation $\left(7\right)$ until finally obtaining $\phi\left(u,0\right)$. This process can be viewed as creating a rooted binary tree beginning with the root node $\phi\left(x,a\right)$; Figure $1$ illustrates this process. In this way, we obtain the following formula:

$$\phi\left(x,a\right)=\sum_{\substack{1\le n\le x\\ P^+\left(n\right)\le y}}{\mu\left(n\right)\left[x/n\right]}$$ $$\begin{gathered} \begin{matrix}&&\phi\left(x,a\right)&&\\ &\swarrow&&\searrow&\\ &\phi\left(x,a-1\right)&&-\phi\left(\frac{x}{p_a},a-1\right)&\\ \swarrow&\downarrow&&\downarrow&\searrow\\ \phi\left(x,a-2\right)&\phi\left(\frac{x}{p_{a-1}},a-2\right)&&-\phi\left(\frac{x}{p_a},a-2\right)&\phi\left(\frac{x}{p_ap_{a-1}},a-2\right)\end{matrix}\\ \vdots\\ \end{gathered}$$

The figure above represents the binary tree for computing $\phi\left(x,a\right)$: the sum of the leaf weights is $\phi\left(x,a\right)$.

However, this requires computing too many things. Since $y\geq x^{1/3}$, merely counting numbers that are products of $3$ primes not exceeding $y$ would produce at least $\dfrac{x}{\log^3 x}$ terms under this method, which does not meet our complexity requirements.

To limit the “growth” of this binary tree, we need to change the original stopping condition. This is the original stopping condition.

Stopping condition $1$: If $b=0$, do not apply the rule to the node $\mu\left(n\right)\phi\left(\dfrac xn,b\right)$ using equation $\left(6\right)$ again.

We replace it with a stronger stopping condition:

Stopping condition $2$: If either of the following $2$ conditions holds, do not apply the rule to the node $\mu\left(n\right)\phi\left(\dfrac xn,b\right)$ using equation $\left(6\right)$ again:

1. $b=0$ and $n\le y$;
2. $n>y$.

According to Stopping condition $2$, we divide the leaves of the original binary tree into two types:

1. If the leaf node $\mu\left(n\right)\phi\left(\dfrac xn,b\right)$ satisfies $n\le y$, it is called an ordinary leaf;
2. If the leaf node $\mu\left(n\right)\phi\left(\dfrac xn,b\right)$ satisfies $n>y$ and $n=mp_b\left(m\le y\right)$, it is called a special leaf.

Therefore, we obtain:

Theorem $5.2$: We have:

$$\phi\left(x,a\right)=S_0+S\tag{7}$$

The contribution of ordinary leaves, $S_0$, is:

$$S_0=\sum_{n\le y}{\mu\left(n\right)\left[\dfrac xn\right]}\tag{8}$$

The contribution of special leaves, $S$, is:

$$S=\sum_{n/\delta\left(n\right)\le y\le n}{\mu\left(n\right)\phi\left(\dfrac{x}{n},\pi\left(\delta\left(n\right)\right)-1 \right)}\tag{9}$$

Computing $S_0$ can clearly be done in $O\left(y\log{\log x}\right)$ time; we now consider how to compute $S$.

Computing S

We have:

$$S=-\sum_{p\le y}{\ \sum_{\substack{\delta\left(m\right)>p\\ m\le y<mp}}{\mu\left(m\right)\phi\left(\dfrac{x}{mp},\pi\left(p\right)-1\right)}}\tag{10}$$

We rewrite this equation as:

$$S=S_1+S_2+S_3$$

where:

$$S_1=-\sum_{x^{1/3}<p\le y}{\ \sum_{\substack{\delta\left(m\right)>p\\ m\le y<mp}}{\mu\left(m\right)\phi\left(\dfrac{x}{mp},\pi\left(p\right)-1\right)}}$$ $$S_2=-\sum_{x^{1/4}<p\le x^{1/3}}{\ \sum_{\substack{\delta\left(m\right)>p\\ m\le y<mp}}{\mu\left(m\right)\phi\left(\dfrac{x}{mp},\pi\left(p\right)-1\right)}}$$ $$S_3=-\sum_{p\le x^{1/4}}{\ \sum_{\substack{\delta\left(m\right)>p\\ m\le y<mp}}{\mu\left(m\right)\phi\left(\dfrac{x}{mp},\pi\left(p\right)-1\right)}}$$

Notice that, in the sums for computing $S_1,S_2$, all occurring $m$ are prime numbers, as follows:

Otherwise, since $\delta\left(m\right)>p>x^{1/4}$, we would have $m>p^2>\sqrt{x}$, contradicting $m\le y$; hence the claim follows.

Moreover, when $mp>x^{1/2}\ge y$, we have $y\le mp$. Therefore:

$$S_1=\sum_{x^{1/3}<p\le y}{\ \sum_{p<q\le y}{\phi\left(\dfrac{x}{pq},\pi\left(p\right)-1\right)}}$$ $$S_2=\sum_{x^{1/4}<p\le x^{1/3}}{\ \sum_{p<q\le y}{\phi\left(\dfrac{x}{pq},\pi\left(p\right)-1\right)}}$$

Computing S₁

Because:

$$\dfrac{x}{pq}<x^{1/3}<p$$

Therefore:

$$\phi\left(\dfrac{x}{pq},\pi\left(p\right)-1\right)=1$$

Thus every term in the sum for $S_1$ is $1$. We therefore only need to count prime pairs $\left(p,q\right)$ satisfying $x^{1/3}<p<q\le y$.

Therefore:

$$S_1=\dfrac{\left(\pi\left(y\right)-\pi\left(x^{1/3}\right)\right)\left(\pi\left(y\right)-\pi\left(x^{1/3}\right)-1\right)}{2}$$

This equation permits an $O\left(1\right)$ computation of $S_1$.

Computing S₂

We have:

$$S_2=\sum_{x^{1/4}<p\le x^{1/3}}{\ \sum_{p<q\le y}{\phi\left(\dfrac{x}{pq},\pi\left(p\right)-1\right)}}$$

We split $S_2$ into the two parts $q>\dfrac x{p^2}$ and $q\le \dfrac x{p^2}$:

$$S_2=U+V$$

where:

$$U=\sum_{x^{1/4}<p\le x^{1/3}}{\ \sum_{\substack{p<q<y\\q>x/p^2}}{\phi\left(\dfrac{x}{pq},\pi\left(p\right)-1 \right)}}$$ $$V=\sum_{x^{1/4}<p\le x^{1/3}}{\ \sum_{\substack{p<q<y\\q\le x/p^2}}{\phi\left(\dfrac{x}{pq},\pi\left(p\right)-1 \right)}}$$

Computing U

Since $q>\dfrac x{p^2}$, we have $p^2>\dfrac xq\le \dfrac xy,p>\sqrt{\dfrac xy}$; therefore:

$$U=\sum_{\sqrt{x/y}<p\le x^{1/3}}{\ \sum_{\substack{p<q\le y\\q>x/p^2}}{\phi\left(\dfrac{x}{pq},\pi\left(p\right)-1 \right)}}$$

Therefore:

$$U=\sum_{\sqrt{x/y}<p\le x^{1/3}}{\#\left\{q\mid \dfrac x{p^2}<q\le y \right\}}$$

Therefore:

$$U=\sum_{\sqrt{x/y}<p\le x^{1/3}}{\left(\pi\left(y\right)-\pi\left(\dfrac{x}{p^2} \right) \right)}$$

Since $\dfrac x{p^2}<y$, precomputing the table $\pi\left(t\right)\left(t\le y\right)$ means that, in $O\left(y\right)$ time, we can compute $U$.

Computing V

For every term in the sum for $V$, we have $p\le \dfrac{x}{pq}<x^{1/2}<p^2$. Therefore:

$$\phi\left(\dfrac{x}{pq},\pi\left(p\right)-1 \right)=1+\pi\left(\dfrac{x}{pq} \right)-\left(\pi\left(p\right)-1\right)=2-\pi\left(p\right)+\pi\left(\dfrac{x}{pq} \right)$$

Thus $V$ can be written as:

$$V=V_1+V_2$$

where:

$$V_1=\sum_{x^{1/4}<p\le x^{1/3}}{\ \sum_{p<q\le \min\left(x/p^2,y\right)}{\left(2-\pi\left(p\right)\right)}}$$ $$V_2=\sum_{x^{1/4}<p\le x^{1/3}}{\ \sum_{p<q\le \min\left(x/p^2,y\right)}{\pi\left(\dfrac{x}{pq} \right)}}$$

After precomputing $\pi\left(t\right)\left(t\le y\right)$, an $O\left(x^{1/3}\right)$ computation gives $V_1$.

Consider how to accelerate the computation of $V_2$. We can split the contribution of $q$ into several intervals on each of which $\pi\left(\dfrac{x}{pq} \right)$ is constant. Then we only need to compute the length of each interval and the change in $\pi\left(\dfrac{x}{pq} \right)$ from one interval to the next.

More precisely, we first split $V_2$ into two parts to simplify the complex condition $q\le \min\left(\dfrac x{p^2},y\right)$:

$$V_2=\sum_{x^{1/4}<p\le \sqrt{x/y}}{\ \sum_{p<q\le y}{\pi\left(\dfrac{x}{pq} \right)}}+\sum_{\sqrt{x/y}<p\le x^{1/3}}{\ \sum_{p<q\le x/p^2}{\pi\left(\dfrac{x}{pq} \right)}}$$

Next, rewrite this expression as:

$$V_2=W_1+W_2+W_3+W_4+W_5$$

where:

$$W_1=\sum_{x^{1/4}<p\le x/y^2}{\ \sum_{p<q\le y}{\pi\left(\dfrac{x}{pq} \right)}}$$ $$W_2=\sum_{x/y^2<p\le \sqrt{x/y}}{\ \sum_{p<q\le \sqrt{x/p}}{\pi\left(\dfrac{x}{pq} \right)}}$$ $$W_3=\sum_{x/y^2<p\le \sqrt{x/y}}{\ \sum_{\sqrt{x/p}<q\le y}{\pi\left(\dfrac{x}{pq} \right)}}$$ $$W_4=\sum_{\sqrt{x/y}<p\le x^{1/3}}{\ \sum_{p<q\le \sqrt{x/p}}{\pi\left(\dfrac{x}{pq} \right)}}$$ $$W_5=\sum_{\sqrt{x/y}<p\le x^{1/3}}{\ \sum_{\sqrt{x/p}<q\le x/p^2}{\pi\left(\dfrac{x}{pq} \right)}}$$

Computing W₁ and W₂

First consider $y<\dfrac{x}{pq}<x^{1/2}$; the required value is $\pi\left(\dfrac{x}{pq} \right)$. It can be obtained by a block sieve on $[1,\sqrt x]$. In each block, for every $(p,q)$ satisfying the conditions, add $\pi\left(\dfrac x{pq}\right)$.

Computing W₃

For each $p$, divide $q$ into several intervals on each of which $\pi\left(\dfrac x{pq}\right)$ is constant; the contribution of each interval can be computed in $O(1)$. When a new $q$ is reached, use the table of $\pi(t)$ for $t\leq y$ to compute $\pi\left(\dfrac x{pq}\right)$. The table of primes up to $y$ gives the $\pi(t)<\pi(t+1)=\pi\left(\dfrac x{pq}\right)$ satisfying $t$. Continuing in this way gives the next $\pi\left(\dfrac x{pq}\right)$ at which $q$ changes.

Computing W₄

Compared with $W_3$, $W_4$ is smaller in $q$, so $\pi\left(\dfrac x{pq}\right)$ changes faster. At this point, using the method for computing $W_3$ to compute $W_4$ offers no advantage, so we directly enumerate the pairs $(p,q)$ to compute $W_4$.

Computing W₅

As for $W_3$, we compute $W_5$ in the same way.

Computing S₃

Using all primes less than $x^{1/4}$, we sieve the interval $\left[1,\dfrac xy\right]$ once. When the sieve reaches $p_k$, we have computed, for all squarefree $m$ with $\delta(m)>p_k$, the values $-\mu(m)\phi\left(\dfrac{x}{mp_k},k-1 \right)$. The sieve is performed in blocks; within each block, we maintain a binary tree to update in real time the intermediate results after sieving by all primes up to a given prime. Thus, we can find the number of values not yet sieved when the sieve reaches a given point in only $O(\log x)$ time.

Time and space complexity of the algorithm

The time and space complexity are determined by the following $3$ processes:

1. Computing $P_2\left(x,a\right)$;
2. Computing $W_1,W_2,W_3,W_4,W_5$;
3. Computing $S_3$.

Complexity of computing P₂(x,y)

We already know that this process takes $O\left(\dfrac{x}{y}\log{\log x}\right)$ time and $O\left(y\right)$ space.

Complexity of computing W₁,W₂,W₃,W₄,W₅

The block sieve used to compute $W_1,W_2$ has length $y$ and takes $O\left(\sqrt{x}\log{\log x}\right)$ time and $O\left(y\right)$ space.

The time complexity required to compute $W_1$ is:

$$\pi\left(\dfrac{x}{y^2} \right)\pi\left(y\right)=O\left(\dfrac{x}{y\log^2 x} \right)$$

The time complexity of computing $W_2$ is:

$$O\left(\sum_{x/y^2<p\le \sqrt{x/y}}{\pi\left(\sqrt{\dfrac xp}\right)} \right)=O\left(\dfrac{x^{3/4}}{y^{1/4}\log^2 x} \right)$$

Therefore, the time complexity of computing $W_3$ is:

$$O\left(\sum_{x/y^2<p\le \sqrt{x/y}}{\pi\left(\sqrt{\dfrac xp}\right)} \right)=O\left(\dfrac{x^{3/4}}{y^{1/4}\log^2 x} \right)$$

The time complexity of computing $W_4$ is:

$$O\left(\sum_{\sqrt{x/y}<p\le x^{1/3}}{\pi\left(\sqrt{\dfrac xp}\right)} \right)=O\left(\dfrac{x^{2/3}}{\log^2 x} \right)$$

The time complexity of computing $W_5$ is:

$$O\left(\sum_{\sqrt{x/y}<p\le x^{1/3}}{\pi\left(\sqrt{\dfrac xp}\right)} \right)=O\left(\dfrac{x^{2/3}}{\log^2 x} \right)$$

Complexity of computing S₃

For preprocessing: because we need to query $\phi(u,b)$ quickly, the ordinary sieve cannot compute it in $O(1)$ time. We must maintain a data structure so that each query takes $O(\log x)$ time, giving a total time complexity of $O\left(\dfrac{x}{y}\log x\log\log x\right)$.

For the summation: for each term in the sum for $S_3$, we query the data structure, for a total of $O\left(\log x\right)$ queries. We also need to count the number of terms in the sum, namely the number of leaves in the binary tree. All leaves have the form $\pm\phi\left(\dfrac{x}{mp_b},b-1\right)$, where $m\le y,b<\pi(x^{1/4})$. Therefore, the number of leaves is on the order of $O\left(y\pi\left(x^{1/4}\right)\right)$. Hence the total time complexity for computing $S_3$ is:

$$O\left(\dfrac{x}{y}\log x\log\log x+yx^{1/4}\right)$$

Overall complexity

The space complexity of this algorithm is $O\left(y\right)$, and its time complexity is:

$$O\left(\dfrac{x}{y}\log{\log x}+\dfrac{x}{y}\log x\log{\log x}+x^{1/4}y+\dfrac{x^{2/3}}{\log^2{x}} \right)$$

Taking $y=x^{1/3}\log^3{x}\log{\log x}$ gives the optimal time complexity $O\left(\dfrac{x^{2/3}}{\log^2 x}\right)$ and space complexity $O\left(x^{1/3}\log^3{x}\log{\log x}\right)$.

Some improvements

Here we give improvements that reduce the constant factors of the algorithm and increase its practical efficiency.

• In Stopping condition $2$, we can use $z$ in place of $y$, where $z$ satisfies $z>y$. It can be proved that the time complexity of computing $S_3$ can then be optimized to:

  $$O\left(\dfrac{x}{z}\log x\log{\log x}+\dfrac{yx^{1/4}}{\log x}+z^{3/2} \right)$$

  This also provides a good way to check the computation by changing the value of $z$.

• For clarity, when describing the algorithm we choose to split at $x^{1/4}$ when computing the total sum $S$, but in fact it suffices to have $p\le \dfrac{x}{pq}<p^2$ to perform the computation. This can be used without changing the asymptotic complexity.

• Preprocessing with the first few primes $2,3,5$ can save additional time.

References and further reading

This article is translated from: Computing $\pi(x)$: the Meissel, Lehmer, Lagarias, Miller, Odlyzko method