おもしろい数え上げの問題

知り合いからおもしろい問題を聞いたので解いてみました。おもしろいかは人によると思いますが、数え上げ方がおもしろかったのでまとめます。

タイトルを解法のキーワードにしようと思ったんですが非常にもったいないネタバレをしてしまいそうなので先に問題を紹介します。自分で考えてみたい人は一段落してからまた見に来てください。

問題

黒い石がN個、白い石がM個あり、これをいくつかの(nonemptyな)山に分けます。2つの分け方が異なるとは山のマルチセットが異なることをいいます(山の中での黒い石と白い石の順序は考えます)。何通りの分け方がありますか。ただし答えは非常に大きくなることがあるので1000003で割ったあまりを求めてください。

制約

  • 0 <= N, M <= 50
  • N + M >= 1

もともと多項式時間で数え上げられるかな〜と思ってたところを多項式時間でできたのでオーダー正直なんでもいいやという感じです。ぼくの解法だとNとMが同じ程度として、 O(N^5 \log N)ですが、本気を出したら縮むかもしれません。最大ケースで0.9秒でした。

入力

N M

サンプル

Input

3 0

Output

3

黒い石を0, 白い石を1で表すことにすると、山の分ける方法は[0,0,0], [0,00], [000]の3通りです。[0,00]と[00,0]はマルチセットとして等しくなることに注意してください。

Input

3 3

Output

131

131通りあります。これらには[001,110]や[010,110]が含まれますが、上述したように[001,110]と[110,001]は区別されないことに注意してください。

解法

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

== ネタバレ回避用スペース ==

解法

まず重複なく数え上げることを目指します。こういう数え上げは、(挿入DPみたいに)大きい順に数えるというのが典型な気がするのでそれに従います。すなわち、高さhの山をhが大きい順に確定していきます。これから扱う山の高さをh-1以下に制限することで、高さhの山の集合を独立に考えることができます。

数え上げる関数をdfs(n, m, h, k) := n個の黒い石とm個の白い石があり、高さhの山をk個つくり、残りの山は高さh-1以下にする場合の数、と定義すると、とりあえず答えはdfs(n, m, n + m, 1) + dfs(n, m, n + m, 0)になります。

関数dfsの擬似コードを以下に示します(再帰の終了条件はちょっと複雑なのであえて省いています。詳細は最後のコードを参照)

dfs(n, m, h, k)
  ret = 0
  N = h * k    // 注目する石の数
  for (i = 0; i <= N; i++)
    if (i <= n && N - i <= m)    // 黒い石をi個使用する
      w = N個から黒い石をi個選んで高さhの山をつくる場合の数
      for (d = 0; d <= (n + m - N) / (h - 1); d++)    // 高さh-1の山をdを作る。dは最大で(n+m-N)/(h-1)個
        ret += w * dfs(n - i, m - (N - i), h - 1, d)
  return ret

さて、あとは途中に書いた「N個から黒い石をi個選んで高さhの山をつくる場合の数」を計算できればおしまいです。

再び、これを重複なく数えることを考えます。このままだとちょっと考えにくいので次のように言い換えてみましょう。

問題

n, m, kはこの問題で新しく定義した変数とします。 n (=h)個の石からなる山がm (=N/h)個ある。 このうちk個の石が黒いです。山をマルチセットとして見たとき、異なるマルチセットは何通りあるか。

これを重複なく数える際にも、先程と同じ考え方ができます。すなわち、ある量の大きい順に確定させていくという考え方です。今の場合、山の中の黒い石の数を大きい順に確定させていきながら考えます。くどいようですが、黒い石がj個含まれている山がs個あるとし、残りの山はすべて黒い石がj-1個以下だとすると、このs個の山は独立に考えることができます。挿入DPですね。

計算の流れは先程と同様なので詳細は省きますが、確定させた部分の計算で次のような問題を考えることになります。

問題

n個の石からなる山がs個あり、それぞれの山にはちょうどj個黒い石がある。山をマルチセットとみたときに何通りのマルチセットがあるか

n個のうちj個が黒い石であるような山はnCj通りあります。これらnCj通りから重複をゆるしてs個選ぶ選び方なのでこの場合の数は(nCj)Hs通りになります。 これらは事前計算することでO(1)で求まるので結局併せて以下のようなプログラムで解くことができます。

コード

constexpr int MODULO = 1000003;
constexpr int MAX = 50;

mint fact[MODULO];
mint factinv[MODULO];
mint dp[MAX + 1][MAX + 1][2 * MAX + 1][MAX + 1], dp2[MAX * MAX + 1][MAX + 1][MAX + 1][MAX + 1];
bool visited[MAX + 1][MAX + 1][2 * MAX + 1][MAX + 1], visited2[MAX * MAX + 1][MAX + 1][MAX + 1][MAX + 1];

void init() {
    fact[0] = 1;
    for (int i = 1; i < MODULO; i++) {
        fact[i] = fact[i - 1] * i;
    }
    factinv[MODULO - 1] = fact[MODULO - 1].inv();
    for (int i = MODULO - 2; i >= 0; i--) {
        factinv[i] = factinv[i + 1] * (i + 1);
    }
    assert(factinv[0] == 1);
}

mint ncr(int n, int r) {
    if (n < r) return 0;
    return fact[n] * factinv[r] * factinv[n - r];
}

mint dfs2(int n, int m, int k, int j) {
    if (j == 0) return k == 0;
    if (j * m < k) return 0;
    assert(n <= MAX * MAX && m <= MAX && k <= MAX && j <= MAX);
    if (visited2[n][m][k][j]) return dp2[n][m][k][j];
    mint ret = 0;
    if (k == 0) return 1;
    for (int s = k / j; s >= 0; s--) {
        mint q = ncr(n, j);
        ret += ncr(q.x + s - 1, s) * dfs2(n, m - s, k - j * s, j - 1);
    }
    visited2[n][m][k][j] = true;
    return dp2[n][m][k][j] = ret;
}

mint cnt(int n, int m, int k) {
    assert(k <= n * m);
    return dfs2(n, m, k, min(n, k));
}

mint dfs(int n, int m, int h, int k) {
    if (n + m == 0) return 1;
    if (h == 1) return k == n + m;
    assert(n <= MAX && m <= MAX && h <= 2 * MAX && k <= MAX);
    if (visited[n][m][h][k]) return dp[n][m][h][k];
    int N = h * k;
    mint ret = 0;
    for (int i = 0; i <= N; i++) {
        if (i <= n && N - i <= m) {
            mint w = cnt(h, N / h, i);
            for (int d = 0; d <= (n + m - N) / (h - 1); d++) {
                ret += w * dfs(n - i, m - (N - i), h - 1, d);
            }
        }
    }
    visited[n][m][h][k] = true;
    return dp[n][m][h][k] = ret;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0);
    init();

    int n, m;
    cin >> n >> m;
    assert(n <= MAX && m <= MAX);
    cout << dfs(n, m, n + m, 1) + dfs(n, m, n + m, 0) << endl;
}

ダブル挿入DP

挿入DPの中で挿入DPを使う問題を初めて見た。うむ、やっぱあまりに汎用性が高いので僕が知らなかっただけで常識なのではないかと思い始めた。

プロたちによる意見

  • 山を黒と白の個数で適当に順序付ける

dp[a][b][c][d] := 黒a白bの山まで考えたときに合計で黒c白dにする通り数

とすると素直に Ο(N^5)になりそう

(by @noshi91)