题目
题目链接：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;
}