[概率 高斯消元 逆矩阵] BZOJ 3640 JC的小苹果

%%%jiry_2

把点根据血量拆开进行高斯消元， O(n3hp3)。
根据 hp 可以把图分成 hp 层，第 i 层对 j(i > j) 层是没有影响的。
每层之间高斯消元，层与层之间递推， O(n3hp)。
每一次高斯消元的系数矩阵都是相同的，可以先高斯消元一次预处
理，之后消元的时候带入就行了， O(n2hp + n3)。

#include<cstdio>#include<cstdlib>#include<algorithm>#include<cmath>using namespace std;inline char nc(){static char buf[100000],*p1=buf,*p2=buf;if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin); if (p1==p2) return EOF; }return *p1++;}inline void read(int &x){char c=nc(),b=1;for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;}const int N=155;int n,m,K,A[N];double a[N][N],num[N][N];double f[15005][N],tmp[N];int d[N][N],D[N];inline void Gauss(){for (int i=1;i<=n;i++){int k=0;for (int j=i;j<=n;j++) if (fabs(a[j][i])>1e-10) { k=j; break; }if (k!=i)for (int j=1;j<=n;j++)swap(a[k][j],a[i][j]),swap(num[k][j],num[i][j]);double tmp=a[i][i];for (int j=1;j<=n;j++)a[i][j]/=tmp,num[i][j]/=tmp;for (int j=1;j<=n;j++){if (j==i) continue;tmp=a[j][i];for (int k=1;k<=n;k++)a[j][k]-=tmp*a[i][k],num[j][k]-=tmp*num[i][k];}}}int main(){int ix,iy;freopen("t.in","r",stdin);freopen("t.out","w",stdout);read(n); read(m); read(K);for (int i=1;i<=n;i++) read(A[i]);for (int i=1;i<=m;i++){read(ix); read(iy);D[ix]++,d[ix][iy]++;if (ix!=iy)D[iy]++,d[iy][ix]++;}for (int i=1;i<n;i++) for (int j=1;j<=n;j++) if (!A[j]) a[j][i]-=(double)d[i][j]/D[i];for(int i=1;i<=n;i++) a[i][i]+=1,num[i][i]=1;Gauss();f[K][1]=1;double ans=0;for (int k=K;k>=1;k--){for (int i=1;i<n;i++) for (int j=1;j<=n;j++) if(A[j] && k+A[j]<=K && D[i])f[k][j]+=f[k+A[j]][i]*d[i][j]/D[i];for (int i=1;i<=n;i++) tmp[i]=f[k][i];for (int i=1;i<=n;i++) if(!A[i]){   f[k][i]=0;   for(int j=1;j<=n;j++)    if (A[j] || j==1)    f[k][i]+=num[i][j]*tmp[j];}ans+=f[k][n];}printf("%.8lf\n",ans);return 0;}