アルゴリズム
最適化
競技プログラミング
動的計画法
ナップザック問題

貰う DP と配る DP、メモ化再帰、個数制限なしナップサック問題

はじめに

ついこないだの AtCoder 上のコンテストで出題された問題 AtCoder Beginner Contest 099 C 問題 - Strange Bank が、初心者向けの 300 点問題としては史上最難ではないかということで話題沸騰になりました。

  • 何も考えずに DP (動的計画法) した
  • 辺の長さが 1 の DAG なので BFS でも DP でも Dijkstra でも OK
  • 個数制限なしナップサック問題を復習しないと
  • 全探索 + Greedy でやった

といった色んな声が飛び交いました。動的計画法アルゴリズムの設計法などを具体的に学べるすごく教育的な問題だと思ったので、上記の事柄をまとめて整理することを試みます。それにしてもこの問題はコイン両替問題の 1 パターンとして位置付けられますが、色んな取り組み方ができますね。

ABC 099 C - Strange Bank 問題へのリンク

問題概要

コインの種類が

  • $1$ 円
  • $6$ 円、$6^2 (= 36)$ 円、$6^3 (= 216)$ 円、...
  • $9$ 円、$9^2 (= 81)$ 円、$9^3 (= 729)$ 円、...

で与えられる。$N$ 円をちょうど払うのに必要なコイン枚数の最小値を求めよ。

制約

  • $1 \le N \le 10^5$

数値例

  • 11 円なら、9 + 1 + 1 = 11 の 3 枚が最適です

普通の感覚だと高い硬貨から順番に払っていくのが良さそうに思えてしまいます。例えば日本円で 7801 円だと、5000 円札 + 1000 円札 2 枚 + 500 円玉 + 100 円玉 3 枚 + 1 円玉で払いたくなりますし、実際日本円ではそれで正しいです。ですが今回はそれだとダメな例があります。

14 円のとき、高い硬貨から払うと 9 + 1 + 1 + 1 + 1 + 1 = 14 の 6 枚になりますが、6 + 6 + 1 + 1 = 14 の 4 枚が最適です。

解法 1-1: トップダウンに N を降下していくメモ化再帰

「何も考えずに DP した」というのは恐らくこの解法か、もしくは次の解法 1-2 のいずれかのことを言っていると思います。この問題を見て、

  • 入力が $N$ の 1 個だけ!!!嬉しい!!!
  • しかも $O(N)$ が間に合うサイズ

という気持ちになります。入力が $N$ だけのときに、$N$ より小さい問題に落とし込んで再帰的に解くのはよくやる感じかなと思います。今回も $N$ についての答えを $f(N)$ とおくと、以下のようにして再帰的に解くことができます:

  • 最初に $1$ 円払った場合は、$f(N-1) + 1$ 枚
  • 最初に $6$ 円払った場合は、$f(N-6) + 1$ 枚
  • 最初に $36$ 円払った場合は、$f(N-36) + 1$ 枚
  • 最初に $216$ 円払った場合は、$f(N-216) + 1$ 枚
  • ...
  • 最初に $9$ 円払った場合は、$f(N-9) + 1$ 枚
  • 最初に $81$ 円払った場合は、$f(N-81) + 1$ 枚
  • ...

と場合分けして、これらの最小値をとる感じです。

$f(N) = \min(f(N-1) + 1, f(N-6) + 1, f(N-36) + 1, ..., f(N-9) + 1, f(N-81) + 1, ...)$

として求まります。注意点として、散々言われていることですが、メモ (メモ化再帰による DP) をしないと大変なことになります。似た例で、フィボナッチ数列 $f(N) = f(N-1) + f(N-2)$ を再帰で求めるときもメモしてあげないと爆発するというのはあちこちで言及されています。

計算量解析

テーブルのノード数は $N$ で、$N \le 10^5$ の場合には、1 ノードあたりの遷移数は多くても

