ICPC 系列でよくありそうな雰囲気の問題！

問題へのリンク

問題概要

二次元平面上で原点から出発し、 $N$ 回にわたって、以下の動きのいずれかを確率 1/4 で選択して実施する。 $N$ 回の動き終了後に座標 $(X, Y)$ にいる確率を求めよ (許容誤差は $10^{-9}$ )。

x 座標を $+D$ する
x 座標を $-D$ する
y 座標を $+D$ する
y 座標を $-D$ する

制約

$1 \le N \le 1000$
$1 \le D \le 10^{9}$
$-10^{9} \le X, Y \le 10^{9}$
入力はすべて整数

考えたこと

まず $X, Y$ が $D$ の倍数でない場合は明らかに到達不可能なので確率は 0.0 となる。 $X, Y$ が $D$ の倍数の場合は、 $X, Y$ をそれぞれ $D$ で割っておけば、以下の 4 択を繰り返して座標 $(X, Y)$ を目指す問題になる。

x 座標を $+1$ する
x 座標を $-1$ する
y 座標を $+1$ する
y 座標を $-1$ する

さて、以下の問題と関連する話なのだが

$|X| + |Y|$ が $N$ より大きい場合は到達不可能
$|X| + |Y|$ が $N$ 以下であっても、 $N - |X| - |Y|$ が奇数の場合は到達不可能

となる。

atcoder.jp

それ以外の場合は到達可能なので、確率は 0.0 より大きな値になる。以降、確率が 0.0 より大きな値になる場合のみ考えていくことにする。

最初に考えたくなる解法として DP があると思う。つまり

dp[ n ][ x ][ y ] := n 回移動後に座標 (x, y) にいる確率

としてあげればできそう。しかしこれでは $O(N^{3})$ の計算量となってしまって間に合わない。何か工夫を考える必要がある。

上下左右に何回ずつ動けばいいのかを考える

$X, Y$ が負の場合は、対称性から、 $X$ を $|X|$ で置き換えて、 $Y$ を $|Y|$ で置き換えても変わらないので、そのようにする。少し考えやすくなる。

$N - X - Y$ は偶数なので、これを $2m$ ( $m$ は整数) とおく。このとき、

右への移動が $X$ 回
上への移動が $Y$ 回

という移動に対して、「左へ 1 回、右へ 1 回」と「上へ 1 回、下へ 1 回」の移動を合計 $m$ 個挟む必要がある。「左へ 1 回、右へ 1 回」をする回数を $i (= 0, 1, \dots, m)$ とすると、次のようになる

右への移動が $X + i$ 回
左への移動が $i$ 回
上への移動が $Y + m-i$ 回
下への移動が $m-i$ 回

これらの合計 $N (= X + Y + 2m)$ 回の移動を並べ替える問題となる。その場合の数は

$C(N, X+Y+m) \times C(X+Y+m, X+i) \times C(m, i)$

で与えられる (先に「右・上」と「左・下」とに分けた)。これを各 $i$ について足し合わせれば OK。計算量は $O(N)$ 。

コード

double 型は $10^{306}$ 程度の値までしか格納できないので、

$C(X+Y+m, X+i) \times C(m, i)$

を合計して、 $4^{N}$ で割ってから、 $C(D, X+Y+m)$ を掛けることにした。それで通った。

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

const int MAX_C = 1100;
long double Com[MAX_C][MAX_C];

void calc_com() {
    for (int i = 0; i < MAX_C; ++i) for (int j = 0; j < MAX_C; ++j) Com[i][j] = 0.0;
    Com[0][0] = 1.0;
    for (int i = 1; i < MAX_C; ++i) {
        Com[i][0] = 1.0;
        for (int j = 1; j < MAX_C; ++j) {
            Com[i][j] = (Com[i-1][j-1] + Com[i-1][j]);
        }
    }
}

long double solve(int N, int D, int X, int Y) {
    calc_com();
    long double all = 1;
    for (int i = 0; i < N; ++i) all *= 4;
    X = abs(X), Y = abs(Y);
    if (X % D != 0) return 0.0;
    if (Y % D != 0) return 0.0;
    X /= D, Y /= D;
    if (N < X + Y || (N - X - Y) % 2 == 1) return 0.0;

    int m = (N - X - Y) / 2;
    long double res = 0;
    for (int i = 0; i <= m; ++i) {
        res += Com[X + Y + m][X + i] * Com[m][i];
    }
    return res / all * Com[N][X + Y + m];
}

int main() {
    int N, D, X, Y;
    cin >> N >> D >> X >> Y;
    cout << setprecision(15) << solve(N, D, X, Y) << endl;
}