えびちゃんの日記

えびちゃん（競プロ）の日記です。

素因数分解を <O(n), O(log(n)/log(log(n)))> で行う

もしかすると全人類知っていたかもしれません。

$\gdef\lpf{\operatorname{lpf}}$ $\gdef\ord{\operatorname{ord}}$

導入

以下、整数 $i$ の最小素因数 (least prime factor; LPF) を $\lpf(i)$ と書きます。 たとえば $240 = 2^4 \times 3 \times 5$ なので $\lpf(240) = 2$ です。

文献や文脈によって、$\lpf(1) = 1$ としたり $\lpf(1) = \infty$ としたりしますが、ここでは $\lpf(1) = 1$ としておきます。$1$ は素数ではないですがよしなに。

$n$ 以下の各正整数 $i$ に対して $\lpf(i)$ が計算できていたとすると、下記のアルゴリズムによって $n$ の素因数分解を worst $\Theta(\log(n))$ time で得ることができます。

/// `n` の素因数（重複あり）の列を返す。
fn factors_dup(mut n: usize, lpf: &[usize]) -> Vec<usize> {
    let mut res = vec![];
    while n > 1 {
        res.push(lpf[n]);
        n /= lpf[n];
    }
    res
}

$n \gt 1$ に対して $\lpf(n) \ge 2$ なのでループごとに $n$ が半分以下になり、ループ回数は $\log_2(n)$ 回で抑えられます。 逆に、たとえば $n = 2^k$ に対して $k = \log(n)$ 回のループが回るため、最悪ケースは $\Theta(\log(n))$ time となります。

これは osa_k 法などと呼ばれています。

何度か言及していますが、えびちゃんは Human Resource Machine の Prime Factory (YEAR 40) を遊んでいるときに自分で思いつけたアルゴリズムなので、結構愛着が強いです。

Eratosthenes の篩を bool を扱う DP として見ると「bool を扱う DP は無駄が多い」というのを思い出し、LPF も求められるじゃんという気持ちになりやすいかもしれません。

前計算

Eratosthenes の篩のようにやると、LPF の配列は $\Theta(n \log\log(n))$ 時間で求められます。 合成数であることがわかったときに LPF がわかります。

let mut lpf: Vec<_> = (0..=n).collect();
for i in 2..=n {
    for j in i..=n / i {
        if lpf[i * j] == i * j {
            lpf[i * j] = i;
        }
    }
}

ここは、線形篩を用いることで $\Theta(n)$ 時間にできます。 えびちゃんのライブラリ や、その参考文献などが役に立つかもしれません。

方針・考察としては、

  • lpf[i] = 1 で初期化しておく。
  • lpf[i] 以下の各素数 j に対して、i * j の LPF は j とわかる。
    • j は素数のリストから順に見ていく。
    • lpf[i * j] = j で更新する。
  • 2 以上の i を順に処理していき、lpf[i] == 1 であれば i は素数とわかる。
    • i を素数のリストに追加する。
    • lpf[i] = i で更新する。

です。lpf[i] の更新は i / lpf[i] から高々一度ずつ行われるだけなので、全体の計算量が $O(n)$ になります。

let mut lpf = vec![1; n + 1];
let mut primes = vec![];
for i in 2..=n {
    if lpf[i] == 1 {
        lpf[i] = i;
        primes.push(i);
    }
    let lpf_i = lpf[i];
    for &j in primes.iter().take_while(|&&j| j <= lpf_i && i * j <= n) {
        lpf[i * j] = j;
    }
}

高速化

ここまでで $\angled{O(n), O(\log(n))}$ になっています。これは多くの人が知っていそうです。 線形篩なしで osa_k 法だけを知っている人だと $\angled{\Theta(n\log(\log(n))), O(\log(n))}$ かもしれません。

さて、$n$ が持つ素因数の種類数は $\log(n)/\log(\log(n))$ 個程度であることが知られています。

素因数分解を (素数, その次数) の形式で出力することにすれば、列の長さは $O(\log(n)/\log(\log(n)))$ なので、そうすることを考えます。

さて、$\lpf(i)$ だけでなく、$\ord_{\lpf(i)}(i)$ と $\lpf(i)^{\ord_{\lpf(i)}(i)}$ の配列も管理することを考えます。 ここで、$\ord_p(a)$ は、素数 $p$ が正整数 $a$ を割り切る回数を表します*1。

たとえば、$1200 = 2^4 \times 3\times 5^2$ なので、$\lpf(1200) = 2$, $\ord_{\lpf(1200)}(1200) = 4$, $\lpf(1200)^{\ord_{\lpf(1200)}(1200)} = 2^4 = 16$ です。

関数呼び出しが多かったり下つきだったり上つきだったりで過剰に複雑に見えますが、やっていることは単純です。

これがあると、下記のアルゴリズムで素因数分解を行えます。

/// `n` の素因数分解を返す。
fn factors(mut n: usize, lpf: &[usize], lpf_ord: &[u32], lpf_pow: &[usize]) -> Vec<(usize, u32)> {
    let mut res = vec![];
    while n > 1 {
        res.push((lpf[n], lpf_ord[n]));
        n /= lpf_pow[n];
    }
    res
}

