D – 大ジャンプ 解説 (AtCoder Beginner Contest 011)

AtCoder二項係数,数え上げ,確率,形式的べき級数・多項式

問題へのリンク

問題概要

座標 \((0,0)\) からスタートして \(N\) 回の移動で \((X,Y)\) に到達する確率を求めたい。
1回の移動では、上下左右それぞれの方向に確率 \(\frac{1}{4}\) で \(D\) だけ移動する。

制約

  • \( 1 \leq N \leq 1000\)
  • \(1 \leq D \leq 10^9\)
  • \(-10^9 \leq X, Y \leq 10^9\)

考え方

まずは少しずつ問題を言い換えて簡単にしていきましょう。

ゴール地点 \((X,Y)\) の座標について、\(X, Y\) がマイナスになる場合があるので面倒です。対称性から符号を変えて常にプラスになるようにしても問題はありません。

また、移動距離が \(D\) のままだと考えにくいです。移動距離が 1 になるように言い換えると以下のようになります。

  • \(X, Y\) が \(D\) で割り切れなければ、ゴールできない
  • \(X, Y\) が \(D\) で割り切れれば、\(\frac{X}{D}, \frac{Y}{D}\) を移動距離を 1 としたときのゴールとする

軸ごとに分離して考える方法と、形式的べき級数による方法の2つについて解説します。

方針1: 軸ごとの移動回数を考える

x座標と y座標 について同時に考えると難しいです。各移動について x軸方向と y軸方向のどちらかにしか移動しないので、それぞれ独立に考えてあげると考えやすくなります。
x軸方向について +1 移動する回数と -1 移動する回数の差はゴール時には \(\frac{X}{D}\) となっているはずです。
y軸方向についても同じで、+1 移動する回数と -1 移動する回数の差はゴール時には \(\frac{Y}{D}\) となっているはずです。

以上を考えるとゴール時の移動回数は、\(\frac{X}{D}+\frac{Y}{D}+2a+2b=N\) とした時に以下のようになっているはずです。

  • x軸方向に +1 : \(\frac{X}{D} + a\)
  • x軸方向に -1 : \(a\)
  • y軸方向に +1 : \(\frac{Y}{D} + b\)
  • y軸方向に -1 : \(b\)

\(a\) を固定してあげると、この移動の仕方が何通りあるかは二項係数を用いて計算ができます。
移動方法は全体で \(4^N\) 通りあるので、これらを用いて以下のように確率を求めることが可能です。

$$\frac{{}_{N}\mathrm{C}_{\frac{X}{D} + a} \times {}_{N – (\frac{X}{D} + a)}\mathrm{C}_{a} \times {}_{N- (\frac{X}{D} + 2a)}\mathrm{C}_{\frac{Y}{D} + b} \times {}_{N- (\frac{X}{D} + \frac{Y}{D} + 2a + b)}\mathrm{C}_{b}} {4^N}$$

\(\frac{X}{D}+\frac{Y}{D}+2a+2b=N\) を変形すると \(a+b=\frac{N-(\frac{X}{D}+\frac{Y}{D})}{2}\) なので、\(N-(\frac{X}{D}+\frac{Y}{D})\) が 2 で割り切れないと偶奇が合わずこれもまたゴールに到達することができません。

方針2: 形式的べき級数で考える

平面上のランダムウォークを形式的べき級数で考えることもできます。

1回の移動を \((\frac{1}{4} x +\frac{1}{4} x^{-1} +\frac{1}{4} y +\frac{1}{4} y^{-1})\) という多項式に対応させると、\(N\) 回目の移動で \(\frac{X}{D}, \frac{Y}{D}\) に到達する確率は

$$ [x^{\frac{X}{D}}, y^{\frac{Y}{D}}] \left(\frac{1}{4} x +\frac{1}{4} x^{-1} +\frac{1}{4} y +\frac{1}{4} y^{-1}\right)^N $$

になります。これを上手く式変形して答えを求めてみます。文字式を因数分解すると、式が簡略化されて計算が簡単になる場合があるのですが、今回はそれを用います。

\(F = \frac{1}{4}(x + x^{-1} + y + y^{-1})\) として、これを因数分解すると

\begin{align}
F &= \frac{1}{4}(x + x^{-1} + y + y^{-1}) \\
&= \frac{1}{4}(xy)^{-1}(x^2y + y + xy^2 + x) \\
&= \frac{1}{4}(xy)^{-1}(xy+1)(x+y)
\end{align}

となります。求めたいのは \([x^{\frac{X}{D}}, y^{\frac{Y}{D}}] F^N\) です。二項定理を用いて展開してみましょう。

\begin{align}
F^N &= \frac{1}{4^N}(xy)^{-N}(xy+1)^N(x+y)^N \\
&= \frac{1}{4^N}(xy)^{-N} \left(\sum_{i}^{N} \binom{N}{i} xy^i 1^{N-i} \right) \left(\sum_{j}^{N} \binom{N}{j} x^j y^{N-j} \right) \\
&= \frac{1}{4^N}(xy)^{-N} \left(\sum_{i,j} \binom{N}{i}\binom{N}{j} x^{i+j} y^{N+i-j} \right) \\
&= \frac{1}{4^N} \left(\sum_{i,j} \binom{N}{i}\binom{N}{j} x^{i+j – N} y^{i-j}\right)
\end{align}

よって係数の計算時には、\(i+j – N = \frac{X}{D} \) かつ \(i-j = \frac{Y}{D} \) を満たすような \(i, j\) について考えればよいのでこれを解くと

  • \( i = \frac{\frac{X}{D} + \frac{Y}{D} + N}{2}\)
  • \( j = \frac{\frac{X}{D} – \frac{Y}{D} + N}{2}\)

