\(N\) 以下の素数を線形時間 \(O(N)\) で列挙する篩（ふるい）を解説します。Atkin の篩という時間 \(O(N/\log\log N)\)、空間 \(O(N^{1/2+o(1)})\) で動作する、より高速な篩も存在しますが、それと比べて実装が簡単で、更に \(N\) 以下の自然数に対する最小の素因数が列挙されるという利点があります。特に後者の性質は、素因数分解や冪乗の列挙に応用できて面白いと思います。

線形篩の解説

素因数分解したときに \(x = p_1^{e_1} p_2^{e_2} \ldots p_n^{e_n} \ (p_1 < p_2 < \ldots < p_n)\) となる合成数 \(x\) について、エラトステネスの篩だと \(p_1, p_2, \ldots, p_n\) の全ての素因数で \(x\) を篩いに掛けて合成数と判断していました。そこで線形篩では最小素因数 \(p_1\) のみで \(x\) を篩に掛け合成数と判断することで計算量を削減します。これは、\(p_1\) で篩うとき、最小素因数が \(p_1\) 以上の自然数 \(d\) との積 \(p_1 d\) だけを篩う対象にすればよいです。Python のコードにすれば下のようになります。最小素因数は各数に対してただ一つですから、各自然数は１回しか篩われません。つまり、計算量は \(O(N)\) です。

N = 10 ** 5
prime_list = [] # 発見した素数のリスト
lpf = [None] * (N + 1) # 発見した最小素因数。
for d in range(2, N + 1):
    if lpf[d] is None:
        lpf[d] = d
        prime_list.append(d)
    for p in prime_list:
        if p * d > N or p > lpf[d]:
            break
        lpf[p * d] = p

応用

ここからは、線形篩を用いて効率的に解ける、面白い問題を２つ紹介します。

線形篩を用いた冪数の列挙

線形篩を用いると、\(m\) を法とする冪数 \(1^k, 2^k ,\ldots, N^k\) の列挙が \(O\left(N\frac{\log m}{\log N}\right)\) でできます。線形篩で \(N\) 以下の数に対する最小素因数と、\(N\) 以下の素数を \(O(N)\) で列挙しておきます。素数に対しては 冪数を繰り返し二乗法で求めます。 素数定理より \(N\) 以下の素数は \(O(N/\log N)\) 個あるので、この列挙は \(O\left(N\frac{\log m}{\log N}\right)\) で行えます。 合成数 \(x\) はどうすれば良いでしょうか。最小素因数 \(p_1\) を使って部分問題に分解しましょう！ 小さい合成数から順に冪乗を計算していきます。 \(x^k=p_1^k\cdot\left(\frac{x}{p_1}\right)^k\) のうち、 \(p_1\) は素数なので冪乗はすでに計算しています。\(\frac{x}{p_1}\) は \(x\) 未満の素数または合成数なのでこちらもすでに冪乗は計算してあります。併せると素数の累乗を計算しておけば \(O(N)\) で合成数の場合は冪乗が列挙できました。以上で冪乗の列挙が \(O\left(N\frac{\log m}{\log N}\right)\) で行えました。

線形篩を用いた合成数を法とする逆元の列挙

線形篩を用いると、合成数 \(m\) を法とする \(N\) 以下の逆元の列挙が \(O\left(N\right)\) でできます。アルゴリズムは冪乗の列挙と似ています。まず線形篩で \(N\) 以下の数の最小素因数を列挙しておきます。素数の逆元については互除法で求めます。一回あたり \(O(\log N)\) 掛かりますが、\(N\) 以下の素数は \(O(N/\log {N})\) 個しかないので全体では \(O\left(N\right)\) です。合成数 \(x\) についてはその最小素因数 \(p_1\) を用いて \(x^{-1}=p_1^{-1}\cdot(x/p_1)^{-1}\) とすれば \(O(1)\) で求まります。以上で逆元の列挙が \(O\left(N\right)\) でできました。

完全乗法的関数

