エラトステネスのふるいをコンパクトにする方法

Squeak Smalltalk OOPL


関連エントリー：

素数を列挙するにあたって、エラトステネスのふるいは比較的高速で有効です。たとえば、404 Blog Not Found:LLR2006 - 1,000,000(番目|まで)の素数 で、いくつかの言語で例示されているのとほぼ同じことをする Squeak を使った Smalltalk の次のコード…

[
   | max primes isPrime |
   max := 1000000.
   primes := ReadWriteStream on: #().
   isPrime := [:n |
      | flag d |
      flag := true.
      primes reset.
      [flag not or: [primes atEnd or: [(d := primes next) * d > n]]]
         whileFalse: [n \\ d = 0 ifTrue: [flag := false]].
      flag ifTrue: [primes setToEnd; nextPut: n] ifFalse: [0]].
   (2 to: max) do: [:i | isPrime value: i].
   primes contents size
]
   timeToRun / 1.0e3

を用い、1000000 までの素数が 78498 個であることを知ろうとしたとき（ブロック内を選択して print it (alt/cmd + p) で値が、全体の選択で算出に要した時間が得られます…）、手元の PowerBook G4 Ti 1GHz（Mac OS X 10.4.6）では 12 秒ほどかかりますが、エラトステネスのふるいを使ったメソッド「Integer class >> #primesUpTo:」ならば、同じ Smalltalk で書かれているにもかかわらず１秒とかかりません。

[(Integer primesUpTo: 1000000) size] timeToRun / 1.0e3   " => 0.884 "


比較のために、手元の同じ環境で さんの Perl と Ruby のコードを走らせたときは、こんな結果が出ています。私の古〜い PowerBook は、小飼さんの最新式の環境（MacBook Pro）の３〜４割程度のパワーしかないこともわかります。ｗ　orz

$ perl -v
This is perl, v5.8.7 built for darwin-2level ...

$ /usr/bin/time perl primesupto.pl 1000000
       27.77 real        24.83 user         0.14 sys
$ ruby -v
ruby 1.8.4 (2005-12-24) [powerpc-darwin8.5.0]

$ /usr/bin/time ruby primesupto.rb 1000000
       85.78 real        77.39 user         0.37 sys


ところで、この #primesUpTo: にはエラトステネスのふるいの「メモリを喰う」という弱点にも対策が施されていて、キミならどう書く 2.0 - ROUND 1 - 00:26 - mputの日記。 (2006-06-17) で提案されているお題である、1000000 個目の素数である 15485863 にも、15 秒強で、メモリを食い尽くすことなく到達できます。

(Integer primesUpTo: 15485863 + 1) size                       " => 1000000 "
[(Integer primesUpTo: 15485863 + 1) size] timeToRun / 1.0e3   " => 15.5 "

これには John Moyer が公開している方法…

が使われており、Integer class >> #largePrimesUpTo:do: にあるコメントによると、計算上、1073741823（16r40000000 - 1。SmallInteger maxVal）までの素数を得るための“ふるい”には、普通、４ギガバイトほどの容量を必要とするが、この方法で圧縮すると 27 メガバイトほどで足りるのだそうです。すげー。


しくみはちょっとややこしいです。簡単にいうと、ふるいをかける前に、周期的に現われる「確実に素数でない数」を落としておくことで、ふるいに“圧縮”をかけるということだそうです。くだんのページに紹介されている 30 個の整数の周期の場合、30 がすでに 2 * 3 * 5 なので、31 から 60 までの数のうち、素数の可能性のある整数は、30 に 1、7、11、13、17、19、23、29 を足してできる８つの整数に限ることができます。このことは、以降の n 番目の周期（30 * n 。n = 2, 3, 4 ...）についても同じです。

換言すれば、８ビット（８個のフラグ）、つまり１バイトあれば、３０個分の整数の素数判定状況を保持することができるわけです。すげー。


なお、くだんのページで公開されている C のプログラムと、Integer class >> #largePrimesUpTo:do: では、64 ビットで 2310（つまり、2 * 3 * 5 * 7 * 11）個の整数の素数判定状況を保存するバージョンが使用されています。


