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

bzoj2707: [SDOI2012]走迷宫

Description

Morenan被困在了一个迷宫里。迷宫可以视为N个点M条边的有向图，其中Morenan处于起点S，迷宫的终点设为T。可惜的是，Morenan非常的脑小，他只会从一个点出发随机沿着一条从该点出发的有向边，到达另一个点。这样，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] = $\sum {f[v] + w[v][u]}$
如果一个环的大小为cnt，那么我们会有cnt个这样的方程。所以我们可以利用高斯消元解决这类问题。

例题分析

这题是一道比较折腾的题。
首先我们要判断会不会走到“死胡同”
先dfs搜索，遍历点和连通块。
如果起点走不到终点，那就肯定不合法。
除了终点所在的连通块之外，只要这个连通块被搜到过且没有出边，那就是死胡同，也不合法。
还有有一个雕虫小技，就是按照拓扑逆序进行DP，因为对于一个拓扑逆序的东西，它的所有出边一定都已经计算完了，所以是正确的。
那么如果是这样，方程也得逆着退咯
f[u]=$\sum (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[N], reto[M], renxt[M], repre[M], top, dtop;
double g[110][110], out[N], f[N];
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[N], low[N], st[N], belong[N], block[N][110], siz[N], num[N], bout[N], cnt;
int U[M], V[M], q[N], n, m, s, t;
bool vis[N], bvis[N];

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;
}