Haskell 中的初筛

Question

我对 Haskell 很陌生，我只是想找出前 200 万个素数的总和。我正在尝试使用筛子生成素数（我认为是 Eratosthenes 的筛子？），但它真的很慢，我不知道为什么。这是我的代码。

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

提前致谢。

Answer 1

这是sieve of Eratosthenes：

P = {3,5, ...} \ ⋃ {{ p², p²+2p, ...} | p in P}

（没有 2）。 :) 或者在“功能”中，即基于列表的 Haskell，

primes = 2 : g (fix g)  where
   g xs = 3 : (gaps 5 $ unionAll [[p*p, p*p+2*p..] | p <- xs])

unionAll ((x:xs):t) = x : union xs (unionAll $ pairs t)  where
  pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t 
fix g = xs where xs = g xs
union (x:xs) (y:ys) = case compare x y of LT -> x : union  xs (y:ys)
                                          EQ -> x : union  xs    ys 
                                          GT -> y : union (x:xs) ys
gaps k s@(x:xs) | k<x  = k:gaps (k+2) s    
                | True =   gaps (k+2) xs 

与answer by augustss 中的试用除法代码相比，它在生成 200k 素数时快 1.9 倍，在 400k 时快 2.1 倍，O(n^1.12..1.15) 与 O(n^1.4) 的 empirical time complexity 在上述范围内。生成 100 万个素数时速度提高了 2.6 倍。

为什么特纳筛子这么慢

因为它为每个素数为时过早打开了多重过滤流，因此以太多结束。在输入中看到素数之前，我们不需要按素数进行过滤。

在stream processing paradigm 下看到，sieve (x:xs) = x:sieve [y|y<-xs, rem y p/=0] 可以看作是在其工作时在其背后创建一个流转换器管道：

[2..] ==> sieve --> 2
[3..] ==> nomult 2 ==> sieve --> 3
[4..] ==> nomult 2 ==> nomult 3 ==> sieve 
[5..] ==> nomult 2 ==> nomult 3 ==> sieve --> 5
[6..] ==> nomult 2 ==> nomult 3 ==> nomult 5 ==> sieve 
[7..] ==> nomult 2 ==> nomult 3 ==> nomult 5 ==> sieve --> 7
[8..] ==> nomult 2 ==> nomult 3 ==> nomult 5 ==> nomult 7 ==> sieve 

在哪里nomult p = filter (\y->rem y p/=0)。但是 8 还不需要检查是否可以被 3 整除，因为它小于3^2 == 9，更不用说被 5 或 7 整除了。

这是该代码中最严重的问题，尽管在每个人都提到的那篇文章的开头，它被认为是无关紧要的。 Fixing it 通过推迟过滤器的创建实现了显着的加速。

Answer 2

就个人而言，我喜欢这种生成素数的方式

primes :: [Integer]
primes = 2 : filter (isPrime primes) [3,5..]
  where isPrime (p:ps) n = p*p > n || n `rem` p /= 0 && isPrime ps n

与此处建议的其他一些方法相比，它的速度也相当快。它仍然是试除法，但它只测试素数。 （不过，这段代码的终止证明有点棘手。）

Answer 3

这很慢，因为该算法是一个不停留在平方根的试除法。

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

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

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

但是素数被所有个较小的素数除，所以直到第k^个素数被发现是素数，它已经被所有k-1个较小的素数除.如果有m 质数低于限制n，则找到所有质数所需的工作是

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

部门。通过Prime number theorem，n 以下的素数数量渐近为n / log n（其中log 表示自然对数）。消除复合材料的工作可以粗略地受到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，因为复合材料是一般便宜）⁽¹⁾ 和更少的列表操作。此外，通过跳过偶数（2 除外）作为候选除数，可以轻松改进此试验除数，从而将所需除数减半。

Eratosthenes 的筛子不需要任何除法，只使用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，它相对于处理素数的成本变得可以忽略不计。但在完全可行的试制范围内，它们的贡献仍然很大。

Answer 4

您使用的算法根本不是筛子，所以就它的速度而言，您应该期望使用试除法。

素数大致出现在平方根函数的频率......即在 1 和 n 之间有大约 n/log(n) 个素数。因此，对于前 200 万个素数，您将需要多达 3200 万个。但是您正在构建一个包含 200 万个元素的数据结构，这些素数将通过这些数据结构。所以你可以开始明白为什么这会这么慢。实际上是 O(n^2)。您可以将其减少到 O(n*(log n)*log(log n))

这里有一个关于各种治疗方法的页面，将引导您了解如何减少这种情况。 http://en.literateprograms.org/Sieve_of_Eratosthenes_(Haskell)

Answer 5

你做的不是埃拉托色尼筛；它是试用部门（注意 mod 运算符）。这是我的埃拉托色尼筛法：

import Control.Monad (forM_, when)
import Control.Monad.ST
import Data.Array.ST
import Data.Array.Unboxed

sieve :: Int -> UArray Int Bool
sieve n = runSTUArray $ do
    let m = (n-1) `div` 2
        r = floor . sqrt $ fromIntegral n
    bits <- newArray (0, m-1) True
    forM_ [0 .. r `div` 2 - 1] $ \i -> do
        isPrime <- readArray bits i
        when isPrime $ do
            forM_ [2*i*i+6*i+3, 2*i*i+8*i+6 .. (m-1)] $ \j -> do
                writeArray bits j False
    return bits

primes :: Int -> [Int]
primes n = 2 : [2*i+3 | (i, True) <- assocs $ sieve n]

您可以在http://ideone.com/mu1RN 运行它。