« 数当て問題 - Hit & Blow - | トップページ | フィボナッチ数を一般項から求めてみる。 »

2012年10月19日 (金)

セクシー素数

先日、「Haskellで見る双子素数, いとこ素数, セクシー素数」というブログを見つけました。

以前から素数に興味があり、プログラミング言語を覚える度に、素数を求めるメソッドや関数を作ったりしてました。
「双子素数」や「三つ子素数」は知ってたんですけどね。「セクシー素数」は知りませんでした。

こちらのブログには「双子素数」「いとこ素数」「セクシー素数」などを求めるHaskell のコードが載っています。
例えば「双子素数」や「いとこ素数」は

twins = [(p1, p2) | (p1, p2) <- zip primes (tail primes), p1 + 2 == p2]
cousins = [(p, p+4) | p <- primes, isPrime (p+4)]

といった具合です。(primes と isPrime は Data.Numbers.Primes というパッケージが使用されているようです)

「双子素数」のコードに異論はないのですが、「いとこ素数」のコードはちょっと気になるところがありました。
せっかく primes という素数列が手元にあるのだから、「いとこ素数」のペアを探すのにわざわざ isPrime 使う必要はないと思うんです。
言葉だけでは分かりにくいので具体的なコードを書いてみます。(以下のコードではすべて、私が考えた primes と isPrime のコードを使用しています)

-- 素数列
-- http://d.hatena.ne.jp/notogawa/20110114/1295006865 を参考に
-- gcd の使い方が秀逸
primes :: Integral a => [a]
primes = map fromIntegral primes'
    where
      primes' = [2, 3, 5] ++ f 5 7 (drop 2 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]]

isPrime :: Integral a => a -> Bool
isPrime n = (n > 1) && isPrime' primes
    where
      isPrime' (p : ps)
          | p * p > n    = True
          | rem n p == 0 = False
          | otherwise    = isPrime' ps

cousinPrimes :: Integral a => [(a, a)]
cousinPrimes = loop primes
  where
    loop (p1 : ps1)
      | p2 - p1 == 4 = (p1, p2) : loop ps1
      | otherwise     = loop ps1
      where (p2 : _) = dropWhile (< p1 + 4) ps1

見てもらいたいのは、 cousinPrimes の

(p2 : _) = dropWhile (< p1 + 4) ps1

の部分です。

dropWhile で (p1 + 4) 未満の素数を落していくと、p2 には当然 (p1 + 4) 以上の素数のうち一番小さい素数が入ります。(p1 + 4) が素数なら、その数が必ず p2 になります。
もし (p1 + 4) が素数でなければ p2 には (p1 + 4) より大きい素数が入るので、p1 と p2 の差が 4 であることが確認できれば、「いとこ素数」のペアを一つ見つけたことになります。

ということで、時間のかかる素数判定は使わなくても「いとこ素数」は見つけられます。
ghci で確認してみると正しく動いているようです。素数判定を行なっていない分、オリジナルよりちょっぴり速くなっています。

*Main> take 10 cousinPrimes
[(3,7),(7,11),(13,17),(19,23),(37,41),(43,47),(67,71),(79,83),(97,101),(103,107)]
(0.00 secs, 514988 bytes)
*Main> take 10 cousins
[(3,7),(7,11),(13,17),(19,23),(37,41),(43,47),(67,71),(79,83),(97,101),(103,107)]
(0.00 secs, 0 bytes)
*Main> cousinPrimes !! 2000
(185869,185873)
(0.41 secs, 39850840 bytes)
*Main> cousins !! 2000
(185869,185873)
(1.51 secs, 137477412 bytes)

 

こちらのブログを読み進めていくと最後の方に「prime k-tuple」という言葉が出てきます。
これを見ていたら自作の cousinPrimes を改良してもっと一般化できそうな気がしてきました。

それで、素数の二つ組、三つ組、四つ組を求めるコードを作ってみました。引数には一つ目の素数との差を渡すことにしました。

-- 素数の二つ組
--
-- 双子素数 : prime2Tuples 2
-- いとこ素数 : prime2Tuples 4
-- セクシー素数 : prime2Tuples 6
--
prime2Tuples :: Integral a => a -> [(a, a)]
prime2Tuples d1 = loop primes
  where
    loop (p1 : ps1)
      | p2 - p1 == d1 = (p1, p2) : loop ps1
      | otherwise     = loop ps1
      where (p2 : _) = dropWhile (< p1 + d1) ps1


-- 素数の三つ組
--
-- 三つ子素数(1) : prime3Tuples 2 6
-- 三つ子素数(2) : prime3Tuples 4 6
-- セクシー素数の三つ組 : 
--   filter (\(p, _, _) -> (not . isPrime) (p + 18)) $ prime3Tuples 6 12
--
prime3Tuples :: Integral a => a -> a -> [(a, a, a)]
prime3Tuples d1 d2 = loop primes
  where
    loop (p1 : ps1)
      | checkDiff = (p1, p2, p3) : loop ps1
      | otherwise = loop ps1
      where
        (p2 : ps2) = dropWhile (< p1 + d1) ps1
        (p3 : _  ) = dropWhile (< p1 + d2) ps2
        checkDiff = map (subtract p1) [p2, p3] == [d1, d2]


