解释这段输出素数流的Haskell代码

Question

解释这段输出素数流的Haskell代码

21

我不理解这段代码的含义：

let
  sieve (p:xs) = p : sieve (filter (\ x -> x `mod` p /= 0) xs)
in sieve [2 .. ]

有人能为我分解一下吗？我知道这里面有递归，但问题是我无法理解这个例子中的递归是如何工作的。

- nubela

9

@所有人：尽管这个算法看起来很优雅，但实际上它并不是非常高效（比试除法表现明显差），而且它肯定不是埃拉托色尼筛法的实现。请参见以下链接：http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf - Juliet

3

@Juliet：嗯...这确实是试除法。 - newacct

4

@newacct：是和不是。一方面，它确实是试除法，但另一方面，这是一个糟糕的实现（上文中的作者称其为“不忠诚的筛子”）。正确的实现只需检查一个数字是否能被之前计算的所有素数整除，直到sqrt(你要检查的数)，时间复杂度大约为theta(n * sqrt(n) / (log n)^2)。而上面的代码实际上会将输入与所有之前计算出的素数进行测试，导致时间复杂度约为theta(n^2 / (log n)^2)。 - Juliet

1

@MihaiMaruseac 这个问题怎么可能是另一个在那时候不存在的问题的重复？ - Patrick Bard

5个回答

17

这实际上非常优雅。

首先，我们定义一个函数sieve，它接受一个元素列表：

sieve (p:xs) =

在sieve的主体中，我们取列表的头部（因为我们传递了无限列表[2..]，并且2被定义为质数），并将其（惰性地！）附加到应用sieve到列表的其余部分的结果中：

p : sieve (filter (\ x -> x 'mod' p /= 0) xs)

那么让我们来看一下处理列表其余部分的代码：

sieve (filter (\ x -> x 'mod' p /= 0) xs)

我们正在对筛选后的列表应用 sieve。让我们分解一下筛选部分的作用：

filter (\ x -> x 'mod' p /= 0) xs

filter函数接受一个函数和一个列表作为参数，该函数被应用于列表中的每个元素，并保留满足函数给定条件的元素。在这种情况下，filter 接受一个匿名函数：

\ x -> x 'mod' p /= 0

这个匿名函数接受一个参数x。它检查x对p（每次调用sieve时的列表头）取模：

 x 'mod' p

如果模数不等于0：

 x 'mod' p /= 0

如果元素x等于0，则过滤掉它。否则将保留它在列表中的位置。这是有道理的：如果x可被p整除，那么x就不仅仅能被1和自身整除，还能被其他数字整除。因此，x不是质数。

- mipadi

“-1”的含义是“我们从列表的开头取出一个元素，并将其附加到...的结果中。”这样不够准确，初学者容易感到困惑。 - Will Ness

1

这与懒惰和时机有关。我们不会在保护递归中“使用”结果（使用递归调用的结果是“递归”）。“结果”将在以后被计算并放置在其位置上。此外，“结果”从未全部出现。真正定义在于逐个计算“结果组成部分”的过程。这只是我的个人看法。 - Will Ness

10

它定义了一个生成器——一个名为“筛子”的流转换器。

Sieve s = 
  while( True ):
      p := s.head
      s := s.tail
      yield p                             -- produce this
      s := Filter (nomultsof p) s         -- go next

primes := Sieve (Nums 2)

使用柯里化形式的匿名函数，相当于：

nomultsof p x = (mod x p) /= 0

Sieve和Filter都是带有内部状态和按值传递参数语义的数据构造操作。

这里我们可以看到，这段代码最显著的问题并不是使用试除法来过滤工作序列中的倍数，而是可以直接通过以p的增量计数来找到它们。如果我们用后者替换前者，结果仍将具有可怕的运行时复杂度。

不，它最显著的问题是在工作序列上过早地放置了一个Filter，当它真正应该在输入中看到质数的平方之后才执行。因此，它创建了比实际需要的二次Filter数量更多的链。它创建的Filter链太长了，其中大部分甚至根本不需要。

经过修正，过滤器的创建被推迟到适当的时刻，如下所示：

Sieve ps s = 
  while( True ):
      x := s.head
      s := s.tail
      yield x                             -- produce this
      p := ps.head
      q := p*p
      while( (s.head) < q ):
          yield (s.head)                  --    and these
          s := s.tail
      ps := ps.tail                       -- go next
      s  := Filter (nomultsof p) s

primes := Sieve primes (Nums 2)

或者在Haskell中，

primes = sieve primes [2..] 
sieve ps (x:xs) = x : h ++ sieve pt [x | x <- t, rem x p /= 0]
     where (p:pt) = ps
           (h,t)  = span (< p*p) xs

这里使用rem而不是mod，因为在某些解释器中它可以更快，而且这里的数字都是正数。

通过在问题大小点n1，n2处以t1，t2的运行时间来测量算法的本地增长次数，如logBase (n2/n1) (t2/t1)，我们得到第一种情况下的O(n^2)，第二种情况略高于O(n^1.4)（在生成的n个质数中）。

为了澄清一下，缺失的部分可以简单地在这种（想象中的）语言中定义为：

Nums x =            -- numbers from x
  while( True ):
      yield x
      x := x+1

Filter pred s =     -- filter a stream by a predicate
  while( True ):
      if pred (s.head) then yield (s.head)
      s := s.tail

见下文。

更新：有趣的是，根据A.J.T. Davie在1992年的Haskell书中所说，这段代码最初出现在David Turner的1976年SASL手册中。

primes = sieve [2..]

