Meissel–Lehmer Algorithm

The "Meissel–Lehmer algorithm" is an algorithm that can compute the number of primes in the range \(1\sim n\) in sublinear time complexity.

Notation Conventions

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

Meissel–Lehmer Algorithm for Computing π(x)

Define \(\phi\left(x,a\right)\) as the number of positive integers less than \(x\) whose all prime factors are greater than \(p_a\), i.e.:

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

Further define \(P_k\left(x,a\right)\) as the number of positive integers less than \(x\) that have exactly \(k\) prime factors (counting multiplicity) and all prime factors greater than \(p_a\), i.e.:

\[ 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} \]

Specifically, we define: \(P_0\left(x,a\right)=1\), so we have:

\[ \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 actually be expressed 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 denote \(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 can deduce:

\[ \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)\) reduces to computing \(\phi\left(x,a\right)\) and \(P_2\left(x,a\right)\).

Computing P₂(x,a)

From equation \(\left(2\right)\), we can deduce that \(P_2\left(x,a\right)\) equals the number of prime pairs \(\left(p,q\right)\) satisfying \(y<p\le q\) and \(pq\le x\).

First, note that \(p\in \left[y+1,\sqrt{x}\right]\). Moreover, 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]\). Therefore, we can sieve the interval \(\left[1,\dfrac{x}{y}\right]\), then for all primes \(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 the above algorithm, we can consider blocking with block length \(L\). If block length \(L=y\), then we can compute \(P_2\left(x,a\right)\) in \(O\left(\dfrac{x}{y}\log{\log{x}}\right)\) time complexity and \(O\left(y\right)\) space complexity.

Computing ϕ(x,a)

For \(b\le a\), consider all positive integers not exceeding \(x\) whose all prime factors are greater than \(p_{b-1}\). These numbers can be divided into two categories:

Those divisible by \(p_b\);
Those not divisible by \(p_b\).

Numbers in category 1 amount to \(\phi\left(\dfrac{x}{p_b},b-1\right)\), and numbers in category 2 amount to \(\phi\left(x,b\right)\).

Thus 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 to compute \(\phi\left(x,a\right)\) can be derived from this theorem: we repeatedly apply equation \(\left(7\right)\) until we obtain \(\phi\left(u,0\right)\). This process can be seen as constructing a rooted binary tree starting from the root node \(\phi\left(x,a\right)\). Figure \(1\) illustrates this process. Through this method, 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 above diagram represents the binary tree for computing \(\phi\left(x,a\right)\): the sum of the leaf node weights is \(\phi\left(x,a\right)\).

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

To limit the "growth" of this binary tree, we change the original termination condition. Here is the original termination condition.

Termination Condition \(1\): If \(b=0\), then do not apply equation \(\left(6\right)\) to node \(\mu\left(n\right)\phi\left(\dfrac xn,b\right)\) anymore.

We change it to a stronger termination condition:

Termination Condition \(2\): Do not apply equation \(\left(6\right)\) to node \(\mu\left(n\right)\phi\left(\dfrac xn,b\right)\) if one of the following two conditions is satisfied:

\(b=0\) and \(n\le y\);

\(n>y\).

According to Termination Condition \(2\), we divide the leaves on the original binary tree into two types:

If leaf node \(\mu\left(n\right)\phi\left(\dfrac xn,b\right)\) satisfies \(n\le y\), it is called a normal leaf;
If 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.

Thus we obtain:

Theorem \(5.2\): We have:

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

where \(S_0\) denotes the contribution of normal leaves:

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

and \(S\) denotes the contribution of special leaves:

\[ 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 complexity. Now we 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)}} \]

Note that the \(m\) involved in the sums for \(S_1,S_2\) are all primes, as proven below:

If not, since \(\delta\left(m\right)>p>x^{1/4}\), we have \(m>p^2>\sqrt{x}\), which contradicts \(m\le y\). Hence the proposition holds.

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

\[ 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₁

Since:

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

we have:

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

So each term in the sum for \(S_1\) is \(1\). Thus we actually need to count the number of 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} \]

With this equation, we can compute \(S_1\) in \(O\left(1\right)\) time.

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 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

From \(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)}} \]

Thus:

\[ 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\), we can precompute all \(\pi\left(t\right)\left(t\le y\right)\), allowing us to compute \(U\) in \(O\left(y\right)\) time complexity.

Computing V

For each 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) \]

So \(V\) can be expressed 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)}} \]

