正文
●BZOJ 3640 JC的小苹果
小程序：扫一扫查出行【扫一扫了解最新限行尾号】
复制小程序
题链：
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;
}