$1, 6, 36, 219, 1296, 7776, 46656, 9, 81, 729, 6561, 59049$

という感じに $12$ 通りになります。一般には $O(\log{N})$ 通りになります。したがって遷移数は多くても $O(N\log{N})$ になります。DP の計算量は通常 $O(ノード数 + 遷移数)$ になるので、結局 $O(N\log{N})$ になります。

#include <iostream>
using namespace std;

const int MAX_N = 110000;

int memo[MAX_N]; // memo[i] が i についての答え
int rec(int n) {
    if (n == 0) return 0; // 終端条件
    if (memo[n] != -1) return memo[n]; // すでに探索済みならリターン

    int res = n; // INF の気持ち

    // 1, 6, 6^2, ..., を試す (n を超えない範囲で)
    for (int pow6 = 1; pow6 <= n; pow6 *= 6) res = min(res, rec(n - pow6) + 1);

    // 1, 9, 9^2, ..., を試す (n を超えない範囲で)
    for (int pow9 = 1; pow9 <= n; pow9 *= 9) res = min(res, rec(n - pow9) + 1);

    // メモしながらリターン
    return memo[n] = res;
}

int main() {
    int N; cin >> N;

    // テーブルを -1 で初期化
    for (int i = 0; i < MAX_N; ++i) memo[i] = -1;
    cout << rec(N) << endl;
}

解法 1-2: ボトムアップにループ回す DP

上の再帰 DP は、図にすると、以下のようなグラフを後ろから戻っていくことをやっています。

DAG.jpg

ここでは、遷移の種類として

  • 1 円払う (青線)
  • 6 円払う (緑線)
  • 9 円払う (赤線)

の 3 種類しか描いてないですが、実際は 36 円の遷移や 216 円の遷移なども含めて、12 種類の遷移が登場し得ることになります。解法 1-1 では後ろから再帰でやってみましたが、前から for 文回してボトムアップにもできます。ボトムアップにやるというのは、前から順番にグラフを辿っていく感じです。貰う DP (集める DP) と配る DP の両方で書いてみます:

貰う DP

dp[n] に遷移を集めます。

dp[n] = min(dp[n], dp[(n より小さい数)] + 1)

みたいな書き方をしてみます。図にすると下のような遷移をやっています (ノード 11 にフォーカス)。

集める DP.jpg

#include <iostream>
using namespace std;

const int MAX_N = 110000;

int dp[MAX_N]; // dp[i] が i についての答え

int main() {
    int N;
    cin >> N;

    // 初期化
    for (int i = 0; i < MAX_N; ++i) dp[i] = N; // INF の気持ち
    dp[0] = 0;

    // 貰う DP --- dp[n] に遷移を集める
    for (int n = 1; n <= N; ++n) {
        for (int pow6 = 1; pow6 <= n; pow6 *= 6) {
            dp[n] = min(dp[n], dp[n - pow6] + 1);
        }
        for (int pow9 = 1; pow9 <= n; pow9 *= 9) {
            dp[n] = min(dp[n], dp[n - pow9] + 1);
        }
    }

    cout << dp[N] << endl;
}

配る DP

dp[(n より大きい数)] = min(dp[(n より大きい数)], dp[n] + 1)

みたいな書き方をしてみます。図にすると下のような遷移をやっています (ノード 2 にフォーカス)。

配る DP.jpg

#include <iostream>
using namespace std;

const int MAX_N = 110000;

int dp[MAX_N]; // dp[i] が i についての答え

int main() {
    int N;
    cin >> N;

    // 初期化
    for (int i = 0; i < MAX_N; ++i) dp[i] = N; // INF の気持ち
    dp[0] = 0;

    // 配る DP --- dp[n] から遷移を出して行く
    for (int n = 0; n < N; ++n) {
        for (int pow6 = 1; n + pow6 <= N; pow6 *= 6) {
            dp[n + pow6] = min(dp[n + pow6], dp[n] + 1);
        }
        for (int pow9 = 1; n + pow9 <= N; pow9 *= 9) {
            dp[n + pow9] = min(dp[n + pow9], dp[n] + 1);
        }
    }

    cout << dp[N] << endl;
}

