haskell - 为什么这个主要测试这么慢？

Question

这段代码取自《Haskell Road to Logic, Math and Programming》一书。它实现了埃拉托色尼筛算法并解决了 Project Euler 问题 10。

sieve :: [Integer] -> [Integer]
sieve (0 : xs) = sieve xs
sieve (n : xs) = n : sieve (mark xs 1 n)
  where
    mark :: [Integer] -> Integer -> Integer -> [Integer]
    mark (y:ys) k m | k == m = 0 : (mark ys 1 m)
                    | otherwise = y : (mark ys (k+1) m)

primes :: [Integer]
primes = sieve [2..]

-- Project Euler #10
main = print $ sum $ takeWhile (< 2000000) primes

实际上，它的运行速度甚至比天真的素数测试还要慢。有人可以解释这种行为吗？

我怀疑这与在标记函数中迭代列表中的每个元素有关。

谢谢。

Answer 1

您正在使用此算法构建二次未评估的 thunk。该算法严重依赖惰性，这也是它无法扩展的原因。

让我们来看看它是如何工作的，希望这能让问题变得明显。简而言之，假设我们想要无穷print大的元素，即我们想要一个接一个地评估列表中的每个单元格。定义为：primesprimes

primes :: [Integer]
primes = sieve [2..]

由于 2 不是 0，因此sieve应用了第二个定义，并且 2 被添加到素数列表中，并且列表的其余部分是未评估的 thunk（我使用而不是fortail中的模式匹配n : xs，所以实际上并不是调用，并且不会在下面的代码中增加任何开销；实际上是唯一的 thunked 函数）：sievexstailmark

primes = 2 : sieve (mark (tail [2..]) 1 2)

现在我们想要第二个primes元素。因此，我们遍历代码（读者练习）并最终得到：

primes = 2 : 3 : sieve (mark (tail (mark (tail [2..]) 1 2)) 1 3)

再次相同的过程，我们要评估下一个素数......

primes = 2 : 3 : 5 : sieve (mark (tail (tail (mark (tail (mark (tail [2..]) 1 2)) 1 3))) 1 5)

这开始看起来像 LISP，但我离题了......开始看到问题了吗？对于列表中的每个元素，必须评估primes越来越多的应用程序堆栈。mark换句话说，对于列表中的每个元素，必须通过评估mark堆栈中的每个应用程序来检查该元素是否被前面的任何素数标记。因此，对于n~=2000000，Haskell 运行时必须调用函数，导致调用堆栈的深度约为……我不知道，137900（let n = 2e6 in n / log n给出下限）？类似的东西。这可能是导致减速的原因；也许vacuum可以告诉你更多（我现在没有一台同时具有 Haskell 和 GUI 的计算机）。

Eratosthenes 的筛子在像 C 这样的语言中起作用的原因是：

1. 您没有使用无限列表。
2. 由于 (1)，您可以在继续 next 之前标记整个数组n，从而完全没有调用堆栈开销。

Answer 2

不仅是 thunk 导致它非常慢，如果在有限位数组上用 C 语言实现，该算法也会非常慢。

sieve :: [Integer] -> [Integer]
sieve (0 : xs) = sieve xs
sieve (n : xs) = n : sieve (mark xs 1 n)
  where
    mark :: [Integer] -> Integer -> Integer -> [Integer]
    mark (y:ys) k m | k == m = 0 : (mark ys 1 m)
                    | otherwise = y : (mark ys (k+1) m)

对于每个素数p，该算法检查从p+1到极限的所有数字是否是 的倍数p。它不是像特纳的筛子那样通过除法来做到这一点，而是通过将计数器与素数进行比较。现在，比较两个数字比除法要快得多，但为此付出的代价n是现在检查每个数字的每个素数< n，而不仅仅是检查直到n的最小素数因子的素数。

结果是该算法的复杂度是 O(N^2 / log N) 与 O( (N/log N)^2 ) 对于 Turner 筛子（和 O(N*log (log N)) 对于一个真实的埃拉托色尼筛）。

这嵌套¹ dflemstr 提到的 thunk 堆叠加剧了问题²，但即使没有这种情况，算法也会比 Turner 的算法更糟糕。我同时感到震惊和着迷。

---

¹“嵌套”可能不是正确的词。尽管每个markthunk 只能通过它上面的一个来访问，但它们不会引用封闭 thunk 范围内的任何内容。

² 不过，thunk 的大小或深度都不是二次方的，而且 thunk 的表现相当好。为了说明，让我们假设mark用相反的参数顺序定义。然后，当发现 7 是素数时，情况是

sieve (mark 5 2 (mark 3 1 (mark 2 1 [7 .. ])))
~> sieve (mark 5 2 (mark 3 1 (7 : mark 2 2 [8 .. ])))
~> sieve (mark 5 2 (7 : mark 3 2 (mark 2 2 [8 .. ])))
~> sieve (7 : mark 5 3 (mark 3 2 (mark 2 2 [8 .. ])))
~> 7 : sieve (mark 7 1 (mark 5 3 (mark 3 2 (mark 2 2 [8 .. ]))))

