Idiom: Cumulative Summation Summary
累積和を実装するとき、いつもインデックスの扱いに悩んでしまいます。これが結構時間を喰うので 、しっかりイディオムとしてまとめようと思いました。
例えば0-indexの配列、\(a_i\)があるとします。
$$ \begin{array}{c|ccccc} i & 0 & 1 & 2 & 3 & 4 & 5\\\hline a_i & 1 & 2 & 3 & -4 & -2 & 1 \end{array} $$
目的はO(1)で\(s_{[i,j)}\)を求めることです(準備にはO(N)かかります)。例えば、\(s_{[0,1)}\)は1、\(s_{[0,3)}\)は6、\(s_{[1,3)}\)は5といった具合です。
これを行うには累積和をテーブルとして用意するのが常套手段です。これを\(\mathit{acc}_i\)とします。\(\mathit{acc}_{i+1} = a_i + \mathit{acc}_i\)とします。また、\(\mathit{acc}_0\)は0です。つまり、
$$ \begin{cases} \mathit{acc}_0 = 0\\ \mathit{acc}_{i+1} = a_i + \mathit{acc}_i\ \text{where}\ 0 \leq i < \mathit{size}\\ \end{cases} $$
作成された累積和テーブルは以下のようになります。
$$ \begin{array}{c|cccccc} i & 0 & 1 & 2 & 3 & 4 & 5 & 6\\\hline \mathit{acc}_i & 0 & 1 & 3 & 6 & 2 & 0 & 1 \end{array} $$
ここで添字がずれるのが混乱の元なんですね。しっかり式を立てておきましょう。
$$ s_{[i, j)} = \mathit{acc}_j - \mathit{acc}_i $$
ポイントは、
となります。結果は以下のようになります。
$$ \begin{array}{c|c} s_{[0,1)} = \mathit{acc}_1 - \mathit{acc}_0 = 1 - 0 = 1\\ s_{[0,3)} = \mathit{acc}_3 - \mathit{acc}_0 = 6 - 0 = 6\\ s_{[1,3)} = \mathit{acc}_3 - \mathit{acc}_1 = 6 - 1 = 5\\ \end{array} $$
さて、1-indexの配列の場合にはどうしたらよいでしょう?
$$ \begin{array}{c|ccccc} i & 1 & 2 & 3 & 4 & 5 & 6\\\hline a_i & 1 & 2 & 3 & -4 & -2 & 1 \end{array} $$
\(a_0\)を0とおけば-0indexの場合と同様に扱うことができそうです。
$$ \begin{array}{c|ccccc} i & 0 & 1 & 2 & 3 & 4 & 5 & 6\\\hline a_i & 0 & 1 & 2 & 3 & -4 & -2 & 1 \end{array} $$
累積和テーブルは一つずれますが、0-indexと同じ考え方で扱うことができそうです。累積和テーブルの先頭の2つの要素を0で初期化しておかなければならないので注意が必要です。
$$ \begin{cases} \mathit{acc}_0 = 0\\ \mathit{acc}_1 = 0\\ \mathit{acc}_{i+1} = a_i + \mathit{acc}_i\ \text{where}\ 1 \leq i \leq \mathit{size}\\ \end{cases} $$
作成された累積和テーブルは以下のようになります。
$$ \begin{array}{c|cccccc} i & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7\\\hline \mathit{acc}_i & 0 & 0 & 1 & 3 & 6 & 2 & 0 & 1 \end{array} $$
\(s_{[i,j)}\)は以下のように算出することができます。
$$ s_{[i, j)} = \mathit{acc}_j - \mathit{acc}_i $$
もう一つ気になることといえば、配列は要素数 + 1、また累積和は要素数 + 2分確保しなければならないことでしょうか。
では累積和を一つずらすことを考えてみましょう。この場合の累積和は要素数 + 1分確保すればよいです。
$$ \begin{cases} \mathit{acc}_0 = 0\\ \mathit{acc}_{i} = a_i + \mathit{acc}_{i-1}\ \text{where}\ 1 \leq i \leq \mathit{size}\\ \end{cases} $$
作成された累積和テーブルは0-indexのときと変わりません。
$$ \begin{array}{c|cccccc} i & 0 & 1 & 2 & 3 & 4 & 5 & 6\\\hline \mathit{acc}_i & 0 & 1 & 3 & 6 & 2 & 0 & 1 \end{array} $$
これは0-indexのときの累積和と全く同じです。\(s_{[i, j)}\)を計算するときは、以下とします。
$$ s_{[i,j)} = \mathit{acc}_{j-1} - \mathit{acc}_{i-1} $$
書き出したことで頭がスッキリしてきました。特に、0-indexのときの悩みはなくなりそうです。1-indexのときはどうしたらいいんでしょうか。もう少し経験を積んでから考えたほうがよさそうです。
Idiom: Check if given string starts with some string
std::stringとして与えられた文字列がある特定の文字列から始まるかどうかを記述するのは簡単なことのように思えます。これが自分にとっては意外に難しく、いつも時間を割いてしまうのです。先日、コードリーディングをしていた際にうまい方法を見つけたため、イディオムとして記載しておこうと思います。
先日(またもや)やってしまった実装は以下のようなものでした。
const std::string s("abc"); // abcから始まる文字を探したい const int size = s.size(); const std::string a[] = { "abcdef", "defabc", "abcabc", "a" }; for (const auto& i : a) if (i.size() >= size) // 文字列の長さが十分かどうかチェックする if (i.substr(0, size) == s) // 先頭の文字列を切り出して比較する // 文字列iはsから始まっている! };
この実装は色々と無駄があるだけではなく、自分にとっては注意すべき点が多く手早く書けないのです。
- 文字列iのサイズを検証する -> この際等号を入れるべきだろうか?
- 文字列iの先頭size文字をstd::stringとして生成する -> この際使うべきsubstrメソッドには引数を1つ与えるべきか。それとも2つ与えるべきか?
そもそもsubstrメソッドは末尾のサイズを適切にチェックしてくれるはずです。
for (const auto& i : a) // if (i.size() >= size) // サイズのチェックはしなくてよい if (i.substr(0, size) == s) // 文字列iはsから始まっている!
これでも悪くはないのですが、やはりsubstrの引数の数は悩んでしまいそうです。どうもあれだなあ、とコードリーディングをしていた際に、以下のイディオムを見つけました。
for (const auto& i : a) if (i.find(s) == 0) // 文字列iはsから始まっている!
なるほど簡潔で覚えやすそうです。次回からはこのイディオムを使おうと思います。
もう少しだけ
findメソッドを使ったこのイディオムは簡潔でよいのですが、先頭から指定して文字列で始まっていない場合、無駄に探索を行ってしまいます(例えば文字列が"xyzxyzxyzxyzxyzxyzxyz"の場合)。これが無駄に思える場合はstrncmpを使った方がいいかもしれません。
for (const auto& i : a) if (strncmp(i.c_str(), s.c_str(), size) == 0) // 文字列iはsから始まっている!
Redirecting Contents of File to Standard Input
1行シェルスクリプトとEmacsで極力物事を済ます方針になって久しいのですが、コンテスト中は時間がないので多少のスクリプトを用意しています。先日コンテスト中にこんなスクリプトがあったらいいなあ、と考えついたので早速作成しました。
- Emacsでテキストを加工する。具体的には「入力」「空行」「出力」「空行」...という形式とする
- コマンドを適用すると、「入力」「出力」がsample1.txt、answer1.txt、...と出力される
例えば以下のようなファイルを用意します。
3 2 1 2 3 1 2 3 2 5 4 3 1 2 3 4 1 2 3 4 3 4 8
これにコマンドを適用すると、以下のようにファイルに格納されるようにしたいのです。
sample1.txt: 3 2 1 2 3 1 2 3 2 answer1.txt: 5 sample2.txt: 4 3 1 2 3 4 1 2 3 4 3 4 answer2.txt: 8
こういう処理はシェルスクリプトでさっさと書いてしまいましょう。
unpack.sh: #!/bin/sh i=1 bodyname="sample" test -f "$bodyname$i.txt" && rm -f "$bodyname$i.txt" while read l; do if ! echo "$l" | grep '^ *$' > /dev/null; then echo $l >> "$bodyname$i.txt" else if [ "$bodyname" == "sample" ]; then bodyname="answer" else bodyname="sample" i=`expr $i + 1` fi test -f "$bodyname$i.txt" && rm -f "$bodyname$i.txt" fi done
早速テストして...いい感じですね。内容もよいようです。
> ls a.txt > unpack.sh < a.txt > ls a.txt answer2.txt sample1.txt sample3.txt answer1.txt answer3.txt sample2.txt
さて、これを自分の環境にリリースしたのですが、どうも自分の操作と相性が合わないのです。どうしてだろう? と観察してみたのですが、Emacs diredから!でファイル操作ができないのが原因でした。Emacs diredから!とするにはシェルスクリプトがファイル名を引数として得て処理できなければなりません。
これを解決するのは簡単で、以下の1行を代入するだけでした。
: test -f "$bodyname$i.txt" && rm -f "$bodyname$i.txt" test -f "$1" && exec <"$1" # この行を挿入 while read l; do :
つまるところ、引数としてファイルが与えられていた場合、ファイルの内容を標準入力にリダイレクトしてしまえばよいのです。
さて、これでEmacs diredから快適にコマンドを使うことができるようになりました。今から使う機会が楽しみです。
少しばかりの続き
さて、このスクリプトであると空行が2つ以上続くとよくない結果となります。
> ls a.txt > cat a.txt 1 # sample1.txtに書き込まれるはず 2 # answer1.txtに書き込まれるはず 3 # sample2.txtに書き込まれるはず 4 # answer2.txtに書き込まれるはず > unpack < a.txt > ls # sample2.txtがない? sample3.txtがある?? a.txt answer1.txt answer2.txt sample1.txt sample3.txt
空白毎に出力するファイルを更新しているのが原因です。これを防ぐために、予めsedを適用して調整しています。
: test -f "$1" && exec <"$1" sed -ne ' s/^ *$// H $ { g s/^\n*//g s/\n*$//g s/\(\n\n\)\n*/\1/g p } ' | while read l; do :
このsedスクリプトも面白いので興味があれば是非処理を追ってみてください。具体的には、「全て読み込んでから連続する改行を取り除く。ただし、各行が空白を含む空行であればこれも処理する」といったことをしています。また、最新のスクリプトはこのGistで管理しています。
Idiom: Calculating Power of 10 using pow Function in Integer Arithmetic
今日参加したコンテストで有効な電話番号の組み合わせ数を求める問題が出題されました。実装ではlong longの大きな値、例えば\(10^{17}\)から小さな値、例えば10や1を引く演算を行う必要がありました。
powを使って10の冪乗を計算する解法も少なくなかったのですが、おそらくそのほとんどはテスト時に落ちてしまいました。どうやら計算誤差が原因のようです。
これは自分にとっても大きな驚きでした。xが整数である場合、pow(10, x)がきっちりと10の冪乗を返すことを知っていたからです。
long long x = 1; for (int i = 0; i < 18; i ++, x *= 10) assert(pow(10, i) == x); // pow(10, i)は問題なし
どこで計算誤差が出ているかを考えるために以下のような実装を用意しました。\(10^{17} - (10^0 + 10^1 + 10^2)\)を計算する実装です。期待している答えは末尾が...,999,889となるはずです。
long long c = pow(10, 17); assert(c == 100000000000000000LL); // ここは問題なし std::vector<int> a = {0, 1, 2}; for (const auto& i : a) c -= pow(10, i); assert(c == 99999999999999889LL); // ここでアサーションエラー!
変数cの内容をチェックしました。期待とは異なり、99,999,999,999,999,888でした。
さて、以下は期待通り...,999,889を返します。
long long c = pow(10, 17); assert(c == 100000000000000000LL); // ここは問題なし astd::vector<int> a = {0, 1, 2}; for (const auto& i : a) c -= static_cast<long long>(pow(10, i)); assert(c == 99999999999999889LL); // ここも問題なし!
どうも自分はこれまで、-=演算子は左側の変数の型で演算される、と思い込んでいたようです。つまり右辺の値がまずlong longにキャストされると考えていました。実際には二項の加減算演算子の振る舞いと同様に、doubleにキャストされてから計算され、改めてlong longにキャストされるようです。
これまで挙動を正確に把握していなかったこともあり、pow関数を使って10の冪乗を計算しないようにしていましたが、次のようなイディオムとして今後積極的に使っていきたいな、と考えています。またこの知見を大好きなハックにも役立てたいです。
long long y = static_cast<long long>(pow(10, x));
最後になりますが、TLで気さくに質問に答えてくださった@n_vipさん、ありがとうございました。
追記
これは実装依存の挙動だそうです。一般的に、pow(x, y)のyが整数のときはpowを使うべきではないとのこと。naoya_tさんにご指摘いただきました。どうもありがとうございます)。
参考: c++ - Definitions of sqrt, sin, cos, pow etc. in cmath - Stack Overflow
Series of Xor Arithmetics
先日参加したコンテストにこんな問題が出題されました。排他的論理和を使う問題です。
nとxが与えられる。n個の異なる整数のビットごとに排他的論理和を取ってxとなるようにしたい。n個の非負整数を出力せよ
ただし、1 ≦ n ≦ 105、0 ≦ x ≦ 105とし、出力する非負整数は106までとする
排他的論理和は融通が効く印象がありました。n - 1個の0から始まる連続した整数を用意してそのビットごとの排他的論理和を取り、最後の一個で辻褄を合わせられるんじゃないか? と考えました(その方針の元で実装した解法は無事、システムテストを通過しました!)
排他的論理和の計算は色々なアルゴリズムに登場します。例えばZobrist HashingやRolling Hashを考えるとき、排他的論理和は欠かせません。今回のエントリでは連続した排他的論理和の計算方法について紹介したいと思います。
さて、{1, 0, 0, 1, 0, 1, 0, 1}という数列の排他的論理和を計算することを考えてみましょう。左から演算を行っていくと以下のように計算することができます。
なかなか時間がかかりますね。工夫をしてもっと速く計算できるようにならないでしょうか? 排他的論理和の演算は順序を入れ替えても結果が変わりません。これは利用できないでしょうか?
試しに、先ほどの数列の隣り合った数を入れ替えてみましょう。
先ほどと同様、左から計算します。
無事同じ結果になりました(まあ答えは0か1にしかならないので、あまり説得力がないようにも思えますが...)。
もっと入れ替えれば数列を整列できそうです。具体的には0の数字の集まりと1の数字の集まりにできそうです。
さて、0と0の排他的論理和は0なのでした。つまり、0 xor 0 = 0。また、0 xor 0 xor 0は(0 xor 0) xor 0と計算すればいいので、0 xor 0となります。この結果も0。つまり、0同士の排他的論理和のグループは0と置き換えることができます。
1のかたまりについても同じように工夫ができるでしょうか? 1と1の排他的論理和は0であることを利用すると先ほどの式は以下のようにまとめることができそうです。1の個数が偶数個であれば、0となることが分かります。
1の個数が奇数個である場合も考えてみましょう。
連続した排他的論理和の計算結果は1が偶数個あれば0、1が奇数個あれば1となることが分かります。1の数を数えるだけで演算結果が分かってしまうのです。
もしこれまで、検算などをする場合に最初の図のように左から計算していたのであれば是非この方法も試してみてください。計算の速さはもちろん、正確さが格段に上がるはずです。
Idiom: Getting the largest number that is less than x
Idiom: Getting the largest number that is less than x
プログラムを書いているとき、思ったことが思った計算量でできることは分かっているのに細かい処理がパッと書けないことってありますよね。今日はこの問題を解いているときに、「std::mapのキーがx未満で最大のものを得たい!」と考え至ったのですが、すぐに繰り出すことができず、悔しい思いをしました。
コンテスト中にこういうことが起こるのはよくあります。よくあることなのであればイディオムとしてまとめられるのではないかな、と思いました。
今回のイディオムは以下です。ただし、要素は整列されているものとします。
- x未満の最大の値を得る
これはO(logN)で行うことができます。具体的には以下のようにします。
- std::lower_boundでxの下限を探す。反復子が示す要素の一つ前の値が求める値である
これを図で表せば以下となります。例では4未満の最大の値を探索しています。
要素群にx未満の値がない場合、std::lower_boundは先頭を指し示す反復子を返します。このとき、一つ前の要素にアクセスすることはできないので注意が必要です。
コンテナ独自のlower_boundメソッドを持つstd::setを使って観察してみましょう。std::setの内容は2、3、5、9、15とします。
std::set<int> set = {2, 3, 5, 9, 15};
さて、STLの反復子は要素の間を指し示すものとするとすっきりと考えることができるものも多いです。ここからは以下の図の左の表現方法を取りたいと思います。
3未満の最大の値を探してみましょう。lower_boundの結果は以下のようになるため、反復子を一つ戻したものが求める値、2になります。
auto it = set.lower_bound(3); // set: 2 3 5 9 15 // ^ // `-- it *(-- it); // 3未満の最大の値は2
次に9未満の値を探してみましょう。
auto it = set.lower_bound(9); // set: 2 3 5 9 15 // ^ // `-- it *(-- it); // 9未満の最大の値は5
7のように要素として存在しない数値の場合はどうでしょうか? その場合も問題なく動作します。
auto it = set.lower_bound(7); // set: 2 3 5 9 15 // ^ // `-- it *(-- it); // 7未満の最大の値は5
もちろん、21のように大きな値を入れてもきちんと動作します。この場合、itはstd::end(set)を指しますが特に問題はありません。
auto it = set.lower_bound(21); // set: 2 3 5 9 15 // ^ // `-- it *(-- it); // 21未満の最大の値は15
最後に先頭の値、例えば、2未満の最大の値を探す場合ですが、これは注意が必要です。itの一つ手間の要素にアクセスすることはできないからです。
auto it = set.lower_bound(2); // set: 2 3 5 9 15 // ^ // `-- it (= std::begin(set)) *(-- it); // これはアクセス違反!
そのため、lower_boundを試したあとに反復子が先頭にあるかどうか調べる必要があります。
auto it = set.lower_bound(x); if (it == std::begin(set)) { // x未満の値は存在しない } else { *(-- it); // x未満の最大の値 }
イディオムがはっきり見えてきました。aをコンテナとすると以下のように書くことができそうです。
auto it = std::lower_bound(std::begin(a), std::end(a), x); // または、auto it = a.lower_bound(x) if (it == std::begin(a)) { // x未満の値は存在しない } else { *(-- it); // x未満の最大の値 }
二分探索で失敗しないために
二分探索で失敗しないために
範囲を求める整数の二分探索を書いていると混乱するときってありますよね。でも全然書けないわけじゃないんです。すっきり書けることも多い。でも、こういうことがよく起こるんです。
- 無限ループに陥る
- 範囲を狭めるとき、下限と上限の何れを更新してよいか迷う
- 答えが一個ずれている
ごく最近幾つかのコツを掴んでようやく迷わず書けるようになってきました。今後も二分探索で失敗しないためにまとめてみることにしました。
探索の型を見極める
範囲を求める整数の二分探索は大きく2つに分けることができます。これを設計段階で見極めておくことが重要であるようです。
本エントリではそれぞれの型をlower_bound型、upper_bound型と呼びます。C++のSTLで実装されているstd::lower_boundとstd::upper_boundを参考としました。日本語にすれば、下界型(下限型)、上界型(上限型)とでも言うのでしょうか。
灰色で図示されている要素は「i番目の要素をaiとしたとき、関数f(ai)が真である」ことを表しています。下限と上限は楔(くさび)で表しています。テキストエディタで言えば、「文字と文字の間にカーソルがある」ような位置の考え方です(余談ですが、実際にテキストエディタを設計する場合はこのような位置の考え方をすると設計が簡潔になるそうです)。
続いて、位置が要素の位置にある考え方も図示します。
上の楔型と要素型の図示の間には不整合性がありますが、後で明らかになりますので今は放っておきましょう。
文章でも言い表しておきましょう。
- lower_bound型
- 関数f(ai)が真となる最小のiを探索する
- upper_bound型
- 関数f(ai)が真となる最大のiを探索する
また、二分探索はその性質上、以下のようなデータには適用できません。
関数f(x)ってなに?
まずは簡単な例で考えてみましょう。整列済みの数値の配列を例にします。
lower_bound型としては以下のような問題を考えます。
aiが3以上である最小のiを求める
図示すると以下となります。0-indexで考えたとき、2が解となります。
このとき、関数f(x)は以下のような関数となります。
bool f(int x) { return x >= 3; }
upper_bound型としては次のような問題を考えます。
aiが3以下である最大のiを求める
図示すると以下となります。0-indexで考えたとき、5が解となります。
このとき、関数f(x)は以下のような関数となります。
bool f(int x) { return x <= 3; }
初期範囲を決める
求める解の範囲が分かっている場合、言い換えれば、ある範囲の中に必ずf(ai)が真となるiがある場合は以下のような範囲を用います。ここで、l*とu*はそれぞれ解が取りうる範囲の下限と上限を示します。
図示すると以下のようになります。
範囲を扱う際は右閉半開区間を用います。例えば、i ∈ [0, 7]の範囲に解があることが分かっているのであれば、初期の範囲を(-1, 7]とします。
C++の実装で表すと次のような感じになります。
int l = -1, u = 7;
右閉半開区間を用いるのは計算機の整数演算を巧みに利用するためで、以下のような特徴があります。
- 無限ループに陥らない
- 同じ要素が複数回評価されることがない
反復の条件
二分探索は範囲をしぼり込むことによって探索を行います。範囲から要素がなくなってしまっては意味がありませんのでその直前で反復を停止します。具体的には範囲がただ一つの整数値を表している場合、つまり範囲が(x - 1, x]となったら反復を終了します。C++での実装は以下のようになります。
while (u - l > 1) { : }
後でまた触れますが、このように実装したときにf(l* - 1)やf(u*)が評価されることはありません。
範囲を更新する
二分探索は反復毎に範囲を約半分に狭めていくことによってO(logN)の計算量を実現しています。各反復ではまず中央の値を計算します。整数演算ですので、端数は切り落とされます。
中央の値を使い、f(am)を評価します。その結果によって範囲を更新します。
- mでuを更新する(更新した範囲は(l, m]となる)
- mでlを更新する(更新した範囲は(m, u]となる)
さて、この節では一先ず細かいことを抜きにしてどのように範囲の更新を行うかを注視したいと思います。そのため、要素を詳細に記した図は一度お休みし、より概念的な図を用いることとします。
範囲を更新する際は、それが全体の範囲のどこにあるか、どの程度の大きさであるか等は意識しないほうがいいです。そのため、左右を切った図とします。lとuに挟まれている区間が注目している区間です(以下の図はuppper_bound型)。
また、更新する範囲は探索の型で異なります。lower_bound型から説明します。
lower_bound型における範囲の更新
- lower_bound型であり、現在注目しているf(am)が真である場合
mは解かもしれませんが、f(ai)が真となる要素iはmよりも左側にまだあるかもしれません。ですので、左側をさらに詳しく探索します。
もちろん次の図のように、中央の値が求める解である場合もありますが、二分探索の途中でこれを意識する必要はありません。また、むしろ意識しないほうがよいようです。慣れてくると、放っておいてもアルゴリズムがこの解を絞り込んだ範囲に放り込んでくれることが分かってくるからです。
- lower_bound型であり、現在注目しているf(am)が偽である場合
求める解はmよりも右側にあります。ですので、右側をさらに詳しく探索します。
C++で実装すると以下のようになります。
if (f(m)) { u = m; } else { l = m; }
upper_bound型における範囲の更新
さて、ここからはupper_bound型を説明します。
- upper_bound型であり、現在注目しているf(am)が真である場合
mは解かもしれませんが、f(ai)が真となる要素iはmよりも右側にまだあるかもしれません。ですので、右側をさらに詳しく探索します。
範囲を右閉半開区間として取ることにしたため、mは絞り込んだ範囲に含まれません。ですが心配することはありません。lower_bound型のときと同様に、最後には辻褄がちゃんと合います。
- upper_bound型であり、現在注目しているf(am)が偽である場合
求める解はmよりも左側にあります。ですので、左側をさらに詳しく探索します。
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; } } /* 解は何? */ }
探索の様子を図示してみましょう。
うまく行きましたね! 範囲は最後に(1, 2]となりました(l = 1、u = 2)。範囲の定義からも、図からもuが答えとなることが分かります。
何故こんなにうまく行くのか考えてみましょう。
先ほどの節で説明した通り、lower_bound型の二分探索では現在注目しているf(am)が真であるとき、f(ai)が真となる要素iはmよりも左側にまだあるかもしれないので、左側をさらに詳しく探索します(図を再掲します)。
範囲を右閉半開区間として取るため、mは絞り込んだ範囲に引き続き含まれます。例えばm未満に解がない場合、それ以降どのように探索されたとしても範囲は最終的に(m - 1, m]となります。
f(am)が偽であるとき、f(ai)が真となる要素iは必ずmよりも右側にあります。そのため、右側をさらに詳しく探索します。
範囲を右閉半開区間として取るため、次の範囲(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; } } /* 解は何? */ }
探索の様子を図示してみましょう。
...こちらは微妙ですね。解がlの位置にあります。どうしてこうなったか考えてみましょう。
まずf(am)が偽であるときのことを考えてみましょう。mより大きいiのaiには解がありませんので、ばっさりと範囲から除くことができます(図は同じく、再掲です)。
さて、現在注目しているf(am)が真であるとき、それは解の候補となります。つまり、現在の要素、amは真に上限であるかもしれません。ですが、f(ai)となる要素iはmよりも右側にまだあるかもしれないのです。右側を詳しく検索し続ける必要があるはそのためです。
範囲は右閉半開区間として取るのでした。次の期間は(m, u]となります。つまり、一時的に解の候補であるmは範囲から取り除かれます。
ではamが真の上限、つまり解であったとしましょう。ai ∈ (m, u]には解がありませんよね。つまり、(m, u]のどのiの要素に対するf(ai)はすべて偽になります。
そのため、これに続く探索は常に範囲の左側を残すように範囲を絞り続けることとなり、結局のところ範囲は(m, m + 1]となります。求める解はmです。
f(am)が解の候補であるけれども真の上限ではない場合、真の上限となるようなm'を有する区間(l', u']が(m, u]の範囲の中に必ずあるはずですので、やはり右側を詳しく探索します。真の上限となるm'が見つかった場合、範囲は(m', u']を経て(m', m' + 1]となります。この場合の解はm'となります。
本当にうまくいくの?
確信を得るために色々な場合の結果を見てましょう。以下はlower_bound型の二分探索を試した結果です。lower_boun型の二分探索の結果はuを見ればよいのでした。それぞれうまく行っているようです。
upper_bound型の二分探索も試してみましょう。
upper_bound型の二分探索の結果はlを見ればよいのでした。upper_bound (1)の戻り値が-1となってしまいました。これはなんでしょうか?
よく考えてみると、何れのiを用いても関数f(ai)が真になることはありません。想定していない入力だったのでした。
また、upper_bound (8)の結果がおかしいですね。この結果はupper_bound (7)の結果と同じようです。試しに比較してみましょう。
これはおかしいですね。どうしたらよいでしょう?
よりよくするために
まずupper_bound (1)について。これはそもそもupper_bound型の二分探索の定義がおかしかったのです。私たちは以下のように定義したのでした。
- upper_bound型
- 関数f(ai)が真となる最大のiを返す
イメージしていたのは以下のような図でした(再掲)。
これを以下のように改訂しましょう。
- upper_bound型
- 関数f(ai)が真となる最大のiを探索し、その次の要素を返す
upper_bound型の二分探索のイメージは以下となります。
数値を入れたときの図は以下のようになります。
これでupper_bound (1)の解釈も通るようになりました。以前のようにupper_bound型の探索でf(ai)が真となる最大のiの要素の位置を知りたい場合はlを使うのではなく、u - 1を使うようにします。
これだとどうも分かり辛い、という方は冒頭に記載した楔型の図示をイメージするとよいと思います。
冒頭で少し触れたように、この方法は範囲を表す際によく使われる概念です。
STLのイテレータも楔型で考えると都合がよいものも多いです。std::begin()、std::end()、std::lower_bound()、std::upper_bound()、std::vector::insert()、std::vector::erase()等々...
例えばstd::vector::insert()は楔の位置に指定した要素群が挿入されると考えるとよいですし、std::vector::erase()も楔と楔の間の要素を削除すると考えるとよいかと思います。
話は少々脱線しましたが、次にupper_bound (8)について考えてみましょう。探索の初期範囲を(l* - 1, u*]ではなく、(l* - 1, u* + 1]としてみましょう。
図示すると以下のようになります。
上限を1つだけ大きくしただけですが、これだけで前の節での問題は不思議と解消されます。
もちろんlower_bound型もきちんと動作します。
それに加え、lower_bound(9)のように解が存在しないも適切に表現できるようになります。
「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(am)が真となるような図を書くのがポイントのようです。
奥深い二分探索
このエントリはCompetitive Programming Advent Calendar 2015の12/16日分のエントリとして記載しました。
私はこのエントリに記載したコツを使うことによって随分迷わなくなりましたが、二分探索の世界は実に奥深く、様々な「必勝法」が存在し、それに関する議論はいつも尽きないようです。
上位コーダーによる記事はシンプルにまとめられています。よく引用されるエントリには必ず目を通しましょう。
- 2012-07-03 - cafelier@SRM - TopCoder部
- @kinabaさんによるエントリ
- 二分探索と三分探索 - komiyamの日記
- @kou_miyamさんによるエントリ
プログラミングコンテストチャレンジブックの存在を忘れてはいけません。本エントリで記載した右閉半開区間を用いるアイデアはこの書籍を読んで得た知識です。プログラミングコンテストチャレンジブックでの二分探索のアプローチは@kinabaさんのアプローチに近いように感じています。
また、私は「ビット逐次決定型」と無理矢理名付けていますが、中央値を求めるような反復を行わず上位ビットから決定していく手法は、二分探索の議論において必ず言及されるようです。二分探索を違った視点から眺めることができますので、ご一読をお勧めします。
- @hirose_golfさんによるツィートとその実装例
- @Mi_Sawaさんによるツィートと実装例
また、@hirose_golfさんによる手法は本エントリで取ったアプローチと同様にlower_bound型とupper_bound型を設計段階で意識したアプローチとなっています。しかしながら、lower_bound型とupper_bound型を明確に区別するように紹介しましたが、そのような説明をする記事は少ないように感じます。少し考えるとupper_bound型はf(ai)の否定(つまり! f(ai))が真となる最小のiを探す問題、つまりlower_bound型に置き換えることができるからかなあと考えています。
