【CF446D】【JZOJ 3748】DZY Loves Games

Description

今天DZY 想要玩一个古老的游戏。他在一个有n 个房间并有m 个走廊互相连接的大迷宫里（每个走廊都允许双向通行）。你可以认为所有房间都被走廊直接或间接连接。

DZY 在迷宫里迷路了。现在他在第一房间并且有k 条命。他将会按如下所述行动：

•首先，他会随机抽取一条从他现在所处房间出发的走廊。每个抽取范围内的走廊选中的机率相等。

•然后他会沿着走廊走到走廊的另一端，并且回到第一步重复这个过程。

迷宫中的一些房间里面埋着陷阱。第一房间明显没有陷阱，第n 号房间明确地有一个陷阱。每次DZY 进入这些有陷阱的房间，他都会失去一条命。现在，DZY 知道如果他恰好有两条命时进入了第n 号房间，那么首先他会失去一条命，但是然后他会开启一个福利关卡。他想要知道他开启福利关卡的机率到底为多少。请帮助他。

Solution

只要把从点1到每个点的概率求出，再求出两两黑点之间的概率，最后再用矩阵乘法即可。
于是这题就是高斯消元；
先列出原方程，
设#d_i#表示点i的度数，$b_i$为点i被经过的期望次数，因为黑点到了就停止，所以对于黑点来说，期望经过次数和到达概率是一样的；
首先每个点都可以从每个与它相邻的点走来，概率是走过来的点的$1/d$；
如果从点1出发，那么与它相邻的点都要加上($1/d_1$)这个常数，有因为点1不是黑点，所以可以经过多次，那么原来1的方程就表示走出去以后有走回来的期望次数；
如果从黑点出发，那么于上面一样，与它相邻的点都要加上一个常数，当然可能从这个黑点出去走了一堆白点又走回来；

这样复杂度是$O(n^4)$

观察方程，发现把未知数移到左边后，每次左边都不会变了，变得是右边的常数，
所以我们可以一开始就求出每个b是由哪几个方程右边的常数的几倍加起来，每次变就直接代入即可，

复杂度：$O(n^3+log(m)*n^3)$

Code

#include <iostream>
#include <cstdio>
#include <cstdlib>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fod(i,a,b) for(int i=a;i>=b;i--)
#define efo(i,q) for(int i=A[q];i;i=B[i][0])
#define abs(q) ((q)>0?(q):-(q))
using namespace std;
typedef double db;
const int N=505,M=100500;
int read(int &n)
{
    char ch=' ';int q=0,w=1;
    for(;(ch!='-')&&((ch<'0')||(ch>'9'));ch=getchar());
    if(ch=='-')w=-1,ch=getchar();
    for(;ch>='0' && ch<='9';ch=getchar())q=q*10+ch-48;n=q*w;return n;
}
int m,n,m1;
int B[M*2][2],A[N],B0,B1[N];
int b1[N],zx[N];
db s[N][N],t[N][N],di[N][N];
void link(int q,int w)
{
    B1[q]++,B1[w]++;
    B[++B0][0]=A[q];A[q]=B0,B[B0][1]=w;
    B[++B0][0]=A[w];A[w]=B0,B[B0][1]=q;
}
void chenans()
{
    fo(i,1,m)fo(j,1,m)t[i][j]=0;
    fo(k,1,m)fo(i,1,m)fo(j,1,m)t[i][j]+=s[i][k]*di[k][j];
    fo(i,1,m)fo(j,1,m)s[i][j]=t[i][j];
}
void chen()
{
    fo(i,1,m)fo(j,1,m)t[i][j]=0;
    fo(k,1,m)fo(i,1,m)fo(j,1,m)t[i][j]+=di[i][k]*di[k][j];
    fo(i,1,m)fo(j,1,m)di[i][j]=t[i][j];
}
db JC()
{
    while(m1)
    {
        if(m1&1)chenans();
        chen();m1>>=1;
    }
    db ans=0;
    fo(i,1,m)ans+=s[i][m];
    return ans;
}
db b[N][N],c[N][N];
int z[N];
void GS()
{
    fo(i,1,n)
    {
        b[i][i]=c[i][i]=1;
        efo(j,i)if(!zx[B[j][1]])b[i][B[j][1]]-=1.0/B1[B[j][1]];
    }
    fo(i,1,n)
    {
        int q=i;
        fo(j,i+1,n)if(abs(b[j][i])>abs(b[q][i]))q=j;
        if(q!=i)
        {
            fo(j,1,n)swap(b[i][j],b[q][j]),swap(c[i][j],c[q][j]);
        }
        fo(j,1,n)if(i!=j)
        {
            db t=b[j][i]/b[i][i];
            fo(k,1,n)b[j][k]-=t*b[i][k],c[j][k]-=t*c[i][k];
        }
    }
    fo(i,1,n)fo(j,1,n)c[i][j]/=b[i][i];
    efo(i,1)z[B[i][1]]++;
    fo(i,1,n)if(zx[i])
    {
        fo(j,1,n)
            s[zx[i]][zx[i]]+=c[i][j]*((db)z[j]/B1[1]);
    }
    fo(I,1,n)if(zx[I]&&B1[I])
    { 
        fo(i,1,n)z[i]=0; 
        efo(i,I)z[B[i][1]]++;
        fo(i,1,n)if(zx[i])
        {
            fo(j,1,n)di[zx[I]][zx[i]]+=c[i][j]*((db)z[j]/B1[I]);
        }
    }
}
int main()
{
    freopen("games.in","r",stdin);
    freopen("games.out","w",stdout);
    int q,w;
    read(n),read(m),read(m1);
    m1-=2;if(m1<0){printf("0\n");return 0;}
    fo(i,1,n)if(read(q))b1[++b1[0]]=i,zx[i]=b1[0];
    fo(i,1,m)read(q),read(w),link(q,w); 
    GS();
    m=b1[0];
    printf("%.12lf\n",JC());
    return 0;
}