Topcoder基本のき：全探索、再帰、メモ化をマスターする！

避けては通れない道

　この前のSRMはレジスターに間に合わなかったので参加できず。

　さて、そろそろ僕も動的計画法を勉強しなくてはいけない。div2 のレベル2でも時々使うからだ。

　動的計画法解説 [TopCoder SRM 468 Div2 250/Div1 500]

　ここで勉強し始めたのだが、Javaで書かれていることもあり、よく分からないところがあった。なんとかメモ化までできたので自分なりの解説を交えながら、コードを晒していく。なんという羞恥プレイ。

やってみるか。

　上記リンク先ではSRM468div2 level1を問題にして、全探索、再帰、メモ化、そして動的計画法を解説している。らきすたで。

王様は長い間出張していたので、できるだけ早く女王のもとに帰りたいと思っています。王様は都市 0 におり、女王は都市 N (1 <= N <= 50)に居ます。全ての都市 i (0 <= i <= N - 1) に対して陸路と空路があります。配列 roadTime, flightTime が与えられ、それぞれの i 番目の要素は都市 i から都市 i + 1 までの陸路、空路の所要時間 (1 以上 1000 以下) を表わします。しかし王国の科学力は低く、空路による移動は墜落の危険を伴います。そのため、王女は王様に空路は K (0 <= K <= N) 回までに留めるように言いました。王様が王女に会うための最短時間を求めなさい。
出題元 SRM 468 Div2 250
TopCoder Statistics - Problem Statement

コードのひな形

これからいくつかコードを載せていくが、それらをコピペしてコンパイルしても何にもならない。煩雑になるのを避けて同じ部分を書いていないためだ。先に進む前にコードのひな形てきなあれを紹介する。

#include <iostream>
#include <sstream>
#include <string>
#include <vector>
#include <algorithm>
#define lengthof(x) (sizeof(x) / sizeof(*(x)))

using namespace std;
 class RoadOrFlightEasy{
	public:
     //中身

	
// BEGIN CUT HERE
	public:
	void run_test(int Case) { if ((Case == -1) || (Case == 0)) test_case_0(); if ((Case == -1) || (Case == 1)) test_case_1(); if ((Case == -1) || (Case == 2)) test_case_2(); if ((Case == -1) || (Case == 3)) test_case_3(); if ((Case == -1) || (Case == 4)) test_case_4(); }
	private:
	template <typename T> string print_array(const vector<T> &V) { ostringstream os; os << "{ "; for (typename vector<T>::const_iterator iter = V.begin(); iter != V.end(); ++iter) os << '\"' << *iter << "\","; os << " }"; return os.str(); }
	void verify_case(int Case, const int &Expected, const int &Received) { cerr << "Test Case #" << Case << "..."; if (Expected == Received) cerr << "PASSED" << endl; else { cerr << "FAILED" << endl; cerr << "\tExpected: \"" << Expected << '\"' << endl; cerr << "\tReceived: \"" << Received << '\"' << endl; } }
	void test_case_0() { int Arg0 = 3; int Arr1[] = {4, 6, 7}; vector <int> Arg1(Arr1, Arr1 + (sizeof(Arr1) / sizeof(Arr1[0]))); int Arr2[] = {1, 2, 3}; vector <int> Arg2(Arr2, Arr2 + (sizeof(Arr2) / sizeof(Arr2[0]))); int Arg3 = 1; int Arg4 = 13; verify_case(0, Arg4, minTime(Arg0, Arg1, Arg2, Arg3)); }
	void test_case_1() { int Arg0 = 3; int Arr1[] = {4, 6, 7}; vector <int> Arg1(Arr1, Arr1 + (sizeof(Arr1) / sizeof(Arr1[0]))); int Arr2[] = {1, 2, 3}; vector <int> Arg2(Arr2, Arr2 + (sizeof(Arr2) / sizeof(Arr2[0]))); int Arg3 = 2; int Arg4 = 9; verify_case(1, Arg4, minTime(Arg0, Arg1, Arg2, Arg3)); }
	void test_case_2() { int Arg0 = 3; int Arr1[] = {4, 6, 7}; vector <int> Arg1(Arr1, Arr1 + (sizeof(Arr1) / sizeof(Arr1[0]))); int Arr2[] = {1, 2, 3}; vector <int> Arg2(Arr2, Arr2 + (sizeof(Arr2) / sizeof(Arr2[0]))); int Arg3 = 3; int Arg4 = 6; verify_case(2, Arg4, minTime(Arg0, Arg1, Arg2, Arg3)); }
	void test_case_3() { int Arg0 = 3; int Arr1[] = {1, 2, 3}; vector <int> Arg1(Arr1, Arr1 + (sizeof(Arr1) / sizeof(Arr1[0]))); int Arr2[] = {2, 3, 4}; vector <int> Arg2(Arr2, Arr2 + (sizeof(Arr2) / sizeof(Arr2[0]))); int Arg3 = 2; int Arg4 = 6; verify_case(3, Arg4, minTime(Arg0, Arg1, Arg2, Arg3)); }
	void test_case_4() { int Arg0 = 7; int Arr1[] = {50, 287, 621, 266, 224, 68, 636}; vector <int> Arg1(Arr1, Arr1 + (sizeof(Arr1) / sizeof(Arr1[0]))); int Arr2[] = {797, 661, 644, 102, 114, 452, 420}; vector <int> Arg2(Arr2, Arr2 + (sizeof(Arr2) / sizeof(Arr2[0]))); int Arg3 = 2; int Arg4 = 1772; verify_case(4, Arg4, minTime(Arg0, Arg1, Arg2, Arg3)); }

// END CUT HERE

};
// BEGIN CUT HERE
int main(){
	RoadOrFlightEasy ___test;
	___test.run_test(-1);
}
// END CUT HERE

