[BZOJ3270]博物館（概率+高斯消元）

2019-11-11 05:02:42

字體：大中小

來源：轉載

供稿：網友

題目描述

傳送門

題解

假設當前在點(i,j)，下一步從這個點走到它某一個相鄰的點的概率即為1?pidi $1-p_i/over d_i$ ，記為goi $go_i$ 設兩個人分別走到i,j的概率為f(i,j)，那么 f(i,j)=f(i,j)pipj+∑(i,x),(j,y)∈Ef(x,j)pjgox+f(i,y)pigoy+f(x,y)goxgoy $f(i,j)=f(i,j)p_ip_j+/sum/limits_{(i,x),(j,y)/in E}f(x,j)p_jgo_x+f(i,y)p_igo_y+f(x,y)go_xgo_y$ 特殊地，f(a,b)的初值為1 這樣得出了n*n個方程，高斯消元即可需要注意的是，方程中等式右邊在同一個f里的兩個點不能相等，因為一旦相等就已經結束，不會再有走到這個點的概率

代碼

#include<algorithm>#include<iostream>#include<cstring>#include<cstdio>#include<cmath>using namespace std;#define N 405const double eps=1e-9;int dcmp(double x){ if (x<=eps&&x>=-eps) return 0; return (x>0)?1:-1;}int n,m,A,B,x,y;int tot,point[N],nxt[N*2],v[N*2];double p[N],d[N],go[N],a[N][N],b[N],ans[N];void add(int x,int y){ ++tot; nxt[tot]=point[x]; point[x]=tot; v[tot]=y;}int id(int x,int y){ return (x-1)*n+y;}void gauss(){ for (int i=1;i<=n*n;++i) { int num=i; for (int j=i+1;j<=n*n;++j) if (dcmp(a[j][i]-a[num][i])>0) num=j; if (num!=i) { for (int j=1;j<=n*n;++j) swap(a[num][j],a[i][j]); swap(b[num],b[i]); } for (int j=i+1;j<=n*n;++j) if (dcmp(a[j][i])) { double t=a[j][i]/a[i][i]; for (int k=1;k<=n*n;++k) a[j][k]-=t*a[i][k]; b[j]-=b[i]*t; } } for (int i=n*n;i>=1;--i) { for (int j=i+1;j<=n*n;++j) b[i]-=a[i][j]*ans[j]; ans[i]=b[i]/a[i][i]; }}int main(){ scanf("%d%d%d%d",&n,&m,&A,&B); for (int i=1;i<=m;++i) { scanf("%d%d",&x,&y); d[x]+=1.0;d[y]+=1.0; add(x,y);add(y,x); } for (int i=1;i<=n;++i) { scanf("%lf",&p[i]); go[i]=(1-p[i])*1/d[i]; } for (int i=1;i<=n;++i) for (int j=1;j<=n;++j) { a[id(i,j)][id(i,j)]=1; if (i!=j) a[id(i,j)][id(i,j)]-=p[i]*p[j]; for (int k=point[i];k;k=nxt[k]) if (v[k]!=i&&v[k]!=j) a[id(i,j)][id(v[k],j)]=-p[j]*go[v[k]]; for (int k=point[j];k;k=nxt[k]) if (v[k]!=i&&v[k]!=j) a[id(i,j)][id(i,v[k])]=-p[i]*go[v[k]]; for (int k=point[i];k;k=nxt[k]) for (int l=point[j];l;l=nxt[l]) if (v[k]!=v[l]) a[id(i,j)][id(v[k],v[l])]=-go[v[k]]*go[v[l]]; } b[id(A,B)]=1.0; gauss(); for (int i=1;i<=n;++i)