feat: etatostenes turner

Author: RyosukeDTomita

math/eratosthenes/eratosthenes.hs

@@ -0,0 +1,65 @@
+# 素数関連のアルゴリズム
+
+## エラトステネスの篩
+
+- 見つけた素数の倍数を消していく方法。
+- 事前に素数かどうかを記録する配列を用意しておく。
+- i=2の時にはn/2回の更新が、i=3の時にはn/3回の更新が行われるため、全体の更新回数はn/2 + n/3 + n/5 + ... = n(1/2 + 1/3 + 1/5 + ...) = O(n log log n)となる。
+
+### 計算量の見積もりの考え方
+
+素数定理より
+π(n) ≈ n / log n
+であるため、素数の個数はn/log n程度であることがわかる。
+そのため、確率で重み付けを行うことで、素数の逆数の和は、
+
+$$\sum_{k=1}^{n} \frac{1}{k} \cdot \frac{1}{\log k}$$
+に近似される。
+
+これを積分近似する。
+
+$$
+\int_{2}^{n} \frac{dx}{x \log x}
+$$
+
+$t = \log x$を使い、置換積分を行う。
+
+$$
+dt = \frac{dx}{x}
+$$
+
+$$
+\int \frac{dx}{x \log x}
+= \int \frac{dt}{t}
+= \log t
+$$
+
+置換を戻して
+
+$$
+\log (\log n)
+$$
+
+よって、素数の逆数の和はO(log log n)であることがわかる。
+
+### 実装例
+
+- [ ] Python
+- [ ] C
+- [x] Haskell
+
+---
+
+## Turnerの篩
+
+1. 2からnまでの数をリストに入れる。
+2. リストの先頭はpを取り出す。pは必ず素数である(2が素数であり、3の処理の結果p以下の素数すべてで割ったあまりが0にならない数がリストに残るため)。pを素数リストに追加する。
+3. pで割ったあまりが0になる数をリストから削除して2に戻る。
+
+洗練された見た目だが、エラトステネスの篩に比べて実際はO(n^2/(log n)^2)
+
+### 実装例
+
+- [ ] Python
+- [ ] C
+- [x] Haskell

math/prime_number/README.md

@@ -2,37 +2,45 @@ import Control.Monad (forM_, when)
 import Data.Vector.Unboxed qualified as VU
 import Data.Vector.Unboxed.Mutable qualified as VUM
 
-type PrimeTable = VU.Vector Bool
+type PrimeTable = VU.Vector Bool -- i番目の数字が素数かどうかをBoolで管理するtable
 
+-- 素数のリストを返す
 primeTableToList :: PrimeTable -> [Int]
 primeTableToList t =
   [ i
     | i <- [0 .. VU.length t - 1],
       t VU.! i
   ]
 
+-- PrimeTableに入っていれば素数であると判定する。
 isInPrimeTable :: PrimeTable -> Int -> Bool
 isInPrimeTable t n
   | 0 <= n && n < VU.length t = t VU.! n
   | otherwise = error "out of range"
 
+-- 先頭の値iはかならす素数なので、iの倍数を生成して素数でないものを弾いていく。
 sieve :: Int -> PrimeTable
 sieve n = VU.create $ do
   vec <- VUM.replicate (n + 1) True
+  -- 0、1は素数でない
   VUM.write vec 0 False
   VUM.write vec 1 False
   forM_ [2 .. n] $ \i -> do
     b <- VUM.read vec i
     when b $ do
       forM_ [2 * i, 3 * i .. n] $ \j -> do
+        -- ghci> [2 * i , 3 * i ..10]=[4,6,8,10]
         VUM.write vec j False
   return vec
 
+main :: IO ()
 main = do
+  -- 100までの素数を作る
   let t = sieve 100
-  print $ primeTableToList t
+  print t
+  let prime = primeTableToList t
+  print prime
 
--- print $ take 20 primes
--- print $ 7 `elem` primes
-
--- print $ 102 `elem` primes -- 素数でない数を探索すると止まらない
+  -- 素数判定
+  print $ isInPrimeTable t 47
+  print $ isInPrimeTable t 57

math/eratosthenes/eratosthenesMonad.hs math/prime_number/eratosthenes.hsmath/eratosthenes/eratosthenesMonad.hs renamed to math/prime_number/eratosthenes.hs

@@ -0,0 +1,20 @@
+-- これはTurnerのふるいであってエラトステネスのふるいではないらしい
+primes :: [Int]
+primes = 2 : sieve [3, 5 ..] -- 2+奇数の無限配列
+  where
+    -- 初期値2は素数であり、pはそれより小さい素数の合成数でないため、リストの先頭は必ず素数である
+    -- pで割ったあまりが0になる数をリストから削除することを繰り返すことで素数のリストができる。
+    sieve (p : xs) =
+      p
+        : sieve
+          [ x
+            | x <- xs,
+              x `mod` p /= 0
+          ]
+
+main :: IO ()
+main = do
+  print $ take 20 primes
+  print $ 7 `elem` primes
+
+-- print $ 102 `elem` primes -- 素数でない数を探索すると止まらない