名古屋で数学するプログラマ(仮)

@antimon2 が趣味兼一部本職の数学で何かするときのブログ。

汎用的で省メモリかつかなり速い素数無限列挙

台風凄かったですね。私の住んでいる地域では、夕方頃に風雨ともかなり強かったのですが、今はもうだいぶ落ち着いてきちゃいました。

さてさて。昨日の記事の続きです。
まずはPrime::A2Generatorの説明、その前に、その元になっているメモ: 上限を決めずに素数を生成の高速化(Ruby版)のPrimeGeneratorの解説から。

log.ttsky.net の PrimeGenerator の解読

大元の記事は、メモ: 素数の逐次生成(Ruby版) v3。たどればもっと元があるようですけれど取り敢えず。ここに書かれているMyPrimeクラスの高速化を図った究極形がこっちの記事のPrimeGeneratorクラス*1
このクラスが素数を列挙するアルゴリズム(succメソッドの中身)は、読み解くと以下の通り*2

  • 取り敢えず最初に2を返す。
  • 2回目以降は2以上の数を2つずつ順に判定:
    • 最初の1個目(=必ず偶数)は、ハッシュにキーとして登録されていれば「処理A」を実施する。
    • 次に2個目(=必ず奇数)は、ハッシュにキーとして登録されていれば「処理A」を実施し、登録されていなければ「処理B」を実施して結果を返して抜ける。
    • 繰り返す。
  • 処理A:
    • ハッシュに登録されている値(=実は素数)を取得 -> b
    • 元の値 n に b を、それがキーとしてハッシュに登録されなていないモノになるまで足す -> m
    • ハッシュに、キー:m、値:b を登録。
    • ハッシュからキー n(元の値)を削除する。
  • 処理B:
    • 元の値 n に n を、それがキーとしてハッシュに登録されていないモノになるまで足す -> m
    • ハッシュに、キー:m、値n を登録。
    • 値 n を返す(=実はこれが新しい素数)。

この処理Aと処理B。これが、実は「エラトステネスの篩」と同じ処理になっています。ただし、『逐次必要な部分だけを生成する』仕組みになっています。
両方に共通する「キーとしてハッシュに登録されていないモノになるまで足す」というのは、「対象の素数の倍数であって、他の素数の倍数として(まだ)マークされていない」数を探す処理になっています。つまり『次に必要な部分』を求める処理です。そんな数が見付かったら「この数は既に見付かっている素数(=b)の倍数だぞ」とマークしている(=ハッシュに登録している)、というわけ。
そして元の値 n が「マークされていない」数ならば、その数よりも小さいどんな素数の倍数でもない⇒新しい素数、というわけです。

拙作のA2Generatorの解説

Prime::A2Generatorは、このPrimeGeneratorのコードを流用しながら、この作者の気付いていない(?)無駄を省くことによるさらなる高速化を図っています。
PrimeGeneratorの無駄は、(本質的にはあまり差異の無い)以下の2点。

  1. エラトステネスの篩をマジメに全部やろうとしている。
  2. 偶数(そして3の倍数)など、見なくても良い部分を省こうとしていない。

例えば、次の n を考えるときに、 n += 1 としていますが、これを n += 2 とするだけで、偶数を省いて奇数だけを対象とすることが出来ます。こうすることで「篩にかける粉の粗さ」を一段階粗くできます。
それに合わせて篩の「目」も同じ粒度にする必要があります。つまり m += b を m += 2*b としないといけません。
そのようにしてこのsuccメソッドを書き直すと、以下のようになります。

def succ()
  if (@n.nil?)
    @n = 1
    return 2
  end
  n = @n
  while true
    m = n += 2
    if (b = @h[n])
      while(@h[m += 2*b]); end
      @h[m] = b
      @h.delete(n)
    else
      while(@h[m += 2*n]); end
      return @n = @h[m] = n
    end
  end
end

見た目もコンパクト。スピードもこれだけで倍近く速くなります。

