【期望DP+高斯消元】BZOJ3270-博物馆

【题目大意】

有m条走廊连接的n间房间，并且满足可以从任何一间房间到任何一间别的房间。两个男孩现在分别处在a，b两个房间，每一分钟有Pi 的概率在这分钟内不去其他地方（即呆在房间不动），有1-Pi 的概率他会在相邻的房间中等可能的选择一间并沿着走廊过去。求两人在每间房间相遇的概率。

【思路】

向用链表把图建出来，建立的时候不要忘记了自己向自己再加一条边(见下面情况1)

我们将男孩A位于x，男孩B位于y的状态定义为id[x][y]。由id[tx][ty]状态转移到id[i][j]状态的概率。如果当前已经知道了抵达id[i][j]状态的期望概率为F[i,j]。

考虑以下几种可能性：

(1)tx=i且ty=j，说明两者都停留在原地不动，概率为Ptx*Pty*F[i,j];

(2)tx=i且ty≠j，说明男孩B移动了而男孩A没有，概率为Ptx*((1-Pty)/Dty)*F[i,ty]，其中Dty为ty的出度;

(3)tx≠i且ty=j，概率为Pty*((1-Ptx)/Dtx)*F[tx,j];

(4)tx≠i且ty≠j，概率为((1-Ptx)/Dtx)*F[i,ty]*((1-Pty)/Dty)*F[tx,j]。

我们把F看成是未知数，它们前面乘上的看作是系数，对于F[i,j]可以列出以下方程：

F[i,j]=情况(1)+情况(2)+情况(3)+情况(4)。

0=(Ptx*Pty-1)*F[i,j]+情况(2)+情况(3)+情况(4)……(*)

如上 (*)式子总共可以列出n*n个。有一个特例是在其实位置，也就是F(a,b)，由于一开始就位于id[a][b]，也就是说初始时概率就为1了，所以在(*)右边还要加上1，移到等式左边就变成了-1，所以有a[id[aa][bb]][n*n+1]=-1。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
#include<cmath>
using namespace std;
const int MAXN=20+5;
double f[MAXN*MAXN][MAXN*MAXN],p[MAXN];
int id[MAXN][MAXN];
vector<int> E[MAXN];
double a[MAXN*MAXN][MAXN*MAXN]; 
int n,m,aa,bb;

void init()
{
    scanf("%d%d%d%d",&n,&m,&aa,&bb);
    int cnt=0;
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++) id[i][j]=++cnt;
    a[id[aa][bb]][n*n+1]=-1;
    for (int i=1;i<=n;i++) E[i].push_back(i);//不要忘记把自己向自己是可以走的
    for (int i=1;i<=m;i++)
    {
        int u,v;
        scanf("%d%d",&u,&v);
        E[u].push_back(v);
        E[v].push_back(u); 
    }
    for (int i=1;i<=n;i++) scanf("%lf",&p[i]);
}

void buildequ()
{
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
        {
            a[id[i][j]][id[i][j]]--;//位于方程右边系数为1，移到左边系数变成-1 
            for (int x=0;x<E[i].size();x++)
                for (int y=0;y<E[j].size();y++)
                {
                    int tx=E[i][x],ty=E[j][y];
                    int fr=id[i][j],to=id[tx][ty];
                    if (tx!=ty)
                    {
                        if (tx==i && ty==j) a[fr][to]+=p[tx]*p[ty];//停留在原地不动
                        else if (tx==i) a[fr][to]+=p[tx]*(1-p[ty])/(E[ty].size()-1);
                        else if (ty==j) a[fr][to]+=p[ty]*(1-p[tx])/(E[tx].size()-1);
                        else a[fr][to]+=(1-p[tx])*(1-p[ty])/(E[ty].size()-1)/(E[tx].size()-1);
                    } 
                } 
        }
}

void guass()
{
    for (int i=1;i<=n*n;i++)
    {
        int t=i;
        for (int j=i+1;j<=n*n;j++) if (fabs(a[j][i])>fabs(a[t][i])) t=j;
        if (t!=i) for (int j=i;j<=n*n+1;j++) swap(a[i][j],a[t][j]);//不要忘记要到n*n+1 
        for (int j=i+1;j<=n*n;j++)
        {
            double x=a[j][i]/a[i][i];
            for (int k=i;k<=n*n+1;k++) a[j][k]-=a[i][k]*x;
        }
    }
    for (int i=n*n;i>=1;i--)
    {
        for (int j=i+1;j<=n*n;j++) a[i][n*n+1]-=a[i][j]*a[j][n*n+1];
        a[i][n*n+1]/=a[i][i]; 
    }
}

void printans()
{
    for (int i=1;i<=n;i++) printf("%.6lf ",a[id[i][i]][n*n+1]);
}

int main()
{
    init();
    buildequ();
    guass();
    printans();
    return 0;
}