Think Stitch
PRINCIPIA  最近の更新


2分探索

2分探索が好きなんです。

「2分探索のコードなんか自分で書くな。ライブラリを使え。」

そうなんですけど、ちょっと次の問題を考えてみませんか?

問題

プロポーショナルフォント(文字によって幅の異なるフォント)を使って文字列を表示するとき、与えられた幅に表示できる最大の文字列長を知りたいことがあります。 文字列の表示幅を計算してくれる関数 width(s, k) があるとしましょう。 ここで s は文字列、k は幅を計算する先頭からの部分文字列の長さです。 文字列 s の長さを s.len とすると、0 <= k <= s.len でなければなりません。 width(s, s.len) が文字列 s 全体の表示幅です。 width(s, 0) = 0 であるとします。

負の幅を持つ文字はないと仮定すると、width(s, k) は k に対して単調に増加することになります。 厳密には、幅 0 の文字があるかもしれないので、非減少になります。 (たまにそういうフォントがありますよね。) そこで、幅 w >= 0 が与えられたとき、表示できる最大の k を2分探索で求められるはずです。 お使いの言語/ライブラリの2分探索関数にはどのような引数を渡せばいいでしょうか。 2分探索関数から返された値と k とのご関係は?

気になること

話を進める前に、ちょっと気になることを2つ書きます。

  • 使っているライブラリの2分探索関数は、非減少の場合をサポートしてますか?。 そもそも非減少の場合をサポートしているかどうか、文書に明言されているでしょうか。
  • 「そもそも幅 w に等しい値を見つけるのが目的ではないのだから、2分探索関数の対象ではないだろう」と考える人もいるかもしれません。 定められた適用範囲という意味では、まあ、そうかなとも思います。 でも、この問題が2分法という考え方で解けないかというと、そんなことはないわけです。 実際、この問題の解と2分探索のコードは同じものです。 だから、もしライブラリの2分探索関数の仕様が、この問題に適用できない形になっているとしたら、ちょっともったいないと思うのです。

問題の精密化

もう少し問題を精密にしてみます。 まず、求める k は、幅 w に収まるのですから、

width(s, k) <= w

を満たさなければなりません。 加えて、これが最大の k だということは、おおざっぱにいうと、

w < width(s, k + 1)

ということです。つまり、1文字でも増やすとはみ出ちゃうってことです。 k は s.len になる可能性がありますから、その場合 k + 1 は条件外になってしまうので、厳密には次のようにする必要があります。

k < s.len --> w < width(s, k + 1)

ここで記号 --> は「ならば」を意味するとします(Isabelle に合わせました)。 記号の左側が成り立たないときには、右側は評価されないと考えておきます(C/C++ 言語の && のように)。

まとめると、求める k が満たすべき条件は次のようになります。

width(s, k) <= w & (k < s.len --> w < width(s, k + 1))

記号 & は「かつ」を意味します。同時に成り立つということです。

2分探索関数の仕様

問題の条件が明確になりました。 これで、2分探索関数をどのように呼び出せばよいか、ライブラリの文書と比較して決められるでしょうか。 手持ちの本やネット上のリファレンスなどをいくつか見てみましたが、日本語で(あるいは英語で)説明が書かれていて、結局上の条件とどう絡めていいのかよくわかりません。(英語だともっとわかりません)。 いつも思うんですけど、リファレンスの文章って、おおまかに何をしてくれるのかを把握するにはもちろんいいんですけど、細かいところとか、極端な状況での振る舞いとか、わからなくて困るんです。 「おれがやりたいのはこういうことなんだよ、なんでそのときのことが書いてないんだ」という感じです。

ここでいきなり飛ばしますが、私がいつも愛用している2分探索のコードは次の2つです。 ライブラリの関数みたいに汎用化(あるいは総称化)すると見にくいと思ったので、簡単に配列の場合で書きました。n は配列 v の大きさです。

int binsearch_A(int x, int v[], int n)
{
    int i = 0, j = n, m;
    while (i != j) {
        m = (i + j) / 2;
        if (v[m] < x)
          i = m + 1;
        else
          j = m;
    }
    return i;
}
int binsearch_B(int x, int v[], int n)
{
    int i = 0, j = n, m;
    while (i != j) {
        m = (i + j) / 2;
        if (v[m] <= x)
          i = m + 1;
        else
          j = m;
    }
    return i;
}

違いは if 文のテストで、不等号に = が入っているかいないかだけです。 どちらの場合も事前条件は配列 v が非減少で整列されていることです。 事後条件はそれぞれ次のようになります。

0 <= i <= n &
(ALL h. 0 <= h < i --> v[h] < x) &
(ALL h. i <= h < n --> x <= v[h])
0 <= i <= n &
(ALL h. 0 <= h < i --> v[h] <= x) &
(ALL h. i <= h < n --> x < v[h])

