Hatena::Grouptopcoder

agwの日記 RSSフィード

2015-12-15二分探索で失敗しないために

二分探索で失敗しないために

06:32 |  二分探索で失敗しないために - agwの日記 を含むブックマーク はてなブックマーク -  二分探索で失敗しないために - agwの日記  二分探索で失敗しないために - agwの日記 のブックマークコメント

範囲を求める整数の二分探索を書いていると混乱するときってありますよね。でも全然書けないわけじゃないんです。すっきり書けることも多い。でも、こういうことがよく起こるんです。

  • 無限ループに陥る
  • 範囲を狭めるとき、下限と上限の何れを更新してよいか迷う
  • 答えが一個ずれている

ごく最近幾つかのコツを掴んでようやく迷わず書けるようになってきました。今後も二分探索で失敗しないためにまとめてみることにしました。

探索の型を見極める

範囲を求める整数の二分探索は大きく2つに分けることができます。これを設計段階で見極めておくことが重要であるようです。

本エントリではそれぞれの型をlower_bound型、upper_bound型と呼びます。C++のSTLで実装されているstd::lower_boundとstd::upper_boundを参考としました。日本語にすれば、下界型(下限型)、上界型(上限型)とでも言うのでしょうか。


f:id:agw:20151210175434p:image


灰色で図示されている要素は「i番目の要素をaiとしたとき、関数f(ai)が真である」ことを表しています。下限と上限は楔(くさび)で表しています。テキストエディタで言えば、「文字と文字の間にカーソルがある」ような位置の考え方です(余談ですが、実際にテキストエディタを設計する場合はこのような位置の考え方をすると設計が簡潔になるそうです)。

続いて、位置が要素の位置にある考え方も図示します。


f:id:agw:20151210175432p:image


上の楔型と要素型の図示の間には不整合性がありますが、後で明らかになりますので今は放っておきましょう。

文章でも言い表しておきましょう。

  • lower_bound型
    • 関数f(ai)が真となる最小のiを探索する
  • upper_bound型
    • 関数f(ai)が真となる最大のiを探索する

また、二分探索はその性質上、以下のようなデータには適用できません。


f:id:agw:20151214232412p:image


関数f(x)ってなに?

まずは簡単な例で考えてみましょう。整列済みの数値の配列を例にします。


f:id:agw:20151211151511p:image

lower_bound型としては以下のような問題を考えます。

aiが3以上である最小のiを求める

図示すると以下となります。0-indexで考えたとき、2が解となります。


f:id:agw:20151210194033p:image


このとき、関数f(x)は以下のような関数となります。

bool f(int x)
{
  return x >= 3;
}

upper_bound型としては次のような問題を考えます。

aiが3以下である最大のiを求める

図示すると以下となります。0-indexで考えたとき、5が解となります。


f:id:agw:20151210194034p:image


このとき、関数f(x)は以下のような関数となります。

bool f(int x)
{
  return x <= 3;
}

初期範囲を決める

求める解の範囲が分かっている場合、言い換えれば、ある範囲の中に必ずf(ai)が真となるiがある場合は以下のような範囲を用います。ここで、l*とu*はそれぞれ解が取りうる範囲の下限と上限を示します。

f:id:agw:20151210175435p:image

図示すると以下のようになります。


f:id:agw:20151210194031p:image


範囲を扱う際は右閉半開区間を用います。例えば、i ∈ [0, 7]の範囲に解があることが分かっているのであれば、初期の範囲を(-1, 7]とします。

f:id:agw:20151210175436p:image

C++の実装で表すと次のような感じになります。

  int l = -1, u = 7;

右閉半開区間を用いるのは計算機の整数演算を巧みに利用するためで、以下のような特徴があります。

  • 無限ループに陥らない
  • 同じ要素が複数回評価されることがない

反復の条件

二分探索は範囲をしぼり込むことによって探索を行います。範囲から要素がなくなってしまっては意味がありませんのでその直前で反復を停止します。具体的には範囲がただ一つの整数値を表している場合、つまり範囲が(x - 1, x]となったら反復を終了します。C++での実装は以下のようになります。

  while (u - l > 1) {
          :
  }

後でまた触れますが、このように実装したときにf(l* - 1)やf(u*)が評価されることはありません。

範囲を更新する

二分探索は反復毎に範囲を約半分に狭めていくことによってO(logN)の計算量を実現しています。各反復ではまず中央の値を計算します。整数演算ですので、端数は切り落とされます。

f:id:agw:20151210194902p:image

中央の値を使い、f(am)を評価します。その結果によって範囲を更新します。

  • mでuを更新する(更新した範囲は(l, m]となる)
  • mでlを更新する(更新した範囲は(m, u]となる)

さて、この節では一先ず細かいことを抜きにしてどのように範囲の更新を行うかを注視したいと思います。そのため、要素を詳細に記した図は一度お休みし、より概念的な図を用いることとします。

範囲を更新する際は、それが全体の範囲のどこにあるか、どの程度の大きさであるか等は意識しないほうがいいです。そのため、左右を切った図とします。lとuに挟まれている区間が注目している区間です(以下の図はuppper_bound型)。


f:id:agw:20151210195552p:image


また、更新する範囲は探索の型で異なります。lower_bound型から説明します。

lower_bound型における範囲の更新
  • lower_bound型であり、現在注目しているf(am)が真である場合

mは解かもしれませんが、f(ai)が真となる要素iはmよりも左側にまだあるかもしれません。ですので、左側をさらに詳しく探索します。


f:id:agw:20151210195555p:image


もちろん次の図のように、中央の値が求める解である場合もありますが、二分探索の途中でこれを意識する必要はありません。また、むしろ意識しないほうがよいようです。慣れてくると、放っておいてもアルゴリズムがこの解を絞り込んだ範囲に放り込んでくれることが分かってくるからです。


f:id:agw:20151214143126p:image


  • lower_bound型であり、現在注目しているf(am)が偽である場合

求める解はmよりも右側にあります。ですので、右側をさらに詳しく探索します。


f:id:agw:20151210195556p:image


C++で実装すると以下のようになります。

  if (f(m)) {
    u = m;
  }
  else {
    l = m;
  }

upper_bound型における範囲の更新

さて、ここからはupper_bound型を説明します。

  • upper_bound型であり、現在注目しているf(am)が真である場合

mは解かもしれませんが、f(ai)が真となる要素iはmよりも右側にまだあるかもしれません。ですので、右側をさらに詳しく探索します。


f:id:agw:20151210195553p:image


範囲を右閉半開区間として取ることにしたため、mは絞り込んだ範囲に含まれません。ですが心配することはありません。lower_bound型のときと同様に、最後には辻褄がちゃんと合います。

  • upper_bound型であり、現在注目しているf(am)が偽である場合

求める解はmよりも左側にあります。ですので、左側をさらに詳しく探索します。


f:id:agw:20151210195554p:image


C++で実装すると以下のようになります。lower_bound型と見比べてみるのもよいでしょう。

  if (f(m)) {
    l = m;
  }
  else {
    u = m;
  }

まずはやってみよう

upper_bound型の場合

まずはやってみましょう。lower_bound型の二分探索の様子を見てみましょう。C++での実装は以下のようになります。

bool f(int x)
{
  return x >= 3;
}

int lower_bound() {
  int l = -1, u = 7;

  while (u - l > 1) {
    int m = (l + u) / 2;

    if (f(m)) {
      u = m;
    }
    else {
      l = m;
    }
  }

  /* 解は何? */
}

探索の様子を図示してみましょう。


f:id:agw:20151210210255p:image


うまく行きましたね! 範囲は最後に(1, 2]となりました(l = 1、u = 2)。範囲の定義からも、図からもuが答えとなることが分かります。

何故こんなにうまく行くのか考えてみましょう。

先ほどの節で説明した通り、lower_bound型の二分探索では現在注目しているf(am)が真であるとき、f(ai)が真となる要素iはmよりも左側にまだあるかもしれないので、左側をさらに詳しく探索します(図を再掲します)。


f:id:agw:20151210195555p:image


範囲を右閉半開区間として取るため、mは絞り込んだ範囲に引き続き含まれます。例えばm未満に解がない場合、それ以降どのように探索されたとしても範囲は最終的に(m - 1, m]となります。

f(am)が偽であるとき、f(ai)が真となる要素iは必ずmよりも右側にあります。そのため、右側をさらに詳しく探索します。