By precomputing \(\pi\left(t\right)\left(t\le y\right)\), we can compute \(V_1\) in \(O\left(x^{1/3}\right)\) time complexity.

Consider how to accelerate the computation of \(V_2\). We can split the contribution of \(q\) into several intervals where \(\pi\left(\dfrac{x}{pq} \right)\) is constant, so 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, simplifying 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)}} \]

Then we rewrite this 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₂

Computing these two values requires computing \(\pi\left(\dfrac{x}{pq} \right)\) for \(y<\dfrac{x}{pq}<x^{1/2}\). This can be done by block sieving over the interval \([1,\sqrt x]\). In each block, we accumulate \(\pi\left(\dfrac x{pq}\right)\) for all eligible \((p,q)\).

Computing W₃

For each \(p\), we divide \(q\) into several intervals such that \(\pi\left(\dfrac x{pq}\right)\) is constant within each interval, and we can compute the contribution of each interval in \(O(1)\). When we obtain a new \(q\), we compute \(\pi\left(\dfrac x{pq}\right)\) using the table of \(\pi(t)\) (for \(t\leq y\)). The prime table up to \(y\) gives the \(t\) such that \(\pi(t)<\pi(t+1)=\pi\left(\dfrac x{pq}\right)\). And so on for the next \(q\) value where \(\pi\left(\dfrac x{pq}\right)\) changes.

Computing W₄

Compared to \(W_3\), in \(W_4\) \(q\) is smaller, so \(\pi\left(\dfrac x{pq}\right)\) changes more rapidly. Using the method for \(W_3\) to compute \(W_4\) offers no advantage. Therefore, we directly brute-force enumerate pairs \((p,q)\) to compute \(W_4\).

Computing W₅

We compute \(W_5\) similarly to \(W_3\).

Computing S₃

We use all primes less than \(x^{1/4}\) to sieve the interval \(\left[1,\dfrac xy\right]\) once. When our sieve reaches \(p_k\), we have computed all \(-\mu(m)\phi\left(\dfrac{x}{mp_k},k-1 \right)\) for \(m\) that are square-free and satisfy \(\delta(m)>p_k\). This sieve is performed in blocks, and we maintain a binary tree during the sieving intervals to maintain intermediate results after sieving up to a given prime in real-time. This allows us to find the count of numbers not yet sieved at a certain stage of the sieve in only \(O(\log x)\) time complexity.

Time and Space Complexity of the Algorithm

The time and space complexity are affected by the following three processes:

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

Complexity of Computing P₂(x,y)

We already know the time complexity of this process is \(O\left(\dfrac{x}{y}\log{\log x}\right)\), and space complexity is \(O\left(y\right)\).

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

The time complexity of sieving with block length \(y\) for computing \(W_1,W_2\) is \(O\left(\sqrt{x}\log{\log x}\right)\), and space complexity is \(O\left(y\right)\).

The time complexity required for computing \(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 for 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 for 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 for 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 for 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: Since we need to quickly query the value of \(\phi(u,b)\), we cannot compute it in \(O(1)\) using ordinary sieving. Instead, we maintain a data structure such that each query takes \(O(\log x)\) time, resulting in a time complexity of \(O\left(\dfrac{x}{y}\log x\log\log x\right)\).

For summation: For each term in the sum for \(S_3\), we query the above data structure, totaling \(O\left(\log x\right)\) queries. We also need to compute the number of terms in the sum, i.e., the number of leaves in the binary tree. All leaves are of 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) \]

Total Complexity

The space complexity of this algorithm is \(O\left(y\right)\), and the 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}\), we achieve the optimal time complexity of \(O\left(\dfrac{x^{2/3}}{\log^2 x}\right)\) and space complexity of \(O\left(x^{1/3}\log^3{x}\log{\log x}\right)\).

Some Improvements

We present some improvements here to reduce the algorithm's constant factors and improve its practical efficiency.

In Termination Condition \(2\), we can use a \(z\) instead of \(y\), where \(z\) satisfies \(z>y\). We can prove that the time complexity for computing \(S_3\) can 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 method for checking the computation by varying \(z\).
For clarity, we chose to split the total sum \(S\) at \(x^{1/4}\) when describing the algorithm, but in fact we only need \(p\le \dfrac{x}{pq}<p^2\) to compute. We can utilize this point while keeping the asymptotic complexity unchanged.
Preprocessing with the first few primes \(2,3,5\) can save more time.

References and Further Reading

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

本页面最近更新：，更新历史
发现错误？想一起完善？在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用