貰う DP と配る DP の比較

どちらも単に書き方の違いで、好みの問題かなと思います。僕自身は「配る DP」を愛用しています。無意識のうちに

  • for 文使ってボトムアップに書きたいときは配る DP
  • メモ化再帰使ってトップダウンに書きたいときは貰う DP

という感じで書いていることが多い気がします (ただし以前に書いたナップサック DP 記事では、貰う DP に統一して書いてみました)。配る DP のメリットとしては

  • 遷移元の配列外参照など気にしなくて良い

というのがあると思います。貰う DP のときは常に遷移先の配列外参照を気にしないといけません。もちろん配る DP でも遷移先の配列外参照は気にしないといけないのですが、多くの場合 DP テーブルに余裕を持たせているので、あまり心配しなくても済むケースも多いです。

さて、ここまで両者の違いは大差ないと書いて来たのですが、インライン DP をはじめとした以下のような DP 高速化テクを駆使する場面では、大きな違いが出るように思います。柔軟に切り替えられるようにしておくことが重要かなと思います。

  • 累積和
  • 行列累乗
  • スライド最小値
  • インライン DP (SegTree 上の DP、実家 DP とも)
  • Convex Hull Trick

DP は (多くの場合) DAG 上の最短路

DP は DAG の最短路 (と同型な問題) を解いているのだという話は、7 年前の tayama さんの記事

によって広く広まりました。今回の問題であれば、まさに

DAG.jpg

のようなグラフの最短路を求める問題とみなすことができます。そしてこうしたグラフには「閉路がない」という際立った特徴があります。DAG は Directed Acyclic Graph の略で、閉路のない有向グラフのことです。閉路がないため、どのノードから順番に見て行けばよいのかが自然に決まるのがミソです。ここで最短路について整理すると以下のようになると思います:

グラフの特性 方法 計算量 備考
一般のグラフ Bellman-Ford 法 $O(VE)$ どのノードから順番に最短路が確定するか自明じゃないので $V$ 回ループ回す
辺の重みがすべて正のグラフ Dijkstra 法 $O(E\log{V})$ どのノードから順番に最短路が確定するかは予め決まらないが、Greedy に決まる
DAG DP $O(V+E)$ どのノードから順番に最短路が確定するかが予め決まる
辺の重みがすべて 1 のグラフ BFS $O(V+E)$ いわゆる「パズルを解く最短手数」などにも有効
辺の重みがすべて 0 か 1 のグラフ 0-1 BFS $O(V+E)$

なお、「ほとんど DAG だけど厳密には DAG じゃないから 2 週ループ回す」という DP も時折見ます。以下の問題はその代表例を言えそうです:

もう 1 つの話として、DP でやっていることは必ずしも「最小値を求める (最短路)」だけでなく

  • 最小値
  • 最大値
  • 解があるか判定
  • 数え上げ
  • 確率

と色々あると思います。これらに共通する構造として、tayama さんも「分配法則」を挙げていますが、それをもっと掘り下げた記事として

があります。この記事の下の方で Dijkstra 法に要請される代数構造が closed semiring で特徴づけられることを述べられています。このことを意識して解くと面白い例として

があります。

DP とメモ化再帰との比較

これについては以下のことが言えそうです:

for 文型の DP のメリット


  1. 「貰う」書き方も「配る」書き方もできる
  2. 配列を再利用するテクニックが使えて、メモリ消費量を抑えられる場合がある
  3. 再帰呼び出しのオーバーヘッドがないので僅かに高速
  4. インライン DP をはじめとした「データ構造を活用した DP 高速化」を実現しやすい

