【问题标题】:Why is this prime test so slow?为什么这个主要测试这么慢?
【发布时间】:2012-07-27 22:02:30
【问题描述】:

这段代码取自《Haskell Road to Logic, Math and Programming》一书。实现了eratosthenes筛分算法,解决了Project Euler Problem 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

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

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

谢谢。

【问题讨论】:

  • 有点相关(但不是答案)-The Genuine Seive of Eratosthenes-评论中已经有一个直接链接到最佳答案,但这是通过 Lambda The Ultimate 页面。

标签: haskell primes sieve-of-eratosthenes


【解决方案1】:

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

让我们来看看它是如何工作的,希望它能让问题变得明显。简单地说,假设我们想要print 无限期地计算primes 的元素,即我们想要一个接一个地评估列表中的每个单元格。 primes 定义为:

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

因为 2 不是 0,所以 sieve 的第二个定义适用,并且 2 被添加到素数列表中,列表的其余部分是未评估的 thunk(我使用 tail 而不是模式匹配n : xs in sieve for xs,所以 tail 实际上并没有被调用,也没有在下面的代码中增加任何开销;mark 实际上是唯一的 thunked 函数):

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),您可以在继续下一个 n 之前标记整个数组,从而完全没有调用堆栈开销。

【讨论】:

  • (替代算法,请阅读this paper。)
  • ...和this one
  • “Eratosthenes 的筛子在 C 之类的语言中起作用的原因” - 你也可以在 Haskell 中进行严格的数组编程,所以我认为在这方面 Haskell 是“一种类似于 C 的语言”。
  • A = "(1) 和 (2) 适用于 C 等语言", B = "(1) 和 (2) 适用于 Haskell 等语言", P = "Sieve of Eratosthenes works" . A→P; B → P。这是否意味着 A = B?
  • 您似乎暗示嵌套mark 应用程序的“越来越大的thunk”正在增长,但每个mark 单独工作,单独工作,实际上产生一串数字,每个都从它的前身,正是因为每个mark 应用程序正在通过模式匹配强制 - 在内存中创建一个实际列表。在 Haskell 中,存储是持久的。因此,那里并没有像您所描述的那样不断增长。
【解决方案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 现在都针对每个素数&lt; n 进行检查,而不仅仅是针对n 的最小素因数的素数.

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

嵌套¹ thunk mentioned by dflemstr 的堆叠加剧了问题²,但即使没有这种情况,算法也会比 Turner 的算法更糟糕。我既震惊又着迷。


¹“嵌套”可能不是正确的词。尽管每个mark thunk 只能通过它上面的那个来访问,但它们不会引用封闭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 .. ]))))

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

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

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

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

【讨论】:

  • 感到震惊?真的很震惊:[n | n&lt;-[2..], not $ elem n [j*k | j&lt;-[1..n-1], k&lt;-[1..n-1]]]。这实际上是 proposed here on SO,虽然在 Python 中。
  • 威尔,我对这种糟糕的算法并不感到震惊。令人震惊的是,代码来自一本书 - 一本据说很好 - 而且显然没有任何警告。
  • 震惊和着迷 - 我应该更全面地引用你的话。 :) 我引用的那个算法也让我着迷,因为它也是定义的真实转录,就像 Turner sieve 一样。毕竟,素数确实都是大于 1 的非复合整数。 Haskell 毕竟不是“数学”。
  • 是的,这个算法很棒。它有一个小缺陷,列表应该从 2 开始以真正反映素数的定义,但除了这个故障之外,它就像特纳的筛子一样,是一个简洁而简洁的声明,隐藏着可怕的复杂性 - 更是如此。
  • 即使使用1s,它也会出于某种原因产生质数。 :) 但它确实是特别,不是吗。值得一看的东西。
【解决方案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 中对 Eratosthenes 的筛子进行了非常详细的处理以及各种优化here

【讨论】:

  • 我没有看到 sieveW 中的任何 mod 计算被记忆,部分原因是每个 y `mod` x 只被要求一次(xs 的元素是唯一的)和部分原因是即使xs 包含重复项,之前计算的y `mod` x 也将无法访问。
  • 是的,你是对的——应该完全考虑到这一点。将编辑答案。
【解决方案4】:

简短回答:计数筛比 Turner 的(又名“朴素”)筛慢,因为它通过顺序计数模拟直接 RAM 访问,这迫使它通过流未筛分 在标记阶段之间。这很讽刺,因为 计数 使它成为 Eratosthenes 的 “真正” 筛子,而不是特纳的试验师筛子。实际上去除倍数,就像特纳的筛子所做的那样,会打乱计数。

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

This answer 解释了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 产生的素数中。但是,当我们可以按增量计数时,为什么要按个数来计数。 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 质数,快速-n-脏测试编译加载到 GHCi 中,产生高达 100k - 150k 质数,在大约 0.5 万个质数时恶化到 O(k^1.3)


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

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

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

对于此答案中的两个版本,v1/W 在 4k43x 时更快 20x8k。 v2/v1 分别是 5.2x20k5.8x40k6.5x 更快地产生 80,000 个素数。

(作为比较,Melissa O'Neill 的优先级队列代码在 k 生成的素数中以大约 O(k^1.2) 的经验复杂度运行。当然,它的可扩展性比这里的代码要好得多。


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

P = {3,5, ...} \ &bigcup; { { p*p, p*p + 2*p, ...} | p in P }

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

【讨论】:

  • “特纳筛子使用nomult p = filter ((/=0).(`rem`p)) ...”。答案正文中的标记有问题。
  • 为您解决了这个问题,双反引号代码-带反引号双反引号也适用于问题正文。
  • @DanielFischer 非常感谢,我在内部尝试了双反引号“rem”,而不是整个事情......如果我们已经在这里,那神话般的未评估 thunk buildup(在 dflemstr 答案中重新“越来越大的 ...” - 另见 cmets)。我在这里说不存在。依赖生产者链,每个都由其消费者强制,不是像foldl 那样的嵌套重击。我们不应该能够以明确的方式使用语言吗?您在回答中也提到了这一点。我认为这是错误的。
  • 好吧,“堆叠的 thunk”可能是更好的表述。但是在第 k 个素数 thunk(切换 arg 顺序)mark p_k x (mark ... (mark 3 x (mark 2 x list))...)、嵌套或堆叠k 层之后,您确实有。所以sieve 中的每个模式匹配都会告诉最外层的 thunk “嘿,给我一个构造函数”,它将消息向下传递,然后 (_:_) 再次向上传递,几个 thunk 在途中将头部替换为 0(通常)。你有一个不断增长的嵌套/堆叠的 thunk(尽管在大小或深度上都不是二次方 - 与找到的素数的数量成正比),这会减慢速度。
  • 链式,不嵌套。每个都以自己的方式存在于内存中。它被强制通过 storage,而不是从它上面的函数内部。它不知道它上面的函数/thunk/producer/generator/transducer,它只知道如何产生它的 next 输出。同样,它只有输入“(y:ys)” - 存储 - 它不知道它来自哪里,来自链下的另一个“标记”或其他任何东西。 no mark (mark (mark ... ))) 在运行时。每个列表生成函数mark 最好被视为生成器,IMO。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 2015-06-09
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2011-03-11
相关资源
最近更新 更多