bzoj3270 博物馆（概率+gauss）

题目链接

分析：
感觉这道题和聪聪可可有点像

实际上我们还是设置状态：$f(x,y)$
表示Petya在$x$，Vasya在$y$的概率（准确人名良心）
设每个结点的度为$deg[i]$
从i到达与ta相邻的任意节点的概率$out[i]$就是$1-p[i] \over deg[i]$

我们枚举$x,y$的前驱结点$x->u,y->v$，则：
$f(x,y)=\sum f(u,v)*out[x]*out[y]+f(u,y)*out[x]*P_y+f(x,v)*out[y]*P_x+f(x,y)*P_x*P_y$

化一下式子（移项）：
$0=f(u,v)*out[x]*out[y]+f(u,y)*out[x]*P_y+f(x,v)*out[y]*P_x+f(x,y)*P_x*P_y-f(x,y)$

其中$f(a,b)$是未知数，上式就可以看成是一个方程
（一开始我以为是单纯型，但是线性规划需要的是不等式）

我们可以得到$n^2$个方程
直接高斯消元即可

需要特别注意的是$f(s,t)$刚开始的概率是1，后来又会有再次经过的概率
但是我们列式子的时候要把左边的常数项赋值成-1

我们在枚举前驱结点的时候，方便起见，每个点添加一条能到达ta自己的边
这样就可以统一处理前进和停留的情况了
（然而这条边不能算在结点的度之中）

tip

列方程的时候还是要想清楚了，系数不要加重了

总结一下这道题的易错点：

• 加自环，便于统一处理前进和停留的情况
• $f(s,t)$概率是1，$a[f(s,t)][f(s,t)]=-1$
• 方程的系数是增加+，不是赋值；注意判断前驱状态（到达不同状态概率不一样）
• 因为gauss和get中都要用到n，所以不要随便改变n的价值
• 最后答案是fabs
• 方程中枚举的是前驱结点（从哪里来）

感觉自己的gauss写的一直很随意，所以从现在开始，把板子固定下来：

void gauss(int n) {
    int now=1,to;
    double t;
    for (int i=1;i<=n;i++) {        //未知量 
        for (to=now;to<=n;to++)
            if (fabs(a[to][i])>eps) break;
        if (to>n) continue;
        if (to!=now) 
            for (int j=1;j<=n+1;j++)
                swap(a[to][j],a[now][j]);
        for (int j=1;j<=n;j++) 
            if (j!=now) {
                t=a[j][i]/a[now][i];
                for (int k=1;k<=n+1;k++)
                    a[j][k]-=t*a[now][k];
            }
        now++;
    }
    for (int i=1;i<=n;i++) a[i][n+1]/=a[i][i];
}

我一开始把系数化为一这一步放在了里面，一直跑不出答案
最后只能把这一步单独拉出来才行
猛然想起蓝书上说过这样的gauss精度损失较大
于是就想开LD，但是不知道为什么精度反而更低了。。。

#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>

using namespace std;

const double eps=1e-10;
const int N=21;
int n,m,s,t,G[N][N],tot,du[N];
double p[N],out[N],a[402][402];

int get(int x,int y) {
    return (x-1)*n+y;
}

void gauss(int n) {
    int now=1,to;
    double t;
    for (int i=1;i<=n;i++) {        //未知量 
        for (to=now;to<=n;to++)
            if (fabs(a[to][i])>eps) break;
        if (to>n) continue;
        if (to!=now) 
            for (int j=1;j<=n+1;j++)
                swap(a[to][j],a[now][j]);
        for (int j=1;j<=n;j++) 
            if (j!=now) {
                t=a[j][i]/a[now][i];
                for (int k=1;k<=n+1;k++)
                    a[j][k]-=t*a[now][k];
            }
        now++;
    }
    for (int i=1;i<=n;i++) a[i][n+1]/=a[i][i];
}

int main()
{
    scanf("%d%d%d%d",&n,&m,&s,&t);
    tot=n*n;
    for (int i=1;i<=m;i++) {
        int x,y;
        scanf("%d%d",&x,&y);
        G[x][++G[x][0]]=y; G[y][++G[y][0]]=x;
    }
    for (int i=1;i<=n;i++) {
        scanf("%lf",&p[i]);
        out[i]=(1.0-p[i])/(double)G[i][0];
        G[i][++G[i][0]]=i;                          //自环 
    }

    a[get(s,t)][tot+1]=-1;                          //概率为1 
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++) {
            a[get(i,j)][get(i,j)]--;                
            for (int s=1;s<=G[i][0];s++)            //枚举前驱结点 
                for (int l=1;l<=G[j][0];l++) {
                    int u=G[i][s];
                    int v=G[j][l];
                    if (u==v) continue;
                    if (u!=i&&v!=j) a[get(i,j)][get(u,v)]+=out[u]*out[v];
                    if (u!=i&&v==j) a[get(i,j)][get(u,v)]+=out[u]*p[v];
                    if (u==i&&v!=j) a[get(i,j)][get(u,v)]+=p[u]*out[v];
                    if (u==i&&v==j) a[get(i,j)][get(u,v)]+=p[u]*p[v];
                } 
        }

    gauss(tot);
    for (int i=1;i<=n;i++) 
        printf("%0.6lf ",fabs(a[get(i,i)][tot+1]));
    return 0;
}