1 についでですが、メモ化再帰では基本的には「貰う」書き方しかできないです。これに関連して 4 のような DP 高速化テクを駆使するときに、メモ化再帰だと柔軟に対応できなくなるケースも多々あります。また、AtCoder だとそのような出題はあまり見ない気がしますが、2 の配列再利用テクを使わないと MLE する問題もあります。

メモ化再帰のメリット


  1. DAG の遷移順を気にしなくて良い
  2. ごく稀にメモ化再帰によって計算量が落ちる場合がある (DP テーブル自体は二次元だが一次元的な部分しか参照しなくてよいケースなど)
  3. 終了条件を簡単に return で返せる
  4. コーナーの処理が楽な場合がある (hirosegolf さんのツイートを参照)

基本的には、for 文型の DP のときには細かく意識しなければならなかったことを意識しないで済むメリットが多いような気がします。DAG が直接与えられて DP する場合や、グラフが与えられてそれを強連結成分分解して DAG にして DP する場合など、遷移順が自明でない場合にはメモ化再帰が有効なように思います。メモ化再帰は DFS とも相性がよく、メモしながら DFS すれば自然にメモ化再帰になるので、ツリー DP なども基本的にはメモ化再帰をしていると言えます。

解法 2: BFS

前述の通り、この問題は以下のようなグラフの最短路問題です。

DAG.jpg

このグラフは DAG ですが、「辺の重みがすべて 1 のグラフ」でもあります。したがって BFS で解くこともできます。

#include <iostream>
#include <queue>
using namespace std;

const int MAX_N = 110000;

int dp[MAX_N]; // dp[i] が i についての答え

int main() {
    int N;
    cin >> N;

    // 初期化
    for (int i = 0; i < MAX_N; ++i) dp[i] = -1;
    dp[0] = 0;
    queue<int> que;
    que.push(0);

    // BFS
    while (!que.empty()) {
        int v = que.front();
        que.pop();

        for (int pow6 = 1; v + pow6 <= N; pow6 *= 6) {
            if (dp[v + pow6] == -1) {
                dp[v + pow6] = dp[v] + 1;
                que.push(v + pow6);
            }
        }
        for (int pow9 = 1; v + pow9 <= N; pow9 *= 9) {
            if (dp[v + pow9] == -1) {
                dp[v + pow9] = dp[v] + 1;
                que.push(v + pow9);
            }
        }
    }

    cout << dp[N] << endl;
}

Dijkstra 法

完全にオーバーキルですが、Dijkstra 法でも解けます。ただし計算量悪化します。自分の実装だと

  • メモ化再帰: 6 ミリ秒 ($O(N\log{N})$)
  • DP: 3 ミリ秒 ($O(N\log{N})$)
  • BFS: 4 ミリ秒 ($O(N\log{N})$)
  • Dijkstra: 13 ミリ秒 ($O(N(\log{N})^{2})$)

でした。

#include <iostream>
#include <queue>
using namespace std;

// よく愛用しています
template <class T> bool chmin(T &a, T b) { if (a > b) {a = b; return true; } return false; }

const int MAX_N = 110000;
int dp[MAX_N]; // dp[i] が i についての答え

int main() {
    int N;
    cin >> N;

    // 初期化
    for (int i = 0; i < MAX_N; ++i) dp[i] = N; // INF
    dp[0] = 0;
    priority_queue<pair<int, int>, vector<pair<int, int> >, greater<pair<int, int> > > que;
    que.push(make_pair(0, 0));

    // Dijkstra 法
    while (!que.empty()) {
        int cur = que.top().first; // 現在の距離
        int v = que.top().second; // 現在のノード
        que.pop();

        // Dijkstra 法でよくやる枝刈り
        if (cur > dp[v]) continue;

        for (int pow6 = 1; v + pow6 <= N; pow6 *= 6) {
            int nv = v + pow6; // 遷移先
            if (chmin(dp[nv], dp[v] + 1)) que.push(make_pair(dp[nv], nv));
        }
        for (int pow9 = 1; v + pow9 <= N; pow9 *= 9) {
            int nv = v + pow9; // 遷移先
            if (chmin(dp[nv], dp[v] + 1)) que.push(make_pair(dp[nv], nv));
        }
    }

    cout << dp[N] << endl;
}