-- 素数の四つ組
--
-- 四つ子素数 : prime4Tuples 2 6 8
--
prime4Tuples :: Integral a => a -> a -> a -> [(a, a, a, a)]
prime4Tuples d1 d2 d3 = loop primes
  where
    loop (p1 : ps1)
      | checkDiff = (p1, p2, p3, p4) : loop ps1
      | otherwise = loop ps1
      where
        (p2 : ps2) = dropWhile (< p1 + d1) ps1
        (p3 : ps3) = dropWhile (< p1 + d2) ps2
        (p4 : _  ) = dropWhile (< p1 + d3) ps3
        checkDiff = map (subtract p1) [p2, p3, p4] == [d1, d2, d3]

と、ここまで作って、「それぞれの関数で同じような処理を繰り返してるな」と気付きました。それで、さらにこれらの関数を一つにまとめてみました。
可変長引数をリストで与えることによって、二つ組でも三つ組でもこの関数一つで求められます。

-- Prime k-tuple の一般化
--
-- 双子素数 : primeKTuples [2]
-- いとこ素数 : primeKTuples [4]
-- セクシー素数 : primeKTuples [6]
--
-- 三つ子素数(1) : primeKTuples [2, 6]
-- 三つ子素数(2) : primeKTuples [4, 6]
-- セクシー素数の三つ組 :
--     filter (\(p : _) -> (not . isPrime) (p + 18)) $ primeKTuples [6, 12]
--
-- 四つ子素数 : primeKTuples [2, 6, 8]
-- セクシー素数の四つ組 : primeKTuples [6, 12, 18]
--
-- 五つ子素数(1) : primeKTuples [2, 6, 8, 12]
-- 五つ子素数(2) : primeKTuples [4, 6, 10, 12]
--
-- 六つ子素数 : primeKTuples [4, 6, 10, 12, 16]
--
primeKTuples :: Integral a => [a] -> [[a]]
primeKTuples ds = loop primes
  where
    loop (p : ps)
      | checkDiff = (p : ps') : loop ps
      | otherwise = loop ps
      where
        checkDiff = map (subtract p) ps' == ds
        ps' = map (head . (\d -> dropWhile (< p + d) ps)) ds

さらに、「三つ子素数」や「五つ子素数」が二つに分れているのが気になったので、マージソートで使われる merge 関数を導入してみました。

-- merge を使用すると「三つ子素数」や「五つ子素数」を一つにまとめられる。
-- 三つ子素数 : merge (primeKTuples [2, 6]) (primeKTuples [4, 6])
-- 五つ子素数 : merge (primeKTuples [2, 6, 8, 12]) (primeKTuples [4, 6, 10, 12])
merge :: Ord a => [a] -> [a] -> [a]
merge xs [] = xs
merge [] ys = ys
merge xs@(x : xs') ys@(y : ys')
  | x <= y = x : merge xs' ys
  | otherwise = y : merge xs ys'

例えば「三つ子素数」ならこんな感じで、まとめて小さい順に表示できます。

*Main> take 10 $ merge (primeKTuples [2, 6]) (primeKTuples [4, 6])
[[5,7,11],[7,11,13],[11,13,17],[13,17,19],[17,19,23],[37,41,43],[41,43,47],[67,71,73],[97,101,103],[101,103,107]]

HTML generated by org-mode 6.33x in emacs 23

« 数当て問題 - Hit & Blow - | トップページ | フィボナッチ数を一般項から求めてみる。 »

Haskell」カテゴリの記事

数学」カテゴリの記事

コメント

素数のペアの場合、まだまだ、たくさんあります。ソフィジェルマン素数と言ってp、2p+1のペアです。さらに、これを拡張すると、いくらでも、面白い素数のペアがつくれます。ちなみに、
ソフィジェルマン素数の分布のCtwは、3/128です。双子素数とおなじですが、そのあとに、(1/logx)*(1/log(2x))隣、実際、120億から120億10000までに、双子素数は29個,ソフィは24個すこし少ない、さらに間隔6のペアだと、61個ある。間隔が偶数であれば、そのペア素数は無限にある、2p+1の替わりに、3p+2としてよい、これを名前がないので、MAY素数とでもつけている。
ちなみに,ソフィの場合横に、4p+3、8p+7、16p+15、とチェインになっている。このぶんぷはとてもすくないが、ある。63419、126839、253679、507359、1014719、2029439などです。
数を高等数学で、見ると、見えない本当の深層世界にある、構造が色々、見せてくれます。その他素因数分解も割らずとも、解に至る方法があります、すべて深層構造が、具体的で、現実と俯瞰することに、よって見えてきます。科学の、カシミール効果のように、色即是空の世界とも、つながってきます。
数は、聖書の書き出しに始めに、言葉ありき、と書かれているのと同じように始めに存在ありきを示唆し、存在が無と存在を生成し、決して無は存在を生成しません。
とても不思議で奥深い世界です。

コメントを書く

コメントは記事投稿者が公開するまで表示されません。

(ウェブ上には掲載しません)

トラックバック

この記事のトラックバックURL:
http://app.cocolog-nifty.com/t/trackback/112020/55926455

この記事へのトラックバック一覧です: セクシー素数:

« 数当て問題 - Hit & Blow - | トップページ | フィボナッチ数を一般項から求めてみる。 »

2016年7月
          1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
31            
フォト

最近のトラックバック

無料ブログはココログ