bzoj 2707: [SDOI2012]走迷宫（Trajan+高斯消元+Dp）

2707: [SDOI2012]走迷宫

Time Limit: 10 Sec   Memory Limit: 128 MB
Submit: 862   Solved: 328
[ Submit][ Status][ Discuss]

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

很经典的题，方法/步骤如下：

①输入边，建图和反相图，无视掉以终点为起点的边

②如果起点终点不连通或者某个点起点能到，但这个点到不了终点就是INF，对原图和反相图分别进行一次DFS就可以判断出来

③如果是个有向无环图，那么很显然可以求出拓扑序然后DP，复杂度O(n+m)

⑤如果有环，那么只能列方程高斯消元，复杂度O(n^3+m)

高斯消元方法：每个点都可以列出一个方程：F[i] = ∑(F[j]/out[i], 存在边(i, j))+1

其中F[i]为从i点走到终点的期望步数，out[i]为i点的出度

这里有个疑问：为什么不能设F[i]为从起到到i点的期望步数？   ------   其实这也是很多概率DP题为什么要倒过来求的原因：其实这样可以，但你会发现情况变得复杂多了，再计算F[i]时，很显然i点的出边全部毫无意义，而在计算非i点时i点出边有有意义，所以这样每次算一个点的时候都要把那个点出边拆掉才能列出方程，这样很显然整个方程都列不出来，处理也很麻烦

而这道题有环，但是特意提了每个强连通分量点数不超过100，而如果将图缩点的话，就变成了有向无环图而可以使用方法③，

而对于每个强连通分量高斯消元即可，复杂度O(n^4)

#include<stdio.h>
#include<stack>
#include<math.h>
#include<vector>
#include<algorithm>
#include<string.h>
using namespace std;
stack<int> st;
vector<int> G[10005], GF[10005], Bel[10005];
int S, T, out[10005], she[10005], ok[10005], vis[10005], cnt, num, scc[10005], time[10005], dfn[10005];
double ave[10005], Jz[105][105], c[105];
void Jud1(int u)
{
	int i, v;
	ok[u]++, vis[u] = 1;
	for(i=0;i<G[u].size();i++)
	{
		v = G[u][i];
		if(vis[v])
			continue;
		Jud1(v);
	}
}
void Jud2(int u)
{
	int i, v;
	ok[u]++, vis[u] = 1;
	for(i=0;i<GF[u].size();i++)
	{
		v = GF[u][i];
		if(vis[v] || ok[v]==0)
			continue;
		Jud2(v);
	}
}
void Trajan(int u)
{
	int i, v;
	vis[u] = 1;
	st.push(u);
	time[u] = dfn[u] = ++cnt;
	for(i=0;i<G[u].size();i++)
	{
		v = G[u][i];
		if(vis[v]==0)
		{
			Trajan(v);
			dfn[u] = min(dfn[u], dfn[v]);
		}
		else if(scc[v]==0)
			dfn[u] = min(dfn[u], time[v]);
	}
	if(dfn[u]==time[u])
	{
		num++;
		while(st.empty()==0)
		{
			v = st.top();
			st.pop();
			scc[v] = num;
			Bel[num].push_back(v);
			if(v==u)
				break;
		}
	}
}
void Gauss(int n, int m)
{
	int i, j, p, q, r;
	p = q = 1;
	while(p<=n && q<=m)
	{
		r = p;
		for(i=p+1;i<=n;i++)
		{
			if(fabs(Jz[i][q])>fabs(Jz[p][q]))
				r = i;
		}
		for(i=q;i<=m;i++)
			swap(Jz[r][i], Jz[p][i]);
		swap(c[r], c[p]);
		c[p] /= Jz[p][q];
		for(i=q+1;i<=m;i++)
			Jz[p][i] /= Jz[p][q];
		for(i=1;i<=n;i++)
		{
			if(Jz[i][q] && i!=p)
			{
				for(j=q+1;j<=m;j++)
					Jz[i][j] -= Jz[p][j]*Jz[i][q];
				c[i] -= c[p]*Jz[i][q];
				Jz[i][q] = 0;
			}
		}
		q++, p++;
	}
}
int main(void)
{
	int i, x, y, n, m, j, k;
	scanf("%d%d%d%d", &n, &m, &S, &T);
	for(i=1;i<=m;i++)
	{
		scanf("%d%d", &x, &y);
		if(x==T)
			continue;
		G[x].push_back(y);
		GF[y].push_back(x);
		out[x]++;
	}
	Jud1(S);
	memset(vis, 0, sizeof(vis));
	Jud2(T);
	for(i=1;i<=n;i++)
	{
		if(ok[i]==1)
		{
			printf("INF\n");
			return 0;
		}
	}
	if(S==T)
	{
		printf("0.000\n");
		return 0;
	}
	memset(vis, 0, sizeof(vis));
	Trajan(S);
	for(i=scc[T]+1;i<=scc[S];i++)
	{
		memset(Jz, 0, sizeof(Jz));
		memset(c, 0, sizeof(c));
		for(j=0;j<Bel[i].size();j++)
			she[Bel[i][j]] = j+1;
		for(j=0;j<Bel[i].size();j++)
		{
			x = Bel[i][j];
			Jz[j+1][j+1] = -1;
			c[j+1] -= 1+ave[x]/out[x];
			for(k=0;k<GF[x].size();k++)
			{
				y = GF[x][k];
				if(scc[y]==scc[x])
					Jz[she[y]][she[x]] += 1.0/out[y];
			}
		}
		Gauss(j, j);
		for(j=0;j<Bel[i].size();j++)
		{
			x = Bel[i][j];
			for(k=0;k<GF[x].size();k++)
			{
				y = GF[x][k];
				if(scc[y]!=scc[x])
					ave[y] += c[she[x]];
			}
		}
	}
	printf("%.3f\n", fabs(c[she[S]]));
	return 0;
}