CF446D. DZY Loves Games

题目大意：

有一个连通的迷宫，你可以理解为一个 $n$ 个点 $m$ 条边的无向连通图。有 些点可能有陷阱，保证 $1$ 号点没有陷阱而 $n$ 号点存在陷阱。一开始你有 $k$  条命，每进入一个有陷阱的点你会损失一条命。

DZY听说这个游戏有一个隐藏关卡：如果进入 $n$ 号点时你恰好剩余 $2$ 条命， 你会先被扣除一条命，然后进入一个隐藏的极限关卡（在触发关卡之前 也可以一次或多次进入 $n$ 号点）。DZY很想进入这个关卡，但是他游戏水平不行，所以他的游戏策略是从 $1$ 号点开始，每次随机走到一个与当前点相邻的点。他想知道他能够触发隐藏关卡的概率，保留 $4$ 位小数。 

算法标签：概率与期望dp，矩阵乘法

思路：

令 $a_u$ 为每个时刻在 $u$ 的概率之和，那么有
$$
a_u=\sum_{v\notin trap,(v,u)\in E}\frac{a_{v}}{deg(v)}+[u为起点的概率]
$$
因为我们这里设定令一个非陷阱的点作为起点，才不会算重，所以要限制只有非陷阱的才能转移到当前点。
$$
a_u-\sum_{v\notin traps,(v,u)\in E}\frac{a_v}{deg(v)}=[u为起点]
$$

根据这个式子就可以进行高斯消元了。

首先我们要计算从 $1$ 出发的方案走到各个陷阱的概率，那么除了 $1$ 为起点的概率为 $1$ 。

从一个陷阱走到另一个陷阱时，我们为了让他变成从一个非陷阱走出来，我们令它的所有出度作为起点的概率为 $\frac{1}{dev(u)}$ 。

发现这两类我们需要的式子除了常数项每一项的系数都相同，所以我们可以考虑合在一起计算。

计算出从一个陷阱走到另一个陷阱的概率之后，我们可以用矩阵乘法优化转移过程。

 以下代码：

#include<bits/stdc++.h>
#define il inline
#define db double
#define _(d) while(d(isdigit(ch=getchar())))
using namespace std;
const int N=605,M=1e5+5;
db f[N][N];
int u[M],v[M];
int n,m,k,d[N],id[N],num[N],tot;
struct node{
    db a[105][105];
}s,t;
il int read(){
   int x,f=1;char ch;
   _(!)ch=='-'?f=-1:f;x=ch^48;
   _()x=(x<<1)+(x<<3)+(ch^48);
   return f*x;
}
il void gauss(){
    for(int i=1;i<=n;i++){
        if(!f[i][i]){
            int mx=i;
            for(int j=i+1;j<=n;j++)if(fabs(f[j][i])>fabs(f[mx][i]))mx=j;
            if(mx^i)for(int j=i;j<=n+tot+1;j++)swap(f[mx][j],f[i][j]);
        }
        if(f[i][i]==0)continue;db tmp=f[i][i];
        for(int j=i;j<=n+tot+1;j++)f[i][j]/=tmp;
        for(int j=1;j<=n;j++)if(i^j&&f[j][i]!=0){
            tmp=f[j][i];
            for(int k=i;k<=n+tot+1;k++)f[j][k]-=tmp*f[i][k];
        }
    }
}
il node C(node x,node y){
    node z;
    for(int i=1;i<=tot;i++)for(int j=1;j<=tot;j++)z.a[i][j]=0;
    for(int i=1;i<=tot;i++)for(int j=1;j<=tot;j++)for(int k=1;k<=tot;k++)
        z.a[i][j]+=x.a[i][k]*y.a[k][j];
    return z;
}
il void ksm(int y){
    while(y){
        if(y&1)s=C(s,t);
        t=C(t,t);y>>=1;
    }
}
int main()
{
    n=read();m=read();k=read();
    for(int i=1;i<=n;i++){
        id[i]=read();
        if(id[i])id[i]=++tot,num[tot]=i;
    }
    for(int i=1;i<=m;i++){
        u[i]=read();v[i]=read();
        d[u[i]]++;d[v[i]]++;
    }
    for(int i=1;i<=m;i++){
        int x=u[i],y=v[i];
        if(!id[x])f[y][x]-=1.0/d[x];
        else f[y][id[x]+n]+=1.0/d[x];
        if(!id[y])f[x][y]-=1.0/d[y];
        else f[x][id[y]+n]+=1.0/d[y];
    }
    f[1][1+tot+n]=1;for(int i=1;i<=n;i++)f[i][i]+=1;
    gauss();
    for(int i=1;i<=tot;i++)for(int j=1;j<=tot;j++)t.a[i][j]=f[num[j]][i+n];
    for(int i=1;i<=tot;i++)s.a[1][i]=f[num[i]][n+tot+1];
    ksm(k-2);printf("%lf\n",s.a[1][tot]);
    return 0;
}