yukicoder No.659 徘徊迷路

yukicoderはとても久しぶりに解きました。最近競プロの進捗を全く生やせていないので、「3月中に毎日1問解いて黄色になろう！」キャンペーンを始めます。ってことで毎日こちらの更新を目指しますが、次の水木金とサークルの大会があるので無理です。

問題概要

ランダムウォークするので、Tターン後に丁度ゴールにいる確率を求めよ。

R、C≦10で、T≦10⁸

No.659 徘徊迷路 - yukicoder

解説

ネックは完全にTで、1ターンずつシミュレーションは無理です。で、こいつが2ⁿの形だった時を考えると、
nxt[i][j][a][b]:(i,j)→(a,b)に移動する確率
とか置いておくと、こいつが2ⁿターン後の遷移状態を示しているとき（最初はn=0で、通常のシミュレーションと同じ）、
nxt[i][j][a][b]=nxt[i][j][x][y]*nxt[x][y][a][b]
とかすると2ⁿ⁺¹ターン後の状態が分かるわけです。

で、これを使うと、1回のループでTの一番上のビットをO(R³C³logT)で消せます。で、これをTのビットが立っている個数≦logT回繰り返すので、全体としてO(R³C³(logT)²)で解けます。最大ケースを当てはめるとこれで間に合いますね。
まあnxtは前計算出来て、それでlogTが一つ落とせるんですが、どうでもいいので残してます。

で、これで解けるはずなんですが、この計算ではものすごく浮動小数点の計算を行っているので、誤差がエグイです。という訳で適切にepsを突っ込む必要があるんですが、これが難しい。足す位置をnxt←nxt2とret←tmp、つまり更新後の配列を移動させるタイミングに指定します（これは完全に適当）。また、このタイミングだと最大T回epsを足していると考えられるので、許容範囲をepsだけで超えないように10^-6×10^-8より小さい10^-15で計算をしています。ここで考えたのは2点で、

• epsがあまりにも小さくなり過ぎないこと

• それなりにepsを足していること

です。普通のdoubleだと有効数字が足りないので、long doubleを使いましょう。

https://yukicoder.me/submissions/241119

int r, c, t;
int sx, sy, gx, gy;
vector<string> mp;
// 確率
ld ret[10][10],tmp[10][10];
// 行先(i,j)->(k,l)の遷移確率（何ターンで移動するかは更新していく）
ld nxt[10][10][10][10], nxt2[10][10][10][10];

bool isInside(int y, int x) {
    if (y < 0 || r <= y)   return false;
    if (x < 0 || c <= x)   return false;
    return true;
}

int main() {
    cin >> r >> c >> t;
    cin >> sy >> sx;
    cin >> gy >> gx;
    mp.resize(r);
    rep(i, r)   cin >> mp[i];
    ret[sy][sx] = 1.0;
    while (t > 0) {
        // nxtの初期化
        rep(i, r)   rep(j, c){
            if (mp[i][j] == '#') continue;
            rep(k, r)   rep(l, c)   nxt[i][j][k][l] = 0.0;
            int cnt = 0;
            rep(k, 4) {
                if (!isInside(i + dy[k], j + dx[k]))  continue;
                if (mp[i + dy[k]][j + dx[k]] == '#') continue;
                cnt++;
            }
            if (cnt == 0)    nxt[i][j][i][j] = 1.0;
            else {
                rep(k, 4) {
                    if (!isInside(i + dy[k], j + dx[k]))  continue;
                    if (mp[i + dy[k]][j + dx[k]] == '#') continue;
                    nxt[i][j][i + dy[k]][j + dx[k]] = 1.0 / cnt;
                }
            }
        }
        // 最初に1ターン進める（nxtの更新はしない）
        int add = 1;
        rep(i, r)   rep(j, c) {
            if (ret[i][j] < eps)   continue;
            rep(k, 5) {
                if (!isInside(i + dy[k], j + dx[k]))  continue;
                tmp[i + dy[k]][j + dx[k]] += ret[i][j] * nxt[i][j][i + dy[k]][j + dx[k]];
            }
        }
        rep(i, r)   rep(j, c) {
            ret[i][j] = tmp[i][j];
            tmp[i][j] = 0.0;
        }

        while (t - add * 2 >= 0) {
            // nxtに沿って確率の更新
            rep(i, r)   rep(j, c){
                if (ret[i][j] < eps)   continue;
                rep(k, r)   rep(l, c) {
                    tmp[k][l] += ret[i][j] * nxt[i][j][k][l];
                }
            }
            rep(i, r)   rep(j, c) {
                ret[i][j] = tmp[i][j] + eps;
                tmp[i][j] = 0.0;
            }
            add <<= 1;
            // nxtの更新
            rep(i, r)   rep(j, c) rep(k, r) rep(l, c) {
                if (nxt[i][j][k][l] < eps) continue;
                rep(a, r)   rep(b, c) {
                    nxt2[i][j][a][b] += nxt[i][j][k][l] * nxt[k][l][a][b];
                }
            }
            rep(i, r)   rep(j, c)   rep(k, r)   rep(l, c) {
                nxt[i][j][k][l] = nxt2[i][j][k][l] + eps;
                nxt2[i][j][k][l] = 0.0;
            }
        }
        t -= add;
    }
    printf("%.15Lf\n", ret[gy][gx]+eps);
    return 0;
}

感想

アルゴリズム自体は簡単だったので、自分にとっては完全な誤差ゲーでしたね…もう1、2問★3をやってみて、簡単そうなら★3.5に行こうと思います。