Bzoj 3640 JC's Little Apple

Chain of title:

http://www.lydsy.com/JudgeOnline/problem.php?id=3640
Exercises

expected DP, Gaussian elimination
set Dp[i][h] in the position of I and the amount of blood is h this state of the expected number of times.
because every time you arrive at N, you stop the game, so the probability of reaching the end is DP[N][1]+DP[N][2]+...+DP[N][HP]
The DP can be divided into several levels according to the amount of blood, and we hope that the problem can be transformed into a DP on the DAG at this level .
However, there is a point of 0 damage, so we may have a ring for the DP formula of N which is listed in each layer.
so we need to use Gauss to eliminate the element. Complexity (HP*N^3)
Notice that the coefficients of these equations are the same at each level of the Gaussian elimination element .
(Specify the form of the equation: a1*x1+a2*x2+a3*x3+...+an*xn=c,a are coefficients, X is n unknowns, C is a constant term)
so we can pre-record how the Gaussian elimination is eliminated, or,
since each unknown must be a linear combination of N constants, let's record how each unknown is made up of a constant term of n equations .
The method is to treat n constant term as a 1xn x matrix, and then consider each unknown XI as a matrix Yi of 1xn respectively .
then, if you know the matrix of a constant term, you can get a 1x1 matrix from the inverted matrix of the yi*x, and the value stored in it is the value of XI.
so Gaussian elimination is to find out the N yi Matrix.
complexity: O (hp*n^2+n^3)
more detailed film here: https://blog.sengxian.com/solutions/bzoj-3640

Code:

#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; 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 &AMP;RTM) Const{matrix 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 &AMP;RTM) Const{return *this+ (-RTM);} Matrix operator * (const double k) Const{matrix 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 &AMP;RTM) const{//The inverted matrix in RTM 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, 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]&GT;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's Little Apple