動的計画法をやってみよう（C++）


## はじめに
本記事は[三重大学 計算研 Advent Calendar 2019](https://qiita.com/advent-calendar/2019/calc-mie) $11$日目です。

## 動的計画法とは
動的計画法とは計算量的に難しい問題などを分割しその分割した問題を解いていき遷移を繰り返すことで答えを出すという手法のことです。

個人的なイメージとして動的計画法は、数列の漸化式のプログラムバージョンです。
実際に動的計画法を用いて数列の計算もできますし、様々な条件分岐を組み合わせることで複雑な問題を解くことができます。

今回は何問か問題を紹介して解説をしていこうと思います。
## 実際に動的計画法が使える問題を見ていきましょう。

本記事ではオーダーに含まれる全変数の最大値の積が$10^9$未満なら解けるとしています。
具体的には$O(N^2M)$で$1\leq N\leq 1000,1\leq M\leq 300$なら$1000^2*300=3*10^8<10^9$

### <span style="font-size: 150%"> **問題１** </span>
ある数列$\{A_n\}$がある。この数列は次の法則に従って構成されています。$$A[1]=1,\ A[2]=1,\ A[i]=A[i-1]+A[i-2]\ (3\leq i)$$
この数列の第$K$項の値を$10^9+7$で割ったあまりを求めてください。

**制約**
$1\leq K\leq 10000000$

**入力**
$K$

**出力**
数列の第$K$項の値を$10^9+7$で割ったあまり

これは時間計算量$O(K)$で解けます。
そのため例え$K$が$10000000$程度でも以下のコードを実行すれば十分高速に求めることが出来る。

```cpp
#include<iostream>
using namespace std;

const long long int mod=1000000007;

long long int K,A[10000001];
int main() {
    cin>>K;
    
    A[1]=1;
    A[2]=1;
    for(int i=3;i<=k;i++){
        A[i]=A[i-1]+A[i-2];//(1)
        A[i]%=mod;
    }
    cout<<A[K];
} 

```

ほとんど問題文をそのまま移しただけのコードです。
しかし$(1)$のような遷移に分割し計算しているから実はこれも動的計画法といえます。
**<参考>**
この数列は有名なフィボナッチ数列であり一般項は
$$A_k=\frac{1}{\sqrt5} \left(\left( \frac{1+\sqrt5}{2} \right)^k- \left( \frac{1-\sqrt5}{2} \right)^k\right\)$$

<br><br>
ここまでなら動的計画法も大して難しくないんじゃないかなんて思うかもしれません。
しかし、次の問題はどうでしょう？
### <span style="font-size: 150%"> **問題２** </span>
あなたは今やりたいことが$N$個あります。
しかし、暇な時間は$T$時間しかありません。
$i$番目のやりたいことは$t_i$時間かかり、やったときの満足度は$x_i$です。
満足度の和が最大になるようにやりたいことをした時の満足度の和はいくつになるでしょうか？

**制約**
$1\leq N\leq 500$
$1\leq T\leq 10000$
$1\leq t_i\leq 10000$
$1\leq x_i\leq 10^9$


**入力**
$N\ T$
$t_1\ x_1$
$t_2\ x_2$
$\ \vdots\ \ \  \vdots$
$t_N\ x_N$

**出力**
満足度の和が最大になるようにやりたいことをした時の満足度の和

どうでしょうか？一気に難易度が上がってしまいました。
これはナップサック問題と言われる問題を言い換えた問題なのですが、動的計画法に慣れていない人ならほぼほぼ解くことは出来ません。

まず愚直に考えてみましょう。
やりたいことの選び方を全部調べて最も満足度の高いものを出力すればいいです。
ただやりたいことは最大で$500$個あるのでその組み合わせは
$2^{500}=32733906078961418700131896968275991522166420460430647894832913680$
$96133796404674554883270092325904157150886684127560071009217256545885393$
$053328527589376$通りあります。
端的に言うと絶望的です。

<span style="font-size: 150%"> **どうにかしましょう。** </span>

ちょっと考えを変えてみましょう。
暇な時間を超えないようにやりたいことを何個かやる時
合計でかかる時間は何通りくらいあるでしょうか？

答えは最大でも$0,1...,T-1,T$の$T+1$通りです。
（$T$は最大でも$10000$なのでこれなら十分少ないです。）

次に今回求めるのは満足度の和の最大値なのでそれ以外の値はどうでもいいです。
この２つの考え方からうまく計算量が減らせないでしょうか？

まず二次元配列$$dp[今何個目のやりたいことを考えているか][合計でかかる時間]=満足度の合計$$
を用意します。
ここで$dp[0][i]\ (0\leq i\leq T)$の値をどうすればよいかちょっと考えてみましょう。

答えは$0$です。
当然ながら、初期値だからそれはそうだろと思う人もいると思いますが次のような考え方をしてみてはどうでしょうか。

**$0$個目のやりたいことは任意の時間かかりその満足度は$0$である。**

この考えに基づいて考えていいことは簡単に証明できます。
満足度が$0$で完全に無駄なやりたいことで満足度の和の最大値には影響を及ぼさないからです。

さてここからが動的計画法の考え方に入ります。
遷移式を先に書いてしまうと
$i$番目のやりたいことをやると決める時
$$dp[i][j]=\max(dp[i][j],dp[i-1][j-t[i]]+x[i])\ \ \  (j-t[i]\geq0)$$
$i$番目のやりたいことをやらないと決める時
$$dp[i][j]=\max(dp[i][j],dp[i-1][j])$$
です。
この式の意味を説明します。
$i$番目のやりたいことをやると決める時、$i-1$番目のやりたいことまで考えたらそれまでにやると決めたやりたいことの合計時間は$j-t[i]$で、合計満足度は$dp[i-1][j-t[i]]$でした。$i$番目のやりたいことをやると決めると合計時間は$(j-t[i])+t[i]=j$ 時間となり合計満足度は$dp[i-1][j-t[i]]+x[i]$と書けます。
もともとの$dp[i][j]$より値が大きい場合は$dp[i][j]$の値を更新します。

$i$番目のやりたいことをやらないと決める時、$i-1$番目のやりたいことまで考えたらそれまでにやると決めたやりたいことの合計時間は$j$で、合計満足度は$dp[i-1][j]$でした。$i$番目のやりたいことをやらないと決めると合計時間は変わらず$j$ 時間となり合計満足度は$dp[i-1][j]$と書けます。
もともとの$dp[i][j]$より値が大きい場合は$dp[i][j]$の値を更新します。

この遷移を$N$番目のやりたいことまで繰り返して$dp[N][T]$に入った値が答えです。

<span style="font-size: 150%"> **ちょっとまってくれ！** </span>
合計時間が暇な時間$T$ぴったりになるとは限らなくないか？そこのところ大丈夫なのか？そう思った人もいると思います。
実はこれで大丈夫です。更に理由もここに至るまでの文章の中に隠れています。
それは$dp[0][i]\ (0\leq i\leq T)$の値を決めるときに用いた
**$0$個目のやりたいことは任意の時間かかりその満足度は$0$である。**
という考えです。
例えば最も満足度が高くやりたいことを選んだ時合計時間が$T-1$だった。この時$0$個目のやりたいことにかかる時間を$1$と捉えれば$dp[0][1]$から遷移がスタートし$dp[N][T]$に最後移ります。
この考え方はコードにはほとんど影響を及ぼしませんが、ナップサック問題を解くときに起きがちな混乱（$dp[0][0]$からの遷移はわかるけどそれ以外の$dp[0][j]$は存在しないからそこからの遷移を考えるのはおかしい）を払拭出来ます。

さて、実際にコードに起こしてみましょう。


```cpp
#include<iostream>
using namespace std;


long long int N,T,t[501],x[501],dp[501][10001]={};
int main() {
    cin>>N>>T;
    
    for(int i=1;i<=N;i++){
        cin>>t[i]>>x[i];
    }
    for(int i=1;i<=N;i++){
        for(int j=0;j<=T;j++){
            //i番目のやりたいことをやると決めた時
            if(j-t[i]>=0){dp[i][j]=max(dp[i][j],dp[i-1][j-t[i]]+x[i]);}

            //i番目のやりたいことをやらないと決めた時
            dp[i][j]=max(dp[i][j],dp[i-1][j]);
        }
    }
    cout<<dp[N][T];
} 

```
時間計算量$O(NT)\ 500\times 10000=5000000$
空間計算量$O(NT)\ 500\times 10000=5000000$


解説の長さ的にもっと長いコードになってもおかしくないような気もしますがこんな短いコードで収まります。
何かちょっと寂しいですね。



これでもまだまだ動的計画法のやさしい問題です。
次の問題に行ってみましょう。
### <span style="font-size: 150%"> **問題３** </span>
あなたの前にボタンが$N$個あり、$i$番目のボタンには数字$A_i$と$C_i$が書かれています。
$i$番目のボタンを押すと得点が$A_i$点入りますがエネルギーを$C_i$消費します。
各ボタンはそれぞれ$1$回押すことが出来、エネルギー$E$を使い切るまではボタンを押すことが出来ます。
ただ最後に押したボタンだけ得点が$10$倍になります。
得点を最大化してください。

**制約**
$1\leq N\leq 3000$
$1\leq E\leq 10000$
$1\leq A_i\leq 10^9$
$1\leq C_i\leq 10000$

**入力**
$N\ \ E$
$C_1\ A_1$
$C_2\ A_2$
$\ \vdots\ \ \ \ \  \vdots$
$C_N\ A_N$

**出力**
獲得できる得点の和の最大値


最後に押したボタンの得点が$10$倍ということを除けばほとんどさっきと同じ問題と言えます。
最後に押したボタンもそのままの点数の時のコードを書いてみてください。

だいたいこうなるはずです。

```cpp
#include<iostream>
using namespace std;

long long int N,E,A[3001],C[3001],dp[3001][10001]={};
int main() {
    cin>>N>>E;
    
    for(int i=1;i<=N;i++){
        cin>>C[i]>>A[i];
    }
    for(int i=1;i<=N;i++){
        for(int j=0;j<=E;j++){
            //i番目のボタンを押す時
            if(j>=C[i]){dp[i][j]=max(dp[i][j],dp[i-1][j-C[i]]+A[i]);}

            //i番目のボタンを押さない時
            dp[i][j]=max(dp[i][j],dp[i-1][j]);
        }
    }
    cout<<dp[N][E];
} 


```

さてこれを基礎に考えていきましょう。


<span style="font-size :140%">では最後に押すボタンをどれにすればいいでしょうか？</span>
### 最も簡単な考え方
最も簡単なのは全部試すことです。
dpの配列を$$dp[今何個目のやりたいことを考えているか][合計でかかる時間][最後に押すボタン]=満足度の合計$$
としてループを回す時今何個目のやりたいことを考えているかと最後に押すボタンの番号が一緒の時は一度遷移を無視します。これは最後に答えを求める際に考えます。
さて最後に出力するものを$res$とすると
$$res=\max(res,dp[N][E-C[fin]][fin]+10*A[i])\ (1\leq fin\leq N)$$
として求めることが出来ます。
ただし、$fin< C[i]$の場合は絶対にそのボタンは押せないため無視します。


```cpp
#include<iostream>
using namespace std;


long long int N,E,A[3001],C[3001],dp[3001][10001][3001]={};
int main() {
    cin>>N>>E;
    
    for(int i=1;i<=N;i++){
        cin>>C[i]>>A[i];
    }
    for(int i=1;i<=N;i++){
        for(int j=0;j<=E;j++){
            for(int fin=1;fin<=N;fin++){
                //一旦最後に押すボタンは無視
                if(i==fin){continue;}   

                //i番目のボタンを押す時
                if(j>=C[i]){dp[i][j][fin]=max(dp[i][j][fin],dp[i-1][j-C[i]][fin]+A[i]);}

                //i番目のボタンを押さない時
                dp[i][j][fin]=max(dp[i][j][fin],dp[i-1][j][fin]);
            }
        }
    }

    long long int res=0;
    for(int fin=1;fin<=N;fin++){
        //EよりC[fin]が大きいときは押せない
        if(E>=C[fin])res=max(res,dp[N][E-C[fin]][fin]+10*A[i]);
    }
    cout<<res;
} 

```
これはかなり厳しいです。
というのは
空間計算量が$O(N^2E)　3000\times3000\times10000=9\times10^{10}$程度
時間計算量が$O(N^2E)　3000\times3000\times10000=9\times10^{10}$程度
であるからです。一応空間計算量は配列を再利用することで１つ目の要素を消すことが出来て$O(NE)$に出来ますが時間計算量はどうしようもありません。

### よりよい計算量になる考え方１

これは明らかだと思うのですが、押す予定のボタンのうち最も得点が大きいものを最後に押せばよいです。
それをどう実装すればいいでしょうか？

累積和の考え方に近いやり方をすればよいです。
まず動的計画法の配列を今回は$2$個使います。
片方は番号小さい順で見ていった時、もう片方は番号大きい順で見ていった時の配列で
$$dpl[今何個目のやりたいことを考えているか(1\rightarrow N)][合計でかかる時間]=満足度の合計$$$$dpr[今何個目のやりたいことを考えているか(N\rightarrow 1)][合計でかかる時間]=満足度の合計$$
といったふうに配列を用意します。
そして最後の処理について考えていきましょう。
左から$i$個の満足度の合計は$dpl[i][j]$です。
そして$i+1$個目のボタンを最後に押すとするとそれ以降の満足度の合計は$dpr[i+2][E-j-C[i+1]]$です。
このとき配列の二つ目が$E-j-C[i+1]$なのは全ての消費エネルギーの和が$E$を超えないようにするためです。
これを式に表すと
$$res=\max(res,dpl[i][j]+10*A[i+1]+dpr[i+2][E-j-C[i+1]])\ (0\leq i\leq N-1,0\leq j\leq E)$$
となります。

```cpp
#include<iostream>
using namespace std;

long long int N,E,A[3001],C[3001],dpl[3001][10001]={},dpr[3002][10001]={};
int main() {
    cin>>N>>E;
    
    for(int i=1;i<=N;i++){
        cin>>C[i]>>A[i];
    }
    for(int i=1;i<=N;i++){
        for(int j=0;j<=E;j++){
            //i番目のボタンを押す時
            if(j>=C[i]){dpl[i][j]=max(dpl[i][j],dpl[i-1][j-C[i]]+A[i]);}

            //i番目のボタンを押さない時
            dpl[i][j]=max(dpl[i][j],dpl[i-1][j]);
        
        }
    }
    for(int i=N;i>=1;i--){
        for(int j=0;j<=E;j++){
            //i番目のボタンを押す時
            if(j>=C[i]){dpr[i][j]=max(dpr[i][j],dpr[i+1][j-C[i]]+A[i]);}

            //i番目のボタンを押さない時
            dpr[i][j]=max(dpr[i][j],dpr[i-1][j]);
        
        }
    }

    long long int res=0;
    for(int i=0;i<=N-1;i++){
        for(int j=0;j<=E;j++){
            //i+1番目のボタンを最後に押したとき
            if(E>=j+C[i+1]){res=max(res,dpl[i][j]+10*A[i+1]+dpr[i+2][E-j-C[i+1]]);}
        }
    }
    cout<<res;
} 

```
このプログラムには二重ループまでしかなく、配列も二次元配列までです。
そのため先程の愚直な実装と比べて計算量が落ちます。
空間計算量が$O(NE)　3000\times10000=3\times10^{7}$程度
時間計算量が$O(NE)　3000\times10000=3\times10^{7}$程度


### よりよい計算量になる考え方２
こちらのほうが簡単だと思います。
よりよい計算量になる考え方１で書いたように得点が高いボタンを最後に押すのが良いです。
そのため、まずボタンを$A_i$が小さい順に並べ替えます。
その後動的計画法を用いて計算するのですが、少し変わったことをします。
まず動的計画法の配列を今回も$2$個使います。
両方とも番号小さい順で見ていきます。
$$dp1[今何個目のやりたいことを考えているか][合計でかかる時間]=満足度の合計$$$$dp2[今何個目のやりたいことを考えているか][合計でかかる時間]=満足度の合計$$
といったふうに配列を用意します。
全く同じ配列のように見えますが遷移が異なります。

$i$番目のボタンを押す時（このあともボタンを押す）
$$dp1[i][j]=\max(dp1[i][j],dp1[i-1][j-C[i]]+A[i])$$$i$番目のボタンを押さない時$$dp1[i][j]=\max(dp1[i][j],dp1[i-1][j])$$$i$番目のボタンを押す時（このあとはボタンを押さない）$$dp2[i][j]=\max(dp2[i][j],dp1[i-1][j-C[i]]+10*A[i])$$$dp2$に遷移したあとはボタンを押さない$$dp2[i][j]=\max(dp2[i][j],dp2[i-1][j])$$

$dp1$から$dp1$への遷移は問題２と変わりません。
$dp1$から$dp2$への遷移は最後のボタンを押す操作で、$dp2$から$dp2$への遷移はボタンを押せないため、ただ値を移すだけの操作です。

```cpp
#include<iostream>
#include<algorithm>
using namespace std;

long long int N,E,dp1[3001][10001]={},dp2[3001][10001]={};
pair<long long int,long long int>AC[3001];
int main() {
    cin>>N>>E;
    
    for(int i=1;i<=N;i++){
        cin>>AC[i].first>>AC[i].first;
    }
    //AとCをAについて小さい順に並び替える。
    sort(AC+1,AC+N+1);
    for(int i=1;i<=N;i++){
        for(int j=0;j<=E;j++){
            //i番目のボタンを押す時(このあともボタンを押す)
            if(j>=AC[i].second){dp1[i][j]=max(dp1[i][j],dp1[i-1][j-AC[i].second]+AC[i].first);}

            //i番目のボタンを押さない時
            dp1[i][j]=max(dp1[i][j],dp1[i-1][j]);

            //i番目のボタンを押す時(このあとはボタンを押さない)
            if(j>=AC[i].second){dp2[i][j]=max(dp2[i][j],dp1[i-1][j-AC[i].second]+10*AC[i].first);}

            //dp2に遷移したあとはボタンを押さない
            dp2[i][j]=max(dp2[i][j],dp2[i-1][j]);
            
        }
    }
    
    cout<<dp2[N][E];
} 

```

<br>









### <span style="font-size: 150%"> **問題４** </span>

あなたは魔王に挑むことになりました。
しかし魔王の謎の力によって攻撃がとどきませんでした。
ただ、あなたたち勇者は諦めません。
それは、魔王の持つ$N$種類のオーブを破壊したら攻撃が通るという伝承を知っていたためです。
さぁ今まで旅の中で身につけた$M$種類の必殺技を使ってオーブを破壊しましょう！
$i$番目の必殺技はエネルギー$E_i$を消費して$K_i$個のオーブを破壊出来ます。
その破壊できる$K_i$個のオーブは$O_{ij}$番目$\ (1\leq j\leq K_i)$のものです。
勇者も人間なのでエネルギーを使いすぎるとその次の戦いで苦労する可能性があるため出来るだけ消費エネルギーを少なく全オーブを破壊したいです。

まず今持っている必殺技で全オーブを破壊できるかを答えてください。
破壊できる場合は"Alive"、破壊できない場合は"Dead"を出力し、破壊できる場合は次の行に消費エネルギーの最低値を出力してください。

**制約**
$1 \leq N\leq 15$
$1 \leq M\leq 200$
$1\leq K_i\leq N$
$1\leq E_i\leq 10^{15}$
$1\leq O_{i1}<O_{i2}...<O_{iK_i}\leq N$

**入力**
$N\ M$
$E_1\ K_1$
$O_{11}\ O_{12}...O_{1K_1}$
$E_2\ K_2$
$O_{21}\ O_{22}...O_{2K_2}$
$\ \vdots\ \ \  \vdots$
$E_M\ K_1M$
$O_{M1}\ O_{M2}...O_{MK_1}$

**出力**
破壊できる場合は"Alive"、破壊できない場合は"Dead"を出力し、破壊できる場合は次の行に消費エネルギーの最低値



さてこの問題ですが$E_i$以外とても小さいことが気になると思います。
しかし小さいとはいえ$M$は最大$200$なので必殺技の選び方は最大$2^{200}=1606938044258990275541962092341162602522202993782792835301376$
通りもあります。

これは、まず無理です。
そこで少し視点を変えます。ある必殺技の組み合わせで破壊されるオーブは最大何通りあるでしょうか？
必殺技とは違いなんと$2^{15}=32768$通りしかありません。
これくらいなら全部列挙できそうですね。
さて、この状態をすべて持つような配列を考えてみましょう。
$$dp[201(必殺技全通り)][2][2][2][2][2][2][2][2][2][2][2][2][2][2][2]$$

冗談です。
こんな配列をいちいちやっていたら頭が爆発してしまいます。
そこで
$$dp[201(必殺技全通り)][32768]$$
のようにまとめましょう。
$2$つ目の配列の番号については$2$進数表示したときの$i$桁目が$0$なら破壊できていない状態、$1$なら破壊できている状態です。

これの遷移はすこし複雑です。
まず最小値を求める問題なので$dp[0][0]$以外巨大な数INFで埋めていこう。（$dp[0][0]$は$0$）
まず、$i$番目の必殺技について値$A[i]$を算出していきます。
$j$番目のオーブを破壊できるなら$2^j$を足す。
少し例を出すと、$1,4,6,7$番目のオーブが破壊できる時値は
$$1101001_{(2)}=105$$
となります。

先程の計算した値を用い、遷移はその必殺技を用いる時$$dp[i][j \ | A[i]]=\min(dp[i][j\ | A[i]],dp[i-1][j]+E[i]) \ (0\leq j \leq2^N-1)$$
その必殺技を用いない時$$dp[i][j]=\min(dp[i][j],dp[i-1][j]) \ (0\leq j \leq2^N-1)$$
と書けます。
ここで$+$でなく$|$を使ったのは、破壊し終わったオーブをもう一度破壊することは出来ないからです。

この問題の最終的な答えは
$$dp[M][2^N-1]$$
この時$2$つ目の配列のインデックスを$2$進数になおすと全て$1$で構成されています。つまり、全部オーブが破壊できていることを意味しています。
ただしこの$dp$の値がINFであればどこからも遷移してきていないため全てのオーブを破壊する方法がないということになるのです。

```cpp
#include<iostream>
using namespace std;

const long long int INF=1000000000000000000;

long long int N,M,E[201],K[201],O[201][15],A[201]={},dp[201][32768]={};
int main() {
    cin>>N>>M;
    
    for(int i=1;i<=M;i++){
        cin>>E[i]>>K[i];
        for(int j=0;j<K[i];j++){
            cin>>O[i][j];
            A[i]+=(1<<(O[i][j]-1));
        }
    }
    
    for(int i=0;i<=M;i++){
        for(int j=0;j<(1<<N);j++){
            dp[i][j]=INF;
        }
    }
    dp[0][0]=0;

    for(int i=1;i<=M;i++){
       for(int j=0;j<(1<<N);j++){
            dp[i][j|A[i]]=min(dp[i][j|A[i]],dp[i-1][j]+E[i]);
            dp[i][j]=min(dp[i][j],dp[i-1][j]);
        }
    }

    if(dp[M][(1<<N)-1]==INF){cout<<"Dead"<<endl;}
    else{cout<<"Alive"<<endl<<dp[M][(1<<N)-1]<<endl;}
} 


```

時間計算量は$O(2^NM)\ 2^{15}\times200=6553600$
空間計算量は$O(2^NM)\ 2^{15}\times200=6553600$

<br>

## お疲れさまでした！
動的計画法は、慣れるまで（慣れてからも）とても難しいためこの記事の内部でもわからないこともあったと思います。
ただ、この記事が少しでも理解の手助けになったとしたら嬉しいです。