HaHaHa!(old)

2006-06-18

誇大宣伝ではないエラトステネスの篩

無限リストの虚仮威し(オッ.)の例としてよくあるのがエラトステネスの篩と

称されている

 primes = [ p | (p:_) <- iterate sieve [2..] ]
 sieve (p:ps) = [ x | x <- ps, mod x p /= 0 ]

とか

 primes = sieve [2..]
 sieve (p:ps) = p : sieve [q | q <- ps, q `mod` p /= 0]

というプログラム.でも

遅延評価によるエラトステネスの篩にも書かれているように「等差数列で倍数が取り除けるということを使って初めて篩と呼べるので、これをエラトステネスの篩と呼ぶのは実は誇大宣伝だ。」

そこで,誇大宣伝ではないエラトステネスの篩の実装をば...

 type Table = ([Int],[Int],[Int]) -- (新規の素数リスト,篩に使う素数列,被篩数列)
 fst3 (x,_,_) = x

 primes :: [Int]
 primes = concatMap fst3 $ iterate sieve ([2],[2],[2..])
 sieve :: Table -> Table
 sieve (ps,q:qs,ss)
  = let (xs,ys) = span (q*q >) ss
    in (xs,qs++xs,filter ((/= 0) . (`mod` q)) ys)

(tさんのコメントを受けて,コードを修正しました: 2006/06/19)

これは 1000000個の素数 - HaHaHa!(old) - haskell

高速版より性能がよいはず.

10,000,000番目の素数は 179,424,671 (179,424,691は10000001番目の素数

した 2006/06/19).primes !! 9999999 を計算するまでに掛った時間は

CPU: Intel(R) Pentium(R) M processor 2.13GHz 
OS : Ubuntu Linux 6.06
GHC option : -O 

で,

1561.57s user 3.05s system 95% cpu 27:13.98 total

という結果でした.

tt2006/06/18 22:29sieveですが、3の倍数を除く必要があるのは9以降、5の倍数を除く必要があるのは25以降なので、以下の定義で良いのではないでしょうか。
sieve (ps,q:qs,ss) = let (xs,ys) = break (q*q <=) ss in (xs,qs++xs,[ s' | s' <- ys, mod s' q /= 0 ])

nobsunnobsun2006/06/19 07:57thanks tさん。なんかずれてましたね。私のコード。

nanasinanasi2006/07/12 17:25primes = [ p | (p:_) <- iterate sieve [2..] ]
sieve (p:ps) = [ x | x <- ps, mod x p /= 0 ]
をGHCでコンパイルすると = のところが parse error
と出てしまいます なぜですか?

RonnieRonnie2012/05/03 06:38It's wonrdfeul to have you on our side, haha!

cstsrgfpmcstsrgfpm2012/05/03 12:26itQeag <a href="http://tuofsavfuqtt.com/">tuofsavfuqtt</a>

wahuvoqwahuvoq2012/05/04 10:08sDxeOL , [url=http://njqrqsfndkxb.com/]njqrqsfndkxb[/url], [link=http://uackowdukbef.com/]uackowdukbef[/link], http://ihicqycyidfm.com/

dtbrggonhadtbrggonha2012/05/05 11:19PgcWlI <a href="http://pwgcmxdgzvmg.com/">pwgcmxdgzvmg</a>

utgqvvutgqvv2012/05/05 16:007O3WdR , [url=http://xbxlgznejfgn.com/]xbxlgznejfgn[/url], [link=http://uydpsdxqmdeu.com/]uydpsdxqmdeu[/link], http://lykfrckungkl.com/