Meissel–Lehmer algorithm

Given ${\displaystyle a\in \mathbb {N} }$, one may define the following functions: Firstly,

${\displaystyle \varphi (x,a):=\left|\left([1,x]\cap \mathbb {N} \right)\setminus \bigcup _{j=1}^{a}p_{j}\mathbb {Z} \right|;}$

this function measures the set ${\displaystyle [1,x]\cap \mathbb {N} }$ (closed interval) when one has sieved out multiples of the first ${\displaystyle a}$ primes ${\displaystyle p_{1},\ldots ,p_{a}}$ (including these primes themselves); that is, ${\displaystyle p_{1},p_{2},p_{3},\ldots =2,3,5,\ldots }$ is the sequence of prime numbers in increasing order.

We may further split that up with the functions

${\displaystyle P_{k}(x,a):=\left|\left[\left([1,x]\cap \mathbb {N} \right)\setminus \bigcup _{j=1}^{a}p_{j}\mathbb {Z} \right]\cap \{n\mid n{\text{ has exactly }}k{\text{ prime factors}}\}\right|,}$

${\displaystyle k\in \mathbb {N} _{0};}$

these measure the set ${\displaystyle \left([1,x]\cap \mathbb {N} \right)\setminus \bigcup _{j=1}^{a}p_{j}\mathbb {Z} }$ but considers only the numbers that have exactly ${\displaystyle k}$ prime factors (this is well-defined by the fundamental theorem of arithmetic). With these, we have

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

the sum on the right is finite, since numbers ${\displaystyle \leq x}$ have, for instance, ${\displaystyle \leq \lfloor x\rfloor }$ prime factors.

The identities

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

prove that one may compute ${\displaystyle \pi (x)}$ by computing ${\displaystyle \varphi (x,a)}$ and ${\displaystyle P_{k}(x,a)}$, ${\displaystyle k\geq 2}$. And this is what the Meissel–Lehmer algorithm does.

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

${\displaystyle {\begin{aligned}P_{2}(x,a)&=\left|\left\{n\mid n\leq x,n=p_{b}p_{c},a<b\leq c\right\}\right|\\&=\sum _{b=a+1}^{\pi (x^{1/2})}\left|\left\{n~{\big |}~n=p_{b}p_{c},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)\\&=\sum _{b=a+1}^{\pi (x^{1/2})}\pi \left({\frac {x}{p_{b}}}\right)+{\binom {a}{2}}-{\binom {\pi (x^{1/2})}{2}}.\end{aligned}}}$

For ${\displaystyle k=3}$, there is a similar expansion.[2]

Using the formula

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

${\displaystyle \varphi (x,a)}$ may be expanded. Each summand, in turn, may be expanded in the same way, so that a very large alternating sum arises.

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