Project Eulerをhaskellで練習していく日記: Problem 10

問題

10未満の素数の和は、2 + 3 + 5 + 7 = 17　である。
同様に2000000未満の素数の和を求めよ。

回答

-- Atkinの篩
-- http://en.wikipedia.org/wiki/Sieve_of_Atkin
-- を参考に実装しました

import qualified List as L
import qualified Data.Set as S

limit=2000000::Integer --ふるいにかける上限
slimit = (truncate $ sqrt $ fromIntegral limit)::Integer -- 上限の平方

-- 条件のテンプレート
gcond f g x y stack
  | x>slimit = stack
  | y>slimit = gcond f g (x+1) 1 stack
  | otherwise = if n<=limit && m
                       then gcond f g x (y+1) (n:stack)
                       else gcond f g x (y+1) stack
                       where n = f x y 
                             m = g n x y

-- 上記テンプレートを使って、個別の条件を実装
cond1 = gcond (\x y -> 4*x^2+y^2) (\n x y -> n `mod` 12==1 || n `mod` 12==5) 1 1 []
cond2 = gcond (\x y -> 3*x^2+y^2) (\n x y -> n `mod` 12==7) 1 1 []
cond3 = gcond (\x y -> 3*x^2-y^2) (\n x y -> x>y && n `mod` 12==11) 1 1 []

-- 上記３条件で奇数個の解を持つものだけを取り出し、候補とする
cond = map head $ 
       filter (\x->length x `mod` 2==1) $ 
       L.group $ L.sort (cond1++cond2++cond3)

-- 素数の平方の倍数の一覧を作る
elim = S.fromList $ concat [[x^2,2*x^2..limit] |x<-cond, x<=slimit]

-- 候補から平方倍数一覧を除外して5以上の素数一覧を作成
prims = S.difference (S.fromList cond) elim

main = do  --以下の5+というのは最初の二つの素数2,3の和の部分
     print $ 5+(foldl1 (+) $ S.toList prims)

感想

エラストテネスの篩は遅いだろうと想像がついたので、wikipediaを参考に、Atkinの篩というのを実装してみた。
まあまあ速いのではないだろうか。

おもいっきり手続き型のシュードコードを参考にしたので、あまりhaskell的ではないかもしれません。

初めてData.Setというものを使ってみた。便利。（集合の差分を計算するため）