[Bzoj3640]JC的小苹果

题意

求$1\to n$点权和$\le hp$的概率

$($有重边和自环$)$

---

题解

因为不知道在终点还剩多少血$,$所以考虑顺推

设$f[x][u]$表示$1\to u$还剩$x$点血的概率$,dg[u]$表示$u$的出度

初值$f[hp][1]=1$

可以得到

$$\forall_{v\to u\in E}f[x][u]=\frac1{dg[u]}\sum f[x+w[u]][v]$$

考虑到血量是单调不增的$;$并且点权可能是$0,$从而导致状态转移有环

所以可以考虑对$hp$分层高斯消元

复杂度$O(hp\times n^3)$显然不行

注意到图是不变的$,$即每层高斯消元$x_i$的系数都是一样的

设$X=\begin{bmatrix}x_1&x_2&\ldots&x_n\end{bmatrix},G$为系数矩阵$,B$为答案矩阵

则$X\times G=B\Rightarrow X=B\times G^{-1}$

对于每一层来说$G$都是不变的$,$并且$B$可以在$O(n^2)$时间内处理出来

所以我们只要求出$G$的逆就可在$O(hp\times n^2)$的时间内完成$dp$

求矩阵的逆要用到高斯消元$,$大家可以好好通过过背代码感性理解一下

总复杂度$O(n^3+hp\times n^2)$

因为有自环所以要注意特判

#include<bits/stdc++.h>
#define fp(i,a,b) for(register int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(register int i=a,I=b-1;i>I;--i)
#define go(u) for(register int i=fi[u],v=e[i].to;i;v=e[i=e[i].nx].to)
#define file(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
template<class T>inline bool cmax(T&a,const T&b){return a<b?a=b,1:0;}
template<class T>inline bool cmin(T&a,const T&b){return a>b?a=b,1:0;}
using namespace std;
char ss[1<<17],*A=ss,*B=ss;
inline char gc(){return A==B&&(B=(A=ss)+fread(ss,1,1<<17,stdin),A==B)?-1:*A++;}
template<class T>inline void sd(T&x){
    char c;T y=1;while(c=gc(),(c<48||57<c)&&c!=-1)if(c==45)y=-1;x=c-48;
    while(c=gc(),47<c&&c<58)x=x*10+c-48;x*=y;
}
char sr[1<<21],z[20];int C=-1,Z;
inline void Ot(){fwrite(sr,1,C+1,stdout),C=-1;}
template<class T>inline void we(T x){
    if(C>1<<20)Ot();if(x<0)sr[++C]=45,x=-x;
    while(z[++Z]=x%10+48,x/=10);
    while(sr[++C]=z[Z],--Z);sr[++C]='\n';
}
const int N=155,M=1e4+5;
const double eps=1e-12;
typedef int arr[N];
typedef double db;
struct eg{int nx,to;}e[M];
struct matrix{
    db a[N][N];
    matrix(int x=0){memset(a,0,sizeof a);if(x)fp(i,1,N-1)a[i][i]=1;}
    inline db*operator[](int x){return a[x];}
}G;
int n,m,hp;arr w,fi;db ans,b[N],dg[N],f[M][N];
inline int cmp(const db x){return fabs(x)<eps?0:x<0?-1:1;}
inline matrix inv(matrix a,int n){
    matrix b(1);db t;int mx;
    fp(i,1,n){mx=i;
        fp(j,i,n)if(cmp(a[mx][i]-a[j][i])<0)mx=j;
        if(mx^i)swap(a.a[i],a.a[mx]),swap(b.a[i],b.a[mx]);
        fp(j,1,n)if(j^i){
            t=a[j][i]/a[i][i];
            fp(k,i,n)a[j][k]-=a[i][k]*t;
            fp(k,1,n)b[j][k]-=b[i][k]*t;
        }
    }
    fp(i,1,n)fp(j,1,n)b[i][j]/=a[i][i];
    return b;
}
inline void add(int u,int v){static int ce=0;e[++ce]={fi[u],v},fi[u]=ce;}
int main(){
    #ifndef ONLINE_JUDGE
        file("s");
    #endif
    sd(n),sd(m),sd(hp);
    fp(i,1,n)sd(w[i]);int u,v;
    while(m--){
        sd(u),sd(v),add(u,v);
        if(v^u)add(v,u),++dg[v];++dg[u];
    }dg[n]=0;
    fp(i,1,n)if(dg[i])dg[i]=1.0/dg[i];
    fp(u,1,n){G[u][u]=1;if(!w[u])go(u)G[u][v]-=dg[v];}
    G=inv(G,n);
    fd(x,hp,1){
        memset(b,0,sizeof b);b[1]=x==hp;
        fp(u,1,n)if(w[u]&&x+w[u]<=hp)
            go(u)b[u]+=dg[v]*f[x+w[u]][v];
        fp(u,1,n)fp(v,1,n)
            f[x][u]+=G[u][v]*b[v];
        ans+=f[x][n];
    }
    printf("%.8lf\n",ans);
return Ot(),0;
}