さらに。
同じ数に同じ偶数を何回も足していくと、3回に1回は必ず3の倍数になります。当たり前ですが、3の倍数は絶対に素数ではありません。
これをうまく省く方法があります。次の数を求めるときに
n += 2 → n += 4 → n += 2 → n += 4 → …
のように、2、4、2、4、… を順に足して行けば良いのです。
倍数を求めるときも同様、
m += 2*p → m += 4*p → m += 2*p → m += 4*p → …
とすれば、うまく3の倍数の生成を防ぎながら次の倍数を求めていくことが出来ます。
これをうまく組み合わせたのが、前の記事のPrime::A2Generatorのsuccメソッドとなっています。

3の倍数をうまく省いたエラトステネスの篩

さて。
昨日の記事で紹介したprimeライブラリの内部クラスPrime::EratosthenesGeneratorを、もう少し掘り下げてみましょう。
昨日も引用した部分をもう一度引用します。

class Prime
  : #中略
  class EratosthenesGenerator < PseudoPrimeGenerator
    : #中略
    def succ
      @last_prime = @last_prime ? EratosthenesSieve.instance.next_to(@last_prime) : 2
    end
    : #中略
  end
  : #中略
end

中でさらにEratosthenesSieveというクラスを参照していることが見て取れます。
"EratosthenesSieve"。そう、つまり日本語訳すると「エラトステネスの篩」です。先の例は「部分的なエラトステネスの篩」だったのに対して、こちらは「ほぼ完全系のエラトステネスの篩」となっています。またエラトステネスの篩なので、素数の列挙はかなり高速です。

でもエラトステネスの篩って、求めたい素数の最大値を決めて、それまでの数のうち合成数を篩にかけて消していく、というイメージがありますよね。でもEratosthenesGeneratorは素数の無限列挙にも対応しています。どういう仕組みになっているのでしょう?
ということで、内容を見てみましょう*3

class Prime
  : #中略
  # Internal use. An implementation of eratosthenes's sieve
  class EratosthenesSieve
    include Singleton

    def initialize # :nodoc:
      # bitmap for odd prime numbers less than 256.
      # For an arbitrary odd number n, @table[i][j] is 1 when n is prime where i,j = n.divmod(32) .
      @table = [0xcb6e, 0x64b4, 0x129a, 0x816d, 0x4c32, 0x864a, 0x820d, 0x2196]
    end

    # returns the least odd prime number which is greater than +n+.
    def next_to(n)
      n = (n-1).div(2)*2+3 # the next odd number of given n
      i,j = n.divmod(32)
      loop do
        extend_table until @table.length > i
        if !@table[i].zero?
          (j...32).step(2) do |k|
            return 32*i+k if !@table[i][k.div(2)].zero?
          end
        end
        i += 1; j = 1
      end
    end

    private
    def extend_table
      orig_len = @table.length
      new_len = [orig_len**2, orig_len+256].min
      lbound = orig_len*32
      ubound = new_len*32
      @table.fill(0xFFFF, orig_len...new_len)
      (3..Integer(Math.sqrt(ubound))).step(2) do |p|
        i, j = p.divmod(32)
        next if @table[i][j.div(2)].zero?

        start = (lbound.div(2*p)*2+1)*p    # odd multiple of p which is greater than or equal to lbound
        (start...ubound).step(2*p) do |n|
          i, j = n.divmod(32)
          @table[i] &= 0xFFFF ^ (1<<(j.div(2)))
        end
      end
    end
  end
  : #中略
end

ポイントは extend_table メソッド。
つまり、次に検索する数が大きくなったら、適宜素数テーブル(篩)を拡大しているのです。これによって無限列挙に対応できているわけですね。
さらにSingletonクラスなので、そのインスタンス変数はずっと保持されます。つまり結果として、今までに見付かった素数はキャッシュされた状態になっています。Prime.eachを2回目に実行したときはこのキャッシュが使用されるので、だから1回目と比べてあんなに速かったというわけです。
あと、英語ですがコメントを読むと分かる通り、こちらは初めから偶数は除外して考えられています。これも高速化の一助となっています。

それじゃもう一工夫して、3の倍数も除外するように改造しちゃいましょう。
ということで、書いてみました。

使用例(というかベンチマーク)。

require 'prime'
require 'benchmark'
require_relative 'prime.a2'
require_relative 'prime.a2e'

n = ARGV.size == 1 ? ARGV[0].to_i : 100
prs = Prime.each(nil, Prime::A2Generator.new).take(n)
puts "The #{n}th prime is #{prs[n-1]}."

