Haskellで有名アルゴリズム実装（4）

第4弾。今日は、素数列挙のアルゴリズム。

ある数が素数かどうかの判定は、擬多項式時間でできます。ある区間[1, n]内の素数を列挙するには、普通にやるとO(n*sqrt(n))かかります。

エラトステネスの篩というアルゴリズムを使用すると、O(n)で計算できます。

では、実装。

Haskellのエラトステネスは美しいです。

view plainprint?

1. eratos []     = []  
2. eratos (x:xs) = x : eratos [y | y <- xs, mod y x /= 0]  

eratos []     = []
eratos (x:xs) = x : eratos [y | y <- xs, mod y x /= 0]

これで十分な場合がほとんどですが、100万までの素数列挙とかすると、上の実装ではたえられません。エラトステネスが速いのは、ある素数pの倍数を2*p, 3*p,..とpおきに消していけるからで、上の方法だとすべての要素をなめているので遅いです。加えてリストは参照にO(n)かかるので遅いです。

ある素数に対して必要な要素だけ参照する+高速に参照することができるデータ構造に変換する。というのが必要です。100万くらいの素数列挙なら、Setを使えば十分でしょう。あとは、thunkがメモリを圧迫するのを防ぐ+遅延評価によるスタックオーバーを防ぐため正格評価を用います。

view plainprint?

1. primes :: [Integer] -> [Integer]  
2. primes = toList . (seive 2) . fromList  
3.   
4. seive :: Integer -> Set Integer -> Set Integer  
5. seive p xs  
6.      | p*p > findMax xs = xs  
7.      | otherwise        = seive (p+1) $ foldl' (\acc x -> delete x acc) xs [2*p,3*p..q]  
8.                           where q = findMax xs  

primes :: [Integer] -> [Integer]
primes = toList . (seive 2) . fromList

seive :: Integer -> Set Integer -> Set Integer
seive p xs
     | p*p > findMax xs = xs
     | otherwise        = seive (p+1) $ foldl' (\acc x -> delete x acc) xs [2*p,3*p..q]
                          where q = findMax xs

ただ、上の実装でも本来のO(n)の速度には届かないので、より大きな範囲（[1, 10000000]など）の素数列挙をしたい場合は、工夫が必要です。Setは参照がO(log n)なので、ランダムアクセスが可能なデータ構造を利用する + なるべく遅延評価の機構をいれないようにするで達成できるかと思います。

Haskellのページに戻る。