f:id:agw:20151210195556p:image


範囲を右閉半開区間として取るため、次の範囲(m, u]にmは含まれません。そもそもmは解の候補ではありませんし、lower_bound型の定義からm未満の要素に解はないはずです。そのため思い切って範囲から除くことができるのです。

upper_bound型の場合

次にupper_bound型の二分探索の様子も見てみましょう。

bool f(int x)
{
  return x <= 3;
}

int upper_bound() {
  int l = -1, u = 7;

  while (u - l > 1) {
    int m = (l + u) / 2;

    if (f(m)) {
      l = m;
    }
    else {
      u = m;
    }
  }

  /* 解は何? */
}

探索の様子を図示してみましょう。


f:id:agw:20151210210257p:image


...こちらは微妙ですね。解がlの位置にあります。どうしてこうなったか考えてみましょう。

まずf(am)が偽であるときのことを考えてみましょう。mより大きいiのaiには解がありませんので、ばっさりと範囲から除くことができます(図は同じく、再掲です)。


f:id:agw:20151210195554p:image


さて、現在注目しているf(am)が真であるとき、それは解の候補となります。つまり、現在の要素、amは真に上限であるかもしれません。ですが、f(ai)となる要素iはmよりも右側にまだあるかもしれないのです。右側を詳しく検索し続ける必要があるはそのためです。

f:id:agw:20151210195553p:image


範囲は右閉半開区間として取るのでした。次の期間は(m, u]となります。つまり、一時的に解の候補であるmは範囲から取り除かれます。

ではamが真の上限、つまり解であったとしましょう。ai ∈ (m, u]には解がありませんよね。つまり、(m, u]のどのiの要素に対するf(ai)はすべて偽になります。

f:id:agw:20151214161832p:image

そのため、これに続く探索は常に範囲の左側を残すように範囲を絞り続けることとなり、結局のところ範囲は(m, m + 1]となります。求める解はmです。


f:id:agw:20151214234525p:image


f(am)が解の候補であるけれども真の上限ではない場合、真の上限となるようなm'を有する区間(l', u']が(m, u]の範囲の中に必ずあるはずですので、やはり右側を詳しく探索します。真の上限となるm'が見つかった場合、範囲は(m', u']を経て(m', m' + 1]となります。この場合の解はm'となります。


本当にうまくいくの?

確信を得るために色々な場合の結果を見てましょう。以下はlower_bound型の二分探索を試した結果です。lower_boun型の二分探索の結果はuを見ればよいのでした。それぞれうまく行っているようです。


f:id:agw:20151210231543p:image


upper_bound型の二分探索も試してみましょう。


f:id:agw:20151211151546p:image


upper_bound型の二分探索の結果はlを見ればよいのでした。upper_bound (1)の戻り値が-1となってしまいました。これはなんでしょうか?

よく考えてみると、何れのiを用いても関数f(ai)が真になることはありません。想定していない入力だったのでした。


また、upper_bound (8)の結果がおかしいですね。この結果はupper_bound (7)の結果と同じようです。試しに比較してみましょう。


f:id:agw:20151210231452p:image


これはおかしいですね。どうしたらよいでしょう?


よりよくするために

まずupper_bound (1)について。これはそもそもupper_bound型の二分探索の定義がおかしかったのです。私たちは以下のように定義したのでした。

  • upper_bound型
    • 関数f(ai)が真となる最大のiを返す

イメージしていたのは以下のような図でした(再掲)。


f:id:agw:20151210175432p:image


これを以下のように改訂しましょう。

  • upper_bound型
    • 関数f(ai)が真となる最大のiを探索し、その次の要素を返す

upper_bound型の二分探索のイメージは以下となります。


f:id:agw:20151210175433p:image


数値を入れたときの図は以下のようになります。


f:id:agw:20151210175442p:image


これでupper_bound (1)の解釈も通るようになりました。以前のようにupper_bound型の探索でf(ai)が真となる最大のiの要素の位置を知りたい場合はlを使うのではなく、u - 1を使うようにします。

これだとどうも分かり辛い、という方は冒頭に記載した楔型の図示をイメージするとよいと思います。


f:id:agw:20151210194036p:image


冒頭で少し触れたように、この方法は範囲を表す際によく使われる概念です。

STLのイテレータも楔型で考えると都合がよいものも多いです。std::begin()、std::end()、std::lower_bound()、std::upper_bound()、std::vector::insert()、std::vector::erase()等々...


f:id:agw:20151210175438p:image


例えばstd::vector::insert()は楔の位置に指定した要素群が挿入されると考えるとよいですし、std::vector::erase()も楔と楔の間の要素を削除すると考えるとよいかと思います。


話は少々脱線しましたが、次にupper_bound (8)について考えてみましょう。探索の初期範囲を(l* - 1, u*]ではなく、(l* - 1, u* + 1]としてみましょう。

f:id:agw:20151210210258p:image

図示すると以下のようになります。


f:id:agw:20151210194032p:image


上限を1つだけ大きくしただけですが、これだけで前の節での問題は不思議と解消されます。


f:id:agw:20151211151601p:image


もちろんlower_bound型もきちんと動作します。


f:id:agw:20151210235038p:image


それに加え、lower_bound(9)のように解が存在しないも適切に表現できるようになります。


f:id:agw:20151210234957p:image


「f(u* + 1)の挙動はそもそも未定義かもしれないじゃないか。もしかしたらアクセス違反するかも」と思うかもしれません。しかし、右閉半開区間を伴った二分探索ではmがu* + 1となることがありません。そのため、アルゴリズム的な曖昧さはなく、関数f(x)を変更する必要もありません。

またそもそも関数f(x)を解が存在する範囲以上の入力に対応させ、ある程度余裕を持たせた初期範囲を設定するのも常套手段であるようです。


まとめ

このエントリでは二分探索で失敗をしないためのコツを考えてきました。

  • 探索の型を見極める
    • lower_bound型か、upper_bound型か
  • 初期範囲を決める
    • 右閉半開区間を用いる
    • 最低限(l* - 1, u* + 1]を設定する
  • 反復の条件はu - l > 1が成立する間とする
  • 中央値mを算出し、f(am)を評価する
  • 評価の結果によって、範囲を更新する
    • lower_bound型でf(am)が真である場合、(l, m]を次の範囲とする
    • lower_bound型でf(am)が偽である場合、(m, u]を次の範囲とする
    • upper_bound型でf(am)が真である場合、(m, u]を次の範囲とする
    • upper_bound型でf(am)が偽である場合、(l, m]を次の範囲とする
  • 探索の型によって解を選ぶ
    • lower_bound型である場合、関数f(ai)が真となる最小のiはuである
    • upper_bound型である場合、関数f(ai)が真となる最大のiはu - 1である

文章にすると意外に長いですが、実際には実装のテンプレートは決まっています。if文の中身を決定することが出来ればelse句は自然と決まりますので、コツを覚えなければならないのは二ヶ所です。

int {lower|upper}_bound() {
  int l = L, u = U;

  while (u - l > 1) {
    int m = (l + u) / 2;

    if (f(m)) {
      /* lower_bound型ではu = m; */
      /* upper_bound型ではl = m; */
    }
    else {
      /* lower_bound型ではl = m; */
      /* upper_bound型ではu = m; */
    }
  }

  /* lower_bound型ではreturn u; */
  /* upper_bound型ではreturn u - 1; */
}

範囲の更新の部分を思い出すのは簡単です。lower_bound型であるなら、以下の図をちょいちょいと書くだけです(図は再々掲)。


f:id:agw:20151210195555p:image


何れの型にしても、f(am)が真となるような図を書くのがポイントのようです。

奥深い二分探索

このエントリはCompetitive Programming Advent Calendar 2015の12/16日分のエントリとして記載しました。

私はこのエントリに記載したコツを使うことによって随分迷わなくなりましたが、二分探索の世界は実に奥深く、様々な「必勝法」が存在し、それに関する議論はいつも尽きないようです。


上位コーダーによる記事はシンプルにまとめられています。よく引用されるエントリには必ず目を通しましょう。


プログラミングコンテストチャレンジブックの存在を忘れてはいけません。本エントリで記載した右閉半開区間を用いるアイデアはこの書籍を読んで得た知識です。プログラミングコンテストチャレンジブックでの二分探索のアプローチは@kinabaさんのアプローチに近いように感じています。


