[C#] エラトステネスの篩で素数列挙と素因数分解を実装する

エラトステネスの篩とは

• エラトステネスの篩 – Wikipedia

エラトステネスの篩 (エラトステネスのふるい、英: Sieve of Eratosthenes) は、指定された整数以下の全ての素数を発見するための単純なアルゴリズムである。古代ギリシアの科学者、エラトステネスが考案したとされるため、この名がある。

指定された整数N以下の素数を列挙できるアルゴリズム、エラトステネスの篩をC#で実装してみます。

アルゴリズム

エラトステネスの篩では、以下の手順で指定された整数N以下の素数をすべて発見します。

1. 2からNまですべての数を素数候補のリストに昇順で入れる。
2. 素数候補のリストの先頭xは素数である。
3. 手順2 で見つけた素数の倍数を素数候補のリストからふるい落とす。
4. 素数候補のリストの先頭xがNの平方根以下であれば、手順2 に戻る。
5. 素数候補のリストに残った数はすべて素数である。

具体例でみてみます。N=10の時は以下のような手順で素数が発見できます。

素数候補のリストを S 発見した素数のリストをPとします。

1. 素数候補のリストSと素数リストPを作成する
   S=[2,3,4,5,6,7,8,9,10], P=[]
2. Sの先頭の 2 は素数である。Pに2を追加して、Sから2の倍数をふるい落とす。
   S=[3,5,7,9], P=[2]
3. Sの先頭の 3 は素数である。Pに3を追加して、Sから3の倍数をふるい落とす。
   S=[5,7], P=[2,3]
4. √10 < 5 なのでふるい落とす手順は終了する。
5. Sに残った数はすべて素数である。よって10以下の素数は以下のとおりである。
   P=[2,3,5,7]

√N より大きい数は確認しなくてもふるい落としの段階ですべての合成数が除かれています。

計算量

• エラトステネスのふるいとその計算量 | 高校数学の美しい物語

2でふるい落とすのに N/2回、3でふるい落とすのに N/3回、以下N/5 … と計算が続きます。

N/2 + N/3 + N/5 + ... 

Nが大きいとき、この計算量は O(NloglogN) くらいになるみたいです。

詳しい計算過程は上のページにありました。

ある数Nを素数かどうか判定するのに試し割を行う実装だと O(√N) の計算量がかかります。N以下のすべての数について素数判定のための試し割を行う場合は O(N√N) かかるので、エラトステネスの篩 O(NloglogN) が非常に高速な素数判定であることがわかります。

実装C#

ではアルゴリズムの紹介が済んだところでC#でエラトステネスの篩を実装してみます。

static List<int> SieveOfEratosthenes(int n)
{
    var p = new List<int>();
    var s = new List<int>();
    if (n < 2) return p;

    // 素数候補を列挙する[2..n]
    for (int i = 2; i <= n; i++) { s.Add(i); };

    // √n以下の素数候補についてふるい落としを行う
    while (s[0] * s[0] <= n)
    {
        var prime = s[0];    // 先頭要素は素数
        p.Add(prime);        // 先頭を素数リストへ
        s.RemoveAt(0);       // 先頭要素削除

        // ふるい落とし
        s.RemoveAll(x => x % prime == 0);
    }
    // 残りの探索リストを素数リストに追加
    p.AddRange(s);
    return p;
}

アルゴリズムの手順を愚直に実装したのが上のコードです。引数n以下の素数リストを返す関数です。

No.713 素数の和 – yukicoder

動作確認に上記問題を使って無事ACできました。

これでも十分動くのですが、リストから要素を削除したりしているのが微妙なのでもう少しスマートにし、できることも増えた実装を紹介します。

素因数分解もできるエラトステネスの篩

エラトステネスの篩で合成数をふるい落とすときに使った素数はその合成数の最小の素因数になっています。

つまりこの数を記録しておくことで、N以下の整数について素因数分解が高速にできるようになります。素因数分解まで対応したエラトステネスの篩の実装は以下のようになります。

class SieveOfEratosthenes
{
    // 最小の素因数を記録した表
    public int[] Table;
    public SieveOfEratosthenes(int max)
    {
        CreateTable(max);
    }
    // n以下の数について、最小素因数を記録した表を生成する
    private void CreateTable(int n)
    {
        // 最初すべての数の最小の素因数はその数自身
        Table = Enumerable.Range(0, n + 1).ToArray();
        if (n <= 3) return;
        // 2の倍数の最小素因数は2
        for (int i = 2; i <= n; i += 2) Table[i] = 2;
        // 3以上の奇数について最小素因数を記録する
        for (int p = 3; p * p <= n; p += 2)
        {
            // 奇数pがp未満の数の倍数である場合は無視
            if (Table[p] < p) continue;
            // 素数pの倍数を記録
            for (int x = p; x <= n; x += p)
            {
                Table[x] = p;
            }
        }
    }
    // 素数判定
    public bool IsPrime(int n)
    {
        return Table[n] == n;
    }
    // 最小素因数表を使って素因数分解を行う
    public List<int> PrimeFactors(int n)
    {
        var res = new List<int>();
        while (Table[n] != 1)
        {
            res.Add(Table[n]);
            n /= Table[n];
        }
        return res;
    }
}

実装の基本はエラトステネスの篩です。

N以下の数の最小の素因数を求める

まずN以下の数について最小の素因数を記録するための表(配列)を作成しておきます。初期状態としてその数自身を最小の素因数としておきます。

次に2の倍数についてはすべて最小の素因数を2にしておきます。

√N以下の奇数pについてはループしてすでに表にp以下の数が素因数として記録されていれば(つまり素数でなければ)無視し、そうでなければ素数なので、倍数すべてに最小の素因数としてpを記録するようにします。この処理がふるい落としの処理に該当します。

private void CreateTable(int n)
{
    // 最初すべての数の最小の素因数はその数自身
    Table = Enumerable.Range(0, n + 1).ToArray();
    if (n <= 3) return;
    // 2の倍数の最小素因数は2
    for (int i = 2; i <= n; i += 2) Table[i] = 2;
    // 3以上の奇数について最小素因数を記録する
    for (int p = 3; p * p <= n; p += 2)
    {
        // 奇数pがp未満の数の倍数である場合は無視
        if (Table[p] < p) continue;
        // 素数pの倍数を記録(ふるい落としの処理)
        for (int x = p; x <= n; x += p)
        {
            Table[x] = p;
        }
    }
}

こうしてすべてのふるい落としが終わった後、最小素因数表には素数の場合その数自身が記録され、それ以外の数については最小の素因数が記録されています。

したがって、素数を判定するには以下のように最小の素因数が自分自身であるかどうかを確認すればよいです。

// 素数判定
public bool IsPrime(int n)
{
    return Table[n] == n;
}

すべての素数が欲しければNまでループして素数の数だけを取り出せばよいです。

素因数分解の実装

最小素因数表を使って素因数分解する関数の実装です。

// 最小素因数表を使って素因数分解を行う
public List<int> PrimeFactors(int n)
{
    var res = new List<int>();
    while (Table[n] != 1)
    {
        res.Add(Table[n]);
        n /= Table[n];
    }
    return res;
}

基本的には現在の数を最小の素因数で割ることを繰り返します。最小の素因数が自分自身になれば終了します。例えば60を素因数分解すると以下のような手順で素因数が得られます。

Table[60] = 2, res=[2]
Table[30] = 2, res=[2,2]
Table[15] = 3, res=[2,2,3]
Table[5] = 5, res=[2,2,3,5]

5の素因数が自分自身なのでここで処理が終了します。

計算量

素因数分解の計算量は、2のべき乗の時に最大になるので O(logN) になります。つまりエラトステネスの篩の前計算 O(NloglogN) を実行しておくことでN以下のすべての数の素因数分解を O(logN) で実行できるようになります。

ある数Nの素因数分解も試し割の要領で O(√N) の計算量がかかりますが、前計算が実行できるならエラトステネスの篩が高速に素因数分解できることがわかります。

ただし前計算が O(NloglogN) がかかるので、実行環境によりますが大きすぎる数について素因数分解したい場合は使えません。その場合は √N 以下の数について割り算を繰り返す方法で素因数分解しましょう。

動作確認

このエラトステネスの篩で素因数分解をできるようにする実装は AtCoder Beginner Contest 152 解説放送 の内容を参考にしました。

E – Flatten

E問題が素因数分解を使った解法の問題です。

10⁶ の素因数分解を 10⁴ 回繰り返す必要があるのですが、エラトステネスの篩を使えば前計算しておくことで高速に素因数分解が実行できました。

ただこの問題では、10⁴回愚直に素因数分解しても間に合います。C#の実装だとエラトステネスの篩で70ミリ秒前後、愚直な素因数分解だと350ミリ秒前後くらいでした。

以上。