AOJ 2336. Spring Tiles

,

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2336

なんだか楽しかったです。

solution

期待値に関する再帰的な等式に帰着し、ある種の全探索により求める。

各座標に居るときのゴールまで移動するのに必要な期待値、バネを踏んだときの期待値を求める。 バネで移動しうる座標の集合$X$、座標$x \in X$に関しそこに飛ぶ確率$p_x$とその座標での期待値$e_x$をおく。 バネ*を踏んだ後の期待値$E$は、$E = \Sigma_{x \in X} (p_x e_x)$である。 バネはどれを踏んでも結果は同じであり、区別はない。

全ての床に等確率で飛ばされるので、$p_x$はどれも$\frac{1}{|X|}$である。 $e_x$に関して、 座標$x$からゴールgまでの距離$g_x \in \mathbb{N}_+ \cup \{ \infty \}$、 座標$x$から最寄りのバネまでの距離$s_x \in \mathbb{N}_+ \cup \{ \infty \}$を使って、 $e_x = \min \{ g_x, s_x + E \} \in \mathbb{R}$と書ける。 ここで$e_x$の値は座標の$g_x, s_x$にのみに依存するので、同じ$(g_x, s_x)$を持つ座標は纏めてしまうと後の計算が速くなる。このとき$p_x$は適当にする。

$E$に関する再帰的な式$E = \Sigma_{x \in X} \min \{ p_x g_x, p_x (s_x + E) \}$が得られた。 $E$以外の全ての変数は求まっている。 代数的に解ける形$E = a E + b$に簡略化したいが、$\min$が問題である。 そこで、式中で$g_x$と$s_x + E$のどちらを使うかの決定に$E$が使われているが、これを$E$の値を適当に仮定することにより除去する。 $E$の値を仮に$e$と固定すると、$E = \Sigma_{x \in X} ( p_x \cdot ( \text{if}\; g_x \le s_x + e \;\text{then}\; g_x \;\text{else}\; s_x + E ) )$とでき、これは一次式$E = a E + b$の形と等価である。 $E = a E + b$であれば$E = \frac{b}{1 - a}$と容易に解け、$\frac{b}{1 - a}$を最小化するような$e$を探せばよい。 この仮に固定する$e$の値は、$g_x, s_x$の値の種類数程度でよい。

$E$の存在性に関する証明等はしていません。

implementation

  • doubleを使うと精度で死ぬ。long doubleを使う。
  • 答えはINT_MAXより大きくなるので、int infを使い回すと死ぬ。
#include <iostream>
#include <vector>
#include <map>
#include <queue>
#include <tuple>
#include <functional>
#include <cstdio>
#include <cmath>
#define repeat(i,n) for (int i = 0; (i) < (n); ++(i))
template <class T> bool setmax(T & l, T const & r) { if (not (l < r)) return false; l = r; return true; }
template <class T> bool setmin(T & l, T const & r) { if (not (r < l)) return false; l = r; return true; }
using namespace std;
const int inf = 1e9+7;
const int dy[] = { -1, 1, 0, 0 };
const int dx[] = { 0, 0, 1, -1 };
int main() {
    // input
    int w, h; cin >> w >> h;
    vector<string> c(h); repeat (y,h) cin >> c[y];
    // modify input
    int sy, sx, gy, gx;
    repeat (y,h) repeat (x,w) {
        if (c[y][x] == 's') {
            sy = y;
            sx = x;
            c[y][x] = '.';
        } else if (c[y][x] == 'g') {
            gy = y;
            gx = x;
        }
    }
    // prepare
    auto on_field = [&](int y, int x) { return 0 <= y and y < h and 0 <= x and x < w; };
    typedef queue<pair<int,int> > points_queue;
    auto bfs = [&](function<void (points_queue &)> init, function<void (points_queue &, int, int, int, int)> update) {
        points_queue que;
        init(que);
        while (not que.empty()) {
            int y, x; tie(y, x) = que.front(); que.pop();
            repeat (i,4) {
                int ny = y + dy[i];
                int nx = x + dx[i];
                if (not on_field(ny, nx)) continue;
                if (c[ny][nx] != '.') continue;
                update(que, y, x, ny, nx);
            }
        }
    };
    vector<vector<int> > goal(h, vector<int>(w, inf));
    bfs([&](points_queue & que) {
        goal[gy][gx] = 0;
        que.push(make_pair(gy, gx));
    }, [&](points_queue & que, int y, int x, int ny, int nx) {
        if (goal[ny][nx] == inf) {
            goal[ny][nx] = goal[y][x] + 1;
            que.push(make_pair(ny, nx));
        }
    });
    vector<vector<int> > jump(h, vector<int>(w, inf));
    bfs([&](points_queue & que) {
        repeat (y,h) repeat (x,w) if (c[y][x] == '*') {
            jump[y][x] = 0;
            que.push(make_pair(y, x));
        }
    }, [&](points_queue & que, int y, int x, int ny, int nx) {
        if (jump[ny][nx] == inf) {
            jump[ny][nx] = jump[y][x] + 1;
            que.push(make_pair(ny, nx));
        }
    });
    map<pair<int,int>,int> freq; // frequency
    int total = 0;
    int max_goal = 0;
    repeat (y,h) repeat (x,w) if (c[y][x] == '.') {
        freq[make_pair(goal[y][x], jump[y][x])] += 1;
        total += 1;
        if (goal[y][x] < inf) setmax(max_goal, goal[y][x]);
    }
    // calc
    long double e = INFINITY;
    repeat (estimate, max_goal + 1) {
        // E = f(E) = aE + b
        long double a = 0;
        long double b = 0;
        for (auto it : freq) {
            int g, s; tie(g, s) = it.first;
            int cnt = it.second;
            long double p = cnt /(long double) total;
            if (g == inf) {
                a += p;
                b += p * s;
            } else if (s == inf) {
                b += p * g;
            } else {
                if (g <= s + estimate) {
                    b += p * g;
                } else {
                    a += p;
                    b += p * s;
                }
            }
        }
        setmin(e, b / (1 - a));
    }
    // output
    long double ans
        = goal[sy][sx] == inf ? jump[sy][sx] + e
        : jump[sy][sx] == inf ? goal[sy][sx]
        : min<long double>(goal[sy][sx], jump[sy][sx] + e); // the answer can become grater than `int inf`, so conditional op. is required.
    printf("%.12Lf\n", ans);
    return 0;
}