ここで ALL h. Q(h) --> R(h) という形の記号は、Q(h) を満たすすべての h について R(h) が成り立つという意味です。 ちょっと複雑ですけど、いっていることはこういうことです。

  • 返される i の値の範囲は 0 <= i <= n です。i = n になる可能性があります。
  • 2番目の項は、i より小さいインデックス h では v[h] < x であるといっています。
  • 3番目の項は、i 以上のインデックス h では x <= v[h] であるといっています。

2つ目の関数 binsearch_B では、事後条件でも = の位置が入れ替わるだけです。

これがそれぞれ上のコードの仕様ということになります。 これで結局何がわかったんでしょうか。 x は見つかったんでしょうか、見つからなかったんでしょうか。 v[i] = x ってこと・・・ではなさそうです。

問題のゴールと関数仕様の比較

ちょっとその結論は後回しにして、文字列幅の問題に戻ります。 (別にもったいぶっているわけではないんですけど。) 求める条件と事後条件を見比べると、すぐにわかることがあります。

width(s, k) <= w & (k < s.len --> w < width(s, k + 1))
(ALL h. 0 <= h < i --> v[h] <= x) & (ALL h. i <= h < n --> x < v[h])

width(s, k) は v[k] に対応することになりますから、類似性が見て取れます。 等号 = の入っている位置から、binsearch_B の方を選ばなければならないことがわかります。

ここで意外なことが2つあるんです。 気づかれたかどうかわかりませんが、事後条件で v[h] <= x が成り立つのは、インデックス h が i よりも小さい場合です。これに対して width(s, k) <= w では、k の時にも不等号が成り立ちます。 もう1つの条件も同じようにずれていますから、結局 k = i - 1 だということです。

もう1つは v[h] ではインデックス h の変域は 0 <= h < n ですけど、k の変域は 0 <= k <= s.len で、k = s.len の場合が含まれています。したがって、n = s.len + 1 としなければなりません。

これにはちょっと驚きました。 もし日本語の解説だけを読んで考えていたら(少なくとも私は)気がつかなかったと思います。 きっと動かしてからおかしいことに気づいたでしょう。

結局、このアルゴリズムを使うとしたら、n = s.len + 1 として呼び出し、k = i - 1 として結果を得ることになります。ここでもう1つ問題があります。事後条件によると、i は 0 になりうるので、その場合はどうなるのかということです。

そこで v[h] と width(s, k) の対応、および n = s.len + 1 の元、事後条件で i = 0 と置いてみると、

(ALL h. 0 <= h < 0 --> width(s, h) <= w) &
(ALL h. 0 <= h < s.len + 1 --> w < width(s, h))

前者は自動的に成り立ちます(前件が偽だからです。理由を知らない人は信じてください:笑)。 後者はすべての h について w < width(s, h) であるといっています。 そうであれば、特に h = 0 としてもよいはずです。 すると w < width(s, 0) となりますが、width(s, 0) = 0 ですから w < 0 となります。 これは前提に反しますから、結局 i = 0 となることはないということです。 したがって、2分探索の関数から返ってきたら、チェックすることなく自信を持って k = i - 1 としてよいわけです。

i > 0 ですから、事後条件から width(s, i - 1) <= w & (i <= s.len --> w < width(s, i)) がいえます。したがって、k = i - 1 とすれば、求める結果が得られたことがわかります。

事前・事後条件による仕様記述の優位性

読むのが面倒な細かい議論になってしまいましたが、詳細はともかく、ポイントは次の点にあります。 事前・事後条件で関数仕様を表しておくと、求める結果との関係を厳密に議論できます。 特に細かいところや極端な場合どうなるのかといったことがわかります。 だから、ライブラリなどでは日本語での大まかな説明に加えて、事前条件、事後条件を書いてくれればいいのにと思うのです。

事前・事後条件を使うと、何が得られるのかということだけを簡潔に表すことができます。 よく、関数の仕様を尋ねると、中で何をどのように計算しているのかを一生懸命説明してくれる人がいますけど、申し訳ないことに、それはどうでもいいんです。欲しいものはそれではないのですよね。

結局のところ探索結果はどうなるの?

保留にしておいた話に戻ります。 ふつうに探索を目的に呼んだ場合、結果はどう判断されるのかということです。 binsearch_A の事後条件を再掲します。

0 <= i <= n &
(ALL h. 0 <= h < i --> v[h] < x) &
(ALL h. i <= h < n --> x <= v[h])

もし i = n であればすべてのインデックス h について v[h] < x となるので、v[j] = x となる j は存在しません。 i < n とすると、x <= v[i] がいえます。しかし x = v[i] かどうかはわかりません。 したがって、x = v[i] かどうかを知るには、関数から返ってきたあとで、次のようにする必要があります。

  if (i < n && v[i] == x)
    /* found */
  else
    /* not found */