-- [Int] -> [Int]
sieve (p:nos) = p : sieve (remove (multsof p) nos)

实际上承认两个对remove和multsof 的实现一起进行--一对用于试除法筛选，如在此问题中，另一对用于直接由计数生成每个质数的倍数有序删除，即伊拉托色尼筛（当然都是非推迟的）。在Haskell中，

-- Int -> (Int -> Bool)                   -- Int -> [Int]
multsof p n = (rem n p)==0                multsof p = [p*p, p*p+p..]

-- (Int -> Bool) -> ([Int] -> [Int])      -- [Int] -> ([Int] -> [Int])
remove m xs = filter (not.m) xs           remove m xs = minus xs m

如果他能够推迟在这里选择实际的实现，那该多好啊...

至于被推迟的代码，在使用“列表模式”的伪代码中可以是：

    primes = [2, ...sieve primes [3..]]

    sieve [p, ...ps] [...h, p*p, ...nos] =
                     [...h, ...sieve ps (remove (multsof p) nos)]

在现代的Haskell中，可以使用ViewPatterns编写如下内容：

{-# LANGUAGE ViewPatterns #-}

primes = 2 : sieve primes [3..]

sieve (p:ps) (span (< p*p) -> (h, _p2 : nos)) 
                             = h ++ sieve ps (remove (multsof p) nos)

- Will Ness

+1. 你确定第二个复杂度是O(n^1.4)吗？你是如何得出这个结果的？ - is7s

log(t2/t1) / log(n2/n1)。经验局部时间复杂度。实际上，在测量的n值较低范围内，它略高于1.40..1.42。我在GHCi中运行了解释代码，对primes!! 1000和primes!! 2000进行了时间统计，并计算了logBase 2(1+t2/t1)（因为在这种情况下计算是累积的）。请参见haskellwiki页面的完整故事。 - Will Ness

当使用 GHC -O2 编译并独立运行时，50-100-200千个素数的时间复杂度为n^1.36、1.46。调用span函数导致非本地空间泄漏。在 inline span函数后，50-100-200-400,000个素数的时间复杂度变为n^1.36、1.43、1.43。 - Will Ness

实际上，我猜它仍然是O(n^2)。据我所知，它仍然是一种试除法算法，并且每次都必须检查与所有先前的质数是否可除。使用STUArrays的命令式可变版本可以立即计算出第一百万个质数，而这个实现需要一分钟。不过，我需要进行更多分析才能准确。 - is7s

1

是的，我的意思是每个质数都对应一个筛选函数调用。我不太知道如何推理这个惰性过滤器链，如果你能更详细地解释一下那部分就更好了。而且我之前已经阅读过Haskellwiki页面，里面有一些不错的想法。 - is7s

显示剩余4条评论

2

它正在实现埃拉托斯特尼筛法

基本上，从一个质数（2）开始，过滤掉其余整数中的所有2的倍数。过滤列表中的下一个数字必须也是质数，因此必须从剩余的数字中过滤出它的所有倍数，依次类推。

- sykora

2

它说，“某个列表的筛子是列表的第一个元素（我们称之为p），过滤掉只允许不被p整除的元素后，其余列表的筛子就是这个过滤后的列表的筛子”。然后，它通过返回从2到无穷大的所有整数的筛子（即2后面跟着所有不能被2整除的整数的筛子，以此类推）来开始处理。

我建议在学习Haskell之前先阅读The Little Schemer。

- Jonathan Feinberg

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Stephan202 · Accepted Answer

与其他人所说的相反，这个函数并没有实现真正的埃拉托斯特尼筛法。它确实返回了质数的初始序列，并且以类似的方式进行，因此可以认为它是埃拉托斯特尼筛法。

我正在解释代码，当 mipadi 发布他的答案时；我可以删除它，但由于我花了一些时间在上面，而且因为我们的答案并不完全相同，所以我会留下它。

首先，请注意你发布的代码中存在一些语法错误。正确的代码应该是：

let sieve (p:xs) = p : sieve (filter (\x -> x `mod` p /= 0) xs) in sieve [2..]

let x in y 定义了 x 并允许在 y 中使用它的定义。这个表达式的结果是 y 的结果。所以在这种情况下，我们定义一个函数 sieve 并返回将 [2..] 应用于 sieve 的结果。
现在让我们更仔细地看一下这个表达式的 let 部分：

sieve (p:xs) = p : sieve (filter (\x -> x `mod` p /= 0) xs)

这定义了一个名为sieve的函数，它以列表作为其第一个参数。
(p:xs)是一个模式，它将p与列表的头匹配，并将xs与尾部（除去头部的所有内容）匹配。
通常，p : xs是一个列表，其第一个元素是p，而xs是包含剩余元素的列表。因此，sieve返回接收到的列表的第一个元素。
不要查看列表的其余部分：

sieve (filter (\x -> x `mod` p /= 0) xs)

sieve 递归调用，因此 filter 表达式将返回一个列表。
filter 接受一个 过滤函数 和一个列表。它仅返回那些过滤函数返回 true 的列表元素。
在这个例子中，xs 是被筛选的列表，
```
(\x -> x `mod` p /= 0)
```
filter函数是一个过滤函数。
filter函数接受一个参数x，如果它不是p的倍数，则返回true。

现在sieve已经被定义，我们将其传递给[2..]，即从2开始的所有自然数的列表。因此，

数字2将被返回。所有其他是2的倍数的自然数都将被丢弃。
第二个数字是3。它将被返回。所有其他3的倍数都将被丢弃。
因此下一个数字将是5。以此类推。