题目

题目链接：https://codeforces.com/contest/446/problem/D
有一张 \(n\) 个点 \(m\) 条边的无向图，有一些点有陷阱（保证点 \(1\) 没有，点 \(n\) 有）。最开始有 \(k\) 条命，每次走到一个有陷阱的点就会减少一条命。每次会从当前点等概率随机选择一条出边走。求从点 \(1\) 开始，某一时刻到达点 \(n\) 后恰好只有 \(1\) 条命的期望。注意点 \(n\) 可以多次到达。
\(n\leq 500\)，\(m\leq 10^5\)，\(k\leq 10^9\)，有陷阱的点数 \(\leq 101\)。

思路

如果所有点都是陷阱点，那么矩阵乘法搞一下就可以了。这启发我们求出从陷阱点 \(i\) 走到陷阱点 \(j\)，中间不经过其他陷阱点的期望。
考虑一个有陷阱的点 \(x\)，记 \(f_i\) 表示点 \(i\) 到达点 \(x\)，中间没有经过除了 \(x\) 以外的陷阱点的期望。
那么如果 \(i\) 不是陷阱点，有

\[f_i=\sum_{(i,j)\in E}\frac{f_j}{\text{deg}[i]} \]

如果 \(i\) 是 \(x\) 以外的陷阱点，则 \(f_i=0\)。
如果 \(i=x\)，则 \(f_i=1\)。
然后跑高斯消元即可。
但是很显然不可能对每一个陷阱点都跑一次高斯消元。注意到对于不同的陷阱点，只有陷阱点的常数项是有变化的，矩阵内其他元素都是不变的。
设有 \(cnt\) 个陷阱点，可以把原来的常数项改为一个 \(cnt\) 维的向量，其中第 \(i\) 维就表示第 \(i\) 个陷阱点的期望。然后只需要跑一次高斯消元就可以了。
最后跑矩阵快速幂即可。
时间复杂度 \(O(n^3+cnt^3\log k+m)\)。

代码

#include <bits/stdc++.h>
using namespace std;

const int N=510;
const double eps=1e-10;
int n,m,t,cnt,typ[N],deg[N],id[N],rk[N],G[N][N];
double a[N][N],b[N][N];

struct Matrix
{
	double a[N][N];
	
	friend Matrix operator *(Matrix a,Matrix b)
	{
		Matrix c;
		for (int i=1;i<=cnt;i++)
			for (int j=1;j<=cnt;j++) c.a[i][j]=0;
		for (int i=1;i<=cnt;i++)
			for (int j=1;j<=cnt;j++)
				for (int k=1;k<=cnt;k++)
					c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j]);
		return c;
	}
}f,c;

Matrix fpow(Matrix x,int k)
{
	Matrix ans;
	for (int i=1;i<=cnt;i++)
		for (int j=1;j<=cnt;j++) ans.a[i][j]=(i==j);
	for (;k;k>>=1,x=x*x)
		if (k&1) ans=ans*x;
	return ans;
}

void gauss()
{
	for (int i=1;i<=n;i++)
	{
		for (int j=i;j<=n;j++)
			if (fabs(a[j][i])>eps)
			{
				for (int k=1;k<=n;k++) swap(a[i][k],a[j][k]);
				for (int k=1;k<=cnt;k++) swap(b[i][k],b[j][k]);
				break;
			}
		for (int j=i+1;j<=n;j++)
		{
			double base=a[j][i]/a[i][i];
			for (int k=1;k<=n;k++) a[j][k]-=a[i][k]*base;
			for (int k=1;k<=cnt;k++) b[j][k]-=b[i][k]*base;
		}
	}
	for (int i=n;i>=1;i--)
		for (int j=1;j<=cnt;j++)
		{
			for (int k=i+1;k<=n;k++)
				b[i][j]-=a[i][k]*b[k][j];
			b[i][j]/=a[i][i];
		}
}

int main()
{
	scanf("%d%d%d",&n,&m,&t);
	for (int i=1;i<=n;i++)
	{
		scanf("%d",&typ[i]);
		if (typ[i]) id[i]=++cnt,rk[cnt]=i;
	}
	for (int i=1,x,y;i<=m;i++)
	{
		scanf("%d%d",&x,&y);
		G[x][y]++; G[y][x]++; deg[x]++; deg[y]++;
	}
	for (int i=1;i<=n;i++)
		if (!typ[i])
		{
			a[i][i]=1;
			for (int j=1;j<=n;j++)
				if (i!=j) a[i][j]=-1.0*G[i][j]/deg[i];
		}
		else a[i][i]=b[i][id[i]]=1;
	gauss();
	for (int i=1;i<=cnt;i++)
	{
		for (int j=1;j<=cnt;j++)
			for (int k=1;k<=n;k++)
				f.a[i][j]+=b[k][j]*G[rk[i]][k]/deg[rk[i]];
		c.a[1][i]=b[1][i];
	}
	c=c*fpow(f,t-2);
	printf("%.10f\n",c.a[1][cnt]);
	return 0;
}