bzoj2707[SDOI2012]走迷宫 关于一类图上有环线性动态规划与高斯消元的问题研究

bzoj2707: [SDOI2012]走迷宫

Description

Morenan被困在了一个迷宫里。迷宫可以视为N个点M条边的有向图，其中Morenan处于起点S，迷宫的终点设为T。可惜的是，Morenan非常的脑小，他只会从一个点出发随机沿着一条从
4000
该点出发的有向边，到达另一个点。这样，Morenan走的步数可能很长，也可能是无限，更可能到不了终点。若到不了终点，则步数视为无穷大。但你必须想方设法求出Morenan所走步数的期望值。

Input

第1行4个整数，N,M,S,T

第[2, M+1]行每行两个整数o1, o2，表示有一条从o1到o2的边。

Output

一个浮点数，保留小数点3位，为步数的期望值。若期望值为无穷大，则输出”INF”。

数据

【样例输入1】

6 6 1 6

1 2

1 3

2 4

3 5

4 6

5 6

【样例输出1】

3.000

【样例输入2】

9 12 1 9

1 2

2 3

3 1

3 4

3 7

4 5

5 6

6 4

6 7

7 8

8 9

9 7

【样例输出2】

9.500

【样例输入3】

2 0 1 2

【样例输出3】

INF

【数据范围】

N<=10000

M<=1000000

保证强连通分量的大小不超过100

另外，均匀分布着40%的数据，图中没有环，也没有自环

模型一般化

第一次发现高斯消元的神奇应用。

由于这个图上的问题并不是DAG，是有环的，所以我们肯定要先tarjan强连通分量缩点。

我们考虑一个环上的点u

f[u] = ∑f[v]+w[v][u]

如果一个环的大小为cnt，那么我们会有cnt个这样的方程。所以我们可以利用高斯消元解决这类问题。

例题分析

这题是一道比较折腾的题。

首先我们要判断会不会走到“死胡同”

先dfs搜索，遍历点和连通块。

如果起点走不到终点，那就肯定不合法。

除了终点所在的连通块之外，只要这个连通块被搜到过且没有出边，那就是死胡同，也不合法。

还有有一个雕虫小技，就是按照拓扑逆序进行DP，因为对于一个拓扑逆序的东西，它的所有出边一定都已经计算完了，所以是正确的。

那么如果是这样，方程也得逆着退咯

f[u]=∑(f[v]+1)∗out[u]

移项高斯消元即可。

剩下就是码农的日常了。这题的确有点小恶心~细节注意下呗

代码

#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const int N = 11000, M = 1100000;
const double eps = 1e-8;
int read() {
char ch = getchar(); int x = 0;
while(ch < '0' || ch > '9') ch = getchar();
while(ch >= '0' && ch <= '9') {x = (x << 1) + (x << 3) + ch - '0'; ch = getchar();}
return x;
}
int to[M], nxt[M], pre
, reto[M], renxt[M], repre[M], top, dtop;
double g[110][110], out
, f
;
void add(int u, int v) {
to[++top] = v; nxt[top] = pre[u]; pre[u] = top;
reto[top] = u; renxt[top] = repre[v]; repre[v] = top;
}
int dfn
, low
, st
, belong
, block
[110], siz
, num
, bout
, cnt;
int U[M], V[M], q
, n, m, s, t;
bool vis
, bvis
;

void tarjan(int u) {
dfn[u] = low[u] = ++dtop; st[++top] = u;
for(int i = pre[u]; i; i = nxt[i]) {
if(!dfn[to[i]]) {
tarjan(to[i]);
low[u] = min(low[u], low[to[i]]);
}
else if(!belong[to[i]]) low[u] = min(low[u], dfn[to[i]]);
}
if(low[u] == dfn[u]) {
++cnt;
for(int x;st[top + 1] != u; block[cnt][num[x] = ++siz[cnt]] = x)
belong[x = st[top--]] = cnt;
}
}

void dfs(int u) {
vis[u] = bvis[belong[u]] = true;
for(int i = pre[u]; i; i = nxt[i])
if(!vis[to[i]]) dfs(to[i]);
}

void gauss(int n) {
for(int i = 1;i <= n; ++i) {
int k = i;
for(int j = i + 1;j <= n; ++j)
if(fabs(g[j][i]) > fabs(g[k][i]))
k = j;
for(int j = 1;j <= n + 1; ++j) swap(g[i][j], g[k][j]);
for(int j = 1;j <= n; ++j)
if(fabs(g[j][i]) > eps && j != i){
double temp = g[j][i] / g[i][i];
for(int k = 1;k <= n + 1; ++k) g[j][k] -= g[i][k] * temp;
}
}
}

int main() {
n = read(); m = read(); s = read(); t = read();
for(int i = 1;i <= m; ++i) {
U[i] = read(), V[i] = read();
add(U[i], V[i]); ++out[U[i]];
}
for(int i = 1;i <= n; ++i) out[i] = 1 / out[i];
top = 0;
for(int i = 1;i <= n; ++i) if(!dfn[i]) tarjan(i);
dfs(s);
if(!vis[t]) return 0 * puts("INF");
for(int i = 1;i<= m; ++i)
if(belong[U[i]] != belong[V[i]])
++bout[belong[U[i]]];
for(int i = 1;i <= cnt; ++i)
if(bvis[i] && i != belong[t] && !bout[i])
return 0 * puts("INF");
int head = 0, tail;
q[tail = 1] = belong[t];
while(head < tail) {
int cur = q[++head];
memset(g, 0, sizeof(g));
for(int i = 1;i <= siz[cur]; ++i) {
int u = block[cur][i]; g[i][i] = 1;
g[i][siz[cur] + 1] = f[u];
if(u == t) continue;
for(int j = pre[u]; j; j = nxt[j])
if(belong[to[j]] == cur) {
g[i][siz[cur] + 1] += out[u];
g[i][num[to[j]]] -= out[u];
}
}
gauss(siz[cur]);
for(int i = 1;i <= siz[cur]; ++i) {
int u = block[cur][i];
f[u] = g[i][siz[cur] + 1] / g[i][i];
for(int j = repre[u]; j; j = renxt[j])
if(belong[reto[j]] != cur) {
--bout[belong[reto[j]]];
if(!bout[belong[reto[j]]]) q[++tail] = belong[reto[j]];
f[reto[j]] += (f[u] + 1) * out[reto[j]];
}
}
}
printf("%.3lf\n", f[s]);
return 0;
}