[SDOI2012]走迷宫 （强连通分量缩点，动态规划，高斯消元）
题面
Morenan被困在了一个迷宫里。迷宫可以视为N个点M条边的有向图，其中Morenan处于起点S，迷宫的终点设为T。可惜的是，Morenan非常的脑小，他只会从一个点出发随机沿着一条从该点出发的有向边，到达另一个点。这样，Morenan走的步数可能很长，也可能是无限，更可能到不了终点。若到不了终点，则步数视为无穷大。但你必须想方设法求出Morenan所走步数的期望值。
N<=10000，M<=1000000，保证强连通分量的大小不超过 200。
题解
做法没什么好细说的，
就是缩点，然后按照拓扑序倒着求每个强连通分量的 DP 值（每个点出发的期望步数），强连通分量内部用高斯消元。如果存在概率走入一个死环（无法走到终点的强连通分量），那么就无穷大。
剩下的就是一些细节：
没必要真的求拓扑序，只需要按照 dfs 序倒着来，也就是 Tarjan 缩点产生强连通分量的正顺序来就行了。
解方程如果无解，说明期望步数无穷大。
解方程之前注意清空（包括 size+1 列），注意各种标号。
CODE
#include<set>
#include<map>
#include<queue>
#include<stack>
#include<vector>
#include<bitset>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
#define MAXN 1000005
#define LL long long
#define DB double
#define lowbit(x) (-(x) & (x))
#define ENDL putchar('\n')
#define FI first
#define SE second
#define eps 1e-9
LL read() {
LL f=1,x=0;int s = getchar();
while(s < '0' || s > '9') {if(s<0) return -1;if(s=='-')f = -f;s = getchar();}
while(s >= '0' && s <= '9') {x=x*10+(s^48);s = getchar();}
return f*x;
}
void putpos(LL x) {if(!x)return ;putpos(x/10);putchar('0'+(x%10));}
void putnum(LL x) {
if(!x) {putchar('0');return ;}
if(x<0) {putchar('-');x = -x;}
return putpos(x);
}
void AIput(LL x,int c) {putnum(x);putchar(c);}
const int MOD = 1000000007;
int n,m,s,o,k;
int od[MAXN],S,T;
int hd[MAXN],v[MAXN<<1],nx[MAXN<<1],cne;
void ins(int x,int y) {
nx[++ cne] = hd[x]; v[cne] = y; hd[x] = cne;
}
int dfn[MAXN],low[MAXN],tim;
bool f[MAXN];
int fail[MAXN];
stack<int> st;
int bel[MAXN],cnt;
vector<int> bu[MAXN];
int rk[MAXN];
void dfs0(int x) {
low[x] = dfn[x] = ++ tim;
st.push(x); f[x] = 1;
for(int i = hd[x];i;i = nx[i]) {
int y = v[i];
if(!dfn[y]) {
dfs0(y);
low[x] = min(low[x],low[y]);
}
else if(f[y]) {
low[x] = min(low[x],dfn[y]);
}
}
if(low[x] == dfn[x]) {
fail[++ cnt] = -1;
while(f[x]) {
int t = st.top(); st.pop();
f[t] = 0;
bel[t] = cnt;
bu[cnt].push_back(t);
rk[t] = bu[cnt].size();
}
}
return ;
}
DB a[205][205],dp[MAXN];
DB Abs(DB x) {return x < 0 ? -x:x;}
int Gauss(int n) {
for(int i = 1;i <= n;i ++) {
for(int j = i+1;j <= n;j ++) {
if(Abs(a[j][i]) > Abs(a[i][i])) {
swap(a[i],a[j]);
}
}
if(Abs(a[i][i]) < eps) continue;
for(int j = i+1;j <= n;j ++) {
DB nm = a[j][i] / a[i][i];
for(int k = n+1;k >= i;k --) {
a[j][k] -= a[i][k] * nm;
}
}
}
for(int i = n;i > 0;i --) {
for(int j = n;j > i;j --) {
a[i][n+1] -= a[i][j] * a[j][n+1];
a[i][j] = 0;
}
if(!a[i][i]) return 1;
a[i][n+1] /= a[i][i];
a[i][i] = 1;
}
return 0;
}
int main() {
n = read(); m = read();
S = read(); T = read();
for(int i = 1;i <= m;i ++) {
s = read();o = read();
ins(s,o); od[s] ++;
}
for(int i = 1;i <= n;i ++) {
if(!dfn[i]) {
dfs0(i);
}
}
for(int i = 1;i <= cnt;i ++) {
int sz = bu[i].size();
for(int j = 1;j <= sz;j ++) {
for(int k = 1;k <= sz+1;k ++) a[j][k] = 0;
}
for(int j = 1;j <= sz;j ++) {
a[j][j] = 1.0;
}
for(int j = 0;j < sz;j ++) {
int x = bu[i][j];
if(x == T) continue;
if(!od[x]) {
fail[i] = 1; break;
}
a[rk[x]][sz+1] += 1.0;
for(int k = hd[x];k;k = nx[k]) {
int y = v[k];
if(y == T) a[rk[x]][sz+1] += 1.0/od[x] * dp[T];
else if(bel[y] != i) {
if(fail[bel[y]] > 0) {fail[i] = 1;break;}
a[rk[x]][sz+1] += 1.0/od[x] * dp[y];
}
else {
a[rk[x]][rk[y]] -= 1.0/od[x];
}
}
if(fail[i] > 0) break;
}
if(fail[i] > 0) continue;
fail[i] = Gauss(sz);
for(int j = 0;j < sz;j ++) {
int x = bu[i][j];
dp[x] = a[rk[x]][sz+1];
}
}
if(fail[bel[S]] > 0) printf("INF\n");
else printf("%.3f\n",dp[S]);
return 0;
}
``
