99 Haskell Problems より、[31] でもその前に

追記

Qiita の記事の中にすでにもっと早くてエレガントなコードを書いていらっしゃる記事を見つけました。
Haskell で「素数列」をつくる。
なお、さらに元ネタはこちらだそうです。
エラトステネスの小型高速遠心分離機

コメントにも有りましたが、これらの実装も実はエラトステネスのふるいではなく（Haskell の遅延評価のせい？）いわゆる試し割りということになるそうですが、コードゴルフにも採用されてる位短くてかつ早いとなると、、、

一見したところ何をしているか分からなかったのでちょっと紐解いてみます、以下本質的にはHaskell で「素数列」をつくる。からのコピペです、ちょっと型とか二乗とか改変していますが（手元の環境では違いが出ませんでしたが、もしかしたら高次の影響があるかも分かりません、掛け算の方が安価なのかも）

p.lhs

Fast primes from
  http://qiita.com/little_Haskeller/items/614a3ae20a517c19bb1f

> primes :: Integral a => [a]
> primes = map fromIntegral primes'
>   where
>     primes' :: Integral a => [a]
>     primes' = [2,3,5] ++ f 5 7 (drop 2 primes')

>     f :: Integral a => a -> a -> [a] -> [a]
>     f m s (p:ps) = [n | n <- ns, gcd m n == 1]
>                    ++ f (m*p) (p^2) ps
>       where -- Primes are either 6n+1 or 6n+5 form.
>         ns = [x + y | x <- [s, s +6 .. p^2 -2], y <- [0, 4]]
 
> main :: IO ()
> main = print $ take 1 $ dropWhile (< 30000000) primes

  $ ghc p.lhs -O2
  [1 of 1] Compiling Main             ( p.lhs, p.o )
  Linking p ...
  $ time ./p
  [30000001]

  real  0m2.524s
  user  0m2.460s
  sys 0m0.058s

primes' 自体はいわゆる遅延リストになっていて、補助関数f がミソとなるわけです。
またgcd がHaskell では安価だというのも効いているみたいです、順序が逆にはなりますが99 Haskell Problems より、[31..41]でやるように引き算だけで実装できますし。

