BZOJ2707 [SDOI2012]走迷宫 【期望DP+高斯消元+tarjan缩点】

题目链接

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%的数据，图中没有环，也没有自环

解法：

INF的情况：
如果一个点无法到达终点，说明这个点的期望步数为 $+\infty$ 我们把它删去
如果一个点的后继无法到达终点，说明这个点的期望步数也为 $+\infty$ ，我们也把它删去
如果最后起点也被删去说明起点期望步数为INF

其他情况：
定义$E\left(u\right)$ 为从点 $u$ 走到终点的期望步数
很明显转移方程为：
$E\left(u\right)=\frac{1}{\mathrm{deg}\left(u\right)}\left({\sum }_{\left(u,v\right)\in E}E\left(v\right)\right)+1$

若图为DAG，则只需要求出逆拓扑序后DP转移即可，但本题可能存在环，所以会有环状转移，于是需要用高斯消元。
但高斯消元整体复杂度$O\left(n^3\right)$ 显然太慢。
题目中提示强连通分量的大小不超过100，所以我们可以先tarjan缩点，再根据缩点后的逆拓扑序依次对每个SCC内用高斯消元来求解。
对于当前SCC内的点，如果它的出点在不在当前SCC内，根据拓扑序，这个点的期望肯定已经在之前被计算出来了，所以放在等号右边当作常量即可

坑点： 可能有重边和自环，所以建立矩阵的时候应该采用 $-$$=$ / $+$$=$
其他细节见代码：

Code：

#include <bits/stdc++.h>
#define pb push_back
#define clr(a,b) memset(a,b,sizeof(a));
using namespace std;
typedef long long ll;
const int INF = 0x3f3f3f3f;
const int MX = 1e6 + 7;
const int K = 103;
const int mod = 1e9 + 7;
double a[K][K], tmp[K],f[MX];
int siz;
void Gauss(){
	for(int i = 0;i < siz;++i){
		int r = i;
		for(int j = i + 1;j < siz;++j)
			if(fabs(a[r][i]) < fabs(a[j][i])) r = j;
		if(fabs(a[r][i]) < eps) return ;
		if(i != r) swap(a[i],a[r]);
		double div = a[i][i];
		for(int j = i;j <= siz;++j) a[i][j] /= div;
		for(int j = i + 1;j < siz;++j){
			div = a[j][i];
			for(int k = i;k <= siz;++k){
				a[j][k] -= a[i][k] * div;
			}
		}
	}

	tmp[siz-1] = a[siz-1][siz];
	for(int i = siz - 2;i >= 0;--i){
		tmp[i] = a[i][siz];
		for(int j = i + 1;j < siz;++j) tmp[i] -= a[i][j] * tmp[j];
	}
	for(int i = 0;i < siz;++i){
		for(int j = 0;j <= siz;++j){
			a[i][j] = 0;
		}
	}
}

int n,m,S,T;
int DFN[MX],LOW[MX],idx = 0,stk[MX],tp = 0,bel[MX],rnk[MX],cnt = 0;
vector<int>scc[MX],inv[MX];
int ecnt = 0,head[MX], indeg[MX], outdeg[MX];
struct Edge{
	int v,next;
}e[MX<<1];

void add(int u,int v){
	outdeg[u]++;//公式中的出度
	e[++ecnt].v = v;
	e[ecnt].next = head[u];
	head[u] = ecnt;
	inv[v].pb(u);
}

bool vis[MX],del[MX],in[MX];
queue<int>q;

void init(){
	vis[T] = true;
	q.push(T);
	//到不了终点的点期望为∞
	while(!q.empty()){
		int u = q.front();q.pop();
		siz = inv[u].size();
		for(int i = 0;i < siz;++i){
			int v = inv[u][i];
			if(!vis[v]) {vis[v] = true;q.push(v);}
		}
	}

	for(int i = 1;i <= n;++i){
		if(!vis[i]) {del[i] = true; q.push(i);}
	}
	//后继到不了终点的点期望也为∞
	while(!q.empty()){
		int u = q.front();q.pop();
		siz = inv[u].size();
		for(int i = 0;i < siz;++i){
			int v = inv[u][i];
			if(!del[v]){
				del[v] = true;q.push(v);
			}
		}
	}
}

void tarjan(int u){
	in[u] = true;stk[++tp] = u;
	DFN[u] = LOW[u] = ++idx;
	for(int i = head[u];i;i = e[i].next){
		int v = e[i].v;
		if(!DFN[v]){
			tarjan(v);LOW[u] = min(LOW[u],LOW[v]);
		} else if(in[v]) LOW[u] = min(LOW[u],DFN[v]);
	}

	if(LOW[u] == DFN[u]){
		cnt++;int x = 0, t = 0;
		while(x != u){
			x = stk[tp--];
			bel[x] = cnt;
			scc[cnt].pb(x);
			rnk[x] = t++;//记录SCC内每个点的序号
			in[x] = false;
		}
	}
}

int p[MX],tot = 0;
vector<int>G[MX];

void topsort(){
	for(int i = 1;i <= cnt;++i) vis[i] = false;//避免SCC之间的重边
	for(int i = 1;i <= cnt;++i){
		siz = scc[i].size();
		for(int j = 0;j < siz;++j){
			int u = scc[i][j];
			for(int k = head[u];k;k = e[k].next){
				int v = e[k].v;
				if(bel[u] != bel[v] && !vis[bel[v]]){
					G[i].pb(bel[v]), vis[bel[v]] = true;
					indeg[bel[v]]++;
				}
			}
		}

		for(int j = 0;j < siz;++j){
			int u = scc[i][j];
			for(int k = head[u];k;k = e[k].next){
				int v = e[k].v;
				if(bel[u] != bel[v]) vis[bel[v]] = false;
			}
		}
	}
	while(!q.empty()) q.pop();
	for(int i = 1;i <= cnt;++i) if(!indeg[i]) q.push(i);
	while(!q.empty()){
		int u = q.front();
		q.pop();
		p[++tot] = u;
		siz = G[u].size();
		for(int i = 0;i < siz;++i){
			int v = G[u][i];
			indeg[v]--;
			if(!indeg[v]) q.push(v);
		}
	}
}

void solve(){
	//拓扑序逆序转移
	while(cnt){
		int x = p[cnt];
		siz = scc[x].size();
		for(int y = 0;y < siz;++y){
			int u = scc[x][y];
			//若当前点已经为终点,则它的出边已经没有意义
			if(u != T){
				for(int i = head[u];i;i = e[i].next){
					int v = e[i].v;
					if(del[v]) continue;
					if(bel[v] == bel[u]) a[y][rnk[v]] -= 1.0 / outdeg[u];//避免重边
					else a[y][siz] += f[v] / outdeg[u];// 拓扑序说明已经计算出来了f[v] 应视为常量
				}
				a[y][siz] += 1.0;
			}
			a[y][y] += 1.0;// 避免自环
		}
		Gauss();
		for(int y = 0;y < siz;++y){
			f[scc[x][y]] = tmp[y];
		}
		cnt--;
	}
}
int main(){
	scanf("%d %d %d %d",&n,&m,&S,&T);
	for(int i = 1;i <= m;++i){
		int u,v;scanf("%d %d",&u,&v);
		add(u,v);
	}
	init();
	if(del[S]) {printf("INF\n");return 0;}
	for(int i = 1;i <= n;++i){
		if(!del[i] && !DFN[i]) tarjan(i);
	}
	topsort();solve();
	printf("%.3f\n", f[S]);
	return 0;
}