题链:
http://www.lydsy.com/JudgeOnline/problem.php?id=3640
题解:
期望dp,高斯消元
设dp[i][h]在i位置且血量为h这个状态的期望经过次数。
因为每当到达n点就停止游戏,所以到达终点的概率就是dp[n][1]+dp[n][2]+...+dp[n][hp]
可以按血量把dp分成若干个层次,我们希望这样分层次后就可以把问题转变为DAG上的dp,
可是存在伤害值为0的点,所以我们对于每一层列出来的n的dp计算式,是可能存在环的。
所以要用高斯消元。复杂度(hp*N^3)
注意到在每个层次做高斯消元时,其实这些方程的系数都相同,
(规定一下方程的形式:a1*x1+a2*x2+a3*x3+...+an*xn=c,a都为系数,x为n个未知数,c为常数项)
所以我们可以预处理记录下高斯消元是怎么消的,或者说,
因为每个未知数一定是由n个常数项线性组合起来的,所以我们记录一下每个未知数是如何由n个方程的常数项构成的
具体做法是把n个常数项看成一个1×n的X矩阵,然后把每个未知数xi分别也看成1×n的一个矩阵Yi
然后如果知道了常数项的矩阵,就可以由Yi*X的倒置矩阵得到一个1×1的矩阵,而里面存的值就是xi的值。
所以高斯消元就是去求出这n个Yi矩阵。
复杂度:O(hp*n^2+n^3)
更详细的膜这里:https://blog.sengxian.com/solutions/bzoj-3640
代码:
#include<bits/stdc++.h> #define MAXN 160 using namespace std; const double eps=1e-12; double ANS,a[MAXN][MAXN],*A[MAXN],dp[MAXN][10005]; int damage[MAXN],cnt[MAXN]; int N,M,HP; int dcmp(double x){ if(fabs(x)<eps) return 0; return x>0?1:-1; } struct Edge{ int ent; int to[MAXN*MAXN],nxt[MAXN*MAXN],head[MAXN]; Edge():ent(2){} void Adde(int u,int v){ to[ent]=v; nxt[ent]=head[u]; head[u]=ent++; } }E; struct Matrix{ int r,c; double a[2][MAXN]; void Reset(int _r,int _c){ r=_r; c=_c; memset(a,0,sizeof(a)); } Matrix operator - () const{ Matrix now; now.Reset(r,c); for(int i=1;i<=now.r;i++) for(int j=1;j<=now.c;j++) now.a[i][j]=-a[i][j]; return now; } Matrix operator + (const Matrix &rtm) const{ Matrix now; now.Reset(r,c); for(int i=1;i<=now.r;i++) for(int j=1;j<=now.c;j++) now.a[i][j]=a[i][j]+rtm.a[i][j]; return now; } Matrix operator - (const Matrix &rtm) const{ return *this+(-rtm); } Matrix operator * (const double k) const{ Matrix now; now.Reset(r,c); for(int i=1;i<=now.r;i++) for(int j=1;j<=now.c;j++) now.a[i][j]=a[i][j]*k; return now; } Matrix operator & (const Matrix &rtm) const{ //乘上rtm的倒置矩阵 Matrix now; now.Reset(r,rtm.r); for(int i=1;i<=now.r;i++) for(int j=1;j<=now.c;j++) for(int k=1;k<=c;k++) now.a[i][j]+=a[i][k]*rtm.a[j][k]; return now; } Matrix operator / (const double k) const{ return *this*(1/k); } }B[MAXN],c[MAXN],X; Matrix *C[MAXN]; void buildequation(){ for(int i=1;i<=N;i++){ a[i][i]=-1; if(damage[i]) continue; for(int e=E.head[i];e;e=E.nxt[e]){ int j=E.to[e]; if(j==N) continue; a[i][j]+=1.0/cnt[j]; } } for(int i=1;i<=N;i++){ B[i].Reset(1,N); A[i]=a[i]; c[i].Reset(1,N); c[i].a[1][i]=1; C[i]=&c[i]; } } void Gausselimination(int pos,int i){ if(pos==N+1||i==N+1) return; for(int j=pos;j<=N;j++) if(dcmp(A[j][i])!=0){ swap(A[pos],A[j]); swap(C[pos],C[j]); break; } if(dcmp(A[pos][i])!=0) for(int j=pos+1;j<=N;j++){ double k=A[j][i]/A[pos][i]; for(int l=i;l<=N;l++) A[j][l]-=A[pos][l]*k; (*C[j])=(*C[j])-(*C[pos])*k; } Gausselimination(pos+(dcmp(A[pos][i])!=0),i+1); if(dcmp(A[pos][i])!=0){ for(int l=i+1;l<=N;l++) B[i]=B[i]+B[l]*A[pos][l]; B[i]=(*C[pos])-B[i]; B[i]=B[i]/A[pos][i]; } } int main(){ ios::sync_with_stdio(0); cin>>N>>M>>HP; for(int i=1;i<=N;i++) cin>>damage[i]; for(int i=1,u,v;i<=M;i++){ cin>>u>>v; E.Adde(u,v),cnt[u]++; if(u!=v) E.Adde(v,u),cnt[v]++; } buildequation(); Gausselimination(1,1); for(int h=HP;h>=1;h--){ X.Reset(1,N); for(int i=1;i<=N;i++){ if(!damage[i]) continue; if(h+damage[i]>HP) continue; for(int e=E.head[i];e;e=E.nxt[e]){ int j=E.to[e]; if(j==N) continue; X.a[1][i]-=dp[j][h+damage[i]]*1.0/cnt[j]; } } if(h==HP) X.a[1][1]-=1; for(int i=1;i<=N;i++) dp[i][h]=(B[i]&X).a[1][1]; ANS+=dp[N][h]; } cout<<fixed<<setprecision(8)<<ANS<<endl; return 0; }