performance - Haskell 中的素筛

Question

我对 Haskell 很陌生，我只是想求前 200 万个素数的和。我正在尝试使用筛子生成素数(我认为埃拉托色尼筛子？)，但它真的很慢，我不知道为什么。这是我的代码。

sieve (x:xs) = x:(sieve $ filter (\a -> a `mod` x /= 0) xs)
ans = sum $ takeWhile (<2000000) (sieve [2..])

提前致谢。

Answer 1

它非常慢，因为该算法是一种试除法，不会停止于平方根。

如果您仔细观察该算法的作用，您会发现对于每个素数 p ，其不具有较小素因数的倍数将从候选列表中删除(具有较小素因数的倍数之前已被删除)。

因此，每个数字都会除以所有素数，直到它作为其最小素数除数的倍数被删除，或者如果它是素数，它就会出现在剩余候选列表的开头。

对于合数来说，这并不是特别糟糕，因为大多数合数都有小的素因数，在最坏的情况下，最小的素因数为 n不超过√n .

但是素数除以所有个较小的素数，因此在发现第 k^th 个素数是素数之前，它已除以所有 k-1较小的素数。如果有m低于极限的素数 n ，找到所有素数所需的工作是

(1-1) + (2-1) + (3-1) + ... + (m-1) = m*(m-1)/2

部门。由Prime number theorem ，低于 n 的素数个数渐近n / log n (其中 log 表示自然对数)。消除复合 Material 的工作可以粗略地限制在 n * √n师，所以不能太小n与花在素数上的工作相比，这是可以忽略不计的。

对于 200 万素数，特纳筛大约需要 10¹⁰ 个分区。此外，它需要解构和重建大量的列表单元。

试除法止于平方根，

isPrime n = go 2
  where
    go d
      | d*d > n        = True
      | n `rem` d == 0 = False
      | otherwise      = go (d+1)

primes = filter isPrime [2 .. ]

需要少于 1.9*10⁹ 个分区(粗略估计，如果每个 isPrime n 检查到 √n - 实际上，只需要 179492732，因为复合 Material 是通常很便宜)⁽¹⁾ 并且列表操作要少得多。此外，通过跳过偶数(2 除外)作为候选除数，可以轻松改进此试除法，从而将所需除法的数量减半。

埃拉托色尼筛不需要任何划分，只使用O(n * log (log n))操作，速度要快得多:

primeSum.hs :

module Main (main) where

import System.Environment (getArgs)
import Math.NumberTheory.Primes

main :: IO ()
main = do
    args <- getArgs
    let lim = case args of
                (a:_) -> read a
                _     -> 1000000
    print . sum $ takeWhile (<= lim) primes

并以 1000 万的限制运行它:

$ ghc -O2 primeSum && time ./primeSum 10000000
[1 of 1] Compiling Main             ( primeSum.hs, primeSum.o )
Linking primeSum ...
3203324994356

real    0m0.085s
user    0m0.084s
sys     0m0.000s

我们让试除法只运行到 100 万(固定类型为 Int ):

$ ghc -O2 tdprimeSum && time ./tdprimeSum 1000000
[1 of 1] Compiling Main             ( tdprimeSum.hs, tdprimeSum.o )
Linking tdprimeSum ...
37550402023

real    0m0.768s
user    0m0.765s
sys     0m0.002s

特纳筛只能达到 100000:

$ ghc -O2 tuprimeSum && time ./tuprimeSum 100000
[1 of 1] Compiling Main             ( tuprimeSum.hs, tuprimeSum.o )
Linking tuprimeSum ...
454396537

real    0m2.712s
user    0m2.703s
sys     0m0.005s

<小时/>

⁽¹⁾ 残酷的估计是

2000000
   ∑ √k ≈ 4/3*√2*10^9
 k = 1

计算为两位有效数字。由于大多数数字都是带有一个小质因数的合数 - 一半的数字是偶数并且只进行一次除法 - 这大大高估了所需的除法数量。

仅考虑素数即可获得所需除法数量的下限:

   ∑ √p ≈ 2/3*N^1.5/log N
 p < N
p prime

其中，对于 N = 2000000给出大约 1.3*10⁸。这是正确的数量级，但低估了一个重要的因素(随着 N 的增长，缓慢减少到 1，对于 N > 10 的增长永远不会大于 2)。

除了素数之外，素数的平方和两个相近素数的乘积也需要试除法达到(接近)√k因此，如果数量足够多，就会对整体工作做出重大贡献。

然而，处理半素数所需的除法数量受到恒定倍数的限制

N^1.5/(log N)^2

所以对于非常大的N相对于处理素数的成本来说，它可以忽略不计。但在试分割完全可行的范围内，它们的贡献仍然很大。