の1通りしかありません。

以上より \(\frac{1}{4^N} \left(\binom{N}{i}\binom{N}{j} x^{i+j – N} y^{i-j}\right)\) が答えです。

解法1: 軸ごとに考える

  1. X, Y が負なら -1 をかけて正にする
  2. \(X, Y\) が \(D\) で割り切れなければ、ゴールできないので 0 を出力して終了
  3. \(N-(X/D+\frac{Y}{D})\) が 2 で割り切れなければ、ゴールできないので 0 を出力して終了
  4. \( 0 \leq a \leq \frac{N-(\frac{X}{D}+\frac{Y}{D})}{2}\) の範囲について以下を計算する
    1. \( b = \frac{N-(\frac{X}{D}+\frac{Y}{D})}{2} – a\) とする
    2. \({}_{N}\mathrm{C}_{\frac{X}{D} + a} \times {}_{N – (\frac{X}{D} + a)}\mathrm{C}_{a} \times {}_{N- (\frac{X}{D} + 2a)}\mathrm{C}_{\frac{Y}{D} + b} \times {}_{N- (\frac{X}{D}+\frac{Y}{D} + 2a + b)}\mathrm{C}_{b} ~ /~ 4^N\) を答えに加える

4 については、ループ回数が \(O(N)\) 程度で、二項係数の計算は愚直に行うと \(O(N)\) なので、計算量は全体で \(O(N^2)\) となり間に合います。

4-2 についての注意ですが、二項係数かけ算を先に行ってしまうと、値が大きくなりすぎてしまい上手く計算できない可能性があります。(double 型は \(10^{308}\) 程度までしか入らないそうです。)
\(4^N\) での割り算を上手く分配して値が大きくなりすぎないように上手く計算しましょう。

C++での実装例

#include <bits/stdc++.h>
using namespace std;

double nCk(int n, int k) {
    double ret = 1.0;
    for (int i = 0; i < k; i++) {
        ret *= (double)(n - i) / (double)(k - i);
    }
    return ret;
}

int main() {
    // 入力
    int N, D;
    int X, Y;
    cin >> N >> D;
    cin >> X >> Y;

    // ゴール地点の座標を正にする
    if (X < 0) X = -X;
    if (Y < 0) Y = -Y;

    if (X % D != 0 || Y % D != 0) {  // D の値によってはピッタリ合わない
        cout << 0 << endl;
        return 0;
    } else {
        X /= D;
        Y /= D;

        if ((N - (X + Y)) % 2 != 0) {  // 偶奇がずれて合わせることができない
            cout << 0 << endl;
            return 0;
        } else {
            double ans = 0.0;
            int res = (N - (X + Y)) / 2;
            for (int a = 0; a <= res; a++) {
                int b = res - a;

                double tmp_ans = 1.0;
                tmp_ans *= nCk(N, X + a) / pow(4.0, X + a);
                tmp_ans *= nCk(N - (X + a), a) / pow(4.0, a);
                tmp_ans *= nCk(N - (X + a + a), Y + b) / pow(4.0, Y + b);
                tmp_ans *= nCk(N - (X + Y + a + a + b), b) / pow(4.0, b);  // nCkは1
                ans += tmp_ans;
            }
            cout << setprecision(11) << ans << endl;
        }
    }
}

解法2:形式的べき級数

移動幅を 1 にするところまでは解法1と同じです。

  1. X, Y が負なら -1 をかけて正にする
  2. \(X, Y\) が \(D\) で割り切れなければ、ゴールできないので 0 を出力して終了
  3. \(\frac{X}{D} + \frac{Y}{D} + N, \frac{X}{D} – \frac{Y}{D} + N\) が 2 で割り切れなければ、ゴールできないので 0 を出力して終了
  4. 以下の計算結果が答え
    1. \( i = \frac{\frac{X}{D} + \frac{Y}{D} + N}{2}\)
    2. \( j = \frac{\frac{X}{D} – \frac{Y}{D} + N}{2}\)
    3. \(\frac{1}{4^N} \left(\binom{N}{i}\binom{N}{j} x^{i+j – N} y^{i-j}\right)\) が答えです。

3. の2で割り切れるかどうかの判定は、2つの偶奇が一致するので片方だけで十分です。

計算量は、二項係数の計算を愚直に行い \(O(N)\) になります。

C++での実装例

#include <bits/stdc++.h>
using namespace std;

double nCk(int n, int k) {
    double ret = 1.0;
    for (int i = 0; i < k; i++) {
        ret *= (double)(n - i) / (double)(k - i);
    }
    return ret;
}

int main() {
    // 入力
    int N, D;
    int X, Y;
    cin >> N >> D;
    cin >> X >> Y;

    // ゴール地点の座標を正にする
    if (X < 0) X = -X;
    if (Y < 0) Y = -Y;

    if (X % D != 0 || Y % D != 0) {  // D の値によってはピッタリ合わない
        cout << 0 << endl;
        return 0;
    } else {
        X /= D;
        Y /= D;

        if ((X + Y + N) % 2 != 0) {  // 偶奇の判定
            cout << 0 << endl;
            return 0;
        } else {
            int i = (X + Y + N) / 2;
            int j = (X - Y + N) / 2;
            double ans = nCk(N, i) / pow(2, N) * nCk(N, j) / pow(2, N);
            cout << setprecision(11) << ans << endl;
        }
    }
}