まず最初には(drop 2 primes') は先頭の[5, しか決まっていません。
処理を追いかけてみますと、f の第二引数はかならず6n+1 型の整数（最初は7 そのあと7^2=49 　）になります。
したがってns の中のx は6n+1 型に、y <- [0,4] との組み合わせでns は6n+1 あるいは6n+5 型の整数になります、この辺は99 Haskell Problems より、[31..41]の31. で扱ったテクニックですね。
p^2 -1 まで考えれば良いのでx の定義域からちょっとオーバーヘッドがあるかもしれませんが、どうでしょうか？

こうして取り出したns をgcd で試し割りをすることで残りカスが素数になります、最初のトライでは[7,11,13,17,19,23] となります、おぉ、25は入ってないのでここではオーバーヘッドないですね。
これでじつはprimes' は[2,3,5,7,11,13,17,19,23 とわかります。
従って(++) の後ろではf 25 25 を[7,11,13,17,19,23, で始まるリストに作用させていることが分かるわけです。
まぁf の定義からこの部分素数列の頭の7 しかいりませんが。
次はgcd 25 でためしわりになります、またns は[25,29,31,35,37,41..　となるので順当に29,31 などが残っていってめでたしですね。

やはり試し割りでgcd を使うことが大きく効きます、おのおのの約数での試し割りも一気に行えるのにコストが安いのはお買い得感が有りますね。

おそらくgcd のコストのやすさと、ns の項の形（6n+1, 6n+5) のチョイスのおかげでこの実装が早いのでしょうね。

はじめに

現在99 問題の31に差し掛かりました。
問題31 では（あるいは問題31 から）素数を扱う必要があり、取り組む前に少々お勉強した足跡です。
今までの途中経過は以下においてあります。

99 Haskell Problems より、[1..10]
99 Haskell Problems より、[11..20]
99 Haskell Problems より、[21..28]

さて、この記事のもとねたは次です。
http://en.literateprograms.org/Sieve_of_Eratosthenes_(Haskell)
http://www.garrisonjensen.com/2015/05/13/haskell-programs-are-lies.html
https://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf

エラトステネスのふるい？

本家にある

primes = filterPrime [2..] 
  where filterPrime (p:xs) = 
          p : filterPrime [x | x <- xs, x `mod` p /= 0]

が本物の篩（ふるい）では無い、と耳かじったところから始まります、どこだったか失念してしまいましたが。
どこが本物の篩と違うかは以下で解説したいと思います。¹

で、まず最初に有限列からなる入力に対して正しい篩を実装しようと思いましたが、せっかく遅延評価なので最初から無限列を扱おうと思ってやめました。
そこで考えあぐねて見つけたのが上記の資料たちです。

アルゴリズムとデータ型

門外漢ながら、まず考えなくてはならないのはアルゴリズムとデータ型は不可分だということです。

数学者²としての観点なら手続きや証明は（極論を言えば）一つあればそれで良い、と言えますが、計算機での実装を考えると例えば時間使用量や空間使用量をきちんと見積もらないとヘタすれば計算出来ない、あるいは間違った値を答えとして返す恐れだってあるわけです。³

素数列をどういう用途で使うかで、例えばそこそこの数が生成できてでもアクセス速い構造がいいとか、なるだけ大きい素数が得られたらそれで良い、とかさらに大胆に原理的に素数が作れれば（アルゴリズムの説明のためなどで）充分だとかいろいろ拘束条件が決まります。
私個人は困ったことにアルゴリズムの紹介にはそれなり出会ってきていますが、データ構造にはそれほどで大変アンバランスな知識になっていまして、困ったものです、すなわちリストと、あとは単純な木ぐらいしか知らない。
ここでは素朴にリストとsize balanced binary trees なるもので実装されたData.Set による実装を紹介します。

素朴な実装

つまりは本家で出てくる奴とほとんど一緒のやつですね。

> primes0 :: [Integer]
> primes0 = naiveSieve (2:[3,5..])
>   where 
>     naiveSieve (p:xs) = p : naiveSieve [x | x <- xs, x `mod` p /= 0]

  *Sieve> take 10000 primes0
   .. ,104729]
  (223.14 secs, 17,790,894,656 bytes)

骨董品のような爆熱MacBookAir の初号機でも10000 素数取ってくるのに200sec のオーダーです、まぁ使い物になりません。
もちろんHaskell で遅延評価によって無限の素数列を生成できますよ、しかもこんなに単純なコードで！という用途にはまぁピッタリでそのために本家でも掲載されているのだと思います、そういう意味ではたしかにシンプルで美しい。

これは素数候補のリストの最初の数2 が素数であることを使って2 の倍数を篩にかけて（ふるいから落として、の方が合ってますかね？）その結果を再帰的にsublist の先頭の要素で篩に掛けていきます。
問題は

That’s 8 operations. When the false algorithm finds 7, it checks every number from 8 to 100, that’s 92 operations!

（稚拙な訳）：8 回でいいのに、この”間違った”アルゴリズムだと7 を見つけたあと、8 から100 までの全数チェックするんだぜ？92 回も！

細かいところはもとねたを見ていただきたいです、シンプルで美しいアニメもあって分かりやすいです！
要は7 を素数として見つけたあと次の素数が11 とわかった瞬間に、121 未満の残った数たちはもう素数と分かるのでもう振るわなくて良いということです、伝わりますかね？

この素朴な”あまりがあるかどうか？”実装ではふるわないで良い数の篩におそらく実行時間の大半を取られているので、時間がかかる訳です。
しかも振るわなくて良い数は素数候補が大きくなっていくほど大きくなるので大きい素数がほしい時こそ使いものにならないアルゴリズムになってしまってるわけです。

無限リストを組み合わせて

ここ のbetter sieve が次になります。

> merge :: (Ord a) => [a] -> [a] -> [a]
> merge xx@(x:xs) yy@(y:ys) 
>   | x <  y = x : (merge xs yy)
>   | x == y = x : (merge xs ys)
>   | x >  y = y : (merge xx ys)
>   | otherwise = error ":merge"

> diff :: (Ord a) => [a] -> [a] -> [a]
> diff xx@(x:xs) yy@(y:ys)
>   | x <  y = x : (diff xs yy)
>   | x == y = diff xs ys
>   | x >  y = diff xx ys
>   | otherwise = error ":diff"

> primes2, nonPrimes2 :: [Integer]
> primes2 = [2,3,5] ++ (diff [7,9..] nonPrimes2)
> nonPrimes2 = foldr1 merge' $ map helper $ tail primes2
>   where
>     merge' (x:xs) ys = x : (merge xs ys)
>     helper p = [n*p | n <- [p,p+2..]]

  *Sieve> take 10000 primes2
  (1.87 secs, 309,541,816 bytes)
  *Sieve> take 100000 primes2
  (49.66 secs, 7,585,213,056 bytes)

まずは（昇順にソートされていると仮定されている）無限リストたちのmerge とdiff を関数として作っておきます。
もとねたの該当箇所ではcase of 構文で書いてましたがガードの方が好きなのでそこだけ書きなおしています。
Haskell のOrd のインスタンスはtotal order なので、いずれの要素もガードの三候補で場合は尽きます。
従って３つめをotherwise で書いても良いと思います、多分。

肝心要のprimes2 たちがやはりヘビーで良いですね、というか劇的に速いコードなのにキーとなるコードはこんなにコンパクトなんですね。
ここでは２つの無限（遅延）リストを相互に見なくてはなりませんが、個人的にこういうのがパズルの答えっぽくて好きです。
まず名前の示す通り、primes2 は明らかな素数たち[2,3,5] に(残りの素数)をつなげています。
この(残りの素数）は上で用意したdiff をつかって7 から始まる奇数列から、非素数nonPrimes2 を引いて作ります。
最後に上手にnonPrimes2 を定義出来れば出来上がりです、とてもHaskell。⁴

関数適用は基本的に”カッコ”の内側⁵ から見ていかなくては行けません、のでtail primes2 からです。
これは[3,5..]という感じで始まる⁶素数列です、もちろん遅延リスト。
これにmap helper が作用します、すなわち

map helper $ tail primes2
= map helper $ [3,5..]
= [[3*3, 5*3, 7*3 ..], [5*5, 7*5, 9*5 ..] ..]

という風な無限リストの無限リストになります。

なんか既視感あるなぁと思ったらこちらで紹介されていた方法そのものですね、モチベーションも似てますがあちらのほうが網羅的で良い記事だと思います。

さて最後に

foldr1 merge'

でこの無限リストの無限リストをただの無限リストに潰します、プライムが付いているのは最初の最初の要素の9 は潰されたリストの先頭に置かれて、あとは各無限リストの先頭から大きさ順に潰されていきます。
この際p^2 という形の項が必ず出てくるわけですね、ここがミソ。
こうやって作った無限リストを[7,9..] という奇数列から取り除きます、このp は素数列から取ってきているので、つまりp^2 より小さくてかつ[7,9..] からなる奇数列に残っている数

(diff [7,9..] nonPrimes2)

は素数であることが分かります。
あとは[2,3,5] と合わせて要らないふるいをふるう手間が省けているというわけです！かしこい。

Data.Set を使って

きちんとプロファイルしたわけではありませんが、おそらく大きな素数を取ってくる際には上記のアルゴリズムより早いと思われるのが以下の実装です、アニメーションの楽しいここから。

> import qualified Data.Set as PQ

> primes3 :: [Integer]
> primes3 = 2: sieve'' [3,5..]
>   where
>     sieve'' (x:xs) = sieve''' xs (insertprime x xs PQ.empty)
>     sieve''' (x:xs) table
>       | nextComposite == x = sieve''' xs (adjust x table)
>       | otherwise          = x : sieve''' xs (insertprime x xs table)
>       where
>         (nextComposite, _) = PQ.findMin table
>     adjust x table
>       | n == x = adjust x (PQ.insert (n', ns) newPQ)
>       | otherwise = table
>       where
>         Just ((n, n':ns), newPQ) = PQ.minView table
>     insertprime p xs = PQ.insert (p^2, map (*p) xs) 

  *Sieve> take 10000 primes3
  (1.86 secs, 459,583,904 bytes)
  *Sieve> take 100000 primes3
  (32.94 secs, 6,977,828,472 bytes)

さて、なにからやるかというとやはり見たこと無い関数達の型の確認からですね

Prelude Data.Set> :type empty 
empty :: Set a
Prelude Data.Set> :type findMin
findMin :: Set a -> a
Prelude Data.Set> :type insert
insert :: Ord a => a -> Set a -> Set a
Prelude Data.Set> :type minView 
minView :: Set a -> Maybe (a, Set a)

なるほど、見たままの名前の通りみたいですね、なおfindMin という名前や、insert などの型を見るとOrd のインスタンスの際は昇順で並び替えてくれるみたいですね、というわけで少々実験。

Prelude Data.Set> fromList [1,3..21]
fromList [1,3,5,7,9,11,13,15,17,19,21]
Prelude Data.Set> :t it
it :: (Enum a, Num a, Ord a) => Set a
Prelude Data.Set> insert 4 it
fromList [1,3,4,5,7,9,11,13,15,17,19,21]
Prelude Data.Set> let aSet = fromList [1,3..21]
Prelude Data.Set> insert 8 aSet 
fromList [1,3,5,7,8,9,11,13,15,17,19,21]
Prelude Data.Set> let bSet = fromList [10,9..1]
Prelude Data.Set> bSet
fromList [1,2,3,4,5,6,7,8,9,10]
Prelude Data.Set> let cSet = fromList "Haskell"
Prelude Data.Set> cSet 
fromList "Haekls"
Prelude Data.Set> insert 'C' cSet
fromList "CHaekls"

そのようですね。
途中のステップで使われている関数達の型が分からなかったのでトップレベルに書き直してghci で読んで型推論してもらい、適当にInteger に書き直して見たのが次になります。

> primes3 :: [Integer]
> primes3 = 2: sieve'' [3,5..]
> --  where
> sieve'' :: [Integer] -> [Integer]
> sieve'' (x:xs) = sieve''' xs (insertprime x xs PQ.empty)

> sieve''' :: [Integer] -> PQ.Set (Integer, [Integer]) -> [Integer]
> sieve''' (x:xs) table
>   | nextComposite == x = sieve''' xs (adjust x table)
>   | otherwise          = x : sieve''' xs (insertprime x xs table)
>   where
>     (nextComposite, _) = PQ.findMin table

> adjust :: Integer -> PQ.Set (Integer, [Integer]) -> PQ.Set (Integer, [Integer])
> adjust x table
>   | n == x = adjust x (PQ.insert (n', ns) newPQ)
>   | otherwise = table
>   where
>     Just ((n, n':ns), newPQ) = PQ.minView table

> insertprime :: Integer -> [Integer] -> PQ.Set (Integer, [Integer]) -> PQ.Set (Integer, [Integer])
> insertprime p xs = PQ.insert (p^2, map (*p) xs) 

ここで謎なのが

table :: PQ.Set (Integer, [Integer])

なる型の人たちです、ともあれフローを追ってみましょうかとりあえず。

sieve'' [3,5,7,9..]
= sieve''' [5,7,9..] $ insertprime 3 [5,7,9..] PQ.empty
= sieve''' [5,7,9..] $ PQ.insert (3^2, map (*3) [5,7,9..]) PQ.empty
= sieve''' [5,7,9..] $ PQ.insert (9, [15, 21, 27..]) PQ.empty
= sieve''' [5,7,9..] PQ.fromList [(9,[15,21,27..])]

さて、これがsieve''' のガードに入ります、where で呼ばれたfindMin は

*Sieve> PQ.findMin $ PQ.fromList [(9,[15,21,27])]
(9,[15,21,27])

となり⁷、nextComposite は9 にバインドされます。
従ってotherwise にヒットするので

sieve'' [3,5,7,9..]
= 3 : sieve''' [5,7,9..] $ insertprime 3 [5,7,9..] PQ.fromList [(9,[15,21,27..])]
= 3 : sieve''' [5,7,9..] $ PQ.insert (9, [15,21,27..]) PQ.fromList [(9,[15,21,27..])]
= 3 : sieve''' [5,7,9..] PQ.fromList [(9,[15,21,27..])]

ここで少々不安なのは無限リストでもinsert がダブリに対しても機能するのかというところですがまぁうまくいくのでしょうね、ちなみに有限リストならば大丈夫でした。

*Sieve> PQ.insert (9, [15,21,27]) $ PQ.fromList [(9,[15,21,27])]
PQ.fromList [(9,[15,21,27])]

このままガードのひとつ目の条件にヒットするまで続きます、つまり9 が落ちるわけですね。
アルゴリズム的には無限リストを組み合わと同じですね。

9 にヒットするとこんどはadjust が呼ばれます。

sieve'' [3,5,7,9..]
= 3:5:7: sieve'' [11,13..] $ adjust 9 PQ.fromList [(9,[15,21,27..])]
= 3:5:7: sieve'' [11,13..] $ adjust 9 PQ.fromList [(9,[15,21,27..])]

ここでadjust のwhere では

Just ((n, n':ns), newPQ) 
= PQ.minView PQ.fromList [(9,[15,21,27..])]
= Just ((9,[15,21,27..]),PQ.fromList [])

にバインドされます、まぁここでも本当は無限列ですね、止まりません、有限リストで試したのが次のになります。

*Sieve> PQ.minView $ PQ.fromList [(9,[15,21,27])]
Just ((9,[15,21,27]), PQ.fromList [])

n もx も9 にバインドされているのでadjust のガードも一つ目にヒットします、これが再帰的にadjust をまた読んで今度は15 にバインドされたn' のおかげでotherwise にヒットしてtable が次のように更新されます

PQ.insert (15, [21,27..]) PQ.fromList []
= PQ.fromList [(15, [21,27..])]

よって次は13 が素数と判定されて、と続いていきます。

さいごに

正直Data.Set を使った実装は始めの幾つかを追うだけで理解と言うには程遠いですが、なかなか楽しかったですね。
個人的にはやはり無限リストを組み合わせた二つ目がバランスが良い気がして好みですね。

コード全景

以下が

Sieve.lhs

> module Sieve where
> import qualified Data.Set as PQ

> primes0 :: [Integer]
> primes0 = naiveSieve (2:[3,5..])
>   where 
>     naiveSieve (p:xs) = p : naiveSieve [x | x <- xs, x `mod` p /= 0]

  *Sieve> take 10000 primes0
  .. ,104729]
  (223.14 secs, 17,790,894,656 bytes)

> primes1 :: [Integer]
> primes1 = 2: sieve' [3] [5,7..]

> sieve' :: [Integer] -> [Integer] -> [Integer]
> sieve' (p:ps) xs = p: sieve' (ps ++ ps') [x | x <- qs, x `mod` p /= 0]
>   where 
>     (ps', qs) = span (<p^2) xs

  *Sieve> take 10000 primes1
  (10.46 secs, 4,485,491,952 bytes)

  From 
  http://en.literateprograms.org/Sieve_of_Eratosthenes_(Haskell)#chunk%20def:primes_better

> merge :: (Ord a) => [a] -> [a] -> [a]
> merge xx@(x:xs) yy@(y:ys) 
>   | x <  y = x : (merge xs yy)
>   | x == y = x : (merge xs ys)
>   | x >  y = y : (merge xx ys)
>   | otherwise = error ":merge"

> diff :: (Ord a) => [a] -> [a] -> [a]
> diff xx@(x:xs) yy@(y:ys)
>   | x <  y = x : (diff xs yy)
>   | x == y = diff xs ys
>   | x >  y = diff xx ys
>   | otherwise = error ":diff"

> primes2, nonPrimes2 :: [Integer]
> primes2 = [2,3,5] ++ (diff [7,9..] nonPrimes2)
> nonPrimes2 = foldr1 merge' $ map helper $ tail primes2
>   where
>     merge' (x:xs) ys = x : (merge xs ys)
>     helper p = [n*p | n <- [p,p+2..]]

  *Sieve> take 10000 primes2
  (1.87 secs, 309,541,816 bytes)
  *Sieve> take 100000 primes2
  (49.66 secs, 7,585,213,056 bytes)

  From
  http://www.garrisonjensen.com/2015/05/13/haskell-programs-are-lies.html

> primes3 :: [Integer]
> primes3 = 2: sieve'' [3,5..]
> --  where
> sieve'' :: [Integer] -> [Integer]
> sieve'' (x:xs) = sieve''' xs (insertprime x xs PQ.empty)
> sieve''' :: [Integer] -> PQ.Set (Integer, [Integer]) -> [Integer]
> sieve''' (x:xs) table
>   | nextComposite == x = sieve''' xs (adjust x table)
>   | otherwise          = x : sieve''' xs (insertprime x xs table)
>   where
>     (nextComposite, _) = PQ.findMin table
> adjust :: Integer -> PQ.Set (Integer, [Integer]) -> PQ.Set (Integer, [Integer])
> adjust x table
>   | n == x = adjust x (PQ.insert (n', ns) newPQ)
>   | otherwise = table
>   where
>     Just ((n, n':ns), newPQ) = PQ.minView table
> insertprime p xs = PQ.insert (p^2, map (*p) xs) 

  *Sieve> take 10000 primes3
  (1.86 secs, 459,583,904 bytes)
  *Sieve> take 100000 primes3
  (32.94 secs, 6,977,828,472 bytes)

1. ここでは本物の篩の定義をちゃんとしてないので、概説だけですが。というかあまりに遅くて使い物にならんというのが実際のところな気がしますが。 ↩

2. 私は数学者ではありません数学者の側、という程度でしょうか？まぁあくまでユーザーです。 ↩

3. Haskell で私のような考え無しが素朴に実装を考えると動くけどリソースを食うコードになりがちです、まぁそれでも動けるコードがかけるからHaskell を勉強しているわけですが。 ↩

4. すなわち手続きではなく定義を書けばよい、という意味でです。 ↩

5. もちろん($) もHaskell では括弧です。 ↩

6. これは本当は正しくありません、というのもHaskell の構文糖衣では[3,5..] == [3,5,7,9..]となって素数列ではありません。が、私の表現力ではここではこう書くしかありませんでした。とりあえず3 と5 は要素に入っています、というくらいの意味です。 ↩

7. もちろん本当は無限列で評価してみたらもちろん止まりませんが、遅延評価のおかげでnextComposite が9 にバインドされたらあとは何も困ることはありません。 ↩