用Ruby生成素数（2009-10-26）

想看看自己写的素数生成代码与专业 coder 的差别。一开始只是想验证自己的一个想法，任何整数都可以分解为素因子的乘积，那么只用素数来测试是否会快一些呢？

这是我写的。用 @primes 来保存找到的素数。

#!ruby
class Prime
  def initialize
    @primes = []
    @n = 1
  end
  
  def isPrime?(n)
    for x in @primes
      return false if n % x == 0
      return true if x > Math.sqrt(n)
    end
    # return [] if @primes is empty
  end
  
  def next
    loop do
      @n = @n + 1
      break if isPrime?(@n)
    end
    @primes << @n
    @n
  end
  
end
p = Prime.new
p p.next # 2
p p.next # 3
p p.next # 5
p p.next # 7
p p.next # 11

测试：

require 'mathn'
require 'benchmark'

n = 5000
class MyPrime
  def initialize
    @primes = []
    @n = 1
  end
  
  def isPrime?(n)
    for x in @primes
      return false if n % x == 0
      return true if x > Math.sqrt(n)
    end
    # return [] if @primes is empty
  end
  
  def next
    loop do
      @n = @n + 1
      break if isPrime?(@n)
    end
    @primes << @n
    @n
  end
  
end

Benchmark.bm do |x|
  x.report {
    list = []
    gen = MyPrime.new
    n.times { list << gen.next }
  }
  x.report {
    list = []
    gen = Prime.new
    n.times { list << gen.succ }
  }
end

结果出乎意料：

当产生500个素数时，我的要慢一些。伤心了！

      user     system      total        real
  0.453000   0.000000   0.453000 (  0.453000)
  0.360000   0.000000   0.360000 (  0.359000)

当产生5000个素数时，我的竟然快了近三倍！！！

      user     system      total        real
 11.141000   0.000000  11.141000 ( 11.312000)
 33.390000   0.000000  33.390000 ( 33.688000)

忍不住测了一下不同数值下的结果。拷了半天机，得出下面这张图（生成 10 到 2000 以内的素数）：（横坐标单位为 10 ，纵坐标单位为秒）

速度对比

狂笑不止……

由于自己的 Prime 类在渐进性能上远远超过 mathn 库的 Prime ，于是很想看看人家是怎么写的（居然可以写得这么慢:)）。但翻出来花了好长时间才看懂，他居然只用了加法，amazing!!!

class Prime
  include Enumerable
  def initialize
    @seed = 1
    @primes = []
    @counts = []
  end
  
  def succ
    i = -1
    size = @primes.size
    while i < size
      if i == -1
    @seed += 1
    i += 1
      else
    while @seed > @counts[i]
      @counts[i] += @primes[i]
    end
    if @seed != @counts[i]
      i += 1
    else
      i = -1
    end
      end
    end
    @primes.push @seed
    @counts.push @seed + @seed
    return @seed
  end
  alias next succ
  def each
    loop do
      yield succ
    end
  end
end

其实总的思路都是只用素数作为测试的因子。他维护了两个数组 @primes 和 @counts (这名字取得太烂俗了)和已经找到的最大素数 @seed （也就是 @primes[-1] ）。显然，@primes 存在已经找到的全部素数，而 @counts 中的第 i 个元素存的是 @seed 与第 i 个素数的最小公倍数。比如，@primse 是 [2,3,5,7] ，@seed 是 7 ，那么 @counts 就是 [8,9,10,14] （暂不考虑最后一个），即 [2×4,3×3,5×2,7×2] 。

我照这个思路验证了一下。只是为了验证，性能完全不行，生成前 100 个素数时，性能差了 130 倍。

require 'mathn'
# 见《The Ruby Way》16.6节

class Integer
  def lcm(other)
    pf1 = self.prime_division.flatten
    pf2 = other.prime_division.flatten
    h1 = Hash[*pf1]
    h2 = Hash[*pf2]
    hash = h2.merge(h1) {|key,old,new| [old,new].max }
    Integer.from_prime_division(hash.to_a)
  end
end

# 主角登场
class MyPrime
  def initialize
    @seed = 1
    @primes = []
  end
  
  def succ
    loop do
      is_composite = false
      @seed += 1
      for x in @primes
        n = x.lcm @seed
        if n == @seed
          is_composite = true
          break
        end
      end
      break if not is_composite
    end
    @primes << @seed
    @seed
  end
  
  alias next succ
end

gen = MyPrime.new
p gen.next #=> 2
p gen.succ #=> 3
p gen.succ #=> 5
p gen.succ #=> 7
p gen.succ #=> 11

据说 Ruby 1.8 的 mathn 库是业余人士写的，1.9 的就快一点了（嗯，很有兴趣）

可见，同一思路的不同实现间的差别还是可以很大的。另外，再次见证了算法的提高比其他一切的优化更经济。

P.S.: 去搜了搜有关素数生成方面的资料，发现了 Atkin 筛法 (Atkin 的论文 )和据说更快的 Zakiya 筛法 。

还可以看看Bob大叔写的性能调优——永远超乎想象（不愧是大神，随手写的都能这么快。当然生成素数的方式不一样）。