Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
5
Help us understand the problem. What is going on with this article?
More than 3 years have passed since last update.

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]]
5
Help us understand the problem. What is going on with this article?
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away

Comments

No comments
Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account Login
5
Help us understand the problem. What is going on with this article?