頭の体操を兼ねて、Ruby で同じものを書いてみました。…いや(^_^;)。書いたっちゅうよりは、Integer class >> #largePrimesUpTo:do: を“変換”しただけですので、ほぼ直訳です。本体は小飼さんのサンプルから拝借いたしました。

#!/usr/bin/ruby

def primes(upto, &block)
  return large_primes(upto, &block) if upto > 25000
  upto -= 1
  flags = Array.new(upto, true)
  (0...upto-1).each do |i|
    if flags[i] then
      prime = i + 2
      k = i + prime
      while k < upto do
        flags[k] = false
        k += prime
      end
      block.call(prime)
    end
  end
end

def large_primes(upto, &block)
  idx_limit = Math.sqrt(upto) + 1

  flags = Array.new((upto + 2309) / 2310 * 60 + 60, 0xFF)
  # flags = "\377" * ((upto + 2309) / 2310 * 60 + 60)

  primes_up_to_2310 = []
  primes(2310) {|prime| primes_up_to_2310 << prime}

  mask_bit_idx = Array.new(2310)
  bit_idx = -1
  mask_bit_idx[0] = (bit_idx += 1)
  mask_bit_idx[1] = (bit_idx += 1)

  (0..4).each {|i| block.call(primes_up_to_2310[i])}

  idx = 5
  (2..2309).each do |n|
    while primes_up_to_2310[idx] < n do idx += 1 end
    if n == primes_up_to_2310[idx] then
      mask_bit_idx[n] = (bit_idx += 1)
    elsif n % 2 == 0 || n % 3 == 0 || n % 5 == 0 || n % 7 == 0 || n % 11 == 0 then
      mask_bit_idx[n] = 0
    else
      mask_bit_idx[n] = (bit_idx += 1)
    end
  end

  (13...upto).step(2) do |n|
    if (mask_bit = mask_bit_idx[n % 2310]) != 0 then
      byte_idx = n / 2310 * 60 + (mask_bit-1 << -3)
      bit_idx = 1 << (mask_bit & 7)
      if flags[byte_idx] & bit_idx != 0 then
        block.call(n)
        if n < idx_limit then
          idx = n * n
          if (idx & 1) == 0 then idx += n end
          while idx < upto do
            if (mask_bit = mask_bit_idx[idx % 2310]) != 0 then
              byte_idx = idx / 2310 * 60 + (mask_bit-1 << -3)
              mask_bit = 255 - (1 << (mask_bit & 7))
              flags[byte_idx] = flags[byte_idx] & mask_bit
            end
            idx += 2 * n
          end
        end
      end
    end
  end
end

max = ARGV.size == 1 ? ARGV[0].to_i : 100
primes = []
primes(max) {|prime| primes << prime}
puts primes.join(" ")


試してみると、速度的には、前述の余りがゼロかを確認するものに比べて、７〜８倍くらい高速のようです。

$ /usr/bin/time ruby primesupto.ar.rb 1000000
       14.07 real        10.76 user         0.10 sys


それでもさすがに Ruby では遅すぎて、私の我慢できる時間では 1000000 個目には到達できませんでした。しかし、top や /proc/self/status でメモリ使用状況を見てみると、large_primes に分岐しないように primes の一行目をコメントアウトした場合と、large_primes を使ったときとでは、断然、large_primes を使ったときの方がメモリの消費は少ないように感じます（UNIX は詳しくないので感じるだけｗ）。


追記：
我慢して待ってみました(^_^;)。３分半くらいしょうか。

$ /usr/bin/time ruby primesupto.ar.rb 15485863
      244.98 real       214.72 user         1.46 sys


ちなみに、nurse さん改良版には負けました（入出力部分は勝手に追加しました）。でも、いい勝負をしていて、遅くてぜんぜんダメってわけでもなさそうです。

$ /usr/bin/time ruby primesupto.nr.rb 15485863
      221.45 real       195.28 user         1.12 sys


追記：
コメント欄で はらさんに ruby-list で紹介された方法を教えていただきました。ありがとうございます。これは確かに速いですね。ぶっちぎりです。

$ /usr/bin/time ruby primesupto.ao.rb 15485863
       80.99 real        69.66 user         0.43 sys


関連：