エラトステネスの篩に基づく高速な素因数分解

概要

• 前処理をした後 $n$ を $O(\log n)$ 時間で素因数分解するよ
• 区間 $[L, R)$ の整数たちを素因数分解するよ
• $R \le 10^{12}$, $R-L \le 10^6$ くらいだよ

導入

Eratosthenes の篩を知っていますか？
$N$ が与えられたとき，$i$ を除く $i$ の倍数にマークをつけていく ($2 \le i \le \sqrt{N}$) ことで，$n \in[2, N)$ に対する素数判定を高速に行う（ための前処理をする）アルゴリズムです．

ところで，このマークをつけていく処理において，true/false よりも多くの情報があるにも関わらずそれを捨てていますよね¹．具体的には，「その整数を割ることができる整数は何か？」という情報です．

実装

素数かどうかを格納した配列ではなく，その整数の持つ最小の素因数を格納した配列を得ることを考えてみましょう．

sieve.cpp

std::vector<int> sieve(int n) {
  std::vector<int> res(n);
  std::iota(res.begin(), res.end(), 0);
  for (int i = 2; i*i < n; ++i) {
    if (res[i] < i) continue;
    for (int j = i*i; j < n; j += i)
      if (res[j] == j) res[j] = i;
  }
  return res;
}

倍数を見ていく部分に関して，$j$ を $i$ で割った結果が $i$ 未満になるなら，その整数でも $j$ を割り切れることになるので，$j \ge i^2$ としてしまってよいです．
ただし，$j \ge i^2$ なら必ずしも更新していいわけではないことに注意してください (e.g. $j = 12$)．

これを用いて以下のように素因数分解を行うことができます．

sieve.cpp

std::vector<int> factor(int n, const std::vector<int>& min_factor) {
  // min_factor は sieve() で得られたものとする
  std::vector<int> res;
  while (n > 1) {
    res.push_back(min_factor[n]);
    n /= min_factor[n];
    // 割った後の値についても素因数を知っているので順次求まる
  }
  return res;
}

もちろん，素数判定もこの篩でできます（sieve[i] < i なら合成数）．

茶番

えびちゃん「このアルゴリズムは先日思いついたばかりでな．なにしろまだ名前すらついておらぬ—————...」
？？？「...いや... 名なら在る」
えびちゃん「...何だと？」
？？？「『osa_k 法』と言う」
えびちゃん「...何を...言っている...？」²

拡張

区間篩と呼ばれる手法がありますね．
ある区間 $[L, R)$ に対する Eratosthenes の篩が欲しいものの $[2, L)$ の部分は不要なときに使えるものです．
$[2, \sqrt{R})$ に関する通常の篩を作りながら構成します．
以下の説明では，$[2, \sqrt{R})$ を小さい篩，$[L, R)$ を大きい篩と呼びます．

sieve.cpp

std::vector<bool> sieve(intmax_t L, intmax_t R) {
  int n = sqrt(R)+1;  // 念のため多めに
  std::vector<bool> aux(n, true);  // 小さい篩
  std::vector<bool> res(R-L, true);  // 大きい篩
  for (intmax_t i = 2; i*i < n; ++i) {
    if (!aux[i]) continue;
    for (intmax_t j = i*i; j < R; j += i)
      aux[j] = false;
    for (intmax_t j = std::max(i, (L+i-1)/i)*i; j < R; j += i)
      res[j] = false;
  }
  return res;
}

これについても素因数を求めながら構築することが可能です．ただし，$k \in [L, R)$ を最小の素因数で割ったものは，必ずしもどちらかの篩がカバーする区間に入っているとは限りませんので，注意が必要です．
たとえば，$(L, R) = (100, 120)$ としたとき，$100$ を $2$ で割って得られた $50$ について $\sqrt{120} \simeq 11 < 50 < 100$ が成り立ちます．

そこで，大きい篩に属する整数については十分多くの素因数を求めておくことにしてみます．
また，いちいち配列を返しているとごちゃごちゃしてくるので，クラスを定義することにします．

sieve.cpp

class smart_sieve {
  intmax_t L, R, M;
  std::vector<int> small;  // 小さい篩
  std::vector<std::vector<intmax_t>> large;  // 大きい篩
  std::vector<intmax_t> aux;  // aux[i] := large[i] の素因数の積

public:
  smart_sieve(intmax_t L, intmax_t R): L(L), R(R), M(sqrt(R)+1) {
    small.resize(M);
    std::iota(small.begin(), small.end(), 0);
    large.resize(R-L);
    aux.assign(R-L, 1);

    for (intmax_t i = 2; i*i < R; ++i) {
      if (small[i] < i) continue;
      small[i] = i;
      for (intmax_t j = i*i; j < M; j += i)
        if (small[j] == j) small[j] = i;

      for (intmax_t j = (L+i-1)/i*i; j < R; j += i) {
        intmax_t k = j;
        do {
          // aux[j-L] > M で判定した方がいいかも？
          // j / aux[j-L] < M の方がいい？（割り算したくない）
          if (aux[j-L] * aux[j-L] > R) break;

          large[j-L].push_back(i);
          aux[j-L] *= i;
          k /= i;
        } while (k % i == 0);
      }
    }
  }

  std::vector<intmax_t> factor(intmax_t n) {
    assert(L <= n && n < R);
    std::vector<intmax_t> res = large[n-L];
    n /= aux[n-L];
    if (n >= M) {
      // この場合，n は素数となることが示せる（はず）
      // n*n >= R だとオーバーフローしそう？
      res.push_back(n);
      return res;
    }
    while (n > 1) {
      res.push_back(small[n]);
      n /= small[n];
    }
    return res;
  }
};

int main() {
  intmax_t L, R;
  scanf("%jd %jd", &L, &R);

  smart_sieve ss(L, R);
  for (intmax_t i = L; i < R; ++i) {
    auto f = ss.factor(i);
    printf("%jd:", i);
    for (auto j: f) printf(" %jd", j);
    printf("\n");
  }
}

計測

しょーぶ！

対戦相手は GNU factor くんです．

$ factor --version
factor (GNU coreutils) 8.29
Copyright (C) 2017 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <https://gnu.org/licenses/gpl.html>.
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.

Written by Paul Rubin, Torbjörn Granlund, and Niels Möller.

$[10^{12}, 10^{12}+10^6)$ を素因数分解するのに掛かる時間を計測します．
factor くん的には閉区間がお望みのようなので，そうします．

$ time echo 100000{0,1}000000 | ./sieve > large.out1

real    0m1.802s
user    0m1.726s
sys     0m0.055s
(*'-')b < Exited successfully

$ time seq 1000000{000000,999999} | factor > large.out2

real    0m4.138s
user    0m3.721s
sys     0m0.440s
(*'-')b < Exited successfully

$ diff large.out{1,2}
(*'-')b < Exited successfully

勝ちました．

区間幅が大きすぎると処理しきれないので別の策を考えましょうね．
クエリが単調増加とかなら，古い篩（激ウマギャグ）を捨てつつ処理することができるかもしれません．

あと，上記コードにおいて main() に出てこない関数たちは特にテストしてないのでちゃんと動かないかもしれません，許してね．

茶番

えびちゃん「これは osa_k 法と区間篩を練り合わせたアルゴリズムでな，私が創り上げたものだ．このアルゴリズムは先日（以下略）」²

1. 蟻本 p. 63「一般に bool 値を求める DP をすることは無駄があることが多く、同じ計算量でもっと多くのことを知ることができます。」 ↩

2. BLEACH 18 巻 pp. 197–198 ↩ ↩²