経緯
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さん,ありがとうございました.