求所有质数——Eratosthenes筛法的Scala实现及证明

Eratosthenes筛法是以筛选(filter)为基本思想的计算质数的一种数学方法。求解n以内的质数步骤如下:

  1. 创建一个[2, n]的自然数列表
  2. val p=2(p是找到的第一个质数)
  3. for (x = 2; x * p <= n; x++) { 从列表中去掉x*p }
  4. p=列表中的下一个数,重复上面一步

当列表中的数遍历完,这个列表就是n以内的所有质数。我们很容易根据上述思路用任何一种语言进行实现。下面我们要利用函数式编程中的延迟求值(lazy evaluation)特性来实现求[2, +Infinity)的所有质数,在这里使用Scala语言中的Stream作为容器。

Stream的功能及操作方法类似与List,但是它的tail部分会在需要tail.head时才对tail.head进行求值,因此我们不担心定义无穷序列会耗尽系统资源,因为只有当取用时才会真正进行求值。

def sieve(s: Stream[Int]): Stream[Int] = 
  s.head #:: sieve(s.tail filter (_ % s.head != 0))

def primes = sieve(from(2))

上面三行就是全部的实现代码了。

其中,sieve就是Eratosthenes筛法中的第三步和第四步的循环,s.head的值就是当前的p,因此留下s.head,在后面所有的元素中,除去能够被p整除的数,使用sieve本身进行筛选。再将s.head prepend到新的tail之前。而from(2)将产生一个从2开始直到正无穷的Stream,作为初始的s,传给sieve方法,就能返回包含所有的质数的Stream。

让我们模拟一下Scala取primes的第一个元素的完整过程:

    primes.head
->  sieve(from(2)).head
->  (from(2).head #:: sieve(from(2).tail filter (_ % from(2).head != 0))).head
->  from(2).head
->  2

可以看到,第三行的#::后面的部分,并不需要进行求值,会等到需要的时候,再进行求值。

让我们再模拟一下取下一个元素的过程,简化的推导过程用->> 表示:

    primes.tail.head
->  sieve(from(2)).tail.head
->> (from(2).head #:: sieve(from(2).tail filter (_ % from(2).head != 0))).tail.head
->> sieve(from(2).tail filter (_ % 2 != 0)).head
->  ((from(2).tail filter (_ % 2 != 0)).head #:: sieve((from(2).tail filter (_ % 2 != 0)).tail) filter (_ % (from(2).tail filter (_ % 2 != 0)).head != 0))).head
->> (3 #:: sieve((from(2).tail filter (_ % 2 != 0)).tail) filter (_ % 3 != 0))).head
->  3

这边filter的语义比较直观,展开又过于繁琐,在这里就进行了简化。

下面我们来简单证明一下这个算法的正确性,即要证明primes中的所有元素都是质数。

首先,我们回顾一下isPrime的逻辑定义:

def isPrime(n: Int) = !((2 to n - 1) exists (n % _ == 0))

首先取s=primes,证明s不包含能被2整除的数(显然)。我们在刚才的模拟过程中已经进行了完整的证明,根据质数的定义,isPrime(2)=true。

然后要证明,如果s中不包含能被从2到s.head-1中的数整除的数(s.head是质数),s中就不包含能被2到s.tail.head-1中的数整除的数(s.tail.head是质数)。

// Condition:
    !(s exists {x => !((2 to s.head - 1) exists {y => x % y == 0})})

    s.tail.head
->> sieve(s.tail filter (_ % s.head != 0)).head

// Therefore (according to definition of "filter"):
->> !((2 to s.head) exists (s.tail.head % _ == 0)

// Hypothesis
    !isPrime(s.tail.head)
->  (2 to s.tail.head - 1) exists (s.tail.head % _ == 0)
->  (s.head + 1 to s.tail.head - 1) exists (s.tail.head % _ == 0)
->> (s.head + 1 to s.tail.head - 1) contains s.tail.head // Impossible, contradiction

    isPrime(s.tail.head)
    !(s exists {x => !((2 to s.tail.head - 1) exists {y => x % y == 0})})