洛谷 P6030 [SDOI2012]走迷宫 概率与期望+高斯消元

题目描述

传送门

分析

首先判掉 (INF) 的情况

第一种情况就是不能从 (s) 走到 (t)

第二种情况就是从 (s) 出发走到了出度为 (0) 的点，这样就再也走不到 (t)

然后我们去考虑 (60) 分的做法

我们设 (dp[u]) 为当前在点 (u) 走到点 (t) 的期望步数,(du[u]) 为 (u) 的出度

那么就有

(dp[u]=sum_{u->v}^v((dp[v]+1) imes frac{1}{du[u]}))

移项之后就变成了

(dp[u]-sum_{u->v}^v(dp[v] imes frac{1}{du[u]})=1)

初始化 (dp[t]=0)

一共有 (n) 个未知数，(n)个方程可以用高斯消元解决

这样做时间复杂度是 (n^3) 的

下面我们去考虑满分做法

因为保证强连通分量的大小不超过 (100)

所以我们可以在强联通分量内使用高斯消元

对于整个 (DAG) 使用拓扑排序去解决

注意拓扑排序要建反图

对于拓扑序小于 (t) 所属强联通分量的强联通分量，我们没有必要去算

因为在到这些点之前就已经到达了 (t)

对于其它的强联通分量

如果分量内的某个点的去边所指向的点在当前的强联通分量中，我们把其作为未知数参与消元

对于不在这个强联通分量中的点，这个点肯定在之前被求出来过

我们把它作为常数

