部分和問題を切る！ Part2

ぶった切れ～～。。
部分和問題なんて楽勝じゃーー。。笑

前回に引き続き部分和問題について。
今回は、動的計画法で解いてみましょう。。

sub(A') == K となるA'を列挙するのではなく、そのようなA'が存在するかどうかを提示するだけでいいというのがポイントであることは前回ちょろっと述べました。。

つまりメモリ上にキャッシュしておくべき情報は、部分集合のすべての要素ではなく、すべてのA'に対するsub(A')の和だけでいいのです！！
これちょっと分かりにくいので、ソースコードみて何を言わんとしているか理解して下さい。

まずは、これ。Aの部分集合A'をすべてキャッシュしています。
割り当てられるメモリサイズはO(2^|A|)です。

bool subsetSumAll(vector<int> set, int K) {
vector<vector<int> > subset;
vector<int> phi;
subset.push_back(phi);

forf(i, set.size()) {
 vector<vector<int> > newSubset;
 forf (j, subset.size()) {
  vector<int> subsetElm = subset[j];
  subsetElm.push_back(set[i]);
  int sum = 0;
  forf (k, subsetElm.size())
   sum += subsetElm[k];
  if (sum == K)
   return true;
  newSubset.push_back(subsetElm);
 }

 forf(j, newSubset.size())
  subset.push_back(newSubset[j]);
}

return false;
}

※一部、実装高速化のためのマクロを使用しています。

これに対し、sum(A')の値をキャッシュする場合。このとき、Aの要素はすべて自然数でかつsum(A') == KとなるA'を探せばいいことから、状態をキャッシュするために用意すべき配列サイズはO(K)です。

bool subsetSumDp(vector<int> set, int K) {
int num[K + 1];
int num_t[K + 1];

memset(num, 0, sizeof(num));
forf(i, set.size()) {

 memcpy(num_t, num, sizeof(num));
 forf(j, K+1) {
  if (!num[j])
   continue;
  if (j + set[i] > K)
   continue;
  num_t[j + set[i]] = 1;
 }
 if (set[i] < K+1)
  num_t[set[i]] = 1;
 memcpy(num, num_t, sizeof(num_t));
 if (num[K])
  return true;
}
return false;
}

上に示した2つのプログラムは同一の解を出力します。なんか騙されたような感じです。
もうちょっと踏み込んで考えてみましょう。
A={1,2,3,4,5,6,7,8,9,10}, K = 5
としましょう。

