はじめに

Meissel–Lehmer algorithm ってな～～～に？

Library-Checker の問題のひとつ, Counting Primes

$N$ 以下の素数の個数を求めてください

$1 \le N \le 10^{11}$

これどうやってやるんですか？

ってことで調べたら出てきました

概要

思考をまとめるためにほとんどwiki写しただけです

素数計数関数を知っていますか？私は知っています

正の実数 $x$ について $\pi(x) = x$ 以下の素数の個数

という関数です

今回の問題は $\pi(N)$ を求めよという問題ですね

このままだと明らかにきついのでちょっと弱めたものを用意します

$p _ i = i$ 番目の素数として

$\pi(x, a) = ( [1, x ] \cap \mathbb{N} ) \setminus \displaystyle\bigcup _ {j = 1}^a p _ j \mathbb{Z}$

とします

？？？？となるかもしれませんが, エラトステネスの篩で $a$ 個目の素数まで篩ったときの素数候補って感じです

つまり, $1$ 以上 $x$ 以下の整数のうち $a$ 番目までの素数の倍数を除いたものです

さらに,

$\phi(x, a) = \left | \pi(x, a) \right |$

$P _ k(x, a) = \pi(x, a)$ のうち素因数をちょうど $k$ 個持っているものの数

とします

ここで, $k$ 個というのは $k$ 種類ということではありません例えば $12 = 2 \times 2 \times 3$ なので $3$ 個です

また, $P _ 0(x, a) = 1$ です ( $1$ は素因数を一つも持たない)

まず, 定義より

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

が成立することはわかると思います

更に,

$\pi(x) = P _ 1(x, a) + a$

が任意の $a$ について成り立ちます

$p _ a$ 以下の素数の倍数を全て篩ったあとに残ってるなかで素因数の数が $1$ のもの $= p _ a$ より大きい素数

よって $x$ 以下の素数の数 $= p _ a$ より大きい素数の数 $+ a$ となります

この二つの式を組み合わせると

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

となります

ということで右辺に出てきたよくわからないのを二つ求められれば終わりです

というのが Meissel–Lehmer algorithm です

後は頑張ってください疲れました

もうひと頑張り

後は簡単ですね！とはいきません

なんか無限級数とか出されても困りますということでまだ続くんじゃ

まずは $\phi(x, a)$ の方から

漸化式のようなものを考えます小さいほうから篩っているので,

$\phi(x, a) = \phi(x, a-1) -$ ( $p _ a$ を最小の素因数とするような整数の数 )

というように考えられます

$p _ a$ を最小の素因数とするような整数の数をどうするかですが,

$p _ a$ の倍数だけを考えてみると, それらを $p _ a$ で割ったもののうち, $a - 1$ までの素数で篩ったものとなります (日本語にすると難しいですね)

ということでこれは $\phi \left( \dfrac{x}{p _ a} , a-1 \right)$ と一致します

したがって

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

という漸化式が得られます

次に $\displaystyle\sum _ {k = 2}^\infty P _ k(x, a)$ です

無限に足していくなんてやってられませんが, よく考えなくても $k$ には上限があります

もっと強くいきましょう

とある $a$ が存在して

$\displaystyle\sum _ {k = 2}^\infty P _ k(x, a) = P _ 2(x, a)$

となってくれたらとてもうれしいですね

素因数を三つ以上持つような合成数を全て篩い落とすには, $a$ をいくつにすればいいでしょうか

ということで $y = \sqrt[3]{x}$ として $a = \pi(y)$ としてみるといいんじゃないでしょうか

そんなわけで $P _ 2(x, a)$ が計算できればいいということになりました

定義からゴリゴリ行きましょう

$P _ 2(x, a) = \Big| \, n \, \big| \, n = p _ b p _ c, \, a \lt b \le c \, \Big |$

と書けますね

このとき $p _ b \le \sqrt{x}$ かつ $p _ c \le \dfrac{x}{p _ b}$ となるはずなので

$P _ 2(x, a) = \displaystyle\sum _ {b = a+1}^{\pi( \sqrt{x} )} \Big| \, n \, \big| \, n = p _ b p _ c, \, b \le c \le \pi \left(\dfrac{x}{p _ b} \right) \, \Big |$

と書けますこれは $c$ の数を数えればいいだけなので

$P _ 2(x, a) = \displaystyle\sum _ {b = a+1}^{\pi( \sqrt{x} )} \left( \pi \left(\dfrac{x}{p _ b} \right) - (b-1) \right)$

と書けて, 等差数列になってる項を括弧から出すと

$P _ 2(x, a) = \dbinom{a}{2} - \dbinom{\pi( \sqrt{x} )}{2} + \displaystyle\sum _ {b = a+1}^{\pi( \sqrt{x} )} \pi \left(\dfrac{x}{p _ b} \right)$

