Meissel–Lehmer algorithm

Let ${\displaystyle p_{1},p_{2},\ldots ,p_{n}}$ be the first n primes. For a natural number a ≥ 1, define

${\displaystyle \varphi (x,a):=\left|\left\{n\leq x:p|n\implies p>p_{a}\right\}\right|,}$

which counts natural numbers no greater than x with all prime factors greater than ${\displaystyle p_{a}}$. Also define for a natural number k,

${\displaystyle P_{k}(x,a):=\left|\left\{n\leq x:n=q_{1}q_{2}\cdots q_{k},~{\text{with}}~q_{1},\ldots ,q_{k}>p_{a}\right\}\right|,}$

which counts natural numbers no greater than x with exactly k prime factors, all greater than ${\displaystyle p_{a}}$. With these, we have

${\displaystyle \varphi (x,a)=\sum _{k=0}^{\infty }P_{k}(x,a),}$

where the sum only has finitely many nonzero terms because ${\displaystyle P_{k}(x,a)=0}$ when ${\displaystyle p_{a}^{k}>x}$. Using the fact that ${\displaystyle P_{0}(x,a)=1}$ and ${\displaystyle P_{1}(x,a)=\pi (x)-a}$, we get

${\displaystyle \pi (x)=\varphi (x,a)+a-1-\sum _{k=2}^{\infty }P_{k}(x,a),}$

which proves that one may compute ${\displaystyle \pi (x)}$ by computing ${\displaystyle \varphi (x,a)}$ and ${\displaystyle P_{k}(x,a)}$ for k ≥ 2. This is what the Meissel–Lehmer algorithm does.

For k = 2, we get the following formula for ${\displaystyle P_{k}(x,a)}$:

${\displaystyle {\begin{aligned}P_{2}(x,a)&=\left|\left\{n:n\leq x,~n=p_{b}p_{c},~{\text{with}}~a<b\leq c\right\}\right|\\&=\sum _{b=a+1}^{\pi (x^{1/2})}\left|\left\{n:n\leq x,~n=p_{b}p_{c},~{\text{with}}~b\leq c\leq \pi \left({\frac {x}{p_{b}}}\right)\right\}\right|\\&=\sum _{b=a+1}^{\pi (x^{1/2})}\left(\pi \left({\frac {x}{p_{b}}}\right)-(b-1)\right)\\&={\binom {a}{2}}-{\binom {\pi (x^{1/2})}{2}}+\sum _{b=a+1}^{\pi (x^{1/2})}\pi \left({\frac {x}{p_{b}}}\right).\end{aligned}}}$

For k ≥ 3, the identities for ${\displaystyle P_{k}(x,a)}$ can be derived similarly.[1]

With the starting condition

${\displaystyle \varphi (x,1)={\Big \lceil }{\frac {x}{2}}{\Big \rceil }=\lfloor x\rfloor -{\Big \lfloor }{\frac {x}{2}}{\Big \rfloor },}$

and the recurrence

${\displaystyle \varphi (x,a)=\varphi (x,a-1)-\varphi \left({\frac {x}{p_{a}}},a-1\right),}$

each value for ${\displaystyle \varphi (x,a)}$ can be calculated recursively.

The only thing that remains to be done is evaluating ${\displaystyle \varphi (x,a)}$ and ${\displaystyle P_{k}(x,a)}$ for k ≥ 2, for certain values of x and a. This can be done by direct sieving and using the above formulas.