Squeak/Pharo次世代高速CogVMで、Project Euler 250に挑戦してみる

OOPL


間があいたわりに似たようなネタが続きますが、き　に　し　な　い。
あと念のため Project Euler にまじめにとりくんでおられる方にはネタバレを含むので以下ご注意ください。


夜中からの Haskell 談義 - Togetter 　より、下のような記事があるのを知って、本題の Haskell とはいっさい関係なくて恐縮ですが、Squeak/Pharo の次世代高速 VM である CogVM でこの手の処理を実行した場合どの程度の効果があるのかな、と試してみました。

もしよろしければhttp://homepage1.nifty.com/herumi/diary/1011.html#10をどうやってHaskellで書いたら綺麗で速いコードになるのか見せていただけますでしょうか．

Twitter / herumi

Project Euler 250
最初見たときは部分集合って2250250個だから無理だろっと思ったけど,やってみたら意外と簡単だった. C++で0.66秒.

http://homepage1.nifty.com/herumi/diary/1011.html#10


ちなみに、掲載されている C++ のコードを手元の環境（MacBook Air(Late 2010)、1.6GHz Core 2 Duo、Cygwin）で実行してみたところ 3.56秒という結果でした。


まず、はずかしながら元のコードを見ても何をやっているかさっぱりだったので、何も考えずそのまま C++ を機械的に Smalltalk に変換したのがこちら。何も考えず―とかいいつつ、Smalltalk の配列は１オリジンなのでちょいちょいいじっていますが―。

| time ans |
time := [
   | n16 N M vec num ipowMod addm |

   n16 := 10000000000000000.
   N := 250250.
   M := 250.
   vec := IntegerArray new: N.
   num := Array new: M withAll: 0.

   ipowMod := [:x :a :m |
      | t out |
      a = 0 ifTrue: [1] ifFalse: [
      a = 1 ifTrue: [x \\ m] ifFalse: [
      t := x.
      out := 1.
      [a > 0] whileTrue: [
         (a bitAnd: 1) = 1 ifTrue: [out := (out * t) \\ m].
         t := (t * t) \\ m.
         a := a >> 1.
      ].
      out]]
   ].

   addm := [:a :b :c |
      | s |
      s := a + b.
      s > c ifTrue: [s - c] ifFalse: [s]
   ].

   1 to: N do: [:i |
      vec at: i put: (ipowMod value: i value: i value: M).
   ].

   num at: M put: 1.
   1 to: N do: [:i |
      | v ad |
      v := vec at: i.
      ad := num copy.
      1 to: M do: [:j |
         | p |
         p := addm value: v value: j value: M.
         num at: p put: (addm value: (ad at: p) value: (ad at: j) value: n16).
      ]
   ].
   ans := (num at: M) - 1.
] timeToRun.
World findATranscript: nil.
Transcript cr; show: ('ans, {1}; time, {2} ms' format: {ans. time})


これを、ついさきごろ公開された Squeak4.2日本語β版でそれぞれの VM で起動し、ワークスペースを開いてそこにコピペして実行（選択後、右クリックメニューから「式を評価 (d)」、もしくは alt+a → alt+d）したところ、旧来のノーマルVM では 210秒ほどかかったのに対して CogVM では 40秒程度と 1/5 の時間で処理を終えられました。すばらしいですね。もちろん C++には遠く及びませんが、10倍程度遅いだけなので許せる範囲です。というか優秀な部類だと思います。

normal ans, 1425480602091519; time, 208176 ms
cogvm  ans, 1425480602091519; time,  41287 ms


ついでに毎度ダシにして申し訳ないのですが比較のため（なにせ Smalltalk からの変換が簡単なもので―）、同じコードを Ruby に直訳して実行したところ 1.8.7 で 315秒ほど、1.9.2 で 95秒ほどという結果になりました。

$ cat euler250.rb
n16 = 10000000000000000
N = 250250
M = 250
num = Array.new(M, 0)

def ipowMod(x, a, m)
  return 1 if a == 0
  return x % m if a == 1
  t = x
  out = 1
  while a > 0
    out = (out * t) % m if a & 1 == 1
    t = (t * t) % m
    a >>= 1
  end
  out
end

def addm(a, b, c)
  s = a + b
  s >= c ? s - c : s
end

num[0] = 1
(1..N).each{ |i|
  v = ipowMod(i, i, M)
  ad = num.dup
  (0...M).each{ |j|
    p = addm(v, j, M)
    num[p] = addm(ad[p], ad[j], n16)
  }
}

p num[0] - 1

$ time ruby -v euler250.rb
ruby 1.8.7 (2011-02-18 patchlevel 334) [i386-cygwin]
1425480602091519

real    5m15.360s
user    5m13.327s
sys     0m0.046s
$ time ruby -v euler250.rb
ruby 1.9.2p0 (2010-08-18 revision 29036) [i386-cygwin]
1425480602091519

real    1m34.505s
user    1m33.600s
sys     0m0.124s

そんなこんなしているうちに処理内容についての理解もそこそこ深まったので、極端な速度低下をおこさない範囲で前述の C++ からの直訳コードを（前後の比較をしやすくするため変数名などはそのままにしつつも―）もうすこし Smalltalk っぽく書き直すこともしてみたのがこちらです。

| time ans |
time := [
   | n16 N M num |
   n16 := 1e16.
   N := 250250.
   M := 250.
   num := Array new: M withAll: 0.
   num at: M put: 1.
   num := (1 to: N) inject: num into: [:ad :i |
      | v |
      v := i raisedTo: i modulo: M.
      (ad flipRotated: v * -2) + ad \\ n16.
   ].
   ans := num last - 1.
] timeToRun.

World findATranscript: nil.
Transcript cr; show: ('ans, {1}; time, {2} ms' format: {ans. time})


実行時間は 1.5倍ほどかかるようになってしまいましたが、CogVM の速度なら許容の範囲内だと思いますし、それより当初のものに比べてかなり簡潔になり、何をやっているかの見通しも（一部 Squeak 独特の APL っぽい謎な配列演算を活用しているところを除けば―）そこそこマシになったような気がします。

ちなみに #filpRotated: は引数が偶数のときにはそのまま、奇数のときには反転してから、(引数+1)//2 分だけローテートする謎メソッド（ｗ。roteteでセレクタ検索 alt+shift+w していて発見）で、#raisedTo:modulo: は、ipowMod の組み込み版（前に見たことあるなーってIntegerを探したら見つかった）です。参考まで。