ここまでの応用例は完全乗法的関数として一般化できます。素因数分解した \(n=\prod p_i^{e_i}\) に対して、\(f(n)=\prod f(p_i)^{e_i}\) が成り立つとき、\(f\) を完全乗法的関数といいます。素数 \(p\) に対して、\(f(p)\) が \(O(\log(N))\) で得られるなら \(f(N)\) が \(O(N)\) で求まります。

乗法的関数

完全乗法的関数の条件を少し弱めた、乗法的関数の場合も同じようなことが言えます（nyaanさんに教えてもらいました！）。素因数分解した \(n=\prod p_i^{e_i}\) に対して、\(f(n)=\prod f(p_i^{e_i})\) が成り立つとき、\(f\) を乗法的関数といいます。素冪 \(p^e\) に対して、\(f(p^e)\) が \(O(1)\) で得られるとします。各数 \(n\) の約数のうち、最小素因数 lpf(n) の冪であってべき指数が最大であるもの \(a_n\) は $$\begin{eqnarray} a_n =\begin{cases}
a_{n/\mathrm{lpf}(n)} \mathrm{lpf}(n) & (\mathrm{lpf}(n)=\mathrm{lpf}(n/\mathrm{lpf}(n)) \\
\mathrm{lpf}(n) & (\text{otherwise})
\end{cases}
\end{eqnarray}$$ に従うと、\(O(N)\) で列挙できます。あとは完全乗法的関数の場合と同様にすると \(f(N)\) が \(O(N)\) で求まります。

例えば、オイラーのトーティエント関数は乗法的関数ですね。

アルゴリズムの弱点：メモリ消費と実速度

実用上、最小素因子が必要なければ、定数倍が大きいため、線形篩の速度・メモリ消費量ともにエラトステネスの篩に及びません。

線形篩のメモリ消費量はエラトステネスの8倍近くあります。エラトステネスの篩は長さ \(n\) の Boolean 型（1byte）の配列が必要だったのに対し、 線形篩では長さ \(n\) の最小素因数（4byte）の配列と長さ \(O(n/\log(n))\) の素数（32bit）のリストが必要になります。長さ \(1.25 \times 10^8\) の32bit整数の配列を２つ取ると、1GB になりますから、この辺りが線形篩の使える \(n\) の上限だと思います。更に、エラトステネスの篩では見つけた素数をハードディスクに書き込むことで、RAM の空間計算量 \(O(\sqrt{N})\) にできました（How?）が、線形篩では難しいと思います。

また、エラトステネスの篩 \(O(n \log\log(n))\) と 線形篩 \(O(n)\) では \(\log \log (n)\) だけ、エラトステネスの篩の方が重いですが、実用上のこの因子は大きくありません。因子の大きさは \(n\) 以下の素数の逆数和に等しく、\(n=10^8\) のとき \(\sum_{p \leq n} \frac{1}{p}=3.17\cdots\) です。

線形篩とエラトステネスの篩（一つ目の実装）で速度を測定しました（LIFEBOOK A550/B（Core i5, 4GB）10回平均）。\(n \leq 10^7\) ではエラトステネスの篩の方が定数倍が小さいことが分かります。最小素因数が必要なければ、エラトステネスの篩を使う方がよいでしょう。

余談：最大素因数で設計

最小素因数の代わりに最大素因数で線形篩を設計することもできます。

\(N\) 以下の各正整数 \(s=2,3,\ldots,N\) に対して、$$ \{sp : \mathrm{mpf}(s) \leq p, p\text{ is prime}\}$$ を合成数として確定させればよいです（mpfは最大素因数）。しかし、\(s\) 未満まで処理を終えた時点では、\(p > s\) なる素数が分かりません。そこで、\(p\) と \(s\) で和を取る順序を交換する（かならいさんのアイディア[1]）と
$$ \bigcup_{s} \{ps : \mathrm{mpf}(s) \leq p , p > s,p\text{ is prime}\} = \bigcup_{p\text{ is prime}}\{ps : p > s\}$$となり、計算できる形になります。計算量は変わらず \(O(N)\) ですが、コードは１行増えます（かならいさんの実装）。

参考文献

1. かならいさんとのやり取り