【Codeforces 446D】DZY Loves Games

Description

一个N个点M条边的无向连通图，每个点是黑点或白点
从1出发，随机抽取一条从当前点连出的边，每条边选中的机率相等，走到该边另一端的点
求从1点到达n点恰好经过k-2个黑点的概率
保证1点白色，n点黑色
有重边无自环
2 <=n <=500; 1 <= m <=10^5; 2 <= k <= 10^9
黑点不超过100个

Analysis

首先你看到k那么大，就知道肯定是矩阵乘法什么的
我们怎样快速算k呢，如果我们能预处理出黑点间两两互相到达，途中不经过其他黑点的概率，然后显然矩乘，这道题就解决了
于是关键变成了预处理那个概率
枚举每个黑点S出发，我们发现对于每个点可以列出一条方程
对于u连出的黑点v，如果不是S，因为走到S就停止了，S不可能走到u
然后对于其他v点，对u点的贡献为$p_u=\sum_{(v,u)\in E} p_v/d[v]$，$d[v]$表示v点的度数
但是注意，这样子算出的$p_v$有可能是大于1的（什么？！概率大于1？！）
其实呢，这个 p <script type="math/tex" id="MathJax-Element-364">p</script>的真正意义是从S出发的经过一个点的期望次数
然后我们知道经过该点的概率=期望经过次数/实际经过次数
对于黑点，它的实际经过次数为1，所以答案不会错QAQ
有了方程组就可以高斯消元了
那么算上枚举的时间，是O(n^4)，会超时
然后发现，枚举每个点，根本没必要重新n^3算（想想如何优化）
提示：只有常数项会改变

Code

#include<cstdio>
#include<cstring>
#include<algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
using namespace std;
typedef double db;
const int N=505,M=105;
int n,m,k,d[N],b[N],bz[N],map[N][N];
db e[N][N],c[N][N],f[N][N],a[M][M],ans[M][M],t[M][M];
void calc(int v)
{
    fo(i,1,n) e[i][0]=map[v][i]*1.0/d[v];
    fo(k,1,n-1)
        if(k!=v)
            fo(j,k+1,n) e[j][0]+=e[k][0]*c[j][k];
    fd(i,n,1)
    {
        e[i][0]/=e[i][i],f[v][i]=e[i][0];
        if(i!=v)
            fo(j,1,i-1) e[j][0]-=e[j][i]*e[i][0];
    }
}
void Gauss()
{
    fo(i,1,n)
        fo(j,i+1,n)
        {
            c[j][i]=-e[j][i]/e[i][i];
            fo(k,i+1,n) e[j][k]+=c[j][i]*e[i][k];
            e[j][i]=0;
        }
    fo(l,1,b[0]) calc(b[l]);
    calc(1);
}
void mul(db C[M][M],db A[M][M],db B[M][M])
{
    memset(t,0,sizeof(t));
    fo(i,1,b[0])
        fo(j,1,b[0])
            fo(k,1,b[0]) t[i][k]+=A[i][j]*B[j][k];
    memcpy(C,t,sizeof(t));
}
void qmi(int k)
{
    for(;k;k>>=1)
    {
        if(k&1) mul(ans,ans,a);
        mul(a,a,a);
    }
}
int main()
{
    freopen("games.in","r",stdin);
    freopen("games.out","w",stdout);
    int u,v;
    scanf("%d %d %d",&n,&m,&k);
    fo(i,1,n)
    {
        scanf("%d",&bz[i]);
        if(bz[i]) b[++b[0]]=i;
    }
    fo(i,1,m)
    {
        scanf("%d %d",&u,&v);
        map[u][v]++,map[v][u]++,d[u]++,d[v]++;
    }
    fo(i,1,n)
    {
        e[i][0]=0,e[i][i]=1;
        fo(j,1,n)
            if(j==1 || !bz[j]) e[i][j]-=map[j][i]*1.0/d[j];
    }
    Gauss();
    fo(i,1,b[0])
        fo(j,1,b[0]) a[i][j]=f[b[i]][b[j]];
    fo(i,1,b[0]) ans[i][i]=1;
    qmi(k-2);
    db t=0;
    fo(i,1,b[0])
        t+=f[1][b[i]]*ans[i][b[0]];
    printf("%.8lf",t);
    return 0;
}