プロジェクト・オイラー15問目

プロジェクト・オイラーの15問目。以下の2x2のグリッドを左上から右下になぞるとき、6つのルートがある。20x20のグリッドではルートは何通りか。

最初にRubyで書いてみた。確認のために経路を配列に残しつつ。

WIDTH, HEIGHT = 2,2
$count = 0

def step_forward (path, x, y)
  if (x < WIDTH && y < HEIGHT) then 
    tmp = path.dup
    tmp2 = path.dup
    step_forward(tmp << [x + 1, y], x + 1, y)
    step_forward(tmp2 << [x, y + 1], x, y + 1)
  elsif (x < WIDTH) then 
    tmp = path.dup
    step_forward(tmp << [x + 1, y], x + 1, y)
  elsif (y < HEIGHT) then 
    tmp = path.dup
    step_forward(tmp << [x, y + 1], x, y + 1)
  end
  if (x == WIDTH && y == HEIGHT) then
    p path
    $count += 1
    return
  end
end

step_forward([[0,0]], 0, 0)
print "total paths : #{$count}\n"

とりあえずちゃんと全部数え上げている。

% ruby p15.rb
[[0, 0], [1, 0], [2, 0], [2, 1], [2, 2]]
[[0, 0], [1, 0], [1, 1], [2, 1], [2, 2]]
[[0, 0], [1, 0], [1, 1], [1, 2], [2, 2]]
[[0, 0], [0, 1], [1, 1], [2, 1], [2, 2]]
[[0, 0], [0, 1], [1, 1], [1, 2], [2, 2]]
[[0, 0], [0, 1], [0, 2], [1, 2], [2, 2]]
total paths : 6

でも、グリッドが少し大きくなると、あっという間に非現実的な時間がかかるようになるので、ひとまずCで書き直してみた。

ルートをいちいち配列とかに記録しないようにしたり、良く見ると判定の順序のせいで無駄があったところも改善。さらに、スタート位置をズルして1つ進めた状態にしてみた。最後にルート数を2倍して正解を出せば所要時間は半分で済む……、って、そんなことやるぐらいなら紙と鉛筆で計算したらどうだという感じもするけど。

#include <stdio.h>
#include <stdlib.h>

int w = 20, h = 20;
double count = 0;

void step_forward (int x, int y) {
  if (x < w) {
    step_forward(x + 1, y);
  }
  if (y < h) {
    step_forward(x, y + 1);
  }
  if (x == w && y == h) {
    count += 1;
    return;
  }
}

int main (int argc, char *argv[]) {
  h = w = atoi(argv[1]);
  //  h = atoi(argv[2]);
  
  step_forward(0, 1);
  printf("total paths : %f\n", count*2);
  return 0;
}

これでもかなり時間がかかる。20x20のグリッドで実行すると2056秒、つまり約30分もかかって結果が出た。コンパイラオプションで-march=nativeはずしてたら1時間かかったかも……。

% gcc -O3 -march=native -o p15c p15c.c
% time ./p15c 20
total paths : 137846528820.000000
./p15c 20  2056.82s user 2.74s system 98% cpu 34:45.32 total

係数的な差を改善するだけでも計算が終わりそうな予測がついたから実行させてみたけど、ひどい富豪っぷり。

で、計算量を減らせないか考えてみて、数え方に激しい重複感があることを理解した。20x20のグリッドで考えたとき、例えば(10,10)にたどり着くルートは18万ほどあるけど、そこから(20,20)に向かうルートも同様に18万。ということで、(10,10)を通るルートだけで18万×18万あるわけだけど、これを調べ上げるのに、(10,10)->(20,20)を毎回調べるのは無駄。それは18万に決まってる。

メモ化? 動的計画法? とか思ったけど、簡単な図を書いてみて、ゴール側から足し算すればいいだけであることに気づいた。ある点からゴールに向かうルートの数は、その右と下のルート数を足したものに等しい。

=begin
* 2 x 2 grid
6 3 1
3 2 1
1 1 1

* 3 x 3 grid
20 10  4  1
10  6  3  1
 4  3  2  1
 1  1  1  1
=end

$w = ARGV.shift.to_i
$grid = Array.new($w + 1) do
  Array.new($w + 1) {0}
end

$w.step(0, -1) do |x|
  $w.step(0, -1) do |y|
    if(x == $w && y == $w) then
      $grid[x][y] = 1;
    else
      count = 0
      if x < $w then
        count += $grid[x + 1][y]
      end
      if y < $w then
        count += $grid[x][y + 1]
      end
      $grid[x][y] = count
    end
  end
end

# p $grid

puts "total paths (grid size: #{$w}) => #{$grid[0][0]}"

これで劇的に高速化できて、20x20のグリッドどころか1000x1000のグリッドでもRubyで5秒で数えられるようになった。

% ruby p15memo.rb 10
total paths (grid size: 10) => 184756
% ruby p15memo.rb 100
total paths (grid size: 100) => 9054851465610328116540417707748416387
4504589675413336841320
% ruby p15memo.rb 1000
total paths (grid size: 1000) =>
204815162698948971433516250298082504439642488798139703382038263767174
818620208375582893299418261020620146476631999802369241548179800452479
201804754976926157856301289663432064714851152395251651227768588611539
546256147907378668464154444533617613770073855673814589630071306510455
959514479888746206368718514551828551173166276253663773084682932255389
049743859481431755030783796444370810085163724827462791417016619883764
840843541430817785947037746565188475514680749694674923803033101818723
298009668567458560252549910118113525353465888794196665367490451130611
0096311906270342502293155911108976733963991149120
%

