Haskellで素数 2 #アルゴリズム

Haskellで素数の続き。
コメントで教えていただいた記事を勉強してみた&自分でもやってみた。
99 Haskell Problems より、[31] でもその前に
 Haskell で「素数列」をつくる。
最新(修正)版を引用。

primes.hs

--
-- 素数列（修正版）
--
-- http://d.hatena.ne.jp/notogawa/20110114/1295006865 を参考に
-- gcd の使い方が秀逸
--
primes :: Integral a => [a]
primes = map fromIntegral $ [2, 3] ++ primes'
    where
      primes' = [5] ++ f 1 7 primes'
      f m s (p : ps) = [n | n <- ns, gcd m n == 1] ++ f (m * p) (p * p) ps
          where ns = [x + y | x <- [s, s + 6 .. p * p - 2], y <- [0, 4]]

　まったく謎なので、ぐるぐる追ってみました。

7から23までの素数候補(6n+1,6n+5型)のリストを作って
これはあらかじめ2と3の倍数は抜いてあるのでgcd 1 n == 1 で全員合格。
25から47までの素数候補のリストを作って
gcd 1*5 n == 1でふるいにかける。
49から119までの素数候補のリストを作って
gcd 1*5*7 n == 1でふるいにかける。
...

なるほど。なんとなくわかってきた。考えた人、すごい!

自分でもやってみた

最初のトライで23までわかっているので、これを使ったらいいんじゃないの?と思いたち、ちょこっと改変してみました。

primes :: Integral a => [a]
primes = map fromIntegral $ [2, 3, 5] ++ primes' 5 5

primes' m p = zonedPrimes m p 
          ++ primes' (m * (product $ zonedPrimes m p)) (last $ zonedPrimes m p)

zonedPrimes m p = [n | n <- croppedPrimables p, gcd m n == 1]

croppedPrimables p 
  | gcd 6 (p + 1) == 6  = [ x + y |
                          x <- [p + 2, p + 8 .. (p + 2) ^ 2 - 2 ]
                          ,  y <- [0, 4]
                          ,  x + y < (p + 2) ^ 2 
                          ]
  | otherwise           = [ x + y |
                          x <- [p + 4, p + 10 .. (p + 4) ^ 2 - 2 ]
                          ,  y <- [0, 2]
                          ,  x + y < (p + 4) ^ 2 
                          ]

main = print . take 1 . dropWhile (<5764799) $ primes

変えたところは、

gcdに使った数と求められた素数を全部かけてgcd用に次回にわたす = 次回トライする素数を素数列の最後の要素にする
pの次の奇数が素数候補( 6n+1 型)とは限らなくなったので場合わけする

です。

paiza.ioでやってみたら dropWhile (< 5764799) までいけました。かなりいけてる?
でも、元のだと 24683983 まで計算できたので、悪くなっている模様。
product とか last とか重たいのかしら...

あと、素数候補の条件として、yで 4(または 2) を足しちゃうんで x + y が (p + 2) ^ 2 (または (p + 4)^2 )に達しちゃうとやだなと思って
x + y < (p + 2) ^ 2
とか
x + y < (p + 4) ^ 2
いれてるんだけど、実際はいらなさそう。
でも、なんで要らないのかうまく説明できないので、そのままいれてます。
なぜなんでしょう?

追記

あとで考えたらわかってきました。

このままだったら x + y < (p + 4) ^ 2 のほうは要るかな。でもxの範囲を適正にすればいらない。

まず、pは素数なので 6n+1 型か 6n+5 型になります。

pが 6n+5 型の場合
次の奇数 p+2 が 6n+1 型で最初の素数候補になります。
それを二乗すると (6n+1)^2 = 6(n^2 + 2n) + 1 で 6n+1 型になります。
最後の素数候補はそれ未満にしたいので、xの最大値は、ひとつ手前の 6n+1 型の素数候補、(p+2)^2-6になります。
このとき y=4 で x+y が最後の素数候補になるのですが、x+y = (p+2)^2 - 2 となり、必ず (p+2)^2 未満になります。

pが 6n+1 型の場合
次の奇数 p+2 が 6n+3 型で素数候補になりません。
その次の奇数 p+4 が最初の素数候補(6n+5 型)になります。(なので、x は 6n+5 型の素数候補をとっていきます。)
それを二乗すると (6n+5)^2 = 6(n^2 + 10n + 4) + 1 で 6n+1 型になります。
最後の素数候補はそれ未満にしたいので、xの最大値は、ふたつ手前の 6n+5 型の素数候補、(p+4)^2-8になります。( ひとつ手前の (p+4)^2-2 にすると、y=2 のとき x+y=(p+4)^2になり不適)
x は 6n+5 型なので、6n+1 型と両方の素数候補を作るための y の取り得る値は、0 または 2 です。
x=(p+4)^2-8 のとき y=2 で x+y が最後の素数候補になるのですが、x+y=(p+4)^2-6 となり、必ず (p+4)^2 未満になります。

追記2

同じアイディアで product や last を使わないで書くとどうなるかやってみました。

primes :: Integral a => [a]
primes = map fromIntegral $ [2, 3, 5] ++ primes' 5 [5]

primes' m [p] = zonedPrimes m p ++ primes' (m * p) (zonedPrimes m p)
primes' m (p:ps) = primes' (m * p) ps

zonedPrimes m p = [n | n <- croppedPrimables p, gcd m n == 1]

croppedPrimables p 
    | gcd 6 (p + 1) == 6  = [x + y | 
                             x <- [p+2,p+8..(p+2)^2-6]
                             ,  y <- [0, 4]
                             ]
    | otherwise           = [x + y |
                             x <- [p+4,p+10..(p+4)^2-2]
                             ,  y <- [0, 2]
                             , x + y < (p+4)^2
                             ]

main = print . take 1 . dropWhile (<5764799) $ primes

同じ桁数まで計算できました。同じアルゴリズムだから当然?

やっぱり元々のやつがいい。
地道に刻んでいくほうが、まとめて一気にいくより無駄なことが少ないんでしょうか...

自分で理解するために、元プログラムを若干修正したものです。
部分的に試せるように関数をトップレベルに持って来て、適当な名前を付けました。

primes :: Integral a => [a]
primes = map fromIntegral $ [2, 3] ++ primes'

primes'  = [5] ++ recursivePrimes 1 7 primes'

recursivePrimes m s (p : ps) = 
  zonedPrimes m s p ++ recursivePrimes (m * p) (p * p) ps

zonedPrimes m s p = 
  [n | n <- croppedPrimables s p, gcd m n == 1] 

croppedPrimables s p = 
  [x + y | x <- [s, s + 6 .. p * p - 2], y <- [0, 4]]