と書けました

はい, ようやく終わりですなんかまた $\pi$ が出てきてしまいましたがもともとの問題よりは小さくなっているので大丈夫です

まとめると

$\pi(x) = a + \phi(x, a) - 1 - P _ 2(x, a) \quad (y = \sqrt[3]{x}, a = \pi(y))$

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

$P _ 2(x, a) = \dbinom{a}{2} - \dbinom{\pi( \sqrt{x} )}{2} + \displaystyle\sum _ {b = a+1}^{\pi( \sqrt{x} )} \pi \left(\dfrac{x}{p _ b} \right)$

となりました

もうひと手間

いい感じに立式することはできましたが計算量はどうなっているでしょうか？

正直よくわからないですが, できればこんなに再帰させないでうまく計算してあげたいです

上で $y = \sqrt[3]{x}$ と書きましたが面倒なので $y$ は $\sqrt[3]{x} \le y \le \sqrt{x}$ である整数とします (最終的に $\sqrt[3]{x} \log ^ ? x$ とします)

まず $P _ 2(x, a)$ ですが, 素因数が二つしかないならなんかもっとうまくやれそうな気分になります

ここで登場するのは尺取り法です素因数が二つということでポインタを二つ使っていい感じにやります

まず, 前計算をします

エラトステネスの篩や線形篩を利用して $\dfrac{x}{y}$ 以下の素数のテーブルを作ります (計算量解析が面倒なので線形篩を使ったことにします)

ついでにそれぞれの整数に対して最小の素因数を求めておきます (これは後で使います)

さて, 定義より $P _ 2(x, a)$ として数えられるような数は, 素因数の小さいほうが $p _ a$ より大きいはずです

そのため素因数の大きいほうは $\dfrac{x}{p _ a}$ より小さいということになりますが, 前計算で $\dfrac{x}{y}$ 以下の素数は全て求めました

ということで, $p _ a \lt p \lt \dfrac{x}{p _ a}$ であるようなそれぞれの素数について, 素の素数がペアになるような素数の数を数えます

これは尺取り法を使えば素数の数に対して線形時間で数えることができますね

ということで $P _ 2(x, a)$ を求めることができました

さて, 残るは $\phi(x, a)$ の計算です, 再帰的に $\phi(m, b)$ を呼び出していくことになりますが

$m \le \dfrac{x}{y}$ のときは後でまとめて処理することにして今は考えないことにします

そうすると, 残りの呼び出されるものは $\phi \left(\dfrac{x}{k}, b\right) \quad k \lt y, b \le a$ です

これは $O(y \pi(y))$ 個ありますが, 素数定理より $\pi(x) \sim \dfrac{x}{\log x}$ という近似ができるので $O\left( \dfrac{y ^ 2}{\log y} \right)$ となります

ということでまぁいい感じに抑えられたということにしておきます

残っているのはさっきためておいた $m \le \dfrac{x}{y}$ の場合です

$\phi(m, b)$ は, 定義より $m$ 以下の自然数のうち最小の素因数が $p _ b$ より大きいものの数です

大きいものより小さいもののほうが数えやすそうですね

最小の素因数が $p _ b$ 以下であるものの数を数えることにします

ここで, だいぶ前に前計算しておいた最小の素因数を使いますうれしいことに, ちょうど $\dfrac{x}{y}$ 以下の整数について求めていました

ということで, $2$ 以上 $m$ 以下の整数の最小の素因数を列挙することはできます

$p _ b$ 以下であるものの数を数えますが, 求めたい $\phi(m, b)$ はたくさんあるので毎回愚直に数えるのではなく効率的に処理したいです

では, $m$ の昇順に見ていくことにしましょう

配列 $d$ を用意して, $d _ i = m$ 以下の整数のうち最小の素因数が $p _ i$ であるようなものの数とすると,

$\phi(m, b) = m - \displaystyle\sum _ {i = 1} ^ {b} d _ i$

と書けます区間和です

区間和なら累積和でできる？いいえ, 更新があるので Binary Indexed Tree (Fenwick Tree) を使いましょう

ということで, $\phi(m, b)$ は $m$ の昇順に見ていくことで合計 $\dfrac{x}{y} \log \dfrac{x}{y}$ で求めることができました

$\phi(x, a)$ 全体ではどうだったかというと $O\left( \dfrac{y ^ 2}{\log y} + \dfrac{x}{y} \log \dfrac{x}{y} \right)$ ですねよくわかりません

しかし, ここで $y = x ^ {1/3} \log ^ {2/3} x$ としてあげます

そうすると

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