はじめに

素数って難しいですよね。今回は(といっても次回があるかはわからないのだが)そんな素数くんと遊んでみた話です。

素因数分解と計算量

任意の自然数が素因数の積として一意にかけるという事実は有名ですが、じゃあ実際にコンピュータで大きな自然数$N$を素因数分解してみようとなると結構大変です。RSA暗号なんかはそれを安全性の根拠として使われているものです。僕は素因数分解アルゴリズムについては全く詳しくないんですが、現在多項式時間アルゴリズムは発見されていないようです。なおここでの多項式時間であるとは$N$の桁数を$d$としたときにオーダーが$d$の多項式で書けることを言います。

エラトステネスのふるい

エラトステネスのふるいとはある自然数$N$以下の素数を列挙する単純なアルゴリズムです。エラトステネスってはるか昔のギリシアにいた賢い人で他にも地球の大きさを計算した話なども有名ですよね。

アルゴリズム

Step 1
サイズ$N+1$の配列$A$を用意し、$0$に初期化。

Step 2
$A[2]$から始めて、値が$0$のまだ見ていないインデックスの中で先頭のものを$i$とし、$A[k \times i] (k=2 ... \lfloor \frac{N}{i} \rfloor)$の値が$0$なら$i$を代入。

Step 3
上記の操作を先頭要素が$\sqrt{N}$に達するまで行う。

Step 4
$A[2]$から初めて、値が0のインデックスをすべて出力して終了。

計算量

アルゴリズム自体は単純ですが、計算量は結構非自明です。素数を昇順に$p_i (i=0, 1 \cdots)$と書くことにすると、$p_i$の倍数をふるい落とす回数がざっくり$\lfloor \frac{N}{p_i} \rfloor$回くらいなので、全体で$N \sum_{p_i\leq\sqrt{N}} \frac{1}{p_i}$で抑えられます。これが$\mathcal{O}(N \log\log N)$らしいです。

素因数分解への応用

実はエラトステネスのふるいは単純に$N$以下の素数の列挙だけでなく、素因数分解も行えます。元々は値が$0$となるインデックスをすべて出力していましたが、インデックスが合成数のときに入っている値は実はその合成数の最小の素因数になっています。よってインデックスをその値で割り切ることができ、割った値を新たなインデックスとして再帰的にアクセスしていけばすべての素因数を列挙できます。例えば$12$だと、$12 \rightarrow 6 \rightarrow 3$の順にアクセスし、最後に素数$3$にたどり着いて終了します。

試し割り

上の議論で最小の素因数さえ求めれば$N$以下の全ての自然数の素因数分解ができることがわかりました。そこで試し割りによって最小の素因数を見つけるようにしたら計算量がどのようになるのか、つまりエラトステネスのふるいと比べて遅いか速いかということが僕の関心事です。ここでの試し割りとは自然数$n$を素数で昇順に割っていって割り切れるものが見つかるか、素数が$\sqrt{n}$を超えたら終了というものです。

計算量解析(偽)

$N$が十分大きいとすると、試し割りが$i+1$回で終わる確率は、大体$\frac{1}{\prod_{j=0}^{i+1} p_j}$になります。なので全体の計算回数は$N \sum_{p_i \leq N} \frac{i+1}{\prod_{j=0}^{i}p_j}$くらいになりそうですよね。しかし総和は等比級数で上から抑えられるので定数になります。よって試し割りの計算量は$\mathcal{O}(n)$だとわかりました。ドヤァ

というのは嘘なんですが、どこが嘘かわかりましたか？

何がおかしかったのか

上のどこが間違っているのか、最初は僕も気付かなかったのですが、結論からいうと、$N$の極限の飛ばし具合を$p_i$によって変えていたことが問題でした。確かに各$p_i$について確率が$\frac{1}{\prod_{j=0}^{i}p_j}$に収束することは正しいのですが、$p_i$が大きくなるほど収束させるのに必要な$N$の大きさが爆発していきます(素数の総積なので階乗よりさらに大きい)。別の言い方をすると$N$をどれだけ大きくしても、出現確率が$\frac{1}{\prod_{j=0}^{i}p_j}$に全然収束していない$p_i$が存在してしまいます。ということで上の議論は間違いです。

論より証拠

結論、よくわからん。
エラトステネスのふるいとの優劣も微妙なんですよね。
例えば素数に関してもエラトステネスのふるいは1回しか調べないので、効率的なんですが合成数については素因数の種類が多いとその種類数だけ調べられるので効率が悪いです。一方、試し割りは素数に関しては$\sqrt{p_i}$まで割らないと素数ということがわからないので効率が悪いんですが、合成数に関しては最小の素因数までしか調べないのでエラトステネスのふるいより効率的です。となると素数定理とかを使って素数の出現頻度とかも考えないといけないのかなあとか考えましたが、いまいちよくわかりませんでした。

てことで、実験しましょう。

エラトステネスのふるいのコード

int A[N];
vector<int> prime;

void era(int n)
{
    for (int i = 2; i * i < n; ++i)
    {
        if (A[i] == 0)
        {
            for (int j = i; i * j < n; ++j)
            {
                if (A[i * j] == 0) A[i * j] = i;
            }
        }
    }

    for (int i = 2; i < n; ++i)
    {
        if (A[i] == 0) prime.push_back(i);
    }
    return;
}

試し割りのコード

int A[N];
vector<int> prime;

void trial(int n)
{
    for (int i = 2; i < n; ++i)
    {
        if (prime.empty()) prime.push_back(i);

        else
        {
            bool isprime = true;
            for (auto j = prime.begin(); j != prime.end() and (*j) * (*j) <= i; ++j)
            {
                int p = *j;
                if (i % p == 0)
                {
                    A[i] = p;
                    isprime = false;
                    break;
                }
            }
            if (isprime) prime.push_back(i);
        }
    }

    return;
}

結果

$N$を変えて実行時間を計測した結果が以下になります。

$N$	エラトステネスのふるい	試し割り
$10^5$	0.02sec	0.02sec
$10^6$	0.02sec	0.06sec
$10^7$	0.14sec	0.79sec
$10^8$	1.36sec	15.7sec

顕著にエラトステネスのふるいの方が高速でした。計算量気になるなあ。わかる方いたら教えてください。

まとめ

エラトステネスさん、強い。早くFGOに実装されねーかな。ギリシアは互除法といい、シンプルかつ高速なアルゴリズムが作られていた印象があります。

おまけ

計算量解析に素数定理とかが使われているらしくて興味深かったです。素数の出現頻度が$\frac{1}{\log N}$で近似できるとかね。
ふるいを$\mathcal{O}(N)$で実行する方法もありますが、それについてはまた今度ということで。

素因数分解小話