下一个模式匹配 bysieve强制mark 7 1thunk，它强制mark 5 3thunk，它强制mark 3 2thunk，它强制mark 2 2thunk，它强制[8 .. ]thunk 并将头部替换为 0，并将尾部包裹在一个mark 2 1thunk 中。这会冒泡到sieve，它会丢弃 0 然后强制下一个 thunk 堆栈。

因此，对于从p_k + 1to p_(k+1)（包括）的每个数字，模式匹配 insieve强制k形成 thunk 的堆栈/链mark p r。这些中的每一个都(y:ys)从封闭的 thunk 中接收（[y ..]对于最里面的mark 2 r），并将尾部包装ys在一个新的 thunk 中，要么保持y不变，要么用 0 替换它，从而在尾部建立一个新的 thunk 堆栈/链列表到达sieve。

对于每个找到的素数，在顶部sieve添加另一个mark p rthunk，因此最后，当找到并完成第一个大于 2000000 的素数时takeWhile (< 2000000)，将有 148933 级 thunk。

这里的 thunk 堆叠不会影响复杂性，它只会影响常数因素。在我们正在处理的情况下，一个延迟生成的无限不可变列表，没有太多可以减少将控制从一个 thunk 转移到下一个的时间。如果我们正在处理一个不是延迟生成的有限可变列表或数组，就像在像 C 或 Java 这样的语言中那样，让每个都mark p完成它的完整工作会更好（这将是一个for开销更少的简单循环比函数调用/控制转移）在检查下一个数字之前，所以永远不会有超过一个标记处于活动状态并且控制传递更少。

Answer 3

好吧，你绝对是对的，它比天真的实现要慢。我从 Wikipedia 中获取了这个，并将其与您使用 GHCI 的代码进行了比较，因此：

-- from Wikipedia
sieveW [] = [] 
sieveW (x:xs) = x : sieveW remaining 
  where 
    remaining = [y | y <- xs, y `mod` x /= 0]

-- your code
sieve :: [Integer] -> [Integer]
sieve (0 : xs) = sieve xs
sieve (n : xs) = n : sieve (mark xs 1 n)
  where
    mark :: [Integer] -> Integer -> Integer -> [Integer]
    mark (y:ys) k m | k == m = 0 : (mark ys 1 m)
                    | otherwise = y : (mark ys (k+1) m)

运行给

[1 of 1] Compiling Main             ( prime.hs, interpreted )
Ok, modules loaded: Main.
*Main> :set +s
*Main> sum $ take 2000 (sieveW [2..])
16274627
(1.54 secs, 351594604 bytes)
*Main> sum $ take 2000 (sieve [2..])
16274627
(12.33 secs, 2903337856 bytes)

为了尝试了解正在发生的事情以及mark代码的工作原理，我尝试手动扩展代码：

  sieve [2..]
= sieve 2 : [3..]
= 2 : sieve (mark [3..] 1 2)
= 2 : sieve (3 : (mark [4..] 2 2))
= 2 : 3 : sieve (mark (mark [4..] 2 2) 1 3)
= 2 : 3 : sieve (mark (0 : (mark [5..] 1 2)) 1 3)
= 2 : 3 : sieve (0 : (mark (mark [5..] 1 2) 1 3))
= 2 : 3 : sieve (mark (mark [5..] 1 2) 1 3)
= 2 : 3 : sieve (mark (5 : (mark [6..] 2 2)) 1 3)
= 2 : 3 : sieve (5 : mark (mark [6..] 2 2) 2 3)
= 2 : 3 : 5 : sieve (mark (mark (mark [6..] 2 2) 2 3) 1 5)
= 2 : 3 : 5 : sieve (mark (mark (0 : (mark [7..] 1 2)) 2 3) 1 5)
= 2 : 3 : 5 : sieve (mark (0 : (mark (mark [7..] 1 2) 3 3)) 1 5)
= 2 : 3 : 5 : sieve (0 : (mark (mark (mark [7..] 1 2) 3 3)) 2 5))
= 2 : 3 : 5 : sieve (mark (mark (mark [7..] 1 2) 3 3)) 2 5)
= 2 : 3 : 5 : sieve (mark (mark (7 : (mark [8..] 2 2)) 3 3)) 2 5)

我想我可能在最后犯了一个小错误，因为 7 看起来要变成 0 并删除，但机制是明确的。这段代码只是创建了一组计数器，计数到每个素数，在正确的时刻发出下一个素数并将其传递到列表中。这相当于在简单实现中仅检查每个先前素数的除法，以及在 thunk 之间传递 0 或素数的额外开销。

这里可能还有一些我遗漏的微妙之处。Haskell 中对埃拉托色尼筛法进行了非常详细的处理，并在此处进行了各种优化。

Answer 4

简短的回答：计数筛比特纳的（又名“幼稚”）筛慢，因为它模拟了具有顺序计数的直接 RAM 访问，这迫使它沿着标记阶段之间未筛分的流传递。这是具有讽刺意味的，因为计数使它成为 Eratosthenes 的“真正”筛子，而不是 Turner 的试验师筛子。实际上去除倍数，就像特纳的筛子正在做的那样，会弄乱计数。

