列表,数组和可变数组之间的Eratosthenes筛选的理想实现是什么？

Question

在Haskell中,我在Rosetta Code页面上找到了三个简单的Eratosthenes Sieve实现.

现在我的问题是,应该在哪些情况下使用哪一个？

纠正我的初步推理也会有所帮助:

我假设列表一是Haskeller中最惯用且易于阅读的.但是这是正确的吗？我想知道它是否遇到了与另一个基于列表的筛子相同的问题,然后我学会了实际上没有实现算法:(
编辑:这里显示的是我知道有问题的基于列表的筛子,而不是来自RosettaCode的筛子,我贴在底部)

primes = sieve [2..] where
         sieve (p:x) = p : sieve [ n | n <- x, n `mod` p > 0 ]

在性能方面,不可变阵列似乎是赢家.随着上限m的2000000,这些时间是约:

• 列表的1.3s
• 数组为0.6s
• 可变阵列的1.8s

所以我选择Array来表现.

当然,Mutable Array也很容易推理,因为我有更迫切的语言背景.我不知道为什么我会选择这个,如果我在Haskell编码,因为它既慢于其他人,也不是非惯用的.

此处复制的代码仅供参考:

列表:

primesTo m = 2 : eratos [3,5..m] where
eratos (p : xs) | p*p>m = p : xs
                | True  = p : eratos (xs `minus` [p*p, p*p+2*p..])

minus a@(x:xs) b@(y:ys) = case compare x y of
     LT -> x : minus  xs b
     EQ ->     minus  xs ys
     GT ->     minus  a  ys
minus a        b        = a 

不可变数组:

import Data.Array.Unboxed

primesToA m = sieve 3 (array (3,m) [(i,odd i) | i<-[3..m]] :: UArray Int Bool)
  where
    sieve p a 
      | p*p > m   = 2 : [i | (i,True) <- assocs a]
      | a!p       = sieve (p+2) $ a//[(i,False) | i <- [p*p, p*p+2*p..m]]
      | otherwise = sieve (p+2) a  

可变阵列:

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

primeSieve :: Integer -> UArray Integer Bool
primeSieve top = runSTUArray $ do
    a <- newArray (2,top) True           -- :: ST s (STUArray s Integer Bool)
    let r = ceiling . sqrt $ fromInteger top
    forM_ [2..r] $ \i -> do
      ai <- readArray a i
      when ai $ do
        forM_ [i*i,i*i+i..top] $ \j -> do
          writeArray a j False
    return a

-- Return primes from sieve as list:  
primesTo :: Integer -> [Integer]
primesTo top = [p | (p,True) <- assocs $ primeSieve top]

编辑

• 我在顶部展示了特纳的筛子,但这不是我与其他两个人比较的基于列表的例子.我只是想知道基于列表的例子是否与特纳的"不正确的Eratosthenes筛子"相同.
• 似乎性能比较是不公平的,因为Mutable Array示例也通过evens Integer而不是Int......

Answer 1

这个

primes = sieve [2..] where
         sieve (p:x) = p : sieve [ n | n <- x, n `mod` p > 0 ]

不是筛子.这是非常低效的试验分工.不要用那个!

我很好奇你是如何度过你的时代的,特纳"筛子"无法在几秒钟内产生不超过2,000,000的素数.让它找到200,000的素数

MUT     time    6.38s  (  6.39s elapsed)
GC      time    9.19s  (  9.20s elapsed)
EXIT    time    0.00s  (  0.00s elapsed)
Total   time   15.57s  ( 15.59s elapsed)

在我的盒子上(64位Linux,ghc-7.6.1,用-O2编译).该算法的复杂性O(N² / log² N)几乎是二次的.让它进入2,000,000大约需要20分钟.

您对阵列版本的时间也是可疑的,但在另一个方向.你测量了解释代码了吗？

筛选到2,000,000,使用优化编译,可变数组代码运行0.35秒,不可变数组代码0.12秒.

现在,仍然有可变数组比不可变数组慢三倍.

但是,这是一个不公平的比较.对于不可变数组,您使用Ints,并使用可变数组Integer.更改可变数组代码以使用Ints - 因为它应该,因为在引擎盖下,数组是Int-indexed,所以使用Integer是不必要的性能牺牲,什么都不买 - 使可变数组代码在0.15秒内运行.接近可变数组代码,但不完全相同.但是,你让可变数组做更多工作,因为在不可变数组代码中你只消除奇数素数的奇数倍,但在可变数组代码中,你标记所有素数的所有倍数.更改可变数组代码以专门处理2,并且仅消除奇数素数的奇数倍也将其降低到0.12秒.

但是,您正在使用范围检查的数组索引,这很慢,并且,因为索引的有效性在代码本身中被检查,所以不必要.将其更改为使用unsafeRead并将unsafeWrite不可变数组的时间减少到0.09秒.

然后你有使用的问题

forM_ [x, y .. z]

使用盒装Int(幸运的是,GHC取消了列表).自己编写循环,以便只使用未装箱Int#的时间,时间减少到0.02秒.

{-# LANGUAGE MonoLocalBinds #-}
import Control.Monad (forM_, when)
import Control.Monad.ST
import Data.Array.ST
import Data.Array.Unboxed
import Data.Array.Base

primeSieve :: Int -> UArray Int Bool
primeSieve top = runSTUArray $ do
    a <- newArray (0,top) True
    unsafeWrite a 0 False
    unsafeWrite a 1 False
    let r = ceiling . sqrt $ fromIntegral top
        mark step idx
            | top < idx = return ()
            | otherwise = do
                unsafeWrite a idx False
                mark step (idx+step)
        sift p
            | r < p     = return a
            | otherwise = do
                prim <- unsafeRead a p
                when prim $ mark (2*p) (p*p)
                sift (p+2)
    mark 2 4
    sift 3

-- Return primes from sieve as list:
primesTo :: Int -> [Int]
primesTo top = [p | (p,True) <- assocs $ primeSieve top]

main :: IO ()
main = print .last $ primesTo 2000000

因此,对于Eratosthenes的Sieve来说,你应该使用一个阵列 - 这并不奇怪,它的效率取决于能够在短时间内从一个多个步进到下一个多个.

你可以使用不可变数组获得非常简单直接的代码,并且该代码可以在不太高的限制下执行(对于更高的限制,它会变得相对更糟,因为数组仍然被复制和垃圾收集,但这并不是太糟糕).

当您需要更好的性能时,您需要可变数组.编写高效的可变数组代码并不是完全无关紧要的,必须知道编译器如何转换各种结构以选择正确的结构,有些人会认为这样的代码是单一的.但是你也可以使用一个库(免责声明:我编写它),它提供了一个相当有效的实现,而不是自己编写.