代码

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<iostream>
#include<queue>
#include<vector>
inline int read(){
	int x=0,fh=1;
	char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*fh;
}
typedef double db;
const db eps=1e-8;
const int maxn=1e6+5;
const int maxk=105;
const int maxm=1e4+5;
int head[maxn],tot=1;
struct asd{
	int to,next;
}b[maxn],b2[maxn];
void ad(int aa,int bb){
	b[tot].to=bb;
	b[tot].next=head[aa];
	head[aa]=tot++;
}
int h2[maxn],t2=1;
void ad2(int aa,int bb){
	b2[t2].to=bb;
	b2[t2].next=h2[aa];
	h2[aa]=t2++;
}
int rd[maxn];
db mp[maxk][maxk];
int jl[maxk][maxk];
int low[maxn],dfn[maxn],dfnc,js,shuyu[maxn],sta[maxn],top,siz[maxn];
void tar(int now){
	dfn[now]=low[now]=++dfnc;
	sta[++top]=now;
	for(int i=head[now];i!=-1;i=b[i].next){
		int u=b[i].to;
		if(!dfn[u]){
			tar(u);
			low[now]=std::min(low[now],low[u]);
		} else if(!shuyu[u]){
			low[now]=std::min(low[now],dfn[u]);
		}
	}
	if(dfn[now]==low[now]){
		js++;
		while(1){
			int y=sta[top--];
			shuyu[y]=js;
			siz[js]++;
			if(y==now) break;
		}
	}
}
int n,m,s,t,rk[maxm],du[maxm],xh[maxm];
std::vector<int> g[maxn];
db dp[maxn];
bool jud;
void gsxy(int id){
	memset(rk,0,sizeof(rk));
	memset(xh,0,sizeof(xh));
	int ncnt=0;
	for(int i=1;i<=n;i++){
		if(shuyu[i]==id){
			rk[i]=++ncnt;
			xh[ncnt]=i;
			g[id].push_back(i);
			if(i==t) jud=1;
		}
	}
	if(jud==0){
		return;
	} else {
		memset(mp,0,sizeof(mp));
		memset(jl,0,sizeof(jl));
		for(int i=1;i<=ncnt;i++){
			mp[i][ncnt+1]=1;
		}
		for(int i=0;i<g[id].size();i++){
			int now=g[id][i];
			for(int i=head[now];i!=-1;i=b[i].next){
				int u=b[i].to;
				if(shuyu[u]==shuyu[now]){
					jl[rk[now]][rk[u]]++;	
				} else {
					mp[rk[now]][ncnt+1]+=(double)dp[u]/du[now];
				}
			}
		}
		for(int i=0;i<g[id].size();i++){
			for(int j=0;j<g[id].size();j++){
				int aa=g[id][i],bb=g[id][j];
				if(i==j){
					mp[rk[aa]][rk[aa]]=(db)(du[aa]-jl[rk[aa]][rk[aa]])/du[aa];
				} else {
					mp[rk[aa]][rk[bb]]=(db)(-1.0)*jl[rk[aa]][rk[bb]]/du[aa];
				}
			}
		}
		for(int i=0;i<g[id].size();i++){
			if(g[id][i]==t){
				int now=rk[g[id][i]];
				for(int j=1;j<=ncnt+1;j++){
					mp[now][j]=0;
				}
				mp[now][now]=1;
			}
		}
		int now=1;
		for(int i=1;i<=ncnt;i++){
			double mmax=0;
			int jl=0;
			for(int j=now;j<=ncnt;j++){
				if(std::fabs(mp[j][i])>std::fabs(mmax)){
					mmax=mp[j][i];
					jl=j;
				}
			}
			if(std::fabs(mmax)<eps) continue;
			if(jl!=now) std::swap(mp[jl],mp[now]);
			for(int j=i+1;j<=ncnt+1;j++){
				mp[now][j]/=mp[now][i];
			}
			mp[now][i]=1.0;
			for(int j=now+1;j<=ncnt;j++){
				db cs=mp[j][i];
				for(int k=i;k<=ncnt+1;k++){
					mp[j][k]-=cs*mp[now][k];
				}
			}
			now++;
		}
		dp[xh[ncnt]]=mp[ncnt][ncnt+1];
		for(int i=ncnt-1;i>=1;i--){
			dp[xh[i]]=mp[i][ncnt+1];
			for(int j=i+1;j<=n;j++){
				dp[xh[i]]-=mp[i][j]*dp[xh[j]];
			}
		}
	}
}
std::queue<int> q;
int huan[maxn];
void tp(){
	for(int i=1;i<=js;i++){
		if(rd[i]==0) q.push(i);
	}
	while(!q.empty()){
		int now=q.front();
		gsxy(now);
		q.pop();
		for(int i=h2[now];i!=-1;i=b2[i].next){
			int u=b2[i].to;
			rd[u]--;
			if(rd[u]==0) q.push(u);
		}
	}
}
bool vis[maxn];
void dfs(int now){
	vis[now]=1;
	for(int i=head[now];i!=-1;i=b[i].next){
		int u=b[i].to;
		if(!vis[u]) dfs(u);
	}
}
int main(){
	memset(h2,-1,sizeof(h2));
	memset(head,-1,sizeof(head));
	n=read(),m=read(),s=read(),t=read();
	for(int i=1;i<=m;i++){
		int aa,bb;
		aa=read(),bb=read();
		ad(aa,bb);
		du[aa]++;
		if(aa==bb) huan[aa]++;
	}
	dfs(s);
	if(vis[t]==0){
		printf("INF
");
		return 0;
	}
	for(int i=1;i<=n;i++){
		if(i==t) continue;
		if(vis[i] && du[i]-huan[i]==0){
			printf("INF
");
			return 0;
		}
	}
	for(int i=1;i<=n;i++){
		if(!dfn[i]) tar(i);
	}
	for(int i=1;i<=n;i++){
		for(int j=head[i];j!=-1;j=b[j].next){
			int u=b[j].to;
			if(shuyu[i]!=shuyu[u]){
				ad2(shuyu[u],shuyu[i]);
				rd[shuyu[i]]++;
			}
		}
	}
	tp();
	printf("%.3f
",dp[s]);
	return 0;
}

更新

补一下判 (INF) 的锅

否则遇到这组数据会挂

5 6 1 5
1 1
4 4 
2 1
1 5 
3 4 
1 3