[高斯消元] BZOJ3640: JC的小苹果

不难想到的题。首先设 f[i][j] 表示血量为 i，位置在 j 的期望概率。

推出方程之后发现，当一个点 a[i]=0 时，不能递推出来。所以要每层进行消元。

发现每次消元的系数矩阵都相同，只是常数不同。可以想到把常数看成一个多项式，先 n3 预处理，然后每层直接 n2 带入就好了。

复杂度 O(n3+n2∗hp)

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int maxn=155,maxe=10005,maxhp=20005;
int n,m,N,hp,a[maxn],fir[maxn],nxt[maxe],son[maxe],tot;
struct data{
double a[maxn];
data(){ memset(a,0,sizeof(a)); }
data operator - (const data &B)const{
data c; for(int i=0;i<=n;i++) c.a[i]=a[i]-B.a[i]; return c;
}
data operator * (const double &val)const{
data c; for(int i=0;i<=n;i++) c.a[i]=a[i]*val; return c;
}
data operator / (const double &val)const{
data c; for(int i=0;i<=n;i++) c.a[i]=a[i]/val; return c;
}
} _K[maxn],Ans[maxn];
double K[maxn][maxn];
void Gs(){
for(int i=1;i<=n;i++){
int where; double _max=0; for(int j=i;j<=n;j++) if(K[j][i]>_max) where=j, _max=K[j][i];
for(int j=1;j<=n;j++) swap(K[i][j],K[where][j]); swap(_K[i],_K[where]);
for(int j=i+1;j<=n;j++){
double t=K[j][i]/K[i][i];
for(int k=1;k<=n;k++) K[j][k]-=t*K[i][k]; _K[j]=_K[j]-_K[i]*t;
}
}
for(int i=n;i>=1;i--){
Ans[i]=_K[i]*(-1);
for(int j=i+1;j<=n;j++) Ans[i]=Ans[i]-Ans[j]*K[i][j];
Ans[i]=Ans[i]/K[i][i];
}
}
int d[maxn];
double ans,f[maxhp][maxn];
void add(int x,int y){
son[++tot]=y; nxt[tot]=fir[x]; fir[x]=tot;
}
int main(){
scanf("%d%d%d",&n,&m,&hp);
if(hp==0) return printf("0.00000000"),0;
for(int i=1;i<=n;i++) scanf("%d",&a[i]);
for(int i=1;i<=m;i++){
int x,y; scanf("%d%d",&x,&y);
add(x,y); d[x]++; if(x!=y) add(y,x), d[y]++;
}
for(int i=1;i<=n;i++){
K[i][i]+=1; if(a[i]) continue;
if(i==1) _K[1].a[0]-=1;
for(int j=fir[i];j;j=nxt[j]) if(son[j]!=n) K[i][son[j]]-=1.0/d[son[j]];
}
Gs();
for(int i=1;i<=n;i++) f[hp][i]=Ans[i].a[0];
ans=f[hp]
;
memset(K,0,sizeof(K)); memset(Ans,0,sizeof(Ans)); memset(_K,0,sizeof(_K));
for(int i=1;i<=n;i++){
K[i][i]+=1;
if(a[i]){ _K[i].a[i]-=1; continue; }
for(int j=fir[i];j;j=nxt[j]) if(son[j]!=n) K[i][son[j]]-=1.0/d[son[j]];
}
Gs();
for(int i=hp-1;i>=1;i--){
for(int j=1;j<=n;j++) if(a[j]){
if(i+a[j]>hp) continue;
for(int k=fir[j];k;k=nxt[k]) if(son[k]!=n) f[i][j]+=f[i+a[j]][son[k]]/d[son[k]];
}
for(int j=1;j<=n;j++) if(!a[j]){
for(int k=1;k<=n;k++) f[i][j]+=Ans[j].a[k]*f[i][k];
}
ans+=f[i]
;
}
printf("%.8lf\n",ans);
//  for(int i=hp;i>=1;i--){
//          for(int j=1;j<=n;j++) printf("%.3lf ",f[i][j]);
//          printf("\n");
//      }
return 0;
}