Haskell中Sundaram的高效On ^ 2筛分

编程问答 2022-06-02

问题描述

SO 上有很多答案解释了如何在 Haskell 中实现 Sundaram 的筛选，但它们都......真的效率低下？

我见过的所有解决方案都是这样工作的：

  1. 找出要排除的数字 <= n
  2. 从[1..n]中过滤这些
  3. 修改剩余的数字* 2 + 1

例如，我的实现是查找 1 和 2n+2 之间的所有素数：

sieveSundaram :: Integer -> [Integer]
sieveSundaram n = map (\x -> 2 * x + 1) $ filter (flip notElem toRemove) [1..n]
  where toRemove = [i + j + 2*i*j | i <- [1..n],j <- [i..n],i + j + 2*i*j <= n]

我遇到的问题是，filter 必须为 toRemove 的每个元素遍历整个 [1..n] 列表，因此这具有复杂度 O(n^3) 而简单的迭代实现具有复杂度 O(n^2)。我如何在 Haskell 中实现这一点？

解决方法

根据评论，base 不应被视为 Haskell 的完整标准库。每个 Haskell 开发人员都知道和使用几个关键包，它们会被视为 Haskell 的事实上标准库的一部分。

通过“直接迭代实现”，我假设您的意思是标记和清除一组标志？通常为此使用 Vector 或 Array。 （两者都被视为“标准”。） O(n^2) Vector 解决方案如下所示。尽管它在内部使用了一个可变向量，但批量更新运算符 (//) 隐藏了这一事实，因此您可以用典型的 Haskell 不可变和无状态风格编写它：

import qualified Data.Vector as V

primesV :: Int -> [Int]
primesV n = V.toList                           -- the primes!
  . V.map (\x -> (x+1)*2+1)                    -- apply transformation
  . V.findIndices id                           -- get remaining indices
  . (V.// [(k - 1,False) | k <- removals n])  -- scratch removals
  $ V.replicate n True                         -- everyone's allowed

removals n = [i + j + 2*i*j | i <- [1..n],j <- [i..n],i + j + 2*i*j <= n]

另一种更简单的可能性是 IntSet，它基本上是一组带有 O(1) 插入/删除和 O(n) 有序遍历的整数。 （这类似于注释中建议的 HashSet，但专门用于整数。）这是在 containers 包中，另一个实际上与 GHC 源捆绑在一起的“标准”包，尽管它不同于base。它给出了一个 O(n^2) 的解决方案，看起来像：

import qualified Data.IntSet as I

primesI :: Int -> [Int]
primesI n = I.toAscList               -- the primes!
  . I.map (\x -> x*2+1)               -- apply transformation
  $ I.fromList [1..n]                 -- integers 1..n ...
    I.\\ I.fromList (removals n)      -- ... except removals

请注意，另一个重要的性能改进是使用更好的 removals 定义，以避免过滤所有 n^2 组合。我相信以下定义会产生相同的删除列表：

removals :: Int -> [Int]
removals n = [i + j + 2*i*j | j <- [1..(n-1) `div` 3],i <- [1..(n-j) `div` (1+2*j)]]

而且我认为是 O(n log(n))。如果你将它与上面的 primesV 或 primesI 一起使用，这是瓶颈，所以我认为最终的整体算法应该是 O(n log(n))。

,

这个问题没有定义低效是什么意思。 OP 似乎一直在使用 Haskell Lazy List 解决方案，这从一开始就是低效的，因为列表上的操作是顺序的，并且在需要为每个包含许多内部“管道”部分的元素分配内存方面具有很高的恒定开销实现可能的懒惰。

正如我在评论中提到的，Sundaram 筛的原始定义是模糊的，并且由于过度订阅了表示奇数的范围而包含许多冗余操作；它可以像那里描述的那样大大简化。

然而，即使在最小化 SoS 的低效率之后，如果使用 List 是人们想要的方式：正如 OP 所确定的那样，对列表的简单重复过滤效率不高，因为每个 List 元素都会有许多重复操作以下修改后的 OP 代码：

sieveSundaram :: Int -> [Int]
sieveSundaram n = map (\x -> 2 * x + 3) $ filter (flip notElem toRemove) [ 0 .. lmt ]
  where lmt = (n - 3) `div` 2
        sqrtlmt = (floor(sqrt(fromIntegral n)) - 3) `div` 2
        mkstrtibp i = ((i + i) * (i + 3) + 3,i + i + 3)
        toRemove = concat [ let (si,bp) = mkstrtibp i in [ si,si + bp .. lmt ]
                            | i <- [ 0 .. sqrtlmt ] ]

main :: IO ()
main = print $ sieveSundaram 1000

由于改进的串联 O(n log n) 列表中有 toRemove 个值，并且必须对所有这些值进行扫描以获得所有奇数到筛分极限，因此其渐近复杂度为 {{1} }，它回答了这个问题，但不是很好。

最快的列表素数过滤技术是懒惰地合并生成的复合剔除列表的树（而不是仅仅连接它），然后生成所有不在合并复合中的奇数的输出列表（以递增的顺序，以避免每次都扫描整个列表）。使用线性合并不是那么有效，但我们可以使用无限树状合并，这只会花费额外的 O(n^2 log n) 因子，当乘以 log n 的复杂度时正确 Sieve of Sundaram 的剔除值给出了 O(n log n) 的综合复杂度，这大大低于之前的实现。

这种合并是有效的，因为每个连续的复合剔除List都是从最后一个奇数的平方开始加2，所以整个List of List的每个基值的剔除序列List的第一个值已经是递增的；因此，List of List 的简单合并排序不会竞争并且很容易实现：

O(n log^2 n)

当然，人们必须问“为什么要筛选 Sundaram？”这个问题。当去除原始 SoS 公式的残骸时 (see the Wikipedia article)，很明显 SoS 和 Eratosthenes 的 Odds-Only Sieve 之间的唯一区别是 SoS 不过滤基数剔除奇数仅适用于像 Odds-Only SoE 那样的质数。以下代码仅对找到的基本素数进行递归反馈：

primesSoS :: () -> [Int]   
primesSoS() = 2 : sel 3 (_U $ map(\n -> [n * n,n * n + n + n..]) [ 3,5.. ]) where
  sel k s@(c:cs) | k < c     = k : sel (k+2) s  -- ~= ([k,k + 2..] \\ s)
                 | otherwise =     sel (k+2) cs --      when null(s\\[k,k + 2..]) 
  _U ((x:xs):t) = x : (merge xs . _U . pairs) t -- tree-shaped folding big union
  pairs (xs:ys:t) = merge xs ys : pairs t
  merge xs@(x:xs') ys@(y:ys') | x < y     = x : merge xs' ys
                              | y < x     = y : merge xs  ys'
                              | otherwise = x : merge xs' ys'

cLIMIT :: Int
cLIMIT = 1000

main :: IO ()
main = print $ takeWhile (<= cLIMIT) $ primesSoS()

fixpoint primesSoE :: () -> [Int] primesSoE() = 2 : _Y ((3:) . sel 5 . _U . map (\n -> [n * n,n * n + n + n..])) where _Y g = g (_Y g) -- = g (g (g ( ... ))) non-sharing multistage fixpoint combinator sel k s@(c:cs) | k < c = k : sel (k+2) s -- ~= ([k,k + 2..]) _U ((x:xs):t) = x : (merge xs . _U . pairs) t -- tree-shaped folding big union pairs (xs:ys:t) = merge xs ys : pairs t merge xs@(x:xs') ys@(y:ys') | x < y = x : merge xs' ys | y < x = y : merge xs ys' | otherwise = x : merge xs' ys' cLIMIT :: Int cLIMIT = 1000 main :: IO () main = print $ takeWhile (<= cLIMIT) $ primesSoE() 组合子负责递归，其余部分相同。此版本将复杂度降低了一个 _Y 因子，因此现在渐近复杂度为 log n。

如果真的想提高效率，就不要使用 List 而是使用可变数组。以下代码使用内置的位打包数组将 SoS 实现到固定范围：

O(n log n log log n)

渐近复杂度为 {-# LANGUAGE FlexibleContexts #-} import Control.Monad.ST ( runST ) import Data.Array.Base ( newArray,unsafeWrite,unsafeFreezeSTUArray,assocs ) primesSoSTo :: Int -> [Int] -- generate a list of primes to given limit... primesSoSTo limit | limit < 2 = [] | otherwise = runST $ do let lmt = (limit - 3) `div` 2 -- limit index! oddcmpsts <- newArray (0,lmt) False -- indexed true is composite let getbpndx i = (i + i + 3,(i + i) * (i + 3) + 3) -- index -> bp,si0 cullcmpst i = unsafeWrite oddcmpsts i True -- cull composite by index cull4bpndx (bp,si0) = mapM_ cullcmpst [ si0,si0 + bp .. lmt ] mapM_ cull4bpndx $ takeWhile ((>=) lmt . snd) -- for bp's <= square root limit [ getbpndx i | i <- [ 0.. ] ] -- all odds! oddcmpstsf <- unsafeFreezeSTUArray oddcmpsts -- frozen in place! return $ 2 : [ i + i + 3 | (i,False) <- assocs oddcmpstsf ] cLIMIT :: Int cLIMIT = 1000 main :: IO () main = print $ primesSoSTo cLIMIT ，以下代码对 Odds-Only SoE 执行相同的操作：

O(n log n)

渐近效率为 {-# LANGUAGE FlexibleContexts #-} import Control.Monad.ST ( runST ) import Data.Array.Base ( newArray,assocs ) primesSoETo :: Int -> [Int] -- generate a list of primes to given limit... primesSoETo limit | limit < 2 = [] | otherwise = runST $ do let lmt = (limit - 3) `div` 2 -- limit index! oddcmpsts <- newArray (0,lmt) False -- when indexed is true is composite oddcmpstsf <- unsafeFreezeSTUArray oddcmpsts -- frozen in place! let getbpndx i = (i + i + 3,si0 + bp .. lmt ] mapM_ cull4bpndx $ takeWhile ((>=) lmt . snd) -- for bp's <= square root limit [ getbpndx i | (i,False) <- assocs oddcmpstsf ] return $ 2 : [ i + i + 3 | (i,False) <- assocs oddcmpstsf ] cLIMIT :: Int cLIMIT = 1000 main :: IO () main = print $ primesSoETo cLIMIT 。

这两个最新版本可能比它们的 List 等价物快一百倍，因为在可变数组操作而不是列表操作中减少了常数因子执行时间，以及减少了渐近线中的 O(n log log n) 因子复杂性。

haskellhaskelloptimizationoptimizationprimessievetime-complexity