解法 3: 個数制限なしナップサック DP

この問題を解く DP 解法はもう 1 つあります。しかしながら冷静に考察を進めて行くと、実は上の DP とほぼ同じになることがわかります。さて、この問題を


【問題】
$12$ 個の品物が与えられる。品物を (重さ, 価値) で表すと

(1, 1), (6, 1), (36, 1), (216, 1), ..., (9, 1), (81, 1), (729, 1), ...

で与えられる。容量 $W$ のナップサックにちょうどピッタリ重さが $W$ になるように品物を入れて行く。ただし同じ品物を何個入れてもよい。品物の価値の総和を最小化せよ。

  • $1 \le W \le 10^5$

という風に考えると個数制限なしナップサック問題 (蟻本の例題 2-3-3 でも登場)です。異なっている部分は

  • 価値の最大化ではなく、最小化である

ですが、DP 遷移式の $\max$ を $\min$ に置き換えるだけなので大きな問題ではないです。品物の個数については今回はとりあえず $n$ 種類あるものと考えておきます。個数制限なしナップサックは $n$ を品物の種類数、$W$ をナップサック容量として、$O(nW)$ で解くことができます。なお、今回は品物の重さと価値は全部与えられていますが、とりあえず $\rm{value}[i]$, $\rm{weight}[i]$ で与えられているものとして漸化式を立ててみます。

普通のナップサック

まずは普通のナップサック DP をおさらいしてみます。


$\rm{dp}[i+1][w]$ := $i$ 番目までの品物の中から重さが $w$ になるように選んだときの、価値の総和の最大値

として、この値を求めるには、以下のうち (今回は) 小さい方をとります:

  • 品物 $(\rm{weight}[i], \rm{value}[i])$ を選ぶ場合 ($w \ge \rm{weight}[i]$ の場合のみ)
    $$\rm{dp}[i+1][w] = \rm{dp}[i][w-\rm{weight}[i]] + \rm{value}[i]$$

  • 品物 $(\rm{weight}[i], \rm{value}[i])$ を選ばない場合
    $$\rm{dp}[i+1][w] = \rm{dp}[i][w]$$


個数制限なしナップサック

個数制限なしナップサックでは、上の場合分けの「品物 i を選ぶ場合」が少し変わります:


  • 品物 $(\rm{weight}[i], \rm{value}[i])$ を選ぶ場合 ($w \ge \rm{weight}[i]$ の場合のみ)
    $$\rm{dp}[i+1][w] = \rm{dp}[i+1][w-\rm{weight}[i]] + \rm{value}[i]$$

  • 品物 $(\rm{weight}[i], \rm{value}[i])$ を選ばない場合
    $$\rm{dp}[i+1][w] = \rm{dp}[i][w]$$


普通のナップサックでは、「品物 i を選ぶ」という遷移によって dp[i+1][w] の状態になったならば、その遷移前は「品物 i-1 までから重さ w-weight[i] 以下で」になります。

個数制限なしナップサックでは、「品物 i を選ぶ」という遷移によって dp[i+1][w] の状態になったとしても、その遷移前も i 番目の品物が入っているかもしれないので「品物 i までから重さ w-weight[i] 以下で」になります。

以上をまとめると、次のような DP になります:


<DP漸化式>

