●BZOJ 3640 JC的小蘋果

標籤：its   處理   his   http   add   [1]   .com   消元   class   

題鏈：

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 160using 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;}

　　

●BZOJ 3640 JC的小蘋果