« Project Euler : Problem 6 | トップページ | Project Euler : Problem 7 ~ 10001 番目の素数 »

2009年12月25日 (金)

Haskell で素数

 "Project Euler" の問題には「素数」に関連したものが時々出てきます。そこで今回は、Haskell でどうやって素数列を扱っていけば良いかを考えてみます。

 Haskell で有名なものに「エラトステネスの篩」があります。(以下のコードは最初から 2 以外の偶数を省いたものです)

primes :: [Integer] primes = 2 : sieve [3, 5 ..] where sieve (p : xs) = p : sieve [x | x <- xs, rem x p /= 0]
 「簡単なコードで素数の無限数列を扱える」ということで Haskell のコードとしては有名(?)なものですが、実は非常に遅いです。
 このコードは、"ghc" で最適化のオプションを付けてコンパイルしても、「20万以下の素数の和」を求めるのに約 40 秒もかかってしまいます。(これ以降、時間に関する記述は、最適化オプションを付けてコンパイルしたコードのものです。ちなみに私の使っているノートパソコンのスペックは "Pentium M, 1500MHz" です)
 Project Euler の Problem 10 のように「200万以下の素数の和」を求めるのにはかなりの時間がかかることが予想されます。

 そこで、自分なりに考えたのが次のコードです。やっていることは、単に素数の候補が本当に素数かどうかを判定しているだけです。
 ただし、素数の候補はただの奇数ではなく、"6 * n ± 1" の形をとる奇数だけを集めることによって、ちょっとだけスピードアップを狙っています。

primes :: [Integer] primes = [2, 3] ++ filter isPrime ns where ns = [x + y | x <- [6, 12 ..], y <- [-1, 1]] isqrt = truncate . sqrt . fromIntegral isPrime n = all ((/= 0) . (rem n)) $ takeWhile (<= isqrt n) primes
 このコードでは、「200万以下の素数の和」が約 4.8 秒で出ます。

 

 自分ではこれ以上いいアイデアを考えつかなかったので、ネット上をいろいろと探してみたところ、こちらのホームページで素晴らしいコードを見つけました。(実際の考え方やアルゴリズムについては、直接こちらのホームページをご覧ください)

primes5 = 2 : primes' where primes' = 3 : sieve 0 5 sieve i x = filter isPrime [x,x+2..p*p-2] ++ sieve (i+1) (p*p+2) where (ps,p:_) = splitAt i primes' isPrime x = all ((/=0).rem x) ps
 たぶん、私の知る限りではもっとも優れたアルゴリズムではないかと思います。実際にコンパイルして「200万以下の素数の和」を求めてみると、約 2.18 秒で答えがでました。

 人のコードをそのまま使うのはちょっと悔しい気がしたので、「このコードを土台にして何とかもっといいコードを書けないか」とかなりの時間考えました。
 まずやったことは、"isPrime" 関数で使われている "ps" を毎回 "primes" から切り出しているのが、少し時間の無駄遣いの用な気がしたので、アルゴリズム自体は変えずに、"sieve" に "ps" に相当する数列を保持させるようにしてみました。

primes :: [Integer] primes = [2, 3] ++ sieve [] (tail primes) 5 where sieve divs (x : xs) n = ps ++ sieve (divs ++ [x]) xs (x * x + 2) where isPrime m = all ((/= 0) . (rem m)) divs ps = filter isPrime [n, n + 2 .. x * x - 2]
 実際にコンパイルしてみると……実行時間はほとんど変わりませんでした。毎回 "divs ++ [x]" としているのが処理に負担をかけているのでしょうか?

 次に、"isPrime" で関数を合成して "all" を使って処理しているところを敢えて再帰関数にしてみました。

primes :: [Integer] primes = [2, 3] ++ sieve [] (tail primes) 5 where sieve divs (x : xs) n = ps ++ sieve (divs ++ [x]) xs (x * x + 2) where isPrime (y : ys) m | rem m y /= 0 = isPrime ys m | otherwise = False isPrime _ _ = True ps = filter (isPrime divs) [n, n + 2 .. x * x - 2]
 すると……なんと「200万以下の素数の和」が約 1.95 秒で出ました。抽象化をするとそれなりに代償がいるということでしょうか?
 もしかしたら、再帰関数の方に末尾再帰の最適化が行われたせいで若干速くなったのでしょうか?

 

 ちなみに冒頭で「遅くて使えない」と言った「エラトステネスの篩」ですが、ちょっと手を加えて「n 以下の素数列を求める」といった使い方に限定するとかなり使えるものに変わります。

primeList :: Integer -> [Integer] primeList n = 2 : sieve [3, 5 .. n] where isqrt = (truncate . sqrt . fromIntegral) n sieve (p : xs) | p > isqrt = p : xs | otherwise = p : sieve [x | x <- xs, rem x p /= 0]
 "main = print $ sum $ primeList 2000000" としてコンパイルすると、約 2.7 秒で答えが出ます。

« Project Euler : Problem 6 | トップページ | Project Euler : Problem 7 ~ 10001 番目の素数 »

Haskell」カテゴリの記事

コメント

おお、確かに速いですね。(特に[Int]とすると爆速です)
Prelude では
all _ [] = True
all p (x:xs) = p x && all p xs
と定義されているので、
isPrime [] _ = True
isPrime (p:ps) x = rem x p /= 0 && isPrime ps x
とすれば同じことだと思うんですが‥。
謎だ。

コメントを書く

コメントは記事投稿者が公開するまで表示されません。

(ウェブ上には掲載しません)

トラックバック

この記事のトラックバックURL:
http://app.cocolog-nifty.com/t/trackback/112020/47109131

この記事へのトラックバック一覧です: Haskell で素数:

» Project Euler : Problem 10 〜 200万以下の素数の和 [ツムジのひとりごと]
 問題はこちらをご覧ください。    問題の解答の説明に入る前に、以前の記事の補足から……。  "Haskell で素数"の記事に rst76 さんが付けてくださったコメントにもあるように、primes の型を [Int] にするとかなり速度が上がることが分かりました。  Pro... [続きを読む]

« Project Euler : Problem 6 | トップページ | Project Euler : Problem 7 ~ 10001 番目の素数 »

2016年7月
          1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
31            
フォト

最近のトラックバック

無料ブログはココログ