[Haskell] 素数判定(試し割り,エラトステネスの篩)

Haskellのお勉強

絶賛Haskell強化週間中です。これまでハノイの塔、ユークリッドの互除法、それぞれのアルゴリズムをHaskellで実装しましたが、今回は素数判定をしてみようと思います。

素数判定にはさまざまな判定法があるらしいのですが、今回は単純な試し割りによる素数判定とエラトステネスの篩をHaskellで書いてみます。 いずれも
Wikipediaに記事が有りますのでそれを参考にしていきます。

試し割りによる素数判定

まずは試し割りによる素数判定です。これはおそらく一番単純な素数判定方法で、約数候補の数で割り切れるまでばしばし割っていくというものです。すべての候補値で割り切れなかった時、晴れて素数と判定されます。

Haskellのプログラム

これは特に難しいこともないのでコードを載せます。

isPrime :: Int -> Bool
isPrime 2 = True
isPrime n
    | n < 2 = False
    | n `mod` 2 == 0 = False
    | otherwise = try n 3
    where
        try :: Int -> Int -> Bool
        try m i
            | m `mod` i == 0 = False
            | i^2 > m = True
            | otherwise =  try m (i+2)

コードの解説ですが、isPrime関数は判定を行う整数を引数に取り、判定結果をBoolで返します。 処理の流れは上から順にだいたいは以下のとおりです。

• 2は素数
• 2未満の数(1とか0)は素数ではない
• 偶数は素数じゃない
• 上記以外の数(2以下の数と偶数以外)について、試し割りを行い確認する(再帰処理)
• 試し割りはnの平方根まででよい(コード上では2乗した値と比較)

なぜnの平方根まででよいか

試し割りはnの平方根以下の値のみの確認でよいのですが、これはなぜでしょうか。調べてみるといろいろ数学の証明が載っているのですが、次のように考えればもっと直感的に理解できるはずです。

仮にnが合成数(素数ではない数)の場合、幾つかの約数を持ちます。n=16の場合、[1,2,4,8,16]が約数となりますが、約数にはかけ合わせると元の数(n)になる組み合わせが存在します。 この場合だと、(2,8)がその組み合わせです。つまり元の数の平方数(この場合4)以下で割り切れない場合、それ以上の値を調べても対になる片割れがなかった以上、もう約数は出てこないため、素数と判定できるのです。

C#のプログラム

ちなみにWikipadiaにはC言語の例がありますが、例によってC#で書き直してみました。冪乗の^(ハット)演算子が使えないことを初めて知りました …

private static bool IsPrime(int n)
{
    if (n < 2)
    {
        return false;
    }
    else if (n == 2)
    {
        return true;
    }
    else if (n % 2 == 0)
    {
        return false;
    }
    for (int i = 3; Math.Pow(i, 2) <= n; i++)
    {
        if (n % i == 0)
        {
            return false;
        }
    }
    return true;
}

エラトステネスの篩

素数判定の方法の一種に
エラトステネスの篩(ふるい) というものがあります。 指定された整数以下の全ての素数を発見するための単純なアルゴリズムである … らしいです。 以下Wikipediaからアルゴリズムの引用です。

アルゴリズム

1. 探索リストに2からxまでの整数を昇順で入れる。
2. 探索リストの先頭の数を素数リストに移動し、その倍数を探索リストから篩い落とす。
3. 上記の篩い落とし操作を探索リストの先頭値がxの平方根に達するまで行う。
4. 探索リストに残った数を素数リストに移動して処理終了。

x = 30 の時の処理の流れ

実際にエラトステネスの篩で30以下の素数を求める際の処理について、流れを追ってみます。

1. 探索リスト=[2..30]
   探索リストの先頭値=2
2. 素数リスト=[2]
   探索リスト[5,7..29]
   探索リストの先頭値=3
3. 素数リスト=[2,3]
   探索リスト[5,7,11,13,17,19,23,25,29]
   探索リストの先頭値=5
4. 素数リスト=[2,3,5]
   探索リスト[7,11,13,17,19,23,29]
   探索リストの先頭値=7
5. √30(=5.4472..) < 7 なのでふるい落とし終了
6. 残った探索リストを素数リストに移動
   素数リスト=[2,3,5,7,11,13,17,19,23,29]

Haskellのプログラム

なんとなくイメージを掴んだのでこれをHaskellのコードで書いてみます。

main = print $ sieveOfEratosthenes 30
-- [2,3,5,7,11,13,17,19,23,29]

sieveOfEratosthenes :: Int -> [Int]
sieveOfEratosthenes n
    | n < 2 = []
    | n == 2 = [2]
    | otherwise = sieve [x | x <- [2..n]] n
    where
        sieve :: [Int] -> Int -> [Int]
        sieve [] _ = []
        sieve (x:xs) n
            | x^2 > n = x:xs
            | otherwise = x : (sieve [a | a <- xs, a `mod` x > 0] n)

sieveOfEratosthenes関数は整数nを引数に取り、その値以下の素数のリストを返します。nが2以下の場合は個別に処理します。n>2の時にふるい(sieve)にかけてやります。

ふるいにかける探索リストと元の数nを引数に取るsieve関数を定義し、これを再帰的に呼び出すことで素数のリストを返します。

最終行のコードが探索リスト(引数のリスト)の先頭要素を素数として取り出し、残りの探索リストをふるいにかけている箇所です。
リスト内包表記を使って、先頭の値で割った余りが0以上、つまり先頭の値の倍数以外のリストを探索リストとして次のふるい落としを呼び出しています。

下から2行目、試し割りのときと同様に2乗値との比較で終了判定を行い、終了条件を満たす場合には残りの探索リストを素数リストとして返しています。
こうすることですでに探索リストの先頭から取り出された素数の後ろに連結され、結果として素数リストが得られます。

確認

wikipediaに1000以下の素数が載っていたのでそれと確認したところ一致していたので、問題無いでしょう。
sieve [] _ = [] にマッチすることはありえないのですが、x:xsには空のリストが渡された時にはマッチしないので一応入れています。

[x]の場合、x:xsには、x:[]としてマッチします。

C#でも

C#でも書いてみます。せっかくなので再帰関数を使わずに。

private static List<int> SieveOfEratosthenes(int n)
{
    var result = new List<int>();
    var list = new List<int>();
    if (n < 2)
    {
        return result;
    }
    else if (n == 2)
    {
        result.Add(2);
        return result;
    }
    for(int i = 2; i <= n; i++) { list.Add(i); };
    while (Math.Pow(list[0], 2) <= n)
    {
        var first = list[0];    // 先頭要素
        result.Add(first);      // 先頭を素数リストへ
        list.Remove(first);     // 先頭要素削除

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

やはり再帰関数を使って書いたほうがわかりやすと思います。