エラトステネスの小型高速遠心分離機

全国の Haskell使いのみなさん Golferのみなさん，縮みましておめでとうございます．今年もよろしくおねがいします．

今回の話題は2011年ということもあり素数問題御用達エラトステネスの篩．正直いろんな人がとり上げてるテーマだから本題のネタがかぶってないかと戦々恐々だ．

よくHaskellの入門記事にも出るみたいなので素数無限リストの基本型はみんな10秒くらいで書けると思われるけど一応紹介すると．

-- 37打基本型
p=f[2..];f(p:x)=p:f[n|n<-x,mod n p>0]

この基本型はとっても遅く，小さい素数までしか使わない問題ならまだしも，大きな素数が必要な問題になると貫禄のTLEくらってﾃｰﾚｯﾃｰFATAL K.O.なのもまた有名な話で，まさにyouはshock．anarchy golf ではそんな大きな素数はあまり出てこないけど，他のゴルフ場とか元々ゴルフ推しじゃないゲーセンでゴルフプレイすると，時間制限の壁が重くのしかかる．

ジャッジシステム「俺はこのままタイムアップでもいいんだが？」

ので，そういった問題を低打数で処理するためには短いながらもそれなりに速い篩が必要となる．

基本型はInt付けたり偶数抜きをしたりしても2倍程度の高速化もできず，そんなんじゃもう全然足りない．

-- 偶数抜き
p=2:f[3,5..];f(p:x)=p:f[n|n<-x,mod n p>0]
-- Int付け
p=f[2..];f(p:x)=p:f[n::Int|n<-x,mod n p>0]

最大公約数と多倍長整数が共に安い言語用のバリエーションとして，こんなのもあったりする

-- 45打変則型
p=1#2;a#b|gcd a b<2=b:(a*b)#(b+1)|0<1=a#(b+1)

が，これにしたって手元の環境では3倍もいかないくらいでしかない．

そんなわけで本日オススメする商品はコチラ．

-- 62打高速型
p=2:3:2#p;m#(a:b:x)=[n|n<-[a^2..b^2-2],gcd m n<2]++(m*b)#(b:x)

最早比較する意味無いレベルの速度差なので，高速型のみ1000000番目を出力する時間を測ると，

$ uname -a
Linux scarlet 2.6.26.8-co-0.7.7.1 #1 PREEMPT Fri Jul 2 17:00:57 UTC 2010 i686 Intel(R) Core(TM)2 Quad CPU Q6600 @ 2.40GHz GenuineIntel GNU/Linux
$ cat prime.hs
main=print$p!!999999
p=2:3:2#p;m#(a:b:x)=[n|n<-[a^2..b^2-2],gcd n m<2]++(m*b)#(b:x)
$ ghc --version
The Glorious Glasgow Haskell Compilation System, version 6.12.3
$ ghc -O --make prime.hs
$ time ./prime
15485863

real    0m11.260s
user    0m10.410s
sys     0m0.840s

と，大体11秒強になる．

まぁ，偶数抜きも実施して，

-- 76打高速型with偶数抜き
p=2:q;q=3:5:7:3#q;m#(a:b:x)=[n|n<-[a^2,a^2+2..b^2-2],gcd n m<2]++(m*b)#(b:x)
-- (追記)68打高速型with偶数抜き(thx koyama41)
p=2:3:1#p;m#(a:b:x)=[n|n<-[a^2..b^2-2],odd n,gcd n m<2]++(m*b)#(b:x)

までやりこめば，

$ time ./prime
15485863

real    0m6.230s
user    0m5.570s
sys     0m0.660s

と，6秒強くらい(追記:68打版は76打版より若干遅く6.7s程度)までいくけど，14(or 6)打差はとても残念な差かな．

(1/16追記 ここから)
いろいろやってたら，こんなのができた．

-- 66打高速型with偶数抜き
p=2:3:5#p;n#x@(m:p:y)=[n|gcd m n<2]++(n+2)#last(x:[m*p:y|p^2-3<n])

koyama41さん案の68打相当の速度で66打といっそうリーズナブルに．
(1/16追記 ここまで)

問題の質によって，まにあう中で最短を選ぶってトコは言うまでもないけど，大抵基本型か今回の高速型の2択でことたりるはず．

入門記事の中でgcdの安さで素数・整数系問題の解き方がうんぬんって書いたネタのひとつとしてこんな使い方アルヨってことで今回はこれまで．