「なんて面倒な!」と思われるかもしれません。 だからといって、これを関数の中に入れてしまうと、文字列の表示問題には適用できなくなってしまいます。

これは関数のインターフェイスをどのように切るかという深い問題であると同時に、コードから得られる収穫を最大にするにはどうしたらいいか、という問題でもあると思うのです。 普段はあまり意識しないかもしれませんが、コードというのは書いた瞬間ある性質を「持ってしまう」ものですよね。 (バグだってそうです :P) 書いた人がどの範囲で使おうと考えていても、それとは関係なく、そのコードがなんらかの結果を生成できるもっとも広い入力の範囲というものがあります。 もちろん結果が有用であるかどうかは問題なのですけど、せっかくならできるだけ適用範囲を広く取れるようにした方がよいと思うのです。

だから2分探索の例でいうと、上の関数はそのままで公開して、値が見つかったかどうかを判定する版は、これに皮をかぶせたものを別に用意すればいいと思うのです。

上の2つの関数は他の用途にも使えます。 例えば、非減少列の中で同じ値がつづく部分を「平原」といったりしますが、この先頭と末尾の位置を求めたりできます。 C++ の標準ライブラリにもあったと思います(C++ はよく知らないのです)。

探索って何をしてるの?

はじめは文字列幅の話だけにするつもりだったのですが、ついだらだらと来てしまったので、ここまで来たらとことんやってしまうことにします。 次のコードは binsearch_A アルゴリズムを Isabelle で検証したものです。 例によって細かい話は抜きにして、ループ不変条件の話だけします。

lemma "ALL i j. i <= j & j < length (v::int list)
         --> v ! i <= v ! j
  ==>
  VARS i j m
  {True}
  i := 0;
  j := length v;
  WHILE j ~= i
  INV {i <= length v & j <= length v & i <= j &
  (ALL k. k < i --> v ! k < x) &
  (ALL k. j <= k & k < length v --> x <= v ! k)}
  DO
    m := (i + j) div 2;
    IF v ! m < x THEN
      i := m + 1
    ELSE
      j := m
    FI
  OD
  {i <= length v &
  (ALL k. k < i --> v ! k < x) &
  (ALL k. i <= k & k < length v --> x <= v ! k)}"

(1点だけ補足をしておくと、Isabelle ではインデックスが自然数(乱暴にいえば unsigned)なので、下限を 0 で押さえる必要はありません。)

ループ不変条件を抜き出すと次のようになっています。

i <= length v & j <= length v & i <= j &
  (ALL k. k < i --> v ! k < x) &
  (ALL k. j <= k & k < length v --> x <= v ! k)

事後条件と違って、ループ変数 i, j が両方使われています。 ループの条件が i ~= j ですから、ループを抜けると i = j になるわけです。 この条件の中で、最後の2つは次のことをいっています。

  • i よりも小さい範囲では v[k] < x
  • j 以上の範囲では x <= v[k]

ループが回るたびに i と j は近づいていきます。 逆に言うと、上の2つの条件でいわれている範囲が広がっているということです。 つまり、この2分探索のコードは、値を探しているというよりむしろ、値がないと確定した部分を広げているという感じなのだと思います。 そういう視点で線形探索なんかを見ると、やはり「ここにはなかった」という範囲を拡大しているように解釈できます。 探索とは「ないところを広げる」ことなのかなと思いました。

補足

コードの次の部分

m = (i + j) / 2;

これはオーバーフロー時のトラブルを避けるために、ほんとは次のように書いてあります。

d = (j - i) / 2;
m = i + d;

2分探索のコードにはいろいろ変種(バリエーション)があります。 変数 j の初期値を n - 1 にするものとか、m = (i + j) / 2 に +/-1 がついたりとか。 配列のインデックスが 1 から始まる言語との輸出入なんか大騒ぎです。 結局、証明しないと、とてもではないですけど確信が持てません。 たいていのものは変数変換をすると同じものであることがわかります。 計算量の異なるものは・・・少しだけ違うものをみたことがあるような、ないような・・・。 実はあり得るすべてのデータを用意して、ループ回数を調べ上げることもできるんですよね。

このコードが好きなもう1つの理由は、ループ変数 i, j が負にならないからです。 unsigned でも使えます。 (負になるループ変数を unsigned で宣言して無限ループになるバグを出したことがあります。まったく間抜けです。)

いろいろと本を調べてみると、あちこちに間違った2分探索のコードが書かれていた(いる?)ようです。 このページは大丈夫かなあ。

「ループの中で v[m] == x の比較はしないのか?/すべきだ!」というご意見があると思います。 これは・・・もうやめにします。今日は疲れました^^;

2013/08/21
© 2013,2014,2015 PRINCIPIA Limited