１．エラトステネスの篩
素数を求める方法で、よく知られているのはエラストテネスの篩（ふるい）である。原理は次のようになっている。
１．2から始まる自然数のリストを用意する。
２．また、素数が求まっていたとする。
３．素数のリストは、自然数のリストからそれぞれの素数の倍数を除いたものである（計算時間を短くするためには、素数をpとした時、pよりも小さい倍数、即ち、p*2, p*3, ..,p * (p-1)まではpよりも小さい素数の倍数を処理するときに既に除いているので、p*p以上の倍数について除けば良い）。
Haskellではこれを素直にプログラムに書けばよい。
２ プログラム化する。
素数をpとするとその倍数の集合はHaskellでは次のように表すことができる。
[p*p, p*p + p,..]
これを全ての素数について行うと「素数の倍数の集合」の集合となる。素数の集合をprimesとすると、Haskellでは
[[ p*p, p*p + p..]|p<-primes]
となる。これは、集合の集合になっているので、集合の中にあるすべての集合の和を取ってあげると一つの集合になる。和集合の演算をunionとすると、foldrを用いて次のように表すことができる。
foldr (\ p r -> union [p*p, p*p + p..] r) [] primes
上記のプログラムで、foldrは、primesの個々の要素pから素数の倍数 [p*p, p*p + p..]を作成し、これ以降に出来上がる素数の倍数のリストrとの和集合を取る。全ての素数の処理が終わった時、最後に（foldrの引数として定義されている）空集合[ ]との和集合を取る。
一見、上記のプログラムで良さそうなのだが、そうはいかない。素数は、2,3,5…であるが、上記のプログラムだと2の倍数がすべて終わってから、3の倍数となる。2の倍数は無限にあるので、3以降の素数の倍数を求めるところには永遠に回ってこない。そこで、素数の倍数が昇順に並ぶようにプログラムを書き換えると次のようになる。
foldr (\p r-> p*p : union [p*p+p, p*p+2*p..] r) [] primes union (x:xs) (y:ys) = case (compare x y) of LT -> x : union xs (y:ys) EQ -> x : union xs ys GT -> y : union (x:xs) ys union xs [] = xs union [] ys = ys
上記のプログラムで、unionは和集合であるが、Data.Listで定義されている和集合ではなく、昇順に並べられた二つのリストから、和集合の昇順のリストを作成する。また、foldrのところで、p*p : unionとなっているのはp*pの値がunionで得られるリストのどの要素よりも小さいためである。
素数は、自然数のリストから、素数の倍の数のリストを除いたもの、即ち、自然数のリストから素数の倍の数のリストの差集合である。差集合をminusとし、この演算で得られる集合も昇順になるようにするとプログラムは次のようになる。
primes :: [Integer] primes = 2 : minus [3..] (foldr (\p r-> p*p : union [p*p+p, p*p+2*p..] r) [] primes) where minus (x:xs) (y:ys) = case (compare x y) of LT -> x : minus xs (y:ys) EQ -> minus xs ys GT -> minus (x:xs) ys minus xs _ = xs union (x:xs) (y:ys) = case (compare x y) of LT -> x : union xs (y:ys) EQ -> x : union xs ys GT -> y : union (x:xs) ys union xs [] = xs union [] ys = ys
プログラムを実行してみる。
take 20 primes
と入力すれば
*Main> take 20 primes [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71] *Main>
が得られる。
素数を求める方法は沢山あるので、計算時間の違いなどを調べながら、他の方法を試してみると面白い。
