BZOJ 3270: 博物馆 [概率DP 高斯消元]

http://www.lydsy.com/JudgeOnline/problem.php?id=3270

---

题意：一张无向图，一开始两人分别在$x$和$y$，每一分钟在点$i$不走的概率为$p[i]$,走的话等概率走到相邻的点，求两人在每个点相遇的概率
对于100%的数据有 n <= 20，n-1 <= m <= n(n-1)/2

---

因为两个人，所以状态肯定要二元组呀
$f(i,j)$表示一人在$i$另一人在$j$的概率，转移方程：
$f(i,j)=f(i,j)p_ip_j-sumlimits_{(x,i) , (y,j) in E}{f(x,i)p_jt_x+f(i,y)p_it_y+f(x,y)t_xt_y}$
$t_i=frac{1-p_i}{d_i}$
然后有环没法$DP$.....
高斯消元解方程...$n^2$个方程$n^2$个变量

注意：
$1.$ $f(i,i)$不能走到其他啦(也不能不走啦)
$2.$ 一开始$f(x,y)$的概率除了走出来的左面还要加上$1$，因为一开始一定在$(x,y)$

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
const int N=405,M=405;
inline int read(){
    char c=getchar();int x=0;
    while(c<'0'||c>'9'){c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x;
}
int n,m,x,y,u,v,nn;
int d[N];
double p[N],a[N][N];
inline int id(int i,int j){return (i-1)*n+j;}
struct edge{
    int v,ne;
}e[M<<1];
int h[N],cnt=0;
inline void ins(int u,int v){
    cnt++;
    e[cnt].v=v;e[cnt].ne=h[u];h[u]=cnt;
    cnt++;
    e[cnt].v=u;e[cnt].ne=h[v];h[v]=cnt;
}
double t[N];
void buildEquation(){
    for(int i=1;i<=n;i++) t[i]=(1-p[i])/d[i];
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++){
            int num=id(i,j);double *g=a[num];
            if(i==j) g[num]=1;
            else g[num]=1-p[i]*p[j];
            int x,y;
            for(int k=h[i];k;k=e[k].ne){
                x=e[k].v;
                if(x!=j) g[id(x,j)]=-p[j]*t[x];
            }
            for(int k=h[j];k;k=e[k].ne){
                y=e[k].v;
                if(y!=i) g[id(i,y)]=-p[i]*t[y];
            }
            for(int k=h[i];k;k=e[k].ne)
                for(int l=h[j];l;l=e[l].ne){
                    x=e[k].v,y=e[l].v;
                    if(x!=y) g[id(x,y)]=-t[x]*t[y];
                }
        }
    a[id(x,y)][nn+1]=1;
}
void GaussElimination(int n){
    for(int i=1;i<=n;i++){
        int r=i;
        for(int j=i+1;j<=n;j++)
            if(abs(a[j][i])>abs(a[r][i])) r=j;
        if(r!=i) for(int k=1;k<=n+1;k++) swap(a[i][k],a[r][k]);

        for(int j=i+1;j<=n;j++){
            double t=a[j][i]/a[i][i];
            for(int k=i;k<=n+1;k++) a[j][k]-=a[i][k]*t;
        }
    }
    for(int i=n;i>=1;i--){
        for(int j=n;j>i;j--) a[i][n+1]-=a[j][n+1]*a[i][j];
        a[i][n+1]/=a[i][i];
    }
}
int main(){
    freopen("in","r",stdin);
    n=read();m=read();nn=n*n;
    x=read();y=read();
    for(int i=1;i<=m;i++) 
        u=read(),v=read(),ins(u,v),d[u]++,d[v]++;
    for(int i=1;i<=n;i++) scanf("%lf",&p[i]);
    buildEquation();
    GaussElimination(nn);
    for(int i=1;i<=n;i++) printf("%.6lf ",a[id(i,i)][nn+1]);
}