読者です 読者をやめる 読者になる 読者になる

エラトステネスの篩の改良

上のソースでは探索範囲1つにつき true/false を作るので、1億まで探索しようと思うとかなりのメモリが必要となる。どうせ真偽値なのだから、BitArray を用いれば候補一つにつき 1bit で済む。また、偶数はどうせ素数ではないので、奇数のみを探索対象にすることで、テーブルサイズを半分にできる。さらに、篩範囲を大きくしていくとどうしてもマシンで確保できるメモリの限界を超えてしまうので、篩範囲を分割できるようにする。といった修正を施したバージョン。NetBSD 3.0 CeleronD 2.6GHz 512MB RAM ruby 1.9 にて 100万は数秒、1千万が1分、1億が10分でした。
参考: エラトステネスの篩い @ IDM

class BitArray
  def initialize(size, bool=false)
    @str = (bool ? "\xFF" : "\x00") * (size / 8 + 1)
  end
  def [](n)
    a = n >> 3
    b = n & 7
    @str[a][b].nonzero?
  end
  def []=(n,val)
    a = n >> 3
    b = n & 7
    if val && val != 0
      @str[a] |= 1 << b if @str[a][b].zero?
    else
      @str[a] &= ~(1 << b)# if @str[a][b].nonzero?
    end
  end
  def step_false(ti,tmax,i)
    ti.step(tmax, i) do |tj|
      a = tj >> 3
      b = tj & 7
      @str[a] &= ~(1 << b)# if @str[a][b].nonzero?
    end
  end
end

def _sieve(primes, min=3, size=1000_0000)
  table = BitArray.new(size,true)
  max = min + size * 2 - 1
  tmax = size
  mid = Math.sqrt(max).floor
  primes.each do |i|
    next if i > mid
    ti = (i - min) / 2
    ti += i while ti < 0
    table.step_false(ti,tmax,i)
  end
#  min.step(mid, 2) do |i|
#    ti = (i - min) / 2
#    next unless table[ti]
  (0..((mid-min)/2)).each do |ti|
    next unless table[ti]
    i = (ti << 1) + min
    primes << i
    table.step_false(ti,tmax,i)
  end
  ((mid > min ? (mid - min) / 2 : 0)...tmax).each do |ti|
    primes << min + ti * 2 if table[ti]
  end
  return primes
end

def sieve(max=100, size=1000_0000)
  max -= 1 unless max % 2 == 0
  primes = []
  3.step(max, size * 2) do |i|
    GC.start
    if max > i + size * 2
      p([:sieve, Time.now.to_f, i, size, primes.length])
      primes = _sieve(primes, i, size)
    else
      size = (max - i) / 2 + 1
      p([:sieve, Time.now.to_f, i, size, primes.length])
      primes = _sieve(primes, i, size)
    end
  end
  return [2]+primes
end

追記: いろいろ高速化していたら、手元で1千万26秒、1億274秒まで短縮(’’
追記2: エラトステネスのふるいをコンパクトにする方法 に偶数以外に3や5の倍数も削る手法が。確かにこうすれば大幅にメモリ容量は減らすことができるのですが、探索する際に ti.step(tmax, i) do |tj| が使えなくなるため、大幅に速度が落ちてしまうような気もします。(偶数だけ3の倍数だけ、ならば使えるのですが)