エラトステネスの篩

$N$ 以下の素数を $O(N\log\log(N))$ で列挙するアルゴリズムである、エラトステネスの篩（ふるい）を解説します。

1 アルゴリズム
2 計算量解析
3 $\sqrt{N}$ で十分：定数倍改善
4 区間篩
5 発展形

アルゴリズム

アルゴリズムは非常にシンプルです。

次の操作を可能な限り繰り返します。
1. $2$ から $N$ のうち、まだ素数とも合成数とも分かっていない、最小の自然数 $x$ を取り出します。iiより、$x$ 未満の素数を因子に持つ合成数は全て確定しているので、$x$ は素数確定です。
2. $x$ 以外の $x$ の正の倍数 $2x, 3x, \ldots$ を合成数として確定させます。

つまり、素数 $2$ の倍数 $4,6,8,\ldots$ を取り除き、次に素数 $3$ の倍数 $6,9,12,\ldots$ を取り除き……というようになります。Python のコードは次の通りです：

is_prime = [True] * (N + 1)
is_prime[0] = False
is_prime[1] = False
for i in range(2, N + 1):
    if not is_prime[i]:
        continue
    for j in range(2 * i, N + 1, i):
        is_prime[j] = False

計算量解析

アルゴリズムのシンプルさに反して、計算量解析は少し難しいです（ですが、計算量が軽いのがこのアルゴリズムの凄いところです）。素数 $p$ の倍数 $2p, 3p, \ldots$ を合成数として取り除くのに必要な計算量は $O(N/p)$ です。各素数 $p$ について足し合わせて、全体の計算量は $O(N\sum_{p \leq N} \frac{1}{p})$ になります。

厳密な $O(\sum_{p \leq N} \frac{1}{p})$ の評価は込み入った計算になるので、フェルミ推定します。カギは素数定理です：

素数定理：
$x$ 以下の素数の個数 $\pi(x)$ について、$$\pi(x)=\Theta\left(\int_2^x \frac{1}{\log{t}} dt\right)$$ が成り立つ。

素数定理より、$t$ が素数である確率密度は $\frac{1}{\log{t}}$ で近似できます。積分で素数の逆数和を近似すると$$\sum_{p \leq N} \frac{1}{p} \approx \int_{2}^{N} \frac{1}{t\log{t}} dt = \log \log(N)-\log \log(2) $$です。従って、厳密な証明ではないですが、エラトステネスの篩の計算量が $O(N \log \log(N))$ だろうということが分かりました。実際、この推定は正しいことが知られています。

$N=10^8$ のとき、$\sum_{p \leq N} \frac{1}{p}=3.17\cdots$ と $O(\log\log(N))$ の因子は実用上小さく抑えられます。この計算量の軽さとアルゴリズムのシンプルさがエラトステネスの篩の魅力だと思います。

$\sqrt{N}$ で十分：定数倍改善

$N$ 以下の合成数は $\sqrt{N}$ 以下の素因数を持つことから、次の二つの計算量改善ができます：

1. 篩に使う素因数 $p$ は $p \leq \sqrt{N}$ で十分。
2. 素数 $p$ によって篩う合成数は $p^2,p(p+1),p(p+2),\ldots$ のみで十分。

$\sum_{p \leq \sqrt{N}} \frac{1}{p} = O(\log\log(\sqrt{N})) = O(\frac{1}{2} \log\log(N))$ ですから、改善1は約半分の定数倍改善です。改善2で削減できた計算量は、積分を長方形で近似して $\sum_{p \leq \sqrt{N}} p \approx O(\sqrt{N}\frac{\sqrt{N}}{\log\sqrt{N}})=O(N/\log{N})$ だと推定できます。どちらもよく用いられる定数倍改善です。

Pythonによるコードは次の通りです：

is_prime = [True]*(N + 1)
is_prime[0] = False
is_prime[1] = False
prime_list = []
for i in range(2, N+1):
    if not is_prime[i]:
        continue
    if i * i >= N:
        break
    for j in range(i * i, N + 1, i):
        is_prime[j] = False

$n=10^7$ で計測（LIFEBOOK A550/B（Core i5, 4GB）10回平均）すると、改善前後で速度が２倍になりました。

$n$	改善前	改善後
$n=10^7$	3.719 sec	1.875 sec

$N^{1/2}$ までの素数の試し割りでよいことから、区間を $N^{1/2}$ ずつ区切れば、空間計算量 $O(N^{1/2})$ で、エラトステネスの篩ができます。ただし、$[iN^{1/2}, (i+1)N^{1/2})$ での計算結果は $i$ をインクリメントすると捨てる必要があるので、使用場面は限定的だと思います。素数カウントでの使用例があります。

区間篩

$[n-x,n]$ の範囲の素数を見つけるエラトステネスの篩を特別に区間篩と呼びます。$p \leq \sqrt{n}$ までの素因数で、$[n-x,n]$ の範囲の素因数を列挙できました。$p \leq \sqrt{n}$ までの素因数を見つけるのに $O(\sqrt{n}\log\log(n))$、$[n-x,n]$ での素数判定に $O(x\log\log(n))$ かかります。

$w \leq \sqrt{N}$ とします。$[1,w),[w,2w),[2w,3w),\ldots $ のそれぞれで区間篩を行い、見つけた素数をハードディスクに書き込むと、$[1,N]$ でのエラトステネスの篩が時間 $O(N \log\log(N))$ 、メモリ（RAM）$O(\sqrt{N})$ でできます。これは、$N$ が大きいときに有効です。

発展形

合成数が一度しか篩に掛けられないようにすることで、計算量を $O(N)$ にできます：線形篩

素数か判断する数を $2$ と奇数に限定すれば、定数倍高速化できます。更に、素数か判断する数を $2,3$ と $6k+1,6k+5$ の形に限定することもできます。この考えを発展させると wheel sieve になり、計算量は $O(N)$ です。wheel sieve と線形篩を合体すると、$O(N/\log\log(N))$ になります。