2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

n番目の素数を表す式をHaskellで書いてみる

Last updated at Posted at 2020-03-20

経緯

ren_265さんのtweetのスレッドより,
n番目の素数を表す式を読み発想を得た.

とりあえずひたすら書いてみる


import Data.List

p :: Integer -> Integer
p n = 1 + (sigma 1 (2^n) (e2 n))
  where
    e2 :: (Integral a) => a -> a -> a
    e2 n m  = floor $ (fromIntegral n / (fromIntegral $ e1 m)) ** (1.0/ fromIntegral n)

    e1 :: (Integral a) => a -> a
    e1 m = sigma 1 m e0

    e0 :: (Integral a) => a -> a
    e0 i = floor $ (cos $ (fromIntegral (product[1..i-1] + 1)) * pi / fromIntegral i )^2 

sigma :: (Enum a, Num a) => a -> a ->  (a -> a) -> a
sigma i j f = foldl1' (\a b -> a + f b) [i..j]

(今考えるとsigmaの引数をIntegralに制約しない方が良さそうですね.)
(sigmaの定義もListを2巡しているので改善の余地がありますね.) 更新しました.

実行してみる

>> map p [1..10]
[2,3,5,7,11,13,129,172,172,172]

6つ目までしか正しく出力されなかった.

考察

172(そもそも素数ではない)以降はどこかで桁溢れか丸め誤差が起こっているのかもしれないが,129(これも素数ではない)が出力されてしまうのが謎.
私の実装もしくは理解に誤りがあるのは確実なのでさらに考察したい.
(ご意見,ご指摘がありましたらコメントにお願いします.)

後半へ続く...

追記(20.04.03)

この件に関してslack(Haskell-jp)の#questionsで質問したところ,

  • cosやpiの精度に問題がある可能性がある.
  • ウィルソンの定理をそのまま書き換えたら動いた.

とのご指摘を受けました.

e0 :: (Integral a) => a -> a
e0 i = if (product[1..i-1] + 1) `mod` i == 0 then 1 else 0

と書き換えると,

>> map p [1..10]
[2,3,5,7,11,13,17,19,23,29]

と正しく計算されます.lotsさん,mizunashi-manaさん,ありがとうございました.

2
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?