Benchmark.bmbm do |b|
  b.report("Prime.each.take(n)") {Prime.each.take(n)[-1]}
  b.report("Prime.each(a2).take(n)") {Prime.each(nil,Prime::A2Generator.new).take(n)[-1]}
  b.report("Prime.each(a2e).take(n)") {Prime.each(nil,Prime::A2EGenerator.new).take(n)[-1]}
end

↑の実行例↓。

$ ruby prime_a2e_benchmark.rb 50000
The 50000th prime is 611953.
Rehearsal -----------------------------------------------------------
Prime.each.take(n)        1.450000   0.010000   1.460000 (  1.446770)
Prime.each(a2).take(n)    0.290000   0.000000   0.290000 (  0.298195)
Prime.each(a2e).take(n)   0.350000   0.000000   0.350000 (  0.345475)
-------------------------------------------------- total: 2.100000sec

                              user     system      total        real
Prime.each.take(n)        0.250000   0.000000   0.250000 (  0.256042)
Prime.each(a2).take(n)    0.270000   0.000000   0.270000 (  0.263810)
Prime.each(a2e).take(n)   0.170000   0.000000   0.170000 (  0.169651)

初回はA2Generatorほどの数字は出ていませんが、それでも標準のEratosthenesGeneratorの約4倍速です。
そして2回目は、EratosthenesGeneratorと同様、今まで見付けた素数はすべてキャッシュされた状態になっていますが、その速さは約1.5倍速。かなり速い♪
また2回の実行時間を足すと、a2が0.56s、a2eが0.52s。3回以上実行すればa2eがa2に比べてどんどん小さくなります。長い目で見れば、こちらの方が実行時間が短くて済むのはこちら、ということが分かると思います。

しかも、メモリ使用量も割と抑えることが出来ています。

$ /usr/bin/time -l ruby -rprime -e "pr=0;Prime.each(600000){|num|pr=num};puts pr"
599999
        1.54 real         1.53 user         0.00 sys
   5111808  maximum resident set size
         0  average shared memory size
         0  average unshared data size
         0  average unshared stack size
      1263  page reclaims
         0  page faults
         0  swaps
         0  block input operations
         0  block output operations
         0  messages sent
         0  messages received
         0  signals received
         1  voluntary context switches
        57  involuntary context switches
$ /usr/bin/time -l ruby -r`pwd`/prime.a2 -e "pr=0;Prime.each(600000,Prime::A2Generator.new){|num|pr=num};puts pr"
599999
        0.30 real         0.28 user         0.00 sys
   7364608  maximum resident set size
         0  average shared memory size
         0  average unshared data size
         0  average unshared stack size
      1813  page reclaims
         0  page faults
         0  swaps
         9  block input operations
        10  block output operations
         0  messages sent
         0  messages received
         0  signals received
        10  voluntary context switches
         0  involuntary context switches
$ /usr/bin/time -l ruby -r`pwd`/prime.a2e -e "pr=0;Prime.each(600000,Prime::A2EGenerator.new){|num|pr=num};puts pr"
599999
        0.40 real         0.39 user         0.00 sys
   5177344  maximum resident set size
         0  average shared memory size
         0  average unshared data size
         0  average unshared stack size
      1279  page reclaims
         0  page faults
         0  swaps
         0  block input operations
         0  block output operations
         0  messages sent
         0  messages received
         0  signals received
         1  voluntary context switches
        18  involuntary context switches

maximum resident set size を見ると、(EratosthenesGenerator)<A2EGenerator<<A2Generator となっています。EratosthenesGeneratorほどメモリ容量を抑えられていないように見えるのは、prime.a2e.rbも内部でprime.rbを読み込んでいるからかもしれません。ただその差は(prime.a2.rbとの差に比べれば)微々たるモノだということが分かると思います。

なかなか良い結果が得られたので、今後個人的に素数列挙するときは、この Prime::A2EGenerator を利用しようと思います♪

*1:記事がRuby1.8全盛の時代のもので、Ruby1.9も出てきたので一応確認してみている、という感じになっているですが、このクラスはRuby1.9のprimeライブラリとは互換性がありません(共存できません)。

*2:ソースを転載する代わりにロジックを日本語で書き下してみました。

*3:実はこのコードはRuby1.9.1のprime.rbの内容です。1.9.3では少し内容が変わっています。