AOJ 2336. Spring Tiles
なんだか楽しかったです。
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;
}