●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;
}