这两种算法都非常慢，因为它们过早地从每个找到的素数而不是它的平方开始多重消除工作，从而创建了太多不需要的流处理阶段（无论是过滤还是标记）——O(n)它们，而不是仅仅~ 2*sqrt n/log n，在产生最多n在价值上。在输入中看到7之前不需要检查 的倍数。49

这个答案解释了如何sieve被视为在其背后构建流处理“转换器”的管道，因为它正在工作：

[2..] ==> sieve --> 2
[3..] ==> mark 1 2 ==> sieve --> 3
[4..] ==> mark 2 2 ==> mark 1 3 ==> sieve 
[5..] ==> mark 1 2 ==> mark 2 3 ==> sieve --> 5
[6..] ==> mark 2 2 ==> mark 3 3 ==> mark 1 5 ==> sieve 
[7..] ==> mark 1 2 ==> mark 1 3 ==> mark 2 5 ==> sieve --> 7
[8..] ==> mark 2 2 ==> mark 2 3 ==> mark 3 5 ==> mark 1 7 ==> sieve
[9..] ==> mark 1 2 ==> mark 3 3 ==> mark 4 5 ==> mark 2 7 ==> sieve
[10..]==> mark 2 2 ==> mark 1 3 ==> mark 5 5 ==> mark 3 7 ==> sieve
[11..]==> mark 1 2 ==> mark 2 3 ==> mark 1 5 ==> mark 4 7 ==> sieve --> 11

特纳筛子用来nomult p = filter ((/=0).(`rem`p))代替mark _ p条目，但在其他方面看起来是一样的：

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

每个这样的转换器都可以实现为闭包框架（也称为“thunk”），或者具有可变状态的生成器，这并不重要。每个这样的生产者的输出直接作为输入进入其链中的继任者。这里没有未评估的 thunk，每个都被其消费者强迫，以产生其下一个输出。

所以，回答你的问题，

我怀疑这与在标记函数中迭代列表中的每个元素有关。

是的，正是。否则，他们都运行非延期计划。

---

因此，可以通过推迟流标记的开始来改进代码：

primes = 2:3:filter (>0) (sieve [5,7..] (tail primes) 9)

sieve (x:xs) ps@ ~(p:t) q 
   | x < q = x:sieve xs ps q
   | x==q  = sieve (mark xs 1 p) t (head t^2)
  where
    mark (y:ys) k p 
       | k == p    = 0 : (mark ys 1 p)      -- mark each p-th number in supply
       | otherwise = y : (mark ys (k+1) p)

O(k^1.5)根据经验，它现在在k产生的素数上运行。但是，当我们可以按增量计数时，为什么要按个数来计数。（每个第 3 个奇数9可以通过添加6, 一次又一次地找到。）然后我们可以立即剔除这些数字，而不是标记，让自己成为一个真正的 Eratosthenes 筛子（即使不是最有效的筛子）：

primes = 2:3:sieve [5,7..] (tail primes) 9

sieve (x:xs) ps@ ~(p:t) q 
   | x < q = x:sieve xs ps q
   | x==q  = sieve (weedOut xs (q+2*p) (2*p)) t (head t^2)
  where
    weedOut i@(y:ys) m s 
       | y < m = y:weedOut ys m s
       | y==m  = weedOut ys (m+s) s
       | y > m = weedOut i (m+s) s

O(k^1.2)这在生成的素数中运行在上面，k编译加载到 GHCi 中的快速 n 脏测试，产生高达 100k - 150k 的素数，退化到O(k^1.3)大约 0.5 万个素数。

---

那么这样可以实现什么样的加速呢？将 OP 代码与“维基百科”的特纳筛进行比较，

primes = sieve [2..] :: [Int]
  where
    sieve (x:xs) = x : sieve [y | y <- xs, rem y x /= 0]

W/OP 有2k的8x加速（即产生 2000 个素数）。但在4k时，这是一个加速。特纳筛似乎在产生素数时以经验复杂度运行，而计数筛在相同范围内运行。15xO(k^1.9 .. 2.3)k = 1000 .. 6000O(k^2.3 .. 2.6)

对于此答案中的两个版本，v1/W在4k和8k时20x更快。v2/v1 为20k ， 40k甚至更快，产生 80,000 个素数。43x5.2x5.8x6.5x

（作为比较，Melissa O'Neill 的优先级队列代码以O(k^1.2)经验复杂度运行，k产生的素数。当然，它的扩展性比这里的代码好得多）。

---

这是埃拉托色尼定义的筛子：

P = { 3 , 5 , ...} \ ⋃ { { p*p , p*p + 2*p , ...} | p中的p }

埃拉托色尼筛法效率的关键是直接生成多个素数，通过从每个素数开始计数（两次）素数值的增量；以及它们的直接消除，这可以通过值和地址的融合来实现，就像整数排序算法一样（仅适用于可变数组）。它是否必须产生预设数量的素数或无限期地工作并不重要，因为它总是可以分段工作。