//BEGIN CUT HERE //END CUT HEREの間の文字列はtopcoderのテンプレを設定すると、問題を開いたときに自動的にテストケースのコードを書いてくれるものだ。これを実行すると

Test Case #0...PASSED
Test Case #1...PASSED
Test Case #2...PASSED

などと書いてくれる。間違っている場合も指摘してくれる。ありがたい。

これ以降はテストケース#0の場合を考えて図を貼っている。先ほど解説サイトとは少し数値が違うので注意してほしい。

全探索

まずは全探索でKの制限がない場合で書いてみる。
ビットが立っている場合を空路とみなす、と書いてあるが、二進法の要領で陸路か空路かを選ぶことだ。C++ではべき乗の演算子がない。代わりに2^Nは「 1 << N 」で表す。「<<」は二進数を左にN個シフトすることで、十進数の1を3個左へシフトすると、00001→01000 つまり1→2^3=8　となる。

※「(i & (1 << j)」の判別。0ならr(road)、そうでないならF（Flight）。

 class RoadOrFlightEasy{
	public:
	int minTime(int N, vector <int> roadTime, vector <int> flightTime, int K) {
        int sum, min = INT_MAX;
        for (int i = 0; i < (1 << N); i++) { // (1 << N ) == 2^N
            sum = 0;
            for (int j = 0; j < N; j++) {
//                printf("i:%d:  2^j:%d",i,(1 <<j));
                if ((i & (1 << j)) != 0) {
                    sum += flightTime[j];
//                    cout << "  Flight" <<endl;
                } else {
                    sum += roadTime[j];
//                    cout << "  Road" <<endl;
                }
            }
            if (min > sum)
                min = sum;
        }
        return min;
	}
};

つぎにKの制限を設けたものを書く。

class RoadOrFlightEasy{
	public:
	int minTime(int N, vector <int> roadTime, vector <int> flightTime, int K) {
        int sum, min = INT_MAX;
        for (int i = 0; i < (1 << N); i++) { // (1 << N ) == 2^N
            if(bitCount(i) > K) continue;
            sum = 0;
            for (int j = 0; j < N; j++) {
//                printf("i:%d:  2^j:%d",i,(1 <<j));
                if ((i & (1 << j)) != 0) {
                    sum += flightTime[j];
//                    cout << "  Flight" <<endl;
                } else {
                    sum += roadTime[j];
//                    cout << "  Road" <<endl;
                }
            }
            if (min > sum)
                min = sum;
        }
        return min;
	}
    int bitCount(int val){
        int ret = 0;
        for (int i=0 ; i < 32; i++) {
            if((val & (1 << i)) != 0) {
                ret ++;
            }
        }
        return ret;
    }
};

bitCount(int val)でビットが立っている数、空路を使った回数をカウントする。それがKより大きかったらcontinue 、つまりそれ以降の命令は実行せず始めのループに戻る。

再帰法/漸化式

　つぎに再帰について説明する。再帰法とは関数がその関数自体を呼び出し、処理を行う方法のことだ。「関数が自分を呼び出したら永遠に処理が終わらないのでは」と思う人がいるが、ちゃんと終了の条件を書いているので大丈夫だ。

a[0] = 1, a[i] = a[i - 1] + 2 の例

public int rec(int N) {
	if (N == 0)
		return 1;
	return rec(N - 1) + 2;
}

　例えばrec(4)と関数を呼び出したとき、
$rec(3)+2=(rec(2)+2)+2=(rec(1)+2)+2+2=(rec(0)+2)+2+2+2=(1+2)+6=9$
となる。

　ややこしいので、再帰法は漸化式に置き換えて考えると良い。

　本題に移る。この問題を漸化式にできるか見てみよう。スタートからi市までの最短所要時間をb_iとする。この場合、初項はb_0=0となる。b_iはi市への陸路と空路のうち、かかる時間が短いほうとb_i-1との和である。つまり、

$b_0 = 0$
$b_i = min ( roadTime_i ~,~flightTime_i ) + b_{i - 1} (i > 0)$

　となる。 $roadTime_i$ はi市へ通ずる陸路の所要時間である。

using namespace std;
vector <int> r;
vector <int> f;
 class RoadOrFlightEasy{
	public:
     
	int minTime(int N, vector <int> roadTime, vector <int> flightTime, int K) {
        f = flightTime;
        r = roadTime;
        return b(N);
    }
     int b(int i){
         if (i == 0) return 0;
         cout << r[i-1] <<":" << f[i-1] <<endl;
         return b(i-1)+min(r[i-1],f[i-1]);
     }
};

※プログラム上は添字は0から始まるのでroad[i-1]となっている。空路も同様だ。注意してほしい。

※rとfはminTime関数の外で宣言する。そうするとb関数でも使うことが出来る。当たり前かもしれないけれど、僕はこれで一時間悩んだ。

　次にKの制限を設けたときを考える。このとき、添字がひとつでは対応できないので $b_{i,j}$ とする。jはあと何回空路が選択できるかを表している。複雑になったがやることは変わらない。まずは初項を考えよう。 $b_{0,j}$ は先と同様に0。 $b_{i,0}$ は空路がない場合だから $b_{i,0}=b_{i-1,j}+r_i$

　 $b_{i,j}$ は空路と陸路のうち所要時間が短いほうと表す。空路を使った場合はjがひとつ減る。
（texでスペースって\hspace{10}とかでできるんだね）

$b_{i, j} = min(b_{i-1,j} + r_i \hspace{10} ,\hspace{10} b_{i - 1, j - 1} + f_i) \hspace{10} (j > 0\hspace{5},\hspace{5} i > 0)$

　以下にコードを載せる。

using namespace std;
vector <int> r;
vector <int> f;
 class RoadOrFlightEasy{
	public:
     
	int minTime(int N, vector <int> roadTime, vector <int> flightTime, int K) {
        f = flightTime;
        r = roadTime;
        return b(N, K);
    }
     int b(int i,int j){
         if (i == 0) return 0;
         int road = b(i-1,j)+r[i-1];
         if (j > 0){
             return min(road,b(i-1,j-1)+f[i-1]);
         }else{
             return road;
         }
    }
};

メモ化

　今までの方法では計算量が2^Nであり、時間がかなり掛かってしまう。その原因は同じ項を何度も計算しているからである。ここでメモ化を行う。

　結局のところ、何度も同じ計算をしてしまうのはもったいないわけで、コンピュータの持つ記憶能力を生かして、前の計算を再利用し、無駄なく力技で計算しよう、というのが本質です。動的計画法は、「ある状態までの計算結果を記憶させる」ものであるのに対し、メモ化再帰は「ある状態以降の計算結果を記憶させる」ものとなります。

最強最速アルゴリズマー養成講座：アルゴリズマーの登竜門、「動的計画法・メモ化再帰」はこんなに簡単だった (2/5) - ITmedia エンタープライズ

　要はプログラム上にチェックリストを作ってしまい、既に計算しているならチェックリスト内にある数値を使ってしまおうというのものだ。

　今回の問題ではmemoという表を用意して、iとjそれぞれに対応する $b_{i,j}$ の値を保存しておく。表の値は-1で初期化して、参照する際に-1なら処理を続け、-1でなければその値を代入する。

　コードは以下の通り。

using namespace std;
vector <int> r;
vector <int> f;
int memo[51][51];

 class RoadOrFlightEasy{
	public:
     
	int minTime(int N, vector <int> roadTime, vector <int> flightTime, int K) {
        f = flightTime;
        r = roadTime;
        fill((int*)memo, (int*)(memo + lengthof(memo)), -1);  // すべての要素を -1 で初期化

        return b(N,K);
    }
     int b(int i,int j){
         if(memo[i][j] != -1) return memo[i][j];
         printf("memo(%d,%d)=%d\n",i,j,memo[i][j]);

         if (i == 0){
             memo[i][j] = 0;
         }else{
             int road = b(i-1,j)+r[i-1];
             if (j > 0){
                 memo[i][j] = min(road,b(i-1,j-1)+f[i-1]);
             }else{
                 memo[i][j] = road;
             }
         }
         printf("memo(%d,%d)=%d\n",i,j,memo[i][j]);
         return memo[i][j];
    }
};

実行結果を示す。

memo(3,1)=-1
memo(2,1)=-1
memo(1,1)=-1
memo(0,1)=-1
memo(0,1)=0
memo(0,0)=-1
memo(0,0)=0
memo(1,1)=1
memo(1,0)=-1
memo(1,0)=4
memo(2,1)=6
memo(2,0)=-1
memo(2,0)=10
memo(3,1)=13
Test Case #0...PASSED

　計算の流れはどうなっているのか。改めて漸化式を書いておく。

$b_{0,j}=0$
$b_{i,0}=b_{i-1,j}+r_i$
$b_{i, j} = min(b_{i-1,j} + r_i \hspace{10} ,\hspace{10} b_{i - 1, j - 1} + f_i) \hspace{10} (j > 0\hspace{5},\hspace{5} i > 0)$