各変数は、

  • lpf[n]: $\lpf(n)$
  • lpf_ord[n]: $\ord_{\lpf(n)}(n)$
  • lpf_pow[n]: $\lpf(n)^{\ord_{\lpf(n)}(n)}$

です。

ループごとに lpf[n] は相異なるので、$n$ の素因数の種類数でループ回数が抑えられます。 最悪ケースは素因数が相異なるときで、worst $\Theta(\log(n)/\log(\log(n)))$ time です。

あとは lpf_ord と lpf_pow を求めるパートですが、これは lpf を利用して $O(n)$ 時間で作ることができます。

let mut lpf_ord = vec![0; n + 1];
let mut lpf_pow = vec![1; n + 1];
for i in 2..=n {
    let j = i / lpf[i];
    lpf_ord[i] = if lpf[i] == lpf[j] { lpf_ord[j] + 1 } else { 1 };
    lpf_pow[i] = if lpf[i] == lpf[j] { lpf_pow[j] * lpf[i] } else { lpf[i] };
}

実装

適当に struct を用意して使いやすくします。

play.rust-lang.org

fn main() {
    let n = 20;
    let ls = LinearSieve::new(n);
    for i in 2..=n {
        eprintln!("{i} -> {:?}", ls.factors(i));
    }
}

pub struct LinearSieve {
    lpf: Vec<usize>,
    lpf_ord: Vec<u32>,
    lpf_pow: Vec<usize>,
}

impl LinearSieve {
    pub fn new(n: usize) -> Self {
        let mut lpf = vec![1; n + 1];
        let mut primes = vec![];
        for i in 2..=n {
            if lpf[i] == 1 {
                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;
            }
        }
        
        let mut lpf_ord = vec![1; n + 1];
        let mut lpf_pow = vec![1; n + 1];
        for i in 2..=n {
            let j = i / lpf[i];
            if lpf[i] == lpf[j] {
                lpf_ord[i] = lpf_ord[j] + 1;
                lpf_pow[i] = lpf_pow[j] * lpf[i];
            } else {
                lpf_ord[i] = 1;
                lpf_pow[i] = lpf[i];
            }
        }
        
        Self { lpf, lpf_ord, lpf_pow }
    }
    
    pub fn factors(&self, mut i: usize) -> Vec<(usize, u32)> {
        let mut res = vec![];
        while i > 1 {
            res.push((self.lpf[i], self.lpf_ord[i]));
            i /= self.lpf_pow[i];
        }
        res
    }
}

各々 C++ や Python など好みの言語で実装してみましょう。

おまけ

Euler の totient 関数

Euler の totientトウシェント 関数というのがあります。 $n$ 以下の正整数のうち $n$ と互いに素なものの個数であり、$\phi(n)$ と書きます。

$n$ の素因数分解が得られていれば $$ \phi{\left(\prod_{p_i\text{: primes}} p_i^{e_i}\right)} = \prod_{p_i\text{: primes}} p_i^{e_i-1}\cdot(p_i-1) $$ から計算できます。線形篩を用いれば、$O(n)$ で一通り計算したり、$\angled{O(n), O(\log(n)/\log(\log(n)))}$ で処理したりできます。

これ以外にも、最小素因数を利用して求められる関数は高速に求められます。$n$ の約数の個数や、$n$ の約数の総和なんかもそうですね。

約数列挙

$n$ の約数の個数を $d(n)$ と書いたりしますが、任意の $\varepsilon\gt 0$ に対して $d(n)\in o(n^{\varepsilon})$ や、 $$ \limsup_{n\to\infty} \frac{\log(d(n))}{\log(n)/\log(\log(n))} = \log(2) $$ が知られています。

en.wikipedia.org

線形篩などで素因数分解を得られていれば、約数列挙は（指数部の辞書順で）$O(d(n))$ time で行えます。

上記の $\limsup$ の式から、最悪ケースの計算量は $$ \begin{aligned} d(n) &\le e^{\log(2)\log(n)/\log(\log(n))} \\ &= 2^{\log(n)/\log(\log(n))} \end{aligned} $$ 程度になりそうです。

$2$ の肩に $\log(\bullet)$ がたくさん乗っている式に対する直感があまりなくて困ります。 pairing heap というデータ構造の decrease-key という操作は、ならし $O(2^{2\sqrt{\log(\log(n))}})$ 時間らしいです。これはなんでしょうか。

おわり

実は意外と知られていない気もするし実はみんな知っているのかもしれないシリーズでした。

線形篩自体や osa_k 法自体も知らない人が多いかもしれません？ みんなにとって愛着の強いアルゴリズムはなんですか？

*1:$\lpf(1)$ の値は使わないので、未定義のままでも問題ないです。各々適当な値を入れましょう。実装の際にわざわざ $-1$ のような値を入れるのではなく、都合のいいように適当に入れていいということです。