線形篩と呼ばれているアルゴリズム・データ構造があります。
「線形時間で前計算できる何々」を「線形何々」と呼ぶのは不適切ですか、もしかして？*1
$\gdef\lpf{\operatorname{lpf}}$
ざっくり紹介
線形篩の構築は、各 $i$ ($2\le i\le n$) に対して 最小素因数 (least prime factor) $\lpf(i)$ を求めるものであり、これを $O(n)$ 時間で行います。 これにより、$i$ の素数判定を $\lpf(i) = i$ でできるわけです。
また、構築の過程で $n$ 以下の素数の列挙も行うので、素数の列挙をしたいときにはそのまま使うこともできます。
構築に関して
下記の考察に基づきます。
整数 $i\ge 2$ に対して、$j\le\lpf(i)$ なる各素数 $j$ に対して $\lpf(i\times j) = j$ が成り立つ。 よって、（配る DP によって）$\lpf(i)$ が求められたら、$j\le\lpf(i)$ なる各素数に対して $\lpf(i\times j)$ が求められる。
DP をしていく際に、$i$ への更新が行われないまま $i$ から配るターンになったら、$i$ が素数であることがわかるので、素数のリストに $i$ を加える。 各 $i$ に対して、更新は $i/{\lpf(i)}$ からしか行われないため、全体として $O(n)$ 時間であることがわかる。
発展編
最小素因数を考えながら DP をできるものがしばしばあります。 たとえば、（重複を含めた）素因数の個数が挙げられます。 $i$ の（重複を含めた）素因数の個数を $\Omega(i)$ とおく*2と、
$$ \Omega(i) = \begin{cases} 0, & \text{if }i = 1; \\ \Omega(i/{\lpf(i)}) + 1, & \text{otherwise} \end{cases} $$ と書けます。
$\gdef\ord{\operatorname{ord}}$
非負整数 $i$ が素数 $p$ で割り切れる回数を $\ord_p(i)$ と書くことにし、$i$ が最小素因数で割り切れる回数を $c(i) = \ord_{\lpf(i)}(i)$ とします。 また、便宜上 $\lpf(1) = 0$ としておきます*3。
$$ c(i) = \begin{cases} 0, & \text{if }i = 1; \\ 1, & \text{if }{\lpf(i)} \ne {\lpf(i/{\lpf(i)})}; \\ c(i/{\lpf(i)}) + 1, & \text{otherwise}. \end{cases} $$
同様の分岐で $+1$ ではなく $\times{\lpf(i)}$ などを考えることで、$\lpf(i)^{c(i)}$ を求めることもできます。
愚直に関して
各 $i$ に対して、愚直に割り続けて回数を調べる方法でも、全体で $O(n)$ 時間を達成できるらしいです。 ${1.77315667\ldots}\times n$ くらいらしいです？ $\eod$
より一般に、乗法的関数 (multiplicative function) $f(n)$ について考えます。乗法的関数とは、下記を満たす関数のことです。
- $f(1) = 1$,
- $\gcd(u, v) = 1 \implies f(uv) = f(u)f(v)$.
$n = p_1^{e_1}\cdots p_k^{e_k}$ と素因数分解できるなら $f(n) = f(p_1^{e_1})\cdots f(p_k^{e_k})$ と書くことができます。
さて、乗法的関数であって以下を満たすものを考えます。
- $f(p)$ は $p$ から計算できる。
- $i\gt 1$ に対して、$f(p^i)$ は $p$, $p^{i-1}$, $f(p^{i-1})$ から計算できる。
このとき、各 $i$ に対する $f(i)$ が DP によって求めることができます。
実装
Euler の $\phi$ 関数や、約数の総和などを求めています。便宜上、$f(0)$ として使う値も渡してみていますが、T::zero() とかにしてもよいかもしれません？
use std::ops::Mul; use num_traits::One; fn main() { let ls = LpfSieve::new(20); println!("{ls:?}"); let phi = ls.dp(0, 1, |&x, p| x * p, |&x, p| x * (p - 1)); println!("{phi:?}"); let phi = ls.mulfn(0, |p| p - 1, |p, _pp, x| x * p); println!("{phi:?}"); let sigma = ls.mulfn(0, |p| p + 1, |p, pp, x| x + p * pp); println!("{sigma:?}"); } #[derive(Debug)] struct LpfSieve { lpf: Vec<usize>, } impl LpfSieve { pub fn new(n: usize) -> Self { let mut lpf: Vec<_> = (0..=n).collect(); let mut primes = vec![]; for i in 2..=n { if lpf[i] == i { primes.push(i); } let lpf_i = lpf[i]; for &j in primes.iter().take_while(|&&j| j <= lpf_i.min(n / i)) { lpf[i * j] = j; } } Self { lpf } } pub fn dp<T>( &self, x0: T, x1: T, mut eq: impl FnMut(&T, usize) -> T, mut gt: impl FnMut(&T, usize) -> T, ) -> Vec<T> { let n = self.lpf.len() - 1; if n == 0 { return vec![x0]; } else if n == 1 { return vec![x0, x1]; } let mut res = vec![x0, x1]; res.reserve(n + 1); for i in 2..=n { let lpf = self.lpf[i]; let tmp = if lpf == self.lpf[i / lpf] { eq(&res[i / lpf], lpf) } else { gt(&res[i / lpf], lpf) }; res.push(tmp); } res } pub fn mulfn<T>( &self, f0: T, mut fp: impl FnMut(usize) -> T, mut fpi: impl FnMut(usize, usize, &T) -> T, // p, p^{i-1}, f(p^{i-1}) ) -> Vec<T> where T: One + Clone, for<'a> &'a T: Mul<&'a T, Output = T>, { let x0 = (f0.clone(), f0.clone(), 1); let x1 = (T::one(), T::one(), 1); let eq = |(lt_f, eq_f, pp): &(T, T, usize), p: usize| { (lt_f.clone(), fpi(p, *pp, eq_f), pp * p) }; let gt = |(lt_f, eq_f, _): &(T, T, usize), p: usize| { (lt_f * eq_f, fp(p), p) }; self.dp(x0, x1, eq, gt) .into_iter().map(|(lt_f, eq_f, _)| lt_f * eq_f).collect() } }
練習問題
Möbius の $\mu$ 関数や約数の個数なども求めてみましょう。
素因数の個数は ABC 368 F など。
おわり
おなかがすきました。