また、私は「ビット逐次決定型」と無理矢理名付けていますが、中央値を求めるような反復を行わず上位ビットから決定していく手法は、二分探索の議論において必ず言及されるようです。二分探索を違った視点から眺めることができますので、ご一読をお勧めします。

また、@hirose_golfさんによる手法は本エントリで取ったアプローチと同様にlower_bound型とupper_bound型を設計段階で意識したアプローチとなっています。しかしながら、lower_bound型とupper_bound型を明確に区別するように紹介しましたが、そのような説明をする記事は少ないように感じます。少し考えるとupper_bound型はf(ai)の否定(つまり! f(ai))が真となる最小のiを探す問題、つまりlower_bound型に置き換えることができるからかなあと考えています。

2015-06-21SRM 535 Div II

SRM 535 Div II

15:50 |  SRM 535 Div II - agwの日記 を含むブックマーク はてなブックマーク -  SRM 535 Div II - agwの日記  SRM 535 Div II - agwの日記 のブックマークコメント

  • 初参加したSRMにて惨敗した問題に挑戦。今でも自分には中々難しく、35分ほど時間を要した

FromAndGCDLCM

  • 未知の数AとBの最大公約数と最小公倍数、GとLが与えられる。ここからAとBを推測し、その最小の和(A + B)を返す問題
  • 最大公約数Gと最小公倍数Lに有効なAとBが存在しなければ-1を返す
  • 問題を見てまずぱっと思い浮かべた式は以下だ

f:id:agw:20150620234031p:image

  • AとBは[1, 1012]を取る値だ。そのため、上述の左辺が1024を取りうる。これは約1,0008、約280の値だ。素直に実装したらlong longから溢れてしまう
  • さて、上述の式を変形すると、以下となる。

f:id:agw:20150620234032p:image

  • 問題に習い、lcm(A, B)L、gcd(A, B)をGとすると、右辺はLGとすることが出来る。
  • AとBを乗してLGとなるようなAとBの組み合わせを求めるような実装は大抵、以下のような形となることが多い。計算量はO(√LG)だ
  for (int i = 1; i * i <= LG; i ++)
                    :
  • さて、この場合、AはGの倍数である。またlong longであることを加味すると、つまり、この実装は以下のような形になると思われる
  for (long long A = G; A * A <= LG; G ++)
                    :
  • LGがlong longに納まりきらないことを承知の上で実装を進めると、以下となる(最大公約数を今一度求めるのは実装をテストして気付いた。サンプルになければWAを出していたところだ)
  long long nm = INF;

  for (long long A = G; A * A <= LG; G ++)
    if (LG % A == 0) {
      long long B = LG / A;

      if (std::__gcd(A, B) == G)
        nm = std::min(A + B, nm);
    }
  • 繰り返すが、A・AやLGはlong longに収まらない。そのため、何らかの回避策を取る必要がある
  • AをiGと表すことにする。iを倍数として、long longから溢れさせないようにする作戦だ

f:id:agw:20150620234033p:image

  • 難しいのは反復の条件だ。A = iGとすると、以下とすることができる

f:id:agw:20150620234034p:image

  • また、iからBを求めるには以下とする

f:id:agw:20150620234035p:image

  • これらを用い、上述の実装をリファクタリングして、以下とすることができる
  long long nm = INF;

  for (long long i = 1; i * i <= L / G; i ++) {
    long long A = i * G;

    if (L % i == 0) {
      long long B = L / i;

      if (std::__gcd(A, B) == G)
        nm = std::min(A + B, nm);
    }
  }

まとめ

  • SRM参戦時に惨敗した問題。解けて素直に嬉しく思うと共に、Div II Mediumとしてはなかなか難しく、またとても良問であったように思えた
  • A、B、lcm(A, B)、gcd(A, B)を考えるとき、以下の関係式は大変重宝する

f:id:agw:20150620234031p:image

  • 積がNとなるようなiとjを組み合わせる際の典型の実装は以下となり、その計算量はO(√N)である
  for (int i = 1; i * i <= N; i ++)
    if (N % i == 0) {
      int j = N / i;
            :
    }

2015-03-18SRM 653 Div I

SRM 653 Div I

15:01 |  SRM 653 Div I - agwの日記 を含むブックマーク はてなブックマーク -  SRM 653 Div I - agwの日記  SRM 653 Div I - agwの日記 のブックマークコメント

http://community.topcoder.com/stat?c=round_overview&er=5&rd=16317

SRM 653 Div Iに参加。制限時間一杯Easyを考えるが、全く手が出ず。SRM後に幾つかの実装を読み、全く考えもしていなかった問題に還元できることを学んだ。

グラフの問題であると捉え作図を行ったため、メモの代わりとしてブログ記事として記載した。


CountryGroupHard

http://community.topcoder.com/stat?c=problem_statement&pm=13688&rd=16317

  • 以下のような数の並びが与えられる

f:id:agw:20150317214445p:image

  • 数値の周りには同じ数値が、数値分だけ集まる
  • 0は数値が隠されていることを示す
  • 1以上の数値は既知であることを示し、必ず正しい
  • 全体が成立するように、数値を推定したい
  • そのとき既知の数値から隠された数値を一意に推定できるかどうかを判断して回答する問題
  • 例えばこの例は数値が一意に求まり、以下のようなものであったと推定できる

f:id:agw:20150317214446p:image

  • 以下も数値が一意に求まる例だ

f:id:agw:20150317214447p:image

f:id:agw:20150317214448p:image

  • 以下は数値が一意に求まらない例

f:id:agw:20150317214449p:image

  • この数列を成立させる組み合わせは複数ある。書き出してみると5種類あった

f:id:agw:20150317214450p:image


  • SRMでは多くの参加者がDPもしくはメモ化再帰の問題として解いていたが、全く理解できなかった
  • SRM終了後に幾つかの実装を読み、ようやく理解をすることができた
  • 自分が理解する上で一番しっくりきたのが、経路の組み合わせを数え上げる問題に還元するものであった
  • 先ほどの例で考えてみよう。数値の列として捉えるのではなく、

f:id:agw:20150317214446p:image

  • 以下のように成立する経路を考え、その経路を数え上げる問題とするような感覚で捉えるのが自分には合っているらしい

f:id:agw:20150317214451p:image


  • さて、実際には上記のように極度に抽象化したグラフではなく、数値群の間に頂点を配置して考えた
    • 初期状態では左端の頂点にのみ経路が設定されていることに注意したい。またその数は1だ

f:id:agw:20150317214452p:image

  • 一番左の数値に1を設定してみよう。結果、左端の頂点から次の頂点への経路が成立する。つまり、左端の頂点から次の頂点への経路は少なくとも1種類存在すると言える
    • そのため、二つ右側への頂点へ左端の頂点に到達する経路の数、1を足し込む

f:id:agw:20150317214453p:image

  • 一番左の数値に2を設定してみよう。このときは左端の頂点から二つ右側への頂点への経路が成立する。先ほどと同様に、左端の頂点から二つ右側への頂点への経路は1種類存在すると言える
    • そのため、二つ右側への頂点へ左端の頂点に到達する経路の数、1を足し込む

f:id:agw:20150317214454p:image

  • 一番左の数値に3を設定してみよう。3の周辺には3つの3が集まるはずだが、隣の数値は2で、条件が成立しない。そのため、経路も成立せず、この場合は左端の頂点から三つ右側への頂点への経路が存在しないことが分かる

f:id:agw:20150317214455p:image

  • 左側の頂点からの経路をまとめてみよう。以下のようになる

f:id:agw:20150317214456p:image

  • 二つ目の頂点から同様に考える。しかし、この頂点から成立する経路がない

f:id:agw:20150317214457p:image

  • 三つ目の頂点からも考えてみよう。以下のように右側の数値群を3で埋めてやると経路が成立することが分かる

f:id:agw:20150317214458p:image

  • 三つ目の頂点からの経路をまとめる

f:id:agw:20150317214459p:image

  • 四つ目の頂点からの経路を考える。四つ目の頂点からの経路は成立するのだが、そもそも左端の頂点から四つ目の頂点に到達する経路がないため、最終的な経路の組み合わせ数に寄与しない

f:id:agw:20150317214500p:image

  • 五つ目の頂点からの経路も同様に寄与しない