sum(A') == K となるA'を探してみましょう。これくらいなら、人間でもできますね（笑）

A' ={2, 3}
A' = {1,4}
A' = {0,5}

はい、1つの値Kに対して、sum(A') == Kを満たすA'が複数存在しますね。。
A -> φ(A)
なる写像φががうまく作用して問題をAの世界ではなく、φ(A)の世界で考えているような感じです。これが何となく騙されたようなトリックじゃないでしょうか。。

さてさて、では、上記2つの実行速度を比較してみましょう。

#define PRB_SIZE 20

int main() {
// init the problem
vector<int> set;
forf(i, PRB_SIZE)
 set.pb(rand() % 100 + 1);

int sum = 0;
forf(i, PRB_SIZE)
 sum += set[i];


// All Possible Combination Search
double start_t = gettimeofday_sec();
forf(k, sum)
 cout << "k = " << k << ": " << subsetSumAll(set, k) << endl;
cout << gettimeofday_sec() - start_t << endl;

// DP
start_t = gettimeofday_sec();
forf(k, sum)
 cout << "k = " << k << ": " << subsetSumDp(set, k) << endl;
cout << gettimeofday_sec() - start_t << endl;

// Proof Calculation
forf(k, sum)
 if(subsetSumDp(set, k) != subsetSumAll(set, k))
  cout << "The two outputs are different!!" << endl;
return 0;
}

全件検索の場合: 524.896(s)
DPの場合 : 0.189(s)

はい（笑）違いは明白です。。
メモリの制限から考えて全件検索の場合は、|A| >= 26くらいになると解くのは不可能でしょう。。
ちなみにwindowsの場合は、メモリ制限は2GBまでです。
|A| = 26の集合Aのすべての部分集合をキャッシュするために必要なリストのサイズは、

部分集合の数 * 部分集合のサイズ * intのサイズ(B)

= 2^26 * 26/2 * 4 = 3.49GB

※部分配列のサイズは、簡単のため|A|/2としています。

こりゃ、解けないわー。。
実際|A| = 26の規模の問題を解くと、以下のエラーが（笑）

This application has requested the Runtime to terminate it in an unusual way.

Please contact the application's support team for more information.

いや、Dynamic Programmingの真髄を見た感じがします。

部分和問題を切る！

DP(Dynamic Programming)の題材として、部分和問題を取り上げよう。

部分和問題はこんな感じ。

与えられた集合Aと整数Kに対して、
sum(A') = KとなるようなAの部分集合A'は存在するか？
※ここでsum(X)は集合Xの全要素の和である。

ちょっと、分かりにくいかもしれないので簡単な例題を。

A = {1, 3, 5, 12, 7}

K = 11

に対して、sum(A') = K となるようなAの部分集合A'は存在するか？

答えは、yes.

Aの部分集合A' = {1,3,7}の全要素の和は11だからである。

ということで、これをコーディングしてみよう。

普通に考えると、Aの部分集合を逐次抽出して行って、sum(A') == Kとなった時点で終了判定してループ抜ければいいのかなって感じだろう。。

しかし、Aの部分集合は全部で、2^|A|となってしまう。（|A|はAの要素数）

ということは、最悪計算量は、2^|A|である。これは、まずい。。。

とりあえず、ソースコードはこんな感じ。。

#define forf(i, n) for(int i=0; i<(int)(n); i++)

vector<int> subsetSum(vector<int> set, int K) {
vector<vector<int> > subset;
vector<int> phi;
subset.push_back(phi);

forf(i, set.size()) {
 vector<vector<int> > newSubset;
 forf (j, subset.size()) {
  vector<int> subsetElm = subset[j];
  subsetElm.push_back(set[i]);
  int sum = 0;
  forf (k, subsetElm.size())
   sum += subsetElm[k];
  if (sum == K)
   return subsetElm;
  newSubset.push_back(subsetElm);
 }

 forf(j, newSubset.size())
  subset.push_back(newSubset[j]);
}

return vector<int> ();
}

さて、ここで気付いて欲しい。

上記のコードでは、sum(A') = KとなるA'をreturnしているが、もともとの問題は、sum(A') = KとなるA'が存在するかどうかだったのではないか。。

つまり、A'を出力しなくても、A'が存在するかどうかを出力すればいいのである。（このような問題を判定問題（または、決定問題）と呼びます。PとかNPとか言ってるやつは、実は判定問題に対するクラスの名称です。なので、決定問題じゃない問題に対してPとかNPとか言うのは間違いです。）

と、PとかNPはちょっと話が脱線してしまったが、

部分和問題は、条件を満たす部分集合の存在を問うているので

その部分集合を出力する必要はない。ということは‥部分集合を生成したり保持しなくてもいいんじゃないの？？

この辺りにトリックが隠されてそうである。

この問題をDPでどのように多項式時間の解法に持ち込むか。

次回に続く・・・・・

フィボナッチ数列を求めよう！

みなさん、フィボナッチ数列をご存知だろうか？
こんな数列。

1, 1, 2, 3, 5, 8, 13, 21, 34, ....

第n項の値が、第n-1項と第n-2項の和になっているような数列だ。

この数列を一般形で書くと、

f(n) = f(n-1) + f(n-2) if n => 2,
f(0) = 1,
f(1) = 1.

である。
さて、このフィボナッチ数列の第n項目を出力する関数 long long fibonacci(int n)を作成してみよう。

long long fibonacci(int n) {
if (n == 0 || n == 1)
 return 1LL;
else
 return fibonacci(n - 1) + fibonacci(n - 2);
}

んー、まあ定義に忠実に実装すれば簡単です。

ではここで、以下のコードを実行してみましょう。

#define forf(i, n) for(int i=0; i<(int)(n); i++)

int main() {
 forf(i, 50)
  cout << i << "," << fibonacci(i) << endl;
 return 0;
}

どうです？？

面白い結果がえられると思います。

私のマシンでは上記処理が終了するまで、543秒かかりました。10分くらいかかっちゃうんですね。。
上記のやり方だと、O(2^n)になっちゃってるんですね。。

では、ちょっとやり方を変えてみましょう。

long long fibonacci2(int n) {
 long long ret = 0;
 long long x0 = 1LL;
 long long x1 = 0;

 forf (i, n) {
  ret = x0 + x1;
  x1 = x0;
  x0 = ret;
 }
 return ret;
}

どうでしょう？どうやら多項式時間で実行できそうですね。私のマシン上で先ほどと同じことをやると0.001秒で済んじゃいました。。

これが、所謂、動的計画法（Dynamic Programming）、通称DPと言われるものです。
ボトムアップ的に必要な情報を積み上げて行って解を出力しています。一方、再帰を使ったものの場合は、トップダウン的な解法になっています。
しかし、これはまだまだ序の口です。
フィボナッチ数列の場合は、初めから何も考えずにDPで解く人が多いと思いますし、パッと見、多項式時間で解けそうな感じがします。
次回は、普通に解くとO(2^n)である問題に対して、少しトリッキーな方法でDPを適用し多項式時間で解く方法を紹介します。

memset()とmemcpy()

恥ずかしながら、Cのmemset()とmemcpy()を今まで使ったことがなかった。

memset()は欠点があるみたいな話を聞いたことがあったので、｢あーじゃ要らない機能なんだー。。｣

くらいに思っていた。

今回、TopCoder editorialに投稿された解答がmemset()、memcpy()を使っていたので改めて勉強しなおした。

これは、なかなか使えそう。。。

まずは、memset()について見てみよう。

ここで問題！あなたは普段、配列を初期化する時どうしていますか？

ちなみに、宣言時にすべての要素を0に初期化する場合は、これでOKです。（disp()は確認用に標準出力へ出力するためのコードです。以下同様。）

#define forf(i, n) for(int i=0; i<(int)(n); i++)
#define disp(x) cout << #x << " : " << x << endl

using namespace std;

int main () {
 int x[10] = {0};

 forf(i, 10)
  disp(x[i]);
 return 0;
}

これは、まだ楽ですね。

じゃあ、宣言した後で、初期化するときは？

int main () {
 int x[10];

 forf(i, 10)
  x[i] = 0;

 forf(i, 10)
  disp(x[i]);
 return 0;
}

上記のようにやってる人。。。ちょっと面倒臭くないですか？

これを解決してくれるのがmemset()です。

int main () {
 int x[10];

 memset(x, 0, sizeof(x));

 forf(i, 10)
  disp(x[i]);
 return 0;
}

さてさて、では0ではなく1に初期化するときは？

int main () {
 int x[10];

 memset(x, 1, sizeof(x));

 forf(i, 10)
  disp(x[i]);
 return 0;
}

と思いきや、これを実行すると、16843009に初期化されてしまいます。ちょっと考えてみて気付きました。

16843009 = 1 + 2^{8} + 2^{16} + 3^{24}

なるほど、配列の要素を初期化するのではなく、メモリを1byte単位で初期化してるのね。。。うちのマシーンは、intが32bit(= 4 byte)なので、0x01010101に初期化されていたということです。

なるほど。。。

つなりこのmemset()という関数は、変数のサイズが1byteであるcharに使うってのが最も基本的な関数なようですね。。そうか！この要件を担保するためにcharのサイズは環境依存なしの1byteなのか！！（個人的な予想。。）

っと、ちょっと話が脱線しましたかね。

ではmemset()のさらに便利な使い方を。。2次元配列の初期化2重ループでやってませんか？そうです。これでイケちゃいます！！

int main () {
 int x[10][10];

 memset(x, 0, sizeof(x));

 forf(i, 10)
  forf(j, 10)
   disp(x[i][j]);
 return 0;
}

これでmemset()の便利さは分かりましたね！

さて、次はmemcpy()。細かい事はここまで読んでくれた人なら、多分分かってくれるはず。

memcpy()を使えば、以下のようにして2次元配列のコピーが可能です！

int main () {
 int x[3][3];
 int y[3][3] = {{1,2,3}, {4,5,6}, {7,8,9}};

 memcpy(x, y, sizeof(y));
 forf(i, 3)
  forf(j, 3)
   disp(x[i][j]);
 return 0;
}

ちなみにmemcpy()の第3引数では、コピー元から何バイト分だけコピー先へコピーするかを指定します。

さー、ここまで来て、memset()、memcpy()はどうやらかなり出来るヤツらしいということが判明しました。では、これらの関数の欠点って何？？

ちょっとネットで調べましたが、よさそうなページは見つけられませんでした。

自分なりに考えてみましたが、この問題点は、先に説明した1で初期化するつもりが0x01010101で初期化されるってことに関係してる気がします。

私が言いたいのは、intのサイズが32bitのマシンで意図的に配列の要素が0x01010101に初期化されるようなコードを書いた時に、これを16bitのマシンに移植すると、0x0101に初期化されることになる。

これが問題なのではないかな～～。と。。

自分ひとりで開発するなら問題ないけど、複数チームで複数環境で開発なんてやってるときに、この辺りを理解せずにむやみやたらにmemset()を使ってたら危険だぞ！ってことですね。

他に何か欠点をご存知の方いらっしゃいましたら、是非是非コメント頂ければ。。