これってパスカルの三角形? ともあれ、再帰とメモ化で書き直してみた。最初からこう書くべきだったんだろうと思う。

$w = ARGV.shift.to_i
$grid = Array.new($w + 1) do
  Array.new($w + 1)
end
$grid[$w][$w] = 1

# $grid = [[nil, nil, nil], [nil, nil, nil], [nil, nil, 1]]

def route (x, y)
  return $grid[x][y] if $grid[x][y]
  count = 0
  if (x < $w) then
    count += route(x + 1, y)
  end
  if (y < $w) then
    count += route(x, y + 1)
  end
  $grid[x][y] = count
  return count
end

puts "total paths (grid size: #{$w}) => #{route(0, 0)}"
% time ruby p15memo2.rb 1000
total paths (grid size: 1000) =>
204815162698948971433516250298082504439642488798139703382038263767174
818620208375582893299418261020620146476631999802369241548179800452479
201804754976926157856301289663432064714851152395251651227768588611539
546256147907378668464154444533617613770073855673814589630071306510455
959514479888746206368718514551828551173166276253663773084682932255389
049743859481431755030783796444370810085163724827462791417016619883764
840843541430817785947037746565188475514680749694674923803033101818723
298009668567458560252549910118113525353465888794196665367490451130611
0096311906270342502293155911108976733963991149120
ruby p15memo2.rb 1000  6.53s user 0.22s system 96% cpu 6.963 total

1000x1000のグリッドで約7秒。計算量はO(n^2)であまり高速といえないけど、ナイーブに再帰で全部数え上げると、たぶんO(2^n)とかなので全然マシ。

もっと速くするには、組み合わせの計算でグリッドサイズがnのときには、ルートの数が、

(2n)!/(n!)(n!)

という関係を使うことで、これを素直に書くと、

def fact(n)
  return 1 if n == 1
  return n * fact(n - 1)
end

def route(n)
  fact(n*2)/(fact(n))**2
end

puts route(ARGV.shift.to_i)

という感じ。これだとサイズ1000のグリッドで0.06秒。

ところで、fact(n)の中でreturnが2つ並ぶのも文字が多くてアレなので、

def fact(n)
  1 if n == 1
  n * fact(n - 1)
end

と書き換えたのだけど、stack level too deep で落ちた。最後の値以外は明示的にreturnしないといけないのだった。かといって、

def fact(n)
  return 1 if n == 1
  n * fact(n - 1)
end

というのもなあ。こういうのがS式万歳のSchemerが嫌がるポイントなのかしらと想像してみたり。

階乗計算を再帰でやっているので、あまり深くなると落ちる。サイズ3000のグリッドはオッケーだけど4000で落ちる。ぼくの環境ではスタックの深さが7000ちょいで落ちるらしい。それで、

def fact(n)
  (1..n).inject(:*)
end

と書き換えたらサイズ2万のグリッドでも17秒ほどで答えが戻ってくるようになった。ただ、これだと1回で済むn!の計算を合計3回やってるので、以下のように、

def route(n)
  fact = (1..n).inject(:*)
  fact_b = fact * ((n+1)..(n*2)).inject(:*)
  fact_b/fact**2
end

puts route(ARGV.shift.to_i)

としたら、17秒→8秒と縮んだ。うーん、あれ、違うなまだ計算に重複がある。

def route(n)
  fact = (1..n).inject(:*)
  fact_b = ((n+1)..(n*2)).inject(:*)
  fact_b/fact
end

puts route(ARGV.shift.to_i)

これで4.6秒にまで縮んだ。

それはそうと、またしてもRubyの多次元配列の初期化でハマった。ふつうの配列は、以下のように要素数と要素を並べて初期化できる。

>> a = Array.new(3, "0")
=> ["0", "0", "0"]
>> a[0] = "1"
=> "1"
>> a
=> ["1", "0", "0"]
>> 

しかし、2次元配列の場合は、以下のような入れ子だと、内側の配列オブジェクトが1度しか生成されず、外側の配列には、それへの参照だけが3つ入るだけになる。

>> b = Array.new(3, Array.new(3, "0"))
=> [["0", "0", "0"], ["0", "0", "0"], ["0", "0", "0"]]
>> b[0][0] = "1"
=> "1"
>> b
=> [["1", "0", "0"], ["1", "0", "0"], ["1", "0", "0"]]
>> 

b[0][0]だけ書き換えたつもりが、b[1][0]もb[2][0]も変わってしまっている。

正しくは、

>> a = Array.new(3) {Array.new(3) {"0"}}
=> [["0", "0", "0"], ["0", "0", "0"], ["0", "0", "0"]]
>> a[0][0] = "1"
=> "1"
>> a
=> [["1", "0", "0"], ["0", "0", "0"], ["0", "0", "0"]]
>> 

のようにブロックで渡すようにする。

いやあ、オイラープロジェクトの問題なんてシンプルすぎていまいちと思ったけど、いい練習素材になるのかもなあと思った。