f:id:agw:20150317214501p:image

  • 結果、左端の頂点から右端の頂点への経路は1種類しかないと言え、問題の答えも一意に推定できるものとすることができる。以下は左端から右端の頂点までの経路をまとめたものだ

f:id:agw:20150318001357p:image


  • 以下は成立する例

f:id:agw:20150317234106p:image

  • その経路をまとめたもの

f:id:agw:20150318001358p:image


  • 以下は成立しない例

f:id:agw:20150317214531p:image

  • 結果、左端の頂点から右端の頂点へ到達しうる経路は5種類であることが分かる

f:id:agw:20150318001359p:image

  • 確かに書き下した組み合わせの数も5種類であった

f:id:agw:20150317214450p:image


まとめ

  • まさか数え上げに還元できる問題であるとは微塵も想像もしておらず、大変驚きであった
  • 次回同様の問題が出たら是非解いて、提出したいところだ

おまけ

SRM後に行った実装(システムテストをパスした)。

class CountryGroupHard {
public:
  std::string solve(std::vector<int> a) {
    int size = a.size();

    long long dp[111] = {0};

    dp[0] = 1;

    For (int i = 0; i < size; i ++)
      for (int j = 1; i + j <= size; j ++) {
        bool valid = true;

        for (int k = 0; k < j; k ++)
          if (a[i + k] && a[i + k] != j)
            valid = false;
        
        if (valid)
          dp[i + j] = (dp[i + j] + dp[i]) % MOD;
      }

    return dp[size] == 1 ? "Sufficient" : "Insufficient";
  };
};

2014-12-05マラソンマッチにおける典型データ構造とアプローチ このエントリーを含むブックマーク このエントリーのブックマークコメント

TopCoderのマラソンマッチは大きく以下の二つの形式に分類することができます。

  • 最適化
  • 機会学習

何れの形式のマッチもそれぞれの魅力的がありますが、今回は特に最適化に主眼を置いたマッチを題材にお話しをします。最適化を主眼においたマッチがどのようなものであるか、2013年の12月に開催されたMM 82を用いて簡単に紹介します。ColorLinkerという問題です。


プログラムには正方形のグリッドが与えられます。グリッドのいくつかのセルは予め彩色がなされています。

我々に与えられた課題は、同じ色のセルが全て繋がるように空のセルを彩色していくことです。セルを彩色するにはコストがかかります。また、同一のセルを複数色に彩色するとコストは指数的に増加します。合計のコストが低ければ低いほどよい彩色であると評価されます。同じ色のセルが連結でない場合は無効となります。

f:id:agw:20141204204722p:image


この問題に対する典型的なアプローチの一つは最短経路探索、特に複数始点から同時にダイクストラ法を適用し焼きなましを行う、というものであったようです。


さて、TopCoderの最適化問題ではこのようにシンプルな二次元グリッドの問題がほどよい頻度で出題されます。そのため、二次元グリッドを効率よく扱うことはTopCoderのマラソンマッチに参加する上で必須のテクニックであると言えます。このエントリでは、TopCoderマラソンマッチにおける典型データ構造とアプローチを主題に以下のトピックに精選し、紹介したいと思います。

  • 盤面の持ちかたと最短経路探索
  • キュー
  • 経路復元
  • コスト付き最短経路探索

盤面の持ちかたと最短経路探索

まず初めに盤面の持ちかたについて考えてみましょう。実際のマラソンマッチの問題は過分に複雑なので、これを模した問題を用意することにします。

迷路が与えられます。迷路は通路と壁で構成されています。その後3秒の間、二点がランダムに、繰り返し与えられます。この二点間の最短経路を出来るだけ速く探索するメソッドを書くのを目標とすることとします。説明やサンプルプログラムを簡単にするために、このプログラムは最短距離の長さも経路も返しません。もはや問題とはなっていませんが、ご容赦ください。

さて、目的のクラス名をsolver_tとしましょう。このクラスは以下のような二つのメソッドを持つものとします。

void solver_t::init(const std::vector<std::string>& a);

void solver_t::bfs(int sx, int sy, int ex, int ey);