\rm{dp}[i+1][w] = \left\{
\begin{array}{ll}
\max (\rm{dp}[i+1][w-\rm{weight}[i]] + \rm{value}[i], \rm{dp}[i][w]) & (w \ge \rm{weight}[i]) \\
\rm{dp}[i][w] & (w < \mathbb{rm}[i])
\end{array}
\right.

<DP初期条件>
$$\rm{dp}[i][0] = 0 (i = 0, 1, \dots, n)$$


計算量

素直に $O(nW)$ なのですが、今回は $n = O(\log{N})$, $W = N$ なので $O(N\log{N})$ になります。

コード

#include <iostream>
#include <vector>
using namespace std;

// DPテーブル
const int MAX_N = 15;
const int MAX_W = 110000;
int dp[MAX_N][MAX_W]; // 12 種類だが余裕をもつ

int main() {
    int W;
    cin >> W;

    // 品物を列挙 (価値はすべて 1)
    vector<int> weight;
    for (int i = 1; i <= W; i *= 6) weight.push_back(i);
    for (int i = 9; i <= W; i *= 9) weight.push_back(i); // 1 が重複しないように (しても問題はない)
    int n = (int)weight.size();

    // DP初期条件: dp[i][0] = 0
    for (int i = 0; i < MAX_N; ++i) for (int w = 0; w < MAX_W; ++w) dp[i][w] = 1 << 29; // INF
    for (int i = 0; i <= n; ++i) dp[i][0] = 0;

    // DPループ
    for (int i = 0; i < n; ++i) {
        for (int w = 0; w <= W; ++w) {
            dp[i + 1][w] = min(dp[i + 1][w], dp[i][w]);
            if (w >= weight[i]) dp[i + 1][w] = min(dp[i + 1][w], dp[i + 1][w - weight[i]] + 1);
        }
    }

    cout << dp[n][W] << endl;
}

配列再利用

配列 dp[i][w] の i の部分をなくして使いまわす実装ができます。メモリも実装量も軽くなります。なお、配列再利用するときの w を回す方向ですが

  • 通常のナップサック: 後ろ向き
  • 個数制限なしナップサック: 前向き

になります。さて、この配列再利用したバージョンの DP を見ると、実は最初の「N について再帰的に解いた DP」とあまり違わない形をしていることに気づきます。for 文ループの i と w の順番を入れ替えると一致します。

#include <iostream>
#include <vector>
using namespace std;

// DPテーブル
const int MAX_W = 110000;
int dp[MAX_W];

int main() {
    int W;
    cin >> W;

    // 品物を列挙 (価値はすべて 1)
    vector<int> weight;
    for (int i = 1; i <= W; i *= 6) weight.push_back(i);
    for (int i = 9; i <= W; i *= 9) weight.push_back(i); // 1 が重複しないように (しても問題はない)
    int n = (int)weight.size();

    // DP初期条件: dp[i][0] = 0
    for (int w = 0; w < MAX_W; ++w) dp[w] = 1 << 29; // INF
    dp[0] = 0;

    // DPループ
    for (int i = 0; i < n; ++i) {
        for (int w = weight[i]; w <= W; ++w) {
            dp[w] = min(dp[w], dp[w - weight[i]] + 1);
        }
    }

    cout << dp[W] << endl;
}

解法 4: 全探索 + Greedy

公式解法です。Kutimoti さんの

を参考にしてもらえたらと思います。ざっくり、

  • もし「$1, 6, 6^2, ... $」だけだったら Greedy でよい (というより $N$ を $6$ 進法表記にすればよい)
  • $N$ を「$6$ 進法側」と「$9$ 進法側」に分ければよくて、分け方を全探索する

という解法です。計算量的には $O(N\log{N})$ になります。

#include <iostream>
using namespace std;

// a 進法で n の桁和
int func(int n, int a) {
    int res = 0;
    while (n > 0) {
        res += n%a;
        n /= a;
    }
    return res;
}

int main() {
    int N;
    cin >> N;

    int res = N; // INF
    for (int i = 0; i <= N; ++i) {
        int tmp = func(i, 6) + func(N - i, 9);
        res = min(res, tmp);
    }

    cout << res << endl;
}

おわりに

この問題を深く学ぶことで様々な手法を整理することができます。教育的な問題で楽しかったです。