別に、おふとんから出たくない朝とかに数えてもいいと思います。
\(n\) 以下の素数の個数を数えたいときどうしますか？ 各整数 \(2\le i\le n\) に対して素数判定して数えます？ 素数判定を試し割り法でやると \(\Theta(n\sqrt{n})\)、線形篩*1を使えば \(\Theta(n)\) とかですね。 これを \(o(n)\) でできたらうれしくないですか？という話です。
今回は \(O(\sqrt{n})\) space, \(O(n^{3/4}/\log(n))\) time の話がメインです。 定数倍があまり軽くない \(O(n^{2/3})\) space, \(O(n^{2/3})\) time の話も少し触れます。
（追記：2021/06/05）少し難しめかもですが、アイディア自体は水色もあれば理解できると思います。
（追記：2022/06/29）勉強会用に作ったスライドを公開しました：https://rsk0315.github.io/slides/prime-counting.pdf
前提知識
以下、\(n\) 以下の個数の素数を \(\pi(n)\) と書きます。
導入
素数の個数を数えるにあたって、素数を実際に列挙する必要はないというのが重要です。 実際に列挙していては、「与えられた素数の次の素数はなに？」というのに \(\Theta(1)\) 時間で答えられても、\(\Theta(n/\log(n))\) 時間かかってしまいます*2。
ここで、Eratosthenes の篩について考え直します。次のようなアルゴリズムです。
- \(1\) から \(n\) の bool 配列 \(a_i\) を用意する。\(a_1\) のみ false とする。
- \(i = 2, 3, \dots, \lfloor\sqrt{n}\rfloor\) に対して、次の処理をする。
- \(a_i\) が false なら次の \(i\) に移る。
- \(j = i, i+1, \dots, \lfloor n/i\rfloor\) に対して、\(a_{ij}\) を false にする。
let mut is_prime = vec![true; n + 1]; is_prime[1] = false; for i in (2..=n).take_while(|&i| i * i <= n) { if !is_prime[i] { continue; } for j in i..=n / i { is_prime[i * j] = false; } }
ここで、各 \(i\) に対して 2.2 のステップでいくつの \(a_{ij}\) が true から false になるかについて考えます。 この個数をうまく管理できれば、篩い終えた後にいくつ true になっているか、すなわち素数の個数がわかることになります。
この手法は、界隈では Lucy DP などと呼ばれることもあるようです*3。
解説
（追記 2021/05/19）画像とかをつけるとわかりやすくなる気がしたので、そのうちつけると思います。
\(S(v, p)\) を、\(i=p\) まで処理を終えた後に \(a_j\) が true である \(2\le j\le v\) の個数と定義します。 \(p\cdot p \gt v\) であれば \(i=p-1\) のときから配列は変化しないので、\(S(v, p) = S(v, p-1)\) とわかります。\(p\) が素数でない場合も同じです。
\(p\) が \(p\cdot p\le v\) なる素数の場合はむずかしいのでちょっと考えます。 \(j = p, p+1, \dots, \lfloor v/p\rfloor\) に対して \(a_{pj}\) が false になりますが、\(j\) が \(p\) 未満の素数で割り切れるものはすでに篩われています。よって、篩われる個数は「\(p, p+1, \dots, \lfloor v/p\rfloor\) のうちで篩われていないもの」であり、これは \(S(\lfloor v/p\rfloor, p-1)-S(p-1, p-1)\) です。すなわち、次が成り立ちます。 \[S(v, p) = S(v, p-1) - (S(\lfloor v/p\rfloor, p-1) - \pi(p-1) ).\]
さて、\(\pi(n)\) は \(S(n, \lfloor\sqrt{n}\rfloor)\) と等しいので、これを求めたいわけですが、\(S(v, \bullet)\) を求めるためには（その時点で見つかっている素数の個数と）\(S(\lfloor v/\bullet\rfloor, \bullet)\) の形で表せるものさえわかっていれば十分とわかります。 \(\lfloor \lfloor n / p\rfloor / q\rfloor = \lfloor n / pq\rfloor\) が重要です*4。
よって、\(v = 1, 2, \dots, \lfloor\sqrt{n}\rfloor\) と \(v = \lfloor n/1\rfloor, \lfloor n/2\rfloor, \dots, \lfloor n/\lfloor\sqrt{n}\rfloor\rfloor\) についてのみ管理します（サイズ \(O(\sqrt{n})\) の配列ふたつで管理できます）。 以下のコード中では \(\texttt{smalls}[i] = S(i, p)\)、\(\texttt{larges}[i] = S(\lfloor n/i\rfloor, p)\) とします。
fn main() { assert_eq!(prime_pi(10), 4); assert_eq!(prime_pi(100), 25); assert_eq!(prime_pi(100000000000), 4118054813); } fn prime_pi(n: usize) -> usize { let n_sqrt = floor_sqrt(n); let mut larges: Vec<_> = (0..=n_sqrt).collect(); for i in 1..=n_sqrt { larges[i] = n / i - 1; } let mut smalls: Vec<_> = (0..n / n_sqrt).map(|x| x.saturating_sub(1)).collect(); for p in 2..=n_sqrt { if p < smalls.len() { if smalls[p] <= smalls[p - 1] { continue; } } else { if larges[n / p] <= smalls[p - 1] { continue; } } let pc = smalls[p - 1]; let q = p * p; for i in (1..=n_sqrt).take_while(|&i| n / i >= q) { // vi = n / i // dp[n / i] -= dp[n / i / p] - pc; let ip = i * p; let cur = *larges.get(ip).unwrap_or_else(|| &smalls[n / ip]) - pc; *larges.get_mut(i).unwrap_or_else(|| &mut smalls[n / i]) -= cur; } for i in (1..smalls.len()).rev().take_while(|&i| i >= q) { // vi = i // dp[i] -= dp[i / p] - pc; let cur = smalls[i / p] - pc; smalls[i] -= cur; } } larges[1] } fn floor_sqrt(n: usize) -> usize { if n <= 1 { return n; } let mut lo = 1; let mut hi = n; while hi - lo > 1 { let mid = lo + (hi - lo) / 2; match mid.overflowing_mul(mid) { (x, false) if x <= n => lo = mid, _ => hi = mid, } } lo }
計算量解析
i に関するループがふたつありますが、ひとつめの計算量は次の値で抑えられます。以下の \(i\) は p に対応します。
\[
\begin{aligned}
\int_2^{\sqrt{n}} \min\{\sqrt{n}, n/i^2\}\, di
&= \int_2^{\sqrt[4]{n}} \sqrt{n}\,di + \int_{\sqrt[4]{n}}^{\sqrt{n}} n/i^2\,di\\
&= \sqrt{n}\cdot\sqrt[4]{n} + O(\sqrt{n}) + (\sqrt[4]{n}-1)\cdot\sqrt{n} \\
&= O(n^{3/4}).
\end{aligned}
\]
ふたつめも。 \[ \begin{aligned} \int_2^{\sqrt{n}} \max\{\sqrt{n}-i^2, 0\}\, di &= \int_2^{\sqrt[4]{n}} (\sqrt{n}-i^2)\, di + \int_{\sqrt[4]{n}}^{\sqrt{n}} 0\,di \\ &= 1/3\cdot (2n^{3/4} - 3\sqrt{n} + 1) \\ &= O(n^{3/4}). \end{aligned} \]
合わせて \(O(n^{3/4})\) です。実際には、素数の分布から \(O(n^{3/4}/\log(n))\) になっている気がします。
工夫
オーダーの改善ではないはずで、定数倍高速化だと思います。
上記のアルゴリズムで、篩う \(i\) を \(\lfloor\sqrt{n}\rfloor\) までではなく \(\lfloor \sqrt[4]{n}\rfloor\) までにしてみます。 このとき、篩われずに残っている合成数は、たかだか 3 つの素数の積で表せることがわかります*5。4 つかけると \(n\) を越えちゃうので。 このことから、全体をふたつの処理に分けて行うことが考えられます。
合成数の個数を数える方について考えます。 篩われずに残っている整数（以前篩うのに使ったのは除く）を昇順に見ていきます。\(i_1, i_2, \dots\) とします。 各 \(k\) に対して、\(i_k\) の倍数である \(n\) 以下の整数の個数を足します。 また、\(k'\lt k\) に対して \(i_k\cdot i_{k'}\) の倍数である整数はすでに数えられているので、引きます。
ここで、\(i_k\) の倍数である \(n\) 以下の整数は、\(S(\lfloor n/i_k\rfloor, \lfloor\sqrt[4]{n}\rfloor) - \pi(i_k -1)\) と表せます。\(i_k i_{k'}\) の方は \(S(\lfloor n/i_k i_{k'} \rfloor, \lfloor\sqrt[4]{n}\rfloor) - \pi(i_{k'} -1)\) です。
実装
アイディアは上記の通りですが、実装上の工夫が少しあるので下で説明します。
変数の概要は次の通りです。
- \(\texttt{smalls}[i] = S(i, p)\)
- \(\texttt{roughs}[i]\)：消されていない \(i\) 番目の整数
- \(1\) は残る、篩うのに使った素数は消える
- サイズは \(\texttt{s}\) で管理される
- \(\texttt{larges}[i] = S(\lfloor n/\texttt{roughs}[i]\rfloor, p)\)
- \(\texttt{pc} = \pi(p-1)\)
- \(\texttt{skip}[i]\)：\(i\) が篩われていれば true
- ただし偶数については更新しない
\(p\) の変化に伴って更新されるので場所によって ±1 がありますが、各ループの最後で上記の定義に合うように管理します。 最初のループの開始前においては、\(p = 2\) での値です。
fn main() { assert_eq!(prime_pi_fast(100000000000), 4118054813); } fn prime_pi_fast(n: usize) -> usize { if n <= 3 { return n.saturating_sub(1); } let v = floor_sqrt(n); let mut smalls: Vec<_> = (0..=v).map(|i| (i + 1) / 2).collect(); let mut s = (v + 1) / 2; let mut roughs: Vec<_> = (0..s).map(|i| 2 * i + 1).collect(); let mut larges: Vec<_> = (0..s).map(|i| (n / (2 * i + 1) + 1) / 2).collect(); let mut skip = vec![false; v + 1]; let mut pc = 0; for p in (3..=v).step_by(2) { if skip[p] { continue; } let q = p * p; pc += 1; if q * q > n { break; } skip[p] = true; for i in (q..=v).step_by(2 * p) { skip[i] = true; } let mut ns = 0; for k in 0..s { let i = roughs[k]; if skip[i] { continue; } let d = i * p; let x = if d <= v { larges[smalls[d] - pc] } else { smalls[n / d] }; larges[ns] = larges[k] + pc - x; roughs[ns] = i; ns += 1; } s = ns; let mut i = v; for j in (p..=v / p).rev() { let c = smalls[j] - pc; let e = j * p; while i >= e { smalls[i] -= c; i -= 1; } } } let roughs = roughs; let pc = pc; let mut res: usize = larges[0] + (s + 2 * (pc - 1)) * (s - 1) / 2 - larges[1..s].iter().sum::<usize>(); for l in 1..s { let q = roughs[l]; let m = n / q; let e = smalls[m / q] - pc; if e <= l { break; } let t: usize = roughs[l + 1..=e].iter().map(|&r| smalls[m / r]).sum(); res += t - (e - l) * (pc + l - 1); } res } fn floor_sqrt(n: usize) -> usize { /* 上と同じなので省略 */ }
実装の説明
冷静になると自明かも。わからなければ、下記の説明を読みつつ、各箇所で値を適宜出力したりしながら確かめるとよいと思います。
roughs[ns] の更新
roughs[..s] が篩われていない整数です。
ns が次のループにおける s です。
ns = 0 で初期化しておき、roughs[ns] に値が入ったら ns += 1 していきます。前から詰め直しているということですね。
larges[smalls[d] - pc] についての説明
smalls[d] - pc は「篩われずに残っている整数のうち d は何番目か？」を意味します*6。
smalls[d] - pc は「d 以下にある篩われていない整数」から「篩うのに使った整数」を除いたものの個数です。
ただし、この段階ではまだ pc は \(\pi(p -1)-1\) であることに注意です。これに、\(1\) のぶんの個数を足して、所望の値であることがわかります。
roughs[それ] が篩われていない d 番目の整数であることと larges の定義から、結局 \(n/d\) 以下で篩われていない整数の個数になります。
smalls[i] の更新
j = i / p を保って更新されるようにします。
i >= j * p である間 smalls[i] -= smalls[j] - pc で更新され、smalls[i] -= smalls[i / p] - pc になっています。
res の初期化
初期化の際に、（上で説明した）\(i_k\) の倍数を先に引いています。等差数列の公式などを思い出すと吉です。
res の更新
q(\(i_{k'}\)) を固定qより大きいr(\(i_k\)) を固定n / (q * r)として表せる篩われていない整数の個数を求める
\(\sqrt[4]{n} \lt i_{k'} \lt i_k\) から、\(\lfloor n / i_k i_{k'}\rfloor \le \lfloor\sqrt{n}\rfloor\) であり、m / r = n / (q * r) で smalls にアクセスしても問題ありません。
3. で求めた個数から、（前半パートですでに数えた）\(\sqrt[4]{n}\) 以下の素数を含む個数を除く処理が -= (e - l) * (pc + l - 1) です。
\(\lfloor n / i_k i_{k'}\rfloor\) として表せるうちで数えたい最小値は \(i_k'\) なので、2. で固定するのは \(\lfloor n / i_{k'}^2\rfloor\) 以下の範囲です。
計算量解析
前半パートでは、各 \(i \le n^{1/4}\) の素数に対して \(O(\sqrt{n})\) 程度の処理をするので \(O(n^{3/4}/\log(n))\) time です。篩の部分も \(O(\sqrt{n}\log(\log(n))) \subset o(n^{3/4}/\log(n))\) です。 後半パートでは、雑に見積もって \(n^{1/4} \lt i \le n^{1/2}\) の各素数に対して \(O(n/i^2)\) 程度の処理をしています。 \[ \int_{n^{1/4}}^{n^{1/2}} \frac{n}{i^2}\, di = n^{3/4} - n^{1/2}. \] であり、素数の分布から \(O(n^{3/4}/\log(n))\) になりそうです。
\(n = 10^{13}\) でも 1 秒前後くらいで終わりました。 \(\log(n)\) で割れる部分のちゃんとした（お気持ちでない）証明ができたら追記します。
その他
\(\sqrt[3]{n}\) とかで区切って \(O(n^{2/3})\) 時間にできるのかな？ 区切り方を変えれば \(O(n^{3/4-\varepsilon})\) にはできそうな気もするけど、定数倍とか諸々のことはわかりません。
Dirichlet 積などを使った方法 とかもあります。軽めの定数倍では実装できませんでした。
参考文献
- Lucy さんの投稿
- math314 さんの記事
- Library Checker の提出たち。Min さん
- この記事で載せた二つ目の実装は、これらの提出コードに基づきます。
おわり
おやすみなさいにゃ。
（記事を書いていたら眠れなくなりました）