initメソッドにはいわゆる「迷路」が渡されます。迷路がプログラムに渡されるのは初めの1回だけとします。データは、std::stringのstd::vectorとします。パウンド(#)が壁、ピリオド(.)が通路を表します。迷路の縦と横の長さは等しくNとします。Nの範囲は毎回異なり、15〜55としましょう。迷路生成の都合から、Nは奇数とします。範囲外には外壁があるものとします。

実際の迷路のデータを覗いてみましょう。競技プログラミングではよく見る感じのデータですので、慣れているかたも多いのではないでしょうか。

.....#.#...#.....#.
.###.#.#.###.#.###.
.#.......#...#.#...
##.###.###.#######.
...#...#...#.#.....
.###.#.#.###.#.###.
.#...#.#.....#.#...
##.###.#.#.###.#.##
...#...#.#.#.#.#.#.
.#####.###.#.#.###.
...#.......#.....#.
.#.#####.###.#####.
.#.#.#.#...........
##.#.#.###.########
.#.#...#...#.......
.#####.#.#.###.####
...#.....#.....#...
.###.#####.###.###.
.......#...#.......

bfsメソッドには開始点と終了点が渡されます。そこから以下のような最短経路を探索したい、ということです。

f:id:agw:20141204204746p:image

この問題を計算機で解く際はグラフの問題として捉えると便利なことがあります。その様子を少し覗いておきましょう。通路を頂点と置き換え、隣接した頂点間に辺を張ります。

f:id:agw:20141204204747p:image


さて、bfsメソッドという名前が諮詢しているように、我々の仕事は速い幅優先探索を書くことです...


では早速実装を開始します。迷路のデータはSTLを使って渡されることですし、そのままデータとして使うことにします。

クラスの定義はこのような感じになるでしょうか。

class solver_t {
pubilc:
  void init(const std::vector<std::string>& a);

  void bfs(int sx, int sy, int ex, int ey);

private:
  std::vector<std::string> a_;

  int N_;

  std::vector<std::vector<bool>> visit_;
};

solver_t::init()では盤面をそのまま保持します。また、盤面の大きさに伴ったvisit_変数を確保しておきましょう。

void solver_t::init(const std::vector<std::string>& a)
{
  a_ = a;

  N_ = a_.size();

  visit_.assign(N_, std::vector<bool>(N_));
};

最後にbfs()メンバ関数を記述します。とてもオーソドックスな幅優先探索の実装です。

void solver_t::bfs(int sx, int sy, int ex, int ey)
{
  for (int i = 0; i < N_; i ++)
    for (int j = 0; j < N_; j ++)
      visit_[i][j] = false;

  visit_[sy][sx] = true;

  std::queue<int> queue;

  queue.push(XY(sx, sy));

  while (! queue.empty()) {
    int p = queue.front(); queue.pop();

    int x = X(p);
    int y = Y(p);

    for (int k = 0; k < 4; k ++) {
      int xx = x + dx[k], yy = y + dy[k];

      if (0 <= xx && xx < N_ && 0 <= yy && yy < N_) {
        if (xx == ex && yy == ey)
          return;

        if (! visit_[yy][xx])
          if (a_[yy][xx] == '.') {
            visit_[yy][xx] = true;

            queue.push(XY(xx, yy));
          }
      }
    }
  }
}

座標値はintにパックするようにしました。パックには以下のマクロを用いています。与えられる迷路のサイズが最大55であるということで、55以上の一番小さなの2の累乗の数(2nの数。2、4、8、16、...等)、64を選びました。これにより、ビットシフトやビットマスクを用いて座標値のパックやアンパックをすることができます。

#define X(x) ((x) & 63)                 // Unpacking
#define Y(x) ((x) >> 6)

#define XY(x, y) (((y) << 6) + (x))     // Packing

下位16ビットの使い方は以下のようになっています。

f:id:agw:20141204204744p:image

また、dx、dyには上左下右の位置に対応する差分値が格納されています。SRMでもよく見る典型的な手法です。

const int dy[] = {-1,  0,  1,  0};
const int dx[] = { 0, -1,  0,  1};

ここで、この実装のベンチマークを取ってみましょう。initメソッド経由で55 × 55の盤面を与え、3秒の間にどれだけbfsメソッドを呼べるかを計測してみました。私の端末では111,410回呼び出すことができました(OS X 10.9.5 + clang + -O2フラグを使用)。


さて、TopCoderのマラソンマッチではアルゴリズムの善し悪しが勝負を決めます。これは決定的で揺るぎない事実です。そのため、闇雲に高速化を行うことは一般的に悪手とされています。しかし、思いついたアルゴリズムがトップ争いを征するポテンシャルを備えていたとしても、アルゴリズムを計算する速度が遅いと勝つことができません。これは勿体ない。

例えばここで説明をしている単純な幅優先探索等はその動作が速ければ速いほど有利になるかもしれません。少しばかり時間をかけ、チューンアップすることを考えてみましょう。


まずSTLのパフォーマンスを疑ってみましょう。盤面をCの二次元配列で持つように直してみます。メンバ変数の宣言と初期化部分の実装を変更するだけですので、変更箇所は少ないです。

class solver_t {
private:
  char a_[64][64];

  int N_;
  
  bool visit_[64][64];
};
void solver_t::init(const std::vector<std::string>& a)
{
  N_ = a.size();

  for (int i = 0; i < N_; i ++)
    for (int j = 0; j < N_; j ++)
      a_[i][j] = a[i][j];
}

二次元配列に置き換える際に考慮すべきこととして、各要素の要素数が上げられます。座標値をパックした際と同様に、まず2の累乗の数を検討してみるのがよいでしょう。ここでも55以上の一番小さな2の累乗の数、つまり64を選択しました。

たったこれだけの変更ですがbfsを呼び出す回数は164,900回に増えました。約1.5倍の速度改善です。STLのコンテナは大変便利ですが、処理速度を重視する必要がある場合、Cの原始的な配列に置き換えることを検討するのは良手の一つです。


ではこの実装をさらに速度改善する方法を考えてみます。一次元配列に直してみましょう。メンバ変数は以下のようになります。盤面を表すa_メンバ変数だけではなく、visit_メンバ変数にも一次元配列を採用していることに注意してください。

class solver_t {
private:
  int N_;

  int LASTROW_;
  
  char a_[4096];		// 64 * 64 = 4096

  bool visit_[4096];
};

LASTROW_は大変便利な変数です。盤面の範囲をチェックする際に用います。これに関しては後ほど説明をします。

盤面の初期化は以下のようになります。

void solver_t::init(const std::vector<std::string>& a)
{
  N_ = a.size();

  LASTROW_ = (N_ - 1) << 6;

  for (int i = 0, p = 0; i < N_; i ++, p += 64)
    for (int j = 0, pp = p; j < N_; j ++, pp ++)
      a_[pp] = a[i][j];
};

また、bfsの実装は以下のようになります。

void solver_t::bfs(int sx, int sy, int ex, int ey)
{
  memset(visit_, 0, sizeof(visit_));

  int s = XY(sx, sy), e = XY(ex, ey);

  std::queue<int> queue;

  queue.push(s);

  visit_[s] = true;

  while (! queue.empty()) {
    int p = queue.front(), pp;

    queue.pop();

    int x = X(p);

    if (p >= 64 && ! visit_[pp = p - 64] && a_[pp] == '.') {
      if (pp == e)
        return;

      visit_[pp] = true;

      queue.push(pp);
    }
    if (x > 0 && ! visit_[pp = p - 1] && a_[pp] == '.') {
      if (pp == e)
        return;

      visit_[pp] = true;

      queue.push(pp);
    }
    if (p < LASTROW_ && ! visit_[pp = p + 64] && a_[pp] == '.') {
      if (pp == e)
        return;

      visit_[pp] = true;

      queue.push(pp);
    }
    if (x < N_ - 1 && ! visit_[pp = p + 1] && a_[pp] == '.') {
      if (pp == e)
        return;

      visit_[pp] = true;

      queue.push(pp);
    }
  }
}

何かすごいことになってしまった感もありますが、これまでと同様に3秒間の呼び出し回数を計測してみると282,420回となりました。Cの二次元配列を用いていたときから2倍弱の速度改善で、STLに比べると3倍弱の速度改善、といったところでしょうか。


何故このような速度改善が達成できるかは後ほど考えてみるとして、まず実装を見て行きましょう。whileによる反復部分です。pはキューから取得した現在の位置です。この値はx座標値とy座標値をintにパックしたものでした。

ここから、以下の処理を上左下右のそれぞれの方向に対して適用しています。

  • 次の位置に進めるかを判定し、
  • ppに次の位置を保持しつつすでに訪れているか否かを判定し、
  • 次の位置が壁でなければ、
    • 目的の位置であれば、探索を終了する
    • 目的の位置でなければ、次の位置をキューに積む。この際、次の頂点に訪れることを記録する

「キューに積む際に次の頂点に訪れることを記録する」のはとても大事なので特に注意しましょう。これにより、同じ頂点が再びキューに積まれることを抑制することができます。


さて、位置pから上の方向に訪れることができるかどうかは、位置pが2行目以降にあればよいことになります。これは、pと64を比較すれば判定できます。

では位置pから下の方向に訪れることができるかどうかはどのように判定すればよいでしょうか? 予め計算しておいたLASTROW_をここに用います。LASTROW_はこのように計算していますので、

  LASTROW_ = (N_ - 1) << 6;

以下の図に示す位置を表すことになります。そのため、現在の位置pがLASTROW_より小さければ、下方向に探索を進めても盤面からはみ出ないことが保証されるのです。

余談ですが、この変数の使い方はTopCoderマラソンマッチでは必ず上位に食い込んでいるwleiteさんの実装を眺めていて見つけました。とてもスマートな方法ですよね。

水平方向に関しては現在のx座標値をデコードしたほうが判定がしやすいのでそのようにします。これをまとめると以下となります。

f:id:agw:20141204225908p:image


さて、何故このような速度改善が達成できたのでしょう? メモリ構造が異なるのでしょうか? いえ、Cの二次元配列と一次元配列は何れも連続して確保されるはずです。例えば、char a[5][10]として配列を確保した場合とchar a[50]として配列を確保した場合、アドレッシングの差はあれど、メモリの確保のされかたは等しくなるはずです(大きな図で見たい場合は画像をクリックしてください)。

f:id:agw:20141204204816p:image


そのため、アドレッシングをする際にかけ算かビットシフトが入ることが性能低下の原因であると考えるのが自然です。

ここで少しだけアセンブラを眺めてみましょう。64 × 64の二次元グリッドにランダムに数値を設定し、ランダムにピックアップした頂点の値のxorを取る関数を書いてみます。繰り返し数は10000回です。

int rand64();		// ランダムに[0, 64)を返す関数が外部定義してあるものとする
int rand4096();		// ランダムに[0, 4096)を返す関数が外部定義してあるものとする


static int a[64][64];

int test_a() {
  int x = 0;

  for (int i = 0; i < 10000; i ++)
    x ^= a[rand64()][rand64()];
  
  return x;
}

int test_b() {
  int* p = reinterpret_cast<int*>(&a[0][0]);

  int x = 0;

  for (int i = 0; i < 10000; i ++)
    x ^= p[rand4096()];
  
  return x;
}

二次元配列を使ったtest_a関数のアセンブラ(の主要な部分)は以下のようになります。

	xorl	%r15d, %r15d		# %r15dレジスタを0にする
	movl	$10000, %ebx		# %ebxレジスタを10000にする
	leaq	__ZL1a(%rip), %r12	# 二次元配列aの先頭のアドレスを%r12レジスタに格納する
LBB1_1:
	callq	__Z2rand64v		# rand64関数をコールする
	movl	%eax, %r14d		# 結果を%r14dレジスタに格納し、
	andl	$63, %r14d		# 63andを取る(これがxとなる)
	callq	__Z2rand64v		# rand64関数をコールする
	andl	$63, %eax		# 63andを取る(これがyとなる)
	shlq	$8, %rax		# y8ビットシフトする
                                        # (int4バイト長であるため、6 + 2ビットシフトしている)
	addq	%r12, %rax		# aの先頭のアドレスを足す
	xorl	(%rax,%r14,4), %r15d	# さらに、x4倍分を足したアドレスの値とxorを取る
	decl	%ebx			# %ebxレジスタから1を引く
	jne	LBB1_1			# %ebxが0でなければ反復する

これに比べ、一次元配列を使ったtest_b関数のアセンブラは以下のようになります。反復の内側で実行されるニモーニックの数がとても少ないことが分かります。

	xorl	%r14d, %r14d		# %r14dレジスタを0にする
	movl	$10000, %ebx		# %ebxレジスタを10000にする
	leaq	__ZL1b(%rip), %r15	# 二次元配列aの先頭のアドレスを%r15レジスタに格納する
LBB2_1:
	callq	__Z2rand4096v		# rand4096関数をコールする
	andl	$4095, %eax             # 4095andを取る(これをpとする)
	xorl	(%r15,%rax,4), %r14d	# b + p × 4のアドレスの値とxorを取る
	decl	%ebx			# %ebxレジスタから1を引く
	jne	LBB2_1			# %ebxが0でなければ反復する

たとえ二次元配列を2の累乗の数で確保していたとしても、アドレッシングの計算の際にビットシフトと加算、また多少のレジスタのコピー操作が入ってしまうのです。たったこれだけなのですが、速度的には相当分が悪いのです。

キュー

さて、さらなる速度改善を考えてみましょう。具体的にはstd::queueもより原始的なデータで置き換えられないでしょうか?

キューの先頭と末尾の位置を保持するインデックスを2つ用意すると、一次元配列でキューを作成することができます。例えばアルゴリズムCに掲載されているキューの実装は以下のようにとてもシンプルです。

#define max 100
static int queue[max+1],head,tail;
put(int v)
  {
    queue[tail++] = v;
    if (tail > max) tail = 0;
  }
int get()
  {
    int t = queue[head++];
    if (head > max) head = 0;
    return t;
  }
queueinit()
  { head = 0; tail = 0; }
int queueempty()
  { return head == tail; }

このデータ構造はstd::queueと比べ、ずっとずっと原始的です。例えば、キューに溜めるべき要素の数がmax個を超えた場合、プログラムは破綻します。そのため、キューに溜まる最大の要素の数を予め考察しておくことが大切なのですが、なかなか難しく、テストを繰り返しながらパラメータを調整する手法に頼りがちです。

問題によっては発想を少し変えることによってキューの実装がより簡単で安全なものになります。今回の迷路の最短経路探索にて使用している幅優先探索はそのような問題の一つです。

幅優先探索が訪れる最大の頂点数がどのくらいか考えてみましょう。すでに訪れている頂点にはもう訪れることがないことを考えると高々N2です(Nは頂点数ではなく、迷路の幅であることを思い出してください)。そのため、要素数がN2の一次元配列を確保しておけばどのような状況であってもキューがはみ出さないことが分かります。ここではアライメントを考慮して、642分の要素を確保することにします。さらに、putやget関数を以下のように簡略化することができます。

#define max 4096
static int queue[max],head,tail;
put(int v)
  { queue[tail++] = v; }

int get()
  { return queue[head++]; }

if文と代入を取り除いているだけですが、速度に大きく影響することにも注意してください(自分の端末だと、10%ほど速度に影響を与えました)。


さて、では実装を見てみましょう。メンバ変数としてqueue_を用意します。

class solver_t {
private:
  int N_;

  int LASTROW_;
  
  char a_[4096];

  bool visit_[4096];

  int queue_[4096];
};

bfsの実装は以下のようになります。

void solver_t::bfs(int sx, int sy, int ex, int ey)
{
  memset(visit_, 0, sizeof(visit_));

  int s = XY(sx, sy), e = XY(ex, ey);

  visit_[s] = true;

  int q1 = 0, q2 = 0;

  queue_[q1 ++] = s;

  while (q1 != q2) {
    int p = queue_[q2 ++], pp;

    int x = X(p);

    if (p >= 64 && ! visit_[pp = p - 64] && a_[pp] == '.') {
      if (pp == e)
        return;

      visit_[queue_[q1 ++] = pp] = true;
    }
    if (x > 0 && ! visit_[pp = p - 1] && a_[pp] == '.') {
      if (pp == e)
        return;

      visit_[queue_[q1 ++] = pp] = true;
    }
    if (p < LASTROW_ && ! visit_[pp = p + 64] && a_[pp] == '.') {
      if (pp == e)
        return;

      visit_[queue_[q1 ++] = pp] = true;
    }
    if (x < N_ - 1 && ! visit_[pp = p + 1] && a_[pp] == '.') {
      if (pp == e)
        return;

      visit_[queue_[q1 ++] = pp] = true;
    }
  }
};

これまでと同様に3秒間の呼び出し回数を計測してみましょう。結果、373,280回となりました。初めに比べると3.5倍以上の速度改善を達成したことになります。

経路復元

さて、次は機能拡張をしてみましょう。bfsメソッドが幅優先探索によって発見した最短経路を返すように変更します。つまり、経路復元を行います(経路復元については@machyさんが去年のアドベントカレンダーで記載したスライドが大変秀逸ですのでご一読することをお勧めいたします)。

経路復元は以下のように行います。

  • 「その頂点にはどの頂点から訪れたか?」という情報を別途記録しておく
  • 終点を発見したら、始点までさかのぼる

「その頂点にはどの頂点から訪れたか?」という情報は相対位置情報(上から来た、下から来た等)でも構いませんが、ここでは絶対位置情報を記録することにします。

探索が終わったときの経路情報は以下のようなものとなります。見方を変えると、これは始点を根とした木であるとも言えます(薄く描画されている辺は、終点が発見されたことにより探索が打ち切られ、調べられなかった辺を表しています)。

f:id:agw:20141204234440p:image

さて、メンバ変数を追加しましょう。木を意識して、parent_という変数名を選択しました。

class solver_t {
private:
  int N_;

  int LASTROW_;
  
  char a_[4096];

  bool visit_[4096];

  int queue_[4096];

  int parent_[4096];
};

bfsの実装は以下のようになります。0はパックされた座標として有効ですので、-1で根を表すことにしました。また、探索を行う毎にparent_の内容を全て初期化する必要がないことに注意しましょう。探索と共に経路情報が木のように育っていることを意識すると分かりやすいかもしれません。

std::vector<int> solver_t::bfs(int sx, int sy, int ex, int ey) {
  memset(visit_, 0, sizeof(visit_));

  visit_[XY(sx, sy)] = true;

  parent_[XY(sx, sy)] = -1;

  int e = XY(ex, ey);

  int q1 = 0, q2 = 0;

  queue_[q1 ++] = XY(sx, sy);

  while (q1 != q2) {
    int p = queue_[q2 ++], pp;

    int x = X(p);

    if (p >= 64 && ! visit_[pp = p - 64] && a_[pp] == '.') {
      parent_[pp] = p;

      if (pp == e)
        break;

      visit_[queue_[q1 ++] = pp] = true;
    }
    if (x > 0 && ! visit_[pp = p - 1] && a_[pp] == '.') {
      parent_[pp] = p;

      if (pp == e)
        break;

      visit_[queue_[q1 ++] = pp] = true;
    }
    if (p < LASTROW_ && ! visit_[pp = p + 64] && a_[pp] == '.') {
      parent_[pp] = p;

      if (pp == e)
        break;

      visit_[queue_[q1 ++] = pp] = true;
    }
    if (x < N_ - 1 && ! visit_[pp = p + 1] && a_[pp] == '.') {
      parent_[pp] = p;

      if (pp == e)
        break;

      visit_[queue_[q1 ++] = pp] = true;
    }
  }

  std::vector<int> path;

  for (int p = e; p != -1; p = parent_[p])
    path.push_back(p);

  std::reverse(std::begin(path), std::end(path));

  return path;
};

さて、今回は絶対位置情報を使いましたが、何れのデータ構造が得であるかは問題と場合に依ります。相対位置情報で記録したほうが省メモリかもしれません。ですが、経路復元をする際に座標値の再計算が必要となるかもしれません。また、絶対位置情報では再計算を行わないかもしれませんが、相対的な方向を解として保持しなけらばならない場合、やはり再計算が必要となるかもしれません。

コスト付き最短経路探索

最後にコスト付き最短経路探索を考えてみましょう。頂点のコストを加味した最短経路探索はこれまで扱ってきたコストを考えない最短経路探索よりもずっと複雑ですが、より現実的な問題設定であると言えると思います。

ここでもこれまでに扱ってきた迷路のデータを使用しましょう。これまで壁は通れないものとして扱ってきましたが、通り抜けられるようにしてしまいましょう。そして、以下のようにコストを設定してみましょう。

  • 通路を一つ進む度にコストが1かかる
  • 壁を一つ乗り越える度にコストが10かかる

この条件下では壁を一つ乗り越えることで、10個以上の通路をスキップしてしまったほうがよい状況を考えなければいけなくなります。例えば、以下の例では始点のすぐ下で一回壁を乗り越えたほうが得であるようです(参考として、壁を乗り越えない場合の最短経路を薄い色の経路で表しています)。

f:id:agw:20141204204751p:image

随分と実際のマラソンマッチの様相を呈してきました。事実、この問題設定は冒頭で紹介したMM82のものにとても類似しています。


先ほども触れましたが、この問題設定はなかなか複雑で難しいものです。前半でやったのと同様に、迷路をグラフとして表現してみましょう。以下のような疎ではありますが、先ほどまでと比べるととても密度の高いグラフとなります(白い頂点は通路を表し、黒い頂点は壁を表しています)。

f:id:agw:20141204204748p:image


さて、この最短経路探索はどのようなアルゴリズムを用いるのが適切でしょうか? アルゴリズムに詳しいかたは迷いなくダイクストラ法を選択するでしょう。

二次元配列を使った実装から見て行くことにしましょう。優先度付きキューを伴ったダイクストラ法の実装は計算量を抑えることが知られています。std::priority_queueを使いましょう。コストが低い頂点から優先して処理を行うため、std::greaterを併用していることに気を付けてください。

コストの上限は64 × 64 × 10として40960ですので、16ビット用意すれば大丈夫そうです。また、コストは初期化時に12ビット左にシフトしてしまうことにしました。位置情報と共にパックする際の実装が簡潔になるからです。そのため、パック、アンパックをするためのマクロも改変します。

#define X(x) ( (x)        &    63)
#define Y(x) (((x) >>  6) &    63)
#define Z(x) ( (x)        & -4096)
#define P(x) ( (x)        &  4095)

#define XY(x, y) (((y) << 6) + (x))

#define PZ(p, z) ((z) + (p))

#define XYZ(x, y, z) PZ(XY((x), (y)), (z))

ビットの使い方は以下のようになりました。

f:id:agw:20141204204745p:image

では実装を見て行きましょう。a_にはもはや文字を収めません。12ビット分左にシフトしたコストを記録します。そのため、char型の配列ではなく、int型の配列として確保します。

class solver_t {
private:
  int a_[64][64];

  int N_;
  
  int cost_[64][64];
};

迷路のデータを受け取り次第、12ビット分左にシフトします。

void solver_t::init(const std::vector<std::string>& a) {
  N_ = a.size();

  for (int i = 0; i < N_; i ++)
    for (int j = 0; j < N_; j ++)
      a_[i][j] = (a[i][j] == '.' ? 1 : 10) << 12;
};

実装は以下のようになります。

int solver_t::dijkstra(int sx, int sy, int ex, int ey) {
  for (int i = 0; i < N_; i ++)
    for (int j = 0; j < N_; j ++)
      cost_[i][j] = INF;

  std::priority_queue<long long, std::vector<long long>, std::greater<long long>> pq;

  pq.push(XYZ(sx, sy, cost_[sy][sx] = a_[sy][sx]));

  while (! pq.empty()) {
    int x = X(pq.top());
    int y = Y(pq.top());
    int z = Z(pq.top());

    pq.pop();

    if (cost_[y][x] < z)
      continue;

    for (int k = 0; k < 4; k ++) {
      int xx = x + dx[k], yy = y + dy[k];

      if (0 <= xx && xx < N_ && 0 <= yy && yy < N_) {
        int zz = z + a_[yy][xx];

        if (xx == ex && yy == ey)
          return zz;

        if (zz < cost_[yy][xx]) {
          cost_[yy][xx] = zz;

          pq.push(XYZ(xx, yy, zz));
        }
      }
    }
  }

  return -1;
};

3秒間の計測を行います。39,490回でした(参考までに、STLを使った場合の実装は30,900回でした)。


さて、こちらを一次元配列を使って書き直してみましょう。std::priority_queueはそのまま使用します。

class solver_t {
private:
  int a_[4096];

  int N_;

  int LASTROW_;
  
  int cost_[4096];
};
int solver_t::dijkstra(int sx, int sy, int ex, int ey) {
  for (int i = 0, p = 0; i < N_; i ++, p += 64)
    for (int j = 0, pp = p; j < N_; j ++, pp ++)
      cost_[pp] = INF;

  int s = XY(sx, sy), e = XY(ex, ey);

  std::priority_queue<int, std::vector<int>, std::greater<int>> pq;

  pq.push(PZ(s, cost_[s] = a_[s]));

  while (! pq.empty()) {
    int x = X(pq.top());
    int p = P(pq.top()), pp;
    int z = Z(pq.top()), zz;

    pq.pop();

    if (cost_[p] < z)
      continue;

    if (p >= 64) {
      pp = p - 64;

      if ((zz = z + a_[pp]) < cost_[pp]) {
        if (pp == e) {
          return zz;
        }
        else {
          pq.push(PZ(pp, cost_[pp] = zz));
        }
      }
    }
    if (x > 0) {
      pp = p - 1;

      if ((zz = z + a_[pp]) < cost_[pp]) {
        if (pp == e) {
          return zz;
        }
        else {
          pq.push(PZ(pp, cost_[pp] = zz));
        }
      }
    }
    if (p < LASTROW_) {
      pp = p + 64;

      if ((zz = z + a_[pp]) < cost_[pp]) {
        if (pp == e) {
          return zz;
        }
        else {
          pq.push(PZ(pp, cost_[pp] = zz));
        }
      }
    }
    if (x < N_ - 1) {
      pp = p + 1;

      if ((zz = z + a_[pp]) < cost_[pp]) {
        if (pp == e) {
          return zz;
        }
        else {
          pq.push(PZ(pp, cost_[pp] = zz));
        }
      }
    }
  }

  return -1;
};

一次元配列を用いた場合、dijkstraメソッドは45,840回呼び出すことができました。STLを用いたものに比べ、1.5倍高速です。


さて、ここまで来ると優先順位付きキューを一次元配列を利用できないのかな? と考えてしまいますね。しかし、優先順位付きキューは単純なキューのように機械的に一次元配列等に変換することはできません。残念。

ここで、優先順位付きキューの代わりに単純なキューを使うと何が起こるか考えてみましょう。この問題では負のコストがないため、優先順位付きキューを使えば最もコストが少ない頂点から訪れることが保証されていました。しかし、単純なキューを使った場合、例えば終点に初めて訪れた際のコストが最小である保証はなされません。

ならば、キューが空になるまで繰り返し、それぞれの頂点に訪れたコストの中で最小のものを記録するようにしてみましょう。言い換えれば、一つ前の頂点からの緩和を繰り返す、ということになります。何かどこかで見たことがありますね...? そうです、これはキューを伴ったベルマンフォードアルゴリズムです。

int bellman_ford(int sx, int sy, int ex, int ey) {
  for (int i = 0, p = 0; i < N_; i ++, p += 64)
    for (int j = 0, pp = p; j < N_; j ++, pp ++)
      cost_[pp] = INF;

  int s = XY(sx, sy), e = XY(ex, ey);

  int q1 = 0, q2 = 0;

  queue_[q1 ++] = PZ(s, cost_[s] = a_[s]);

  int zm = INF;

  while (q1 != q2) {
    int q = queue_[q2 ++];

    int x = X(q);
    int p = P(q), pp;
    int z = Z(q), zz;

    if (cost_[p] < z || z > zm)
      continue;

    if (p >= 64) {
      pp = p - 64;

      if ((zz = z + a_[pp]) < cost_[pp]) {
        if (pp == e) {
          zm = std::min(cost_[pp] = zz, zm);
        }
        else {
          queue_[q1 ++] = PZ(pp, cost_[pp] = zz);
        }
      }
    }
    if (x > 0) {
      pp = p - 1;

      if ((zz = z + a_[pp]) < cost_[pp]) {
        if (pp == e) {
          zm = std::min(cost_[pp] = zz, zm);
        }
        else {
          queue_[q1 ++] = PZ(pp, cost_[pp] = zz);
        }
      }
    }
    if (p < LASTROW_) {
      pp = p + 64;

      if ((zz = z + a_[pp]) < cost_[pp]) {
        if (pp == e) {
          zm = std::min(cost_[pp] = zz, zm);
        }
        else {
          queue_[q1 ++] = PZ(pp, cost_[pp] = zz);
        }
      }
    }
    if (x < N_ - 1) {
      pp = p + 1;

      if ((zz = z + a_[pp]) < cost_[pp]) {
        if (pp == e) {
          zm = std::min(cost_[pp] = zz, zm);
        }
        else {
          queue_[q1 ++] = PZ(pp, cost_[pp] = zz);
        }
      }
    }
  }

  return zm;
};

実装部分はまあこれでよいとして、とても注意すべきことがあります。もはやキューの最大要素数はN2ではありません。緩和が終わるまでそれぞれの頂点に複数回訪れることを考えると不思議ではありません。

この問題でデータのアライメントを64と設定した場合、4,096の16倍、65,536個分を用意する必要がありました。理論的に計算できればカッコいいのですが、私はこれを実験的に決定しました。実際には50,000個ほどでも十分だと思います。

  int queue_[65536];

優先順位付きキューを伴ったダイクストラ法では45,840回の探索でしたが、一次元配列によるキューを伴ったベルマンフォード法では59,600回の探索を行うことが出来ました。ほぼ1.3倍の速度改善です。STLを使ったものに比べると約2倍のパフォーマンス改善を行ったことになります。


この実装のポイントとして、以下の3種類の枝刈りを行っていることに触れておきましょう。

  • 注目する頂点により低いコストで到達する経路がある場合
  • 注目する頂点にたどり着いたコストが、現在までに見つかっている終点までのコストよりも大きい場合
  • 隣接する次の候補の頂点により低いコストで到達する経路がある場合

2番目の枝刈りの条件は忘れやすいと思いますが(「z > zm」の条件)、特に重要です。これを忘れると、ダイクストラ法のパフォーマンスに遥かに及びません。


「ベルマンフォード法の計算量はダイクストラ法と比べ随分大きいんじゃないのか? 遅くなるんじゃないのか?」と考えたかたも多いのではないのでしょうか? 確かにその通りで、ベルマンフォード法を適用できるか否かは問題の性質に大きく依存します。例えば大部分が同じ設定であったとしても、盤面の大きさが999 × 999である場合は以下のことが起こります。

  • キューの要素数を相当大きく確保しておく必要がある(例えば、67,108,864要素(= 1024 × 1024 × 64))
  • 優先順位付きキューを伴ったダイクストラ法と比べ、半分以下の速度となってしまう。

最後に、この条件下での経路復元について簡単に触れておきます。キューを伴ったベルマンフォード法の場合も幅優先探索と変わらず、「その頂点にはどの頂点から訪れたか?」という情報を記録すれば十分に経路復元をすることが可能です。

こちらの問題設定では先ほどと比べ、多少複雑なグラフを扱うことになります(以下の図は再掲)。

f:id:agw:20141204204748p:image

しかし経路の情報は、やはり始点を根とした木の様相を成します。以下の図をよく見ると、薄い辺があります。これらの辺は緩和の課程で一度は接続されたものの他の最短な経路によって切断された辺です。ベルマンフォード法では同じ頂点に何度も訪れるため、親の頂点情報も何度も入れ替わるのです。

f:id:agw:20141204204753p:image

まとめ(または、回顧録)

今回はマラソンマッチにおける典型データ構造とアプローチという題材で多少のトピックを解説しました。それぞれのトピックを見直すと、以下の共通点があることが挙げられます。

  • 原始的なデータ構造およびロジックは速い
  • 二次元配列は意外に遅い
  • 一次元配列は十分に単純であるため、速い
  • メモリを動的に確保するのは遅い
  • 問題に適したミニマムかつ効率的な設計を行う

「ひゃー、マラソンマッチではこんなに雑多な実装を行わなきゃいけないのか。つまらないじゃないか」と思われた方もいらっしゃるかもしれません。再掲となりますが、マラソンマッチで真に大切なのは問題にフィットしたロジックであり、アルゴリズムです。このエントリに記載されていることはある意味「おまけ」なのですが、トップ陣はそのおまけの部分も異様に強い。まるで魔術師を見ているようです。しかし、一つ一つのテクニックを丁寧に観察しひも解くと、実装の大部分に明確な理屈があり、ここで記載したような小さな工夫やチューニングが巧みに全体を構成しているのが見えてきます。このエントリがそんな「面倒臭い部分」に向き合うときの少しばかりの助けになれば大変嬉しく思います。


さて、この記事はCompetitive Programming Advent Calendar 2014 - PARTAKEのマラソン部門、12月6日分のエントリとして記載しました。これまでも、またこれからも続々とマラソンに関する記事が執筆されます。大変貴重なイベントだと思います。

アドベントカレンダーが始まってから実力者のエントリが続き、「問題をどう解くか」という核心の部分を熱く熱く執筆くださっています。素晴らしい記事を読むにつれ、自分の選んだ題材はあまりにも稚拙だよなあ? と迷いながらも、特に日本のマラソンマッチコンペティターの底上げに少しでも貢献できればと思い、なんとか書き上げました。


さて、そもそも今年の競技プログラミングアドベントカレンダーにマラソン部門が出来たことに関して少し触れておきたいと思います。私が日本に一時帰国した際にhotpepsiさんが開催してくださった飲み会の席でのことでした。突然、simezi_tanが熱く、そして激しく、「マラソンは何をしていいか分からない。競技プログラミングのアドベントカレンダーにマラソン部門が欲しい!」と声を上げたのです。tomerunさんやmachyさん達の援護が加わったこと、また、_tanzaku_さんが柔軟な企画、そして対応をしてくださったこともあり、気がついたときには見たこともない素晴らしいイベントとして運用されることになっていました。日本の競技プログラミングコミュニティの柔軟さには舌を巻くばかりです。


そんなsimezi_tanに伝えたいことがあります。この場を借りて記したいと思います。

「おいしめたん。そもそもの言い出しっぺなんだから参加しなはれ(何が分からないか書いてみるとか)」


明日、12月7日はyowaさんが「マラソンマッチで1位になる方法(最終順位とは言ってない)」という題材で執筆なさる予定です。お楽しみに。

2014-04-07Past, Present and Future

Past, Present and Future

13:40 |  Past, Present and Future - agwの日記 を含むブックマーク はてなブックマーク -  Past, Present and Future - agwの日記  Past, Present and Future - agwの日記 のブックマークコメント

コンテスト中に自分のブログから該当するエントリを探す作業をよくしていることに気付いたので、インデックスを作成しておくことにした。

Google Code Jam 2012 Round 1A
  • Kingdom Rush
    • std::sort()のコンパレータオブジェクトについて考察
SRM 549 Div II

f:id:agw:20140406193454p:image

SRM 553 Div II

f:id:agw:20140406173301p:image

  • Suminator
    • 二分探索について考察
    • 特に、閉区間を用いた探索を考察した
SRM 556 Div I

f:id:agw:20140406174217p:image

SRM 564 Div I

f:id:agw:20140406173302p:image

TCO13 Round 1C

f:id:agw:20140406173306p:image

  • TheOlympiadInInformatics
    • 二分探索について考察
    • 特に、右閉半開区間を用いた範囲の探索について考察を行い、閉区間を用いた範囲の探索との比較をした
SRM 340 Div II

f:id:agw:20140406175725p:image

Codeforces Round # 185 Div II

f:id:agw:20140406173303p:image

SRM 594 Div I(1)

f:id:agw:20140406213943p:image

  • FoxAndGo3
    • 最大安定集合問題について考察
Codeforces # 223 Div. 2

f:id:agw:20140406173304p:image

Telling If a Number is Prime

f:id:agw:20140406181333p:image

Implementing Sieve of Eratosthenes

f:id:agw:20140406181334p:image

How to Find Palindromic Substrings?

f:id:agw:20140406181336p:image

Understanding How std::string::substr works

f:id:agw:20140406181335p:image

SRM 611 Div II

f:id:agw:20140406173305p:image

  • LCMSetEasy
    • gcdおよびlcmについて考察。また、gcd(x1, x2, ..., xn)、lcm(x1, x2, ..., xn)についても考察
Precalculating a Combination

f:id:agw:20140624004013p:image