MeCab Ruby と CaboCha Ruby をインストールしてみた
機械学習も統計もあまり勉強進んでいませんゴメンナサイ。
ちょっとそのきっかけ作り（という名目の個人的な実験）のために、形態素解析とか日本語係り受け解析とかを試してみることにしました。
個人的な理由のために Ruby からの利用を前提として。
MeCab と MeCab Ruby のインストール
最初は「形態素解析」というものに触れてみようということで、 MeCab をインストール。
技術ブログでは「失敗したログも含めて全部記録する」のが定石なんですが、めんどくさいし個人的には成功例と「なぜそれで成功したのか（なぜ失敗を回避できたか）」が分かれば良いと思うのでその方針で書きます。
ちなみに環境は、Mac OSX Lion 10.7.4 です。
インストールしたのは以下の通り：
まずは MeCab 本体。最新のtar ballをDLして以下の要領でインストール。
$ tar xzvf mecab-0.994.tar.gz $ cd mecab-0.994 $ ./configure --enable-utf8-only $ make $ sudo make install
ポイントは以下の通り：
./configure時に「--enable-utf8-only」オプションを付加。Ruby 1.9 で UTF-8 での使用しか想定していないので。またこうするとeucやsjisの変換テーブルを埋め込まず実行バイナリを小さく出来ると言うことなので*1。
これはすんなり完了。
続いて辞書のインストール。
$ tar xzvf mecab-ipadic-2.7.0-20070801.tar.gz $ cd mecab-ipadic-2.7.0-20070801 $ ./configure --with-charset=utf8 $ make $ sudo make install
ポイント：
こちらもすんなり。なんだMeCabてばインストール簡単じゃん。
早速実験。sumomo.txt は utf-8 のテキストファイルです。
$ cat sumomo.txt すもももももももものうち $ mecab sumomo.txt すもも 名詞,一般,*,*,*,*,すもも,スモモ,スモモ も 助詞,係助詞,*,*,*,*,も,モ,モ もも 名詞,一般,*,*,*,*,もも,モモ,モモ も 助詞,係助詞,*,*,*,*,も,モ,モ もも 名詞,一般,*,*,*,*,もも,モモ,モモ の 助詞,連体化,*,*,*,*,の,ノ,ノ うち 名詞,非自立,副詞可能,*,*,*,うち,ウチ,ウチ EOS
OK。
では続いて Ruby で直接いじれるようにするためのモジュールをインストール。Downloads - mecab - Japanese morphological analyzer - Google Project Hostingから最新のモジュールのtar ballをDLして、以下の要領でインストール。
$ tar xzvf mecab-ruby-0.994.tar.gz $ cd mecab-ruby-0.994 $ ruby extconf.rb $ make $ sudo make install $ ruby test.rb 《実行結果省略》
こちらは何も難しいことはないですね。
当然のことながら、今実行している ruby に新しいモジュールとしてインストールされるので、rvmやrbenv等で複数バージョンを管理しているなら、必要なバージョンのRubyごとにインストールしないといけません。rvmなら、専用のgemsetを作ってからインストール、っていう方法も良いかも。私はめんどくさいので素の（普段使用している）Ruby だけにインストールしました。
簡単に実行確認。
[1] 1.9.3-p125(main)> require 'MeCab' => true [2] 1.9.3-p125(main)> sentence = "\u592a\u90ce\u306f\u3053\u306e\u672c\u3092\u4e8c\u90ce\u3092\u898b\u305f\u5973\u6027\u306b\u6e21\u3057\u305f\u3002" => "太郎はこの本を二郎を見た女性に渡した。" [3] 1.9.3-p125(main)> mecab = MeCab::Tagger.new => #<MeCab::Tagger:0x007f885c18ae20 @__swigtype__="_p_MeCab__Tagger"> [4] 1.9.3-p125(main)> puts mecab.parse sentence 太郎 名詞,固有名詞,人名,名,*,*,太郎,タロウ,タロー は 助詞,係助詞,*,*,*,*,は,ハ,ワ この 連体詞,*,*,*,*,*,この,コノ,コノ 本 名詞,一般,*,*,*,*,本,ホン,ホン を 助詞,格助詞,一般,*,*,*,を,ヲ,ヲ 二 名詞,数,*,*,*,*,二,ニ,ニ 郎 名詞,一般,*,*,*,*,郎,ロウ,ロー を 助詞,格助詞,一般,*,*,*,を,ヲ,ヲ 見 動詞,自立,*,*,一段,連用形,見る,ミ,ミ た 助動詞,*,*,*,特殊・タ,基本形,た,タ,タ 女性 名詞,一般,*,*,*,*,女性,ジョセイ,ジョセイ に 助詞,格助詞,一般,*,*,*,に,ニ,ニ 渡し 動詞,自立,*,*,五段・サ行,連用形,渡す,ワタシ,ワタシ た 助動詞,*,*,*,特殊・タ,基本形,た,タ,タ 。 記号,句点,*,*,*,*,。,。,。 EOS => nil
内容については、ここでは深く突っ込みません。
CaboCha のインストール
つづいて CaboCha 。こちらは「日本語係り受け解析器」。
実は MeCab の存在は以前から知っていて、それでちょっと試したんですが、「結局これ、単語（正確には形態素）の分割は出来るけれど、これだけだと係り受けは分かんないじゃん」てことに気付いて（私がやりたかったことがそっちだったことに気付いて）、後から探して見付けてインストールしたのが CaboCha 、だったりします。
インストールしたのは以下の通り：
なお他に、MeCab がインストールされていることが前提です。
まずは CRF++。
$ tar xzvf CRF++-0.57.tar.gz $ cd CRF++-0.57 $ ./configure $ make $ sudo make install
インストールそのものは問題なし。一応、ポイント：
- CaboCha のインストールをまとめたサイトで、よく「CRF++ を一番最初に（MeCab より前に）インストールする」、と記述があります（ 日本語係り受け解析器 CaboCha Ruby 拡張の基本的な使い方とちょっとした応用 - デー等）が、CaboChaの動作に必要なモノなので、今回のように MeCab の後 CaboCha の前でも大丈夫なようです。
続いて CaboCha。MeCab のインストールの時に言及した通り、成功例だけ示します。
$ tar xzvf cabocha-0.64.tar.gz $ cd cabocha-0.64 $ LIBS=-liconv ./configure --with-charset=utf8 $ make $ sudo make install
ポイントは以下の通り：
./configure時に「--with-charset=utf8」オプションを付加。Ruby 1.9 で UTF-8 での使用しか想定していないので。- そのまま ./configure すると、iconv 周りで「
ld: symbol(s) not found for architecture x86_64」というエラーになって落ちてしまいます（実際落ちました）。対処法として調べてみたところ：- 「homebrewでlibiconvをインストール（Mac OS X 10.7 (Lion)へのcabochaのインストール - 犬ぶよツールズ制作記録等）」→homebrew 以外でインストールしているパッケージが色々あったので移行がめんどくさかったので今回は不採用。
- 「LIBS=-liconv（MacOSX10.6 にCaboCha0.60pre4をインストールとPythonバインディング - Do not forget this.）」→つまりliconvを無効に。手っ取り早いのでこちらを採用。今のところ動作に問題なし。
取り敢えず動作確認。
$ cat sumomo.txt
すもももももももものうち
$ cabocha sumomo.txt
すももも-D
ももも---D
ももの-D
うち
EOS
…例が悪かったですねｗ これじゃ何が起きているかよく分からない。別の例。
$ cabocha
私が「あなたは右の道から来ましたか？」と質問したら、あなたは「はい」と答えますか？
私が---------D
「あなたは-----D |
右の-D | |
道から-D |
来ましたか？」と-D
質問したら、-----D
あなたは---D
「はい」と-D
答えますか？
EOS
複文、係り受け、きちんと解析できてますね。素晴らしい(^-^)
CaboCha Ruby のインストール
続いて CaboCha Ruby。
MeCab には Ruby で直接利用するためのラッパーモジュールが用意されていましたが、CaboCha は自分で用意する必要があります。といってもSWIG (Simplified Wrapper and Interface Generator)というツールを使って簡単にラッパーモジュールを作れるようになっている*2ので、ありがたくそれを利用します。
ということで、さらに以下をインストール：
$ tar xjvf pcre-8.31.tar.bz2 $ cd pcre-8.31 $ ./configure $ make $ sudo make install
$ tar xzvf swig-2.0.7.tar.gz $ cd swig-2.0.7 $ ./configure $ make $ sudo make install
ポイント：
- そのまま SWIG をインストールしようとすると、./configure で pcre がないと怒られたので、PCRE の最新を探してきて先にインストール。予め何かの目的で PCRE を入れてあれば不要かも。
さて、では CaboCha の Ruby 用のモジュールを作ってインストール。
$ cd cabocha-0.64/swig $ make ruby $ cd ../ruby/ $ ruby extconf.rb $ make $ sudo make install
ポイント：
- SWIG が入ってないと、
make rubyの時点で「swig: No such file or directory」のエラーになります。 ruby extconf.rb以降は、MaCab Ruby と同じ手順ですね。
さてテストですが、用意されている test.rb が Ruby 1.9 だと 正常に動作しません。文字コード指定のマジックコメントが必要です。
$ vi test.rb #!/usr/bin/ruby # -*- coding: utf-8 -*- # ↑この1行を追記 require 'CaboCha' sentence = "太郎はこの本を二郎を見た女性に渡した。" 《以下略》
で、test.rb 実行。
$ ruby test.rb <PERSON>太郎</PERSON>は-----------D この-D | 本を---D | 二郎を-D | 見た-D | 女性に-D 渡した。 EOS <PERSON>太郎</PERSON>は-----------D この-D | 本を---D | 二郎を-D | 見た-D | 女性に-D 渡した。 EOS * 0 6D 0/1 2.909358 太郎 名詞,固有名詞,人名,名,*,*,太郎,タロウ,タロー B-PERSON は 助詞,係助詞,*,*,*,*,は,ハ,ワ O * 1 2D 0/0 1.257926 この 連体詞,*,*,*,*,*,この,コノ,コノ O * 2 4D 0/1 0.638994 本 名詞,一般,*,*,*,*,本,ホン,ホン O を 助詞,格助詞,一般,*,*,*,を,ヲ,ヲ O * 3 4D 1/2 1.696047 二 名詞,数,*,*,*,*,二,ニ,ニ O 郎 名詞,一般,*,*,*,*,郎,ロウ,ロー O を 助詞,格助詞,一般,*,*,*,を,ヲ,ヲ O * 4 5D 0/1 0.000000 見 動詞,自立,*,*,一段,連用形,見る,ミ,ミ O た 助動詞,*,*,*,特殊・タ,基本形,た,タ,タ O * 5 6D 0/1 0.000000 女性 名詞,一般,*,*,*,*,女性,ジョセイ,ジョセイ O に 助詞,格助詞,一般,*,*,*,に,ニ,ニ O * 6 -1D 0/1 0.000000 渡し 動詞,自立,*,*,五段・サ行,連用形,渡す,ワタシ,ワタシ O た 助動詞,*,*,*,特殊・タ,基本形,た,タ,タ O 。 記号,句点,*,*,*,*,。,。,。 O EOS
Ruby の実行環境上で実験。
[1] 1.9.3-p125(main)> require 'CaboCha' => true [2] 1.9.3-p125(main)> sentence = "\u592a\u90ce\u306f\u3053\u306e\u672c\u3092\u4e8c\u90ce\u3092\u898b\u305f\u5973\u6027\u306b\u6e21\u3057\u305f\u3002" => "太郎はこの本を二郎を見た女性に渡した。" [3] 1.9.3-p125(main)> cabocha = CaboCha::Parser.new => #<CaboCha::Parser:0x007f9b690d1fe0 @__swigtype__="_p_CaboCha__Parser"> => true [4] 1.9.3-p125(main)> tree = cabocha.parse sentence => #<CaboCha::Tree:0x007f9b6990d948 @__swigtype__="_p_CaboCha__Tree"> [5] 1.9.3-p125(main)> puts tree.toString(CaboCha::OUTPUT_RAW_SENTENCE) <PERSON>太郎</PERSON>は-----------D この-D | 本を---D | 二郎を-D | 見た-D | 女性に-D 渡した。 EOS => nil [6] 1.9.3-p125(main)> chunk_to_s = ->(ch){(0...ch.token_size).map{|j|tree.token(ch.token_pos+j).normalized_surface}.join("-")} => #<Proc:0x007f9b699611b0@(pry):38 (lambda)> [7] 1.9.3-p125(main)> chunk_link = ->(ch){if ch.link>=0;ch_from=chunk_to_s[ch];ch_to=chunk_to_s[tree.chunk(ch.link)];"#{ch_from} => #{ch_to}";end;} => #<Proc:0x007f9b6903e0b0@(pry):41 (lambda)> [8] 1.9.3-p125(main)> puts chunk_link[tree.chunk(0)] => 太郎-は => 渡し-た-。 [9] 1.9.3-p125(main)> puts chunk_link[tree.chunk(1)] => この => 本-を
つまり、treeは文の解析ツリー、chunkが文節、tokenが単語（形態素）、ですね。
そしてCaboCha::Chunk#link が係り受けを表している（自分が「係り」、linkで表されるインデックスのChunkが「受け」）、と。
でもこれ、本当に CaboCha を単純にラップしたもので、Ruby 的に使いにくいところとか気になるところがいくつかあるので、余裕があったらより Ruby っぽくなるようにさらにラップしたいところ。
取り敢えず今日はここまで。
汎用的で省メモリかつかなり速い素数無限列挙
台風凄かったですね。私の住んでいる地域では、夕方頃に風雨ともかなり強かったのですが、今はもうだいぶ落ち着いてきちゃいました。
さてさて。昨日の記事の続きです。
まずは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点。
- エラトステネスの篩をマジメに全部やろうとしている。
- 偶数（そして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 switchesmaximum resident set size を見ると、(EratosthenesGenerator)<A2EGenerator<<A2Generator となっています。EratosthenesGeneratorほどメモリ容量を抑えられていないように見えるのは、prime.a2e.rbも内部でprime.rbを読み込んでいるからかもしれません。ただその差は（prime.a2.rbとの差に比べれば）微々たるモノだということが分かると思います。
なかなか良い結果が得られたので、今後個人的に素数列挙するときは、この Prime::A2EGenerator を利用しようと思います♪
汎用的でけっこう速い素数無限列挙
二週間ぶりです。風邪引いてました。二週間かかってようやく落ち着きました。この時期の夏風邪ってば以下略。
さてさて、今日から暫く「素数列挙」の話をします。
Ruby1.9 の素数列挙メソッド
Ruby1.9には、素数列挙・素数判定・素因数分解を一手に扱う素数ライブラリが用意されています。
このブログでも素数列挙する時とかに何の説明もなく使ってました。
require 'prime' Prime.each {|pr|puts pr} # => 2 # => 3 # => 5 # => 7 # => 11 # => …（以降無限に列挙）
これ、割と速いです。というか結構速いです。
python+素数列挙でググるとよく出てくる、こんなサンプル
from itertools import ifilter, count def prime_generator(): g = count(2) while True: prime = g.next() yield prime g = ifilter(lambda x, prime=prime: x%prime, g) for pr in prime_generator(): print pr # => 2 # => 3 # => 5 # => 7 # => 11 # => …（以降無限に列挙）
よりもよっぽど速いです*1。その理由もちゃんとあるのですが、ここでは触れる程度にしておきます。
結構速いのですが、色々調べてみると、もっと速くできるらしいことが分かります。
特に「指定した数までの素数を列挙する」ことに特化すれば、かなり速くできます。例えば
エラトステネスのふるいをコンパクトにする方法 - Smalltalkのtは小文字ですとか、[ruby-list:16645] (http://blade.nagaokaut.ac.jp/cgi-bin/scat.rb/ruby/ruby-list/16645)とかを参照のこと。
でも、今回は「無限列挙」。しかも、Ruby1.9のprimeライブラリよりも高速なもの。しかもできれば、Prime#eachによる列挙と互換性を持たせたい*2。
Ruby1.9 の素数列挙メソッドの仕組み
Rubyのライブラリリファレンスを見れば分かるのですが、RubyのprimeライブラリのPrime#eachメソッドには、第2引数もあります。
$RUBY_HOME/lib/ruby/1.9.1/prime.rb の中を覗いてみると、こんな感じになっています（かなり抜粋）。
class Prime : #中略 def each(ubound = nil, generator = EratosthenesGenerator.new, &block) generator.upper_bound = ubound generator.each(&block) end : #中略 end
EratosthenesGeneratorというのが肝。Primeクラスの内部クラスになっています。やはりprime.rbから定義を見てみます。
class Prime : #中略 class PseudoPrimeGenerator include Enumerable : #中略 def succ raise NotImplementedError, "need to define `succ'" end : #中略 end class EratosthenesGenerator < PseudoPrimeGenerator : #中略 def succ @last_prime = @last_prime ? EratosthenesSieve.instance.next_to(@last_prime) : 2 end : #中略 end : #中略 end
つまりPseudoPrimeGeneratorというクラスがあって、所謂抽象クラス(Abstract Class)で、EratosthenesGeneratorはそれを継承したクラスで、succメソッドを実装してやっている（そして上の抜粋コードでは省略しましたが、PseudoPrimeGenerator#eachメソッドはPseudoPrimeGenerator#succメソッドを利用して実装されている）、という仕組みになっています。
なので。
オリジナルの素数ジェネレータをPrime::PseudoPrimeGeneratorクラスを継承して作って*3*4、そのインスタンスをPrime#eachメソッドの第2引数に渡せば、独自実装による素数列挙であってしかもprimeライブラリのPrimeクラスを拡張した（＝同じAPI仕様で動く）ものができる、というわけですね。
五月雨
ということで、そのようなクラスを1つ作ってみました。
メモ: 上限を決めずに素数を生成の高速化(Ruby版)を参考に、それを改良してさらに高速かつコンパクトに。
使用例。
require 'prime' require 'benchmark' require_relative 'prime.a2' max = ARGV.size == 1 ? ARGV[0].to_i : 100 pr = 0 Prime.each(max, Prime::A2Generator.new) {|num| pr = num} puts "Max prime num(<=#{max}) is #{pr}." block = Proc.new {} Benchmark.bmbm do |b| b.report("Prime.each(max)") {Prime.each(max, &block)} b.report("Prime.each(max,a2)") {Prime.each(max,Prime::A2Generator.new,&block)} end
↑の実行例↓。
$ ruby prime_a2_benchmark.rb 600000
Max prime num(<=600000) is 599999.
Rehearsal ------------------------------------------------------
Prime.each(max) 1.680000 0.010000 1.690000 ( 1.688865)
Prime.each(max,a2) 0.340000 0.000000 0.340000 ( 0.347410)
--------------------------------------------- total: 2.530000sec
user system total real
Prime.each(max) 0.330000 0.000000 0.330000 ( 0.331884)
Prime.each(max,a2) 0.340000 0.000000 0.340000 ( 0.341556)うん、かなり速いです。少なくとも1回目は約5倍速です。
あれ、でも2回目はほとんど同じ数字をたたき出しています。というか数字の上では微妙に標準のPrime#eachの方が速いですね。これはなぜでしょう？
実は標準で使用されるPrime::EratosthenesGeneratorは、ある仕組みで生成した素数をキャッシュしているのです。2回目以降に呼ばれたときにはそのキャッシュを利用するので1回目よりもかなり速くなるのです。その結果、今回作ったPrime::A2Generatorと同じくらい速くなってしまっているのです。
じゃ、A2Generatorもキャッシュを持たせれば速くなる？
それも一つのアプローチですが、ここではそっちには進まず、違うアプローチに進みたいと思います。
キーワードは、"EratosthenesSieve"。
ということで今回はここまで。続きは次回に。
*1:なぜRubyの話なのにpythonの例を出したのかというと、この例そのままをRubyで同じような動きをするモノを書くのが意外と面倒だったためです。言語レベルで「どっちが優れているか」という話ではないことを言及しておきます。
*2:最後のはもともと、著者によるワガママですけれど。
*3:先の抜粋コードで省略してしまっているのですが、このときsuccメソッドだけでなく、next、rewindの2つのメソッドも実装する必要があります。
*4:というかRubyなので、upper_bound=メソッドとeachメソッドの2つを実装したクラスであればPrime::PseudoPrimeGeneratorクラスを継承したクラスでなくても実はOKです。
セクシー素数
なにやら数学嫌いでも興味を持たずにはいられない数学用語ランキング - 学びランキング
- goo ランキングで1位だということで。
記念に「セクシー素数（http://ja.wikipedia.org/wiki/セクシー素数）」を列挙するメソッドをRuby*1で書いてみました。
実行例↓
#セクシー素数の組（最初の10組） Prime.sexy_pair.take(10) # => [[5, 11], [7, 13], [11, 17], [13, 19], [17, 23], [23, 29], [31, 37], [37, 43], [41, 47], [47, 53]] #セクシー素数の組（＜500） Prime.sexy_pair(500).to_a # => [[5, 11], [7, 13], [11, 17], [13, 19], [17, 23], [23, 29], [31, 37], [37, 43], [41, 47], [47, 53], [53, 59], [61, 67], [67, 73], [73, 79], [83, 89], [97, 103], [101, 107], [103, 109], [107, 113], [131, 137], [151, 157], [157, 163], [167, 173], [173, 179], [191, 197], [193, 199], [223, 229], [227, 233], [233, 239], [251, 257], [257, 263], [263, 269], [271, 277], [277, 283], [307, 313], [311, 317], [331, 337], [347, 353], [353, 359], [367, 373], [373, 379], [383, 389], [433, 439], [443, 449], [457, 463], [461, 467]] #セクシー素数の三つ組（＜1000） Prime.sexy_triplets(1000).to_a # => [[7, 13, 19], [17, 23, 29], [31, 37, 43], [47, 53, 59], [67, 73, 79], [97, 103, 109], [101, 107, 113], [151, 157, 163], [167, 173, 179], [227, 233, 239], [257, 263, 269], [271, 277, 283], [347, 353, 359], [367, 373, 379], [557, 563, 569], [587, 593, 599], [607, 613, 619], [647, 653, 659], [727, 733, 739], [941, 947, 953], [971, 977, 983]] #セクシー素数の四つ組（＜1000） Prime.sexy_quadruplets(1000).to_a # => [[5, 11, 17, 23], [11, 17, 23, 29], [41, 47, 53, 59], [61, 67, 73, 79], [251, 257, 263, 269], [601, 607, 613, 619], [641, 647, 653, 659]]
フィボナッチ数
今日は、みんな大好きフィボナッチ数列。
柏餅おいしかったとかちまき上品な甘さだったとか実家のタケノコごはんがとかそういうGW後半ネタはすぱっと省略。
フィボナッチ数列。
「1, 1,」から始まって、それ以降は「前の2つの数の和」が続く数列、のこと。数式（漸化式）による定義は以下の通り（0番目を0とした定義）。
最初の方を列挙すると以下の通り。
(0, )1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, …
できるだけ短いコードで列挙する、という所謂CodeGolfのネタにもなっています。
私の最高記録はRubyで27バイト…このあたりまでは誰でも行き着くんですよね、あと4バイトどうやって短くするんだろう(+_+)
でも今日は「短いコード」の話ではなく、一般項の話。つまり、
「整数nを与えたときに、フィボナッチ数列のn番目の数（＝n番目のフィボナッチ数）を返す関数」
を実装する話。
それをアルゴリズムの話だけでなく、もうちょっと数学のスパイスを加えた話をしてみます。
手続き指向の解答例
プログラムを少しでも組んだことのある人なら、きっと誰でも思いつくのは、「定義に従って順番に計算して答えを出す方法」。
Rubyで書くと、こんな感じ。
def fib n return n.even? ? -fib(-n) : fib(-n) if n < 0 # for negative n a, b = 0, 1 n.times { a = b + b = a } a end
関数*1内1行目は、n<0だった場合の対応。今回はあまり気にしないでください。本質はここじゃないので。
さて、もちろんこれは正しいフィボナッチ数を返しますが、効率の方はどうでしょう？
ベンチマークを取ってみました。
n=50,100,500,1000,5000,10000 の場合の各計算を1000回ずつ実施してその実行時間を計測（Rubyの benchmark ライブラリを使用）
user system total real fib(50)x1000 0.010000 0.000000 0.010000 ( 0.006343) fib(100)x1000 0.010000 0.000000 0.010000 ( 0.011449) fib(500)x1000 0.200000 0.000000 0.200000 ( 0.198729) fib(1000)x1000 0.480000 0.000000 0.480000 ( 0.475711) fib(5000)x1000 3.460000 0.160000 3.620000 ( 3.633362) fib(10000)x1000 8.080000 0.860000 8.940000 ( 8.946042)
うん、ほぼ当然のごとく、nにだいたい比例した数値が出ていますね。つまり計算量はO(n)です。
十分小さなnに対しては問題ありませんが、大きくなるに従ってどんどん時間がかかります。
関数指向の解答例
続いて、少しプログラムやアルゴリズムを勉強して「再帰関数」を覚えた人、もしくは「名古屋人なら関数型言語」な人なら、こんなコードを思いつくでしょう。
def fib n return n.even? ? -fib(-n) : fib(-n) if n < 0 # for negative n return n if n < 2 fib(n - 2) + fib(n - 1) end
つまり「定義をそのまま実装」した例です。
これを実行してみると、もちろん正しい答えは出るのですが…。
user system total real fib(10) 0.000000 0.000000 0.000000 ( 0.000037) fib(20) 0.000000 0.000000 0.000000 ( 0.004549) fib(30) 0.400000 0.000000 0.400000 ( 0.399396) fib(40) 48.120000 0.040000 48.160000 ( 48.225261) # fib(50) takes about 1:40' # fib(60) takes about 8.6days…
はい。非常に時間がかかります。上のベンチマーク、先ほどと違って1000回繰り返しているのではなく、1回だけしか実行していない処理時間であることに注意。
n=50で約1時間40分、n=60で約8日と14時間かかります。n=100,500,1000,… なんて、恐ろしくて実行できません。計算量はO(en)。
これは関数が、同じ引数を与えると毎回同じ処理をするから。
関数指向＋メモ化
同じ引数が与えられたら同じ結果を返すようにするために、前に計算した値を覚えておくという方法があります。これを「メモ化」と言います。
先ほどの例に、メモ化の要素を加えた例を示します。
def fib n return n.even? ? -fib(-n) : fib(-n) if n < 0 # for negative n fib_hash = Hash.new {|h, n| h[n] = case n when 0, 1 n else h[n - 2] + h[n - 1] end } fib_hash[n] end
Ruby特有の書き方になっているので関数型言語のメモ化とはちょっと違うし分かりにくい実装ですが、やってることは同じ。
「一度計算した（他の）フィボナッチ数は、ハッシュに記憶しておいて、後で再利用できるようにする」コードになっています。
これだけで、実行時間は劇的に短くなります。、が…。
user system total real fib(50)x1000 0.050000 0.000000 0.050000 ( 0.050571) fib(100)x1000 0.100000 0.000000 0.100000 ( 0.100457) fib(500)x1000 0.740000 0.000000 0.740000 ( 0.742134) fib(1000)x1000 1.760000 0.010000 1.770000 ( 1.774478) #fib(5000)x1000 #@> stack level too deep (SystemStackError)
n=1000 を超えたあたりになると、再帰呼び出しの深さが深くなりすぎて、コンピュータが処理できずにエラーになったりします。
それに、計算量そのものは同じO(n)なのですが、ベンチマークの結果を見てみると、最初の「手続き指向」よりもわずかながら時間がかかっています。
つまり関数型指向の実装は、必ずしも効率が良い方法ではない、ということなのでしょうか？
そうとも言えますし、実はそうとばかりは言えない部分もあったりします。
メモ化＋倍数公式
さぁ、ここからがこの記事のメインです。
フィボナッチ数は、その定義となる漸化式から、色々な公式が導き出せます。
そして余り知られていないのですが、2n、2n+1番目のフィボナッチ数について、こんな面白い公式があります。
これをFn=Fn-1+Fn-2の代わりに使ってみましょう。
def fib n return n.even? ? -fib(-n) : fib(-n) if n < 0 # for negative n fib_hash = Hash.new {|h, n| h[n] = case n when 0, 1 n else n.even? ? (h[n/2-1]*2 + h[n/2]) * h[n/2] : h[n/2]**2 + h[n/2+1]**2 end } fib_hash[n] end
ベンチマークは以下の通り。
user system total real fib(50)x1000 0.020000 0.000000 0.020000 ( 0.019347) fib(100)x1000 0.020000 0.000000 0.020000 ( 0.019735) fib(500)x1000 0.030000 0.000000 0.030000 ( 0.032989) fib(1000)x1000 0.050000 0.000000 0.050000 ( 0.043373) fib(5000)x1000 0.080000 0.000000 0.080000 ( 0.079305) fib(10000)x1000 0.130000 0.000000 0.130000 ( 0.135395)
n=10000でもSystemStackErrorが発生せず、しかもnが100を超えたあたりからは手続き指向の解答例よりも明らかに高速になっています。n=10000だと、「メモ化＋倍数公式」は「手続き指向」の約60倍速！
これは、今までの方法は結局「n番目以下の全てのフィボナッチ数も算出」していたのに対して、今回の方法が「計算に必要なフィボナッチ数だけを算出」しているから、というのが最も大きな理由。
算出する個数は、最初にn/2番目あたりを探し、次にそのさらにn/2あたりを探し…ということで、つまりlog2n に大体比例することになります。つまり計算量はO(logn)。かなり理想的です。
ただ、nが十分に小さい場合は2乗や掛け算の分だけ余分に計算することになるので、実時間は手続き指向の場合に敵わないこともある、というわけです。
おまけ：フィボナッチ数（＝フィボナッチ数列の一般項）の公式
n番目のフィボナッチ数は、nのみに依る式で以下のように表せます。
証明はめんどくさいしブログ1ページじゃHeavyなのでしませんが、興味のある方は、「数学ガール」（ISBN: 4797341378）でも読んでみてください。
これをうまく利用して、「途中のフィボナッチ数を一切計算せず直接算出する」方法ももちろん考えられます。
ただ、そのまま実数（浮動小数点数）を用いてやろうとすると誤差の問題が、それを無理矢理整数範囲で収めようとすると計算量の問題が、結局ついて回ります。
世の中うまくいかないものですね。
まとめ
計算量の考え方は大事。
あと、アルゴリズムの工夫も大事だけど、数学的な裏付けとか、「数学的な考察による数式変形（公式導出）によるさらなる効率化」が、とっても大事。
ちなみに「倍数公式」という呼び方は正式な用語かどうか分かりませんが、同様にF3nやF4nに対する公式も存在し（というか探せばいくらでも作り出すことができ）ますが、それらを利用してプログラムを書き換えてみても、必ずしもこれ以上効率的な結果にはなりません。
なぜでしょう？興味があればご自分で考察してみてください。
*1:Rubyの場合正式には「メソッド」と言いますが、言語はここでは重要ではないので一般的名称として「関数」という用語で通します。
コラッツの予想
勉強は明日からにして、今日はまたちょと趣味に走ります。
コラッツの予想（wikipedia:コラッツの問題）というものがあります。
任意の0より大きい整数 n をとり、
- n が偶数の場合、n を 2 で割る
- n が奇数の場合、n に 3 をかけて 1 を足す
文章で書くとこんな感じ（Wikipediaの記載をほぼそのまま引用）。簡単ですね。
確かめてみましょう。最初の数を 3 とすると、
3 → 10 → 5 → 16 → 8 → 4 → 2 → 1
はい、確かに最後に 1 になりました。
Rubyで整数からコラッツの数列を列挙するコードを書いてみました。
Integer#collatzメソッドは、ブロックを指定しなければEnumerable(Enumerator)を返すので注意。配列化するにはコメントにあるように.to_aを、列の長さを取得するには.countを使用してください。
てか 27 から始めると 111 手（最初の数 27 も含めて長さが 112 なので）もかかるんですね。意外と長いですね。
さて。なぜ始めから配列で返さず、Enumerable（言い換えると「ストリーム」）で返したのかというと。
「必ず有限回で 1 に到達する」かどうか、は、まだ証明できていないからです。
つまり、「ひょっとしたら無限大に発散するかもしれない」し、「ひょっとしたらどこかで永久ループに入っちゃって、抜け出せないかもしれない」のです。それで安全性を考えて、ストリームにしているのです。
実際、「3をかけて」「1を足す」の部分（ソース中の「n*3+1 # <- *1」の部分）を色々変えてみると、ループになったり、数がどんどん大きくなるパターンが出てきます。
例えば「n*5+1」にすると、9.collatz{|n|p n}が表示する数はどんどんどんどん大きくなっていき、永久に終わりません*1。
「n*3-1」にすると、5.collatz{|n|p n}はループします（5→14→7→20→10→5→…）*2。
本当にn*3+1のときは必ず 1 になるのか。
どういうパターンの時に、発散やループに入らずに1に到達するのか。
所謂「数学の未解決問題」の一つです。
もし。
上のコードを適当な整数から始めて、「発散」したり、「ループ」したりしたら。
『反例発見！』
とTwitterで叫びましょう。
きっと世界中の数学者（主に整数論とか代数学者）がなんらかの反応をしますｗ
【追記：2012/05/05】勉強のために、コラッツの数列を列挙するコード（関数）を Python でも書いてみました。
参考までに。
