Hatena::ブログ(Diary)

Smalltalkのtは小文字です

id:sumim:about

Smalltalk を本格的に勉強する気はないけれど、うんちく程度に知っておきたいなら→Smalltalkをちょっとかじってみたい人のための、チュートリアルまとめ - Qiita

オブジェクト指向の“モヤッと”の正体を知りたくなったらこちらの記事が役に立つかも→id:sumim:20080415:p1 とか id:sumim:20040525:p1


 | 

2006-06-19

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


関連エントリー:



素数を列挙するにあたって、エラトステネスのふるいは比較的高速で有効です。たとえば、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 で書かれているにもかかわらず1秒とかかりません。

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

比較のために、手元の同じ環境で さんの PerlRuby のコードを走らせたときは、こんな結果が出ています。私の古〜い PowerBook は、小飼さんの最新式の環境(MacBook Pro)の3〜4割程度のパワーしかないこともわかります。w 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)までの素数を得るための“ふるい”には、普通、4ギガバイトほどの容量を必要とするが、この方法で圧縮すると 27 メガバイトほどで足りるのだそうです。すげー。


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

換言すれば、8ビット(8個のフラグ)、つまり1バイトあれば、30個分の整数の素数判定状況を保持することができるわけです。すげー。


なお、くだんのページで公開されている 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(" ")

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

$ /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 は詳しくないので感じるだけw)。


追記:

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

$ /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

関連:

はらはら 2006/06/21 11:32 スピードだけなら、こちらの方がかなり速いみたいです。
[ruby-list:16645] (http://blade.nagaokaut.ac.jp/cgi-bin/scat.rb/ruby/ruby-list/16645)

nursenurse 2006/06/22 18:15 うーむ、そこまで速くなるのですか。ループの内側が簡素なのもかなり効いてるようですねぇ。

nursenurse 2006/06/25 01:25 たしかSmalltalkの実装って、内部でバイトコードに変換しているんじゃないですっけ。なので、Rubyのコードを呼び出しているが故の速度の低下をもともと受けないのでは、と予想。log.ttsky.net による YARV によって実行した例で、わたしやsuminさんのコードが大幅に速くなっているのも、その辺が原因ではないかなぁと。

スパム対策のためのダミーです。もし見えても何も入力しないでください
ゲスト

コメントを書くには、なぞなぞ認証に回答する必要があります。

 | 
2004 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 |
2005 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 |
2006 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 |
2007 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 |
2008 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 12 |
2009 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 |
2010 | 01 | 02 | 03 | 04 | 06 | 07 | 08 | 09 | 10 | 11 | 12 |
2011 | 01 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 10 | 11 | 12 |
2012 | 01 | 02 | 03 | 04 | 05 | 06 | 08 | 10 | 12 |
2013 | 02 | 03 | 04 | 05 | 06 | 07 | 08 | 09 | 11 | 12 |
2014 | 01 | 02 | 05 | 07 | 08 | 09 | 10 | 11 |
2015 | 04 | 07 | 08 | 11 | 12 |
2016 | 02 | 03 | 06 | 07 | 08 |
2017 | 04 | 05 | 06 |

最近のコメント

1. 06/28 sumim
2. 06/28 squeaker
3. 06/28 sumim
4. 06/28 squeaker
5. 06/25 sumim

最近のトラックバック

1. 05/25 プラグインレスでSVGを表示する「SIE」開発ブログ - メッセージをや...
2. 01/30 no_orz_no_life - Erlangとジャンケン
3. 12/31 檜山正幸のキマイラ飼育記 - J言語、だってぇー?
4. 09/04 Twitter / @atsushifx
5. 07/06 みねこあ - オブジェクト指向 と FizzBuzz

この日記のはてなブックマーク数
1734667
Connection: close