bzoj3640 JC的小苹果（期望dp+Gauss+矩阵求逆）

我们有比较显然的期望dp，f[x][i]表示到x点血量为i的期望次数。
我们有转移$f[x][i]=\sum_y\frac{f[y][i+w[x]]}{du[y]}$
直接高斯消元复杂度$O((nhp)^3)$gg
由于血量不增，因此我们可以按剩余血量分层dp。
可是存在w[x]=0的点，每层内部还是需要一个高斯消元。复杂度$O(hpn^3)$
还是gg
我们发现每层内的高斯消元其实系数矩阵是相同的，只有常数矩阵在变。
即每次都是根据$Ax=B$求解列向量x，那么我们有$x=A^{-1}B$，即我们求出A矩阵的逆，每次只需把逆矩阵与常数矩阵乘起来即可得到x。因为常数矩阵也是一个列向量，所以每次的复杂度是$O(n^2)$的。
这样复杂度就优化到了$O(n^3+hpn^2)$，可以通过。

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define inf 0x3f3f3f3f
#define N 160
inline char gc(){
    static char buf[1<<16],*S,*T;
    if(T==S){T=(S=buf)+fread(buf,1,1<<16,stdin);if(T==S) return EOF;}
    return *S++;
}
inline int read(){
    int x=0,f=1;char ch=gc();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=gc();}
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=gc();
    return x*f;
}
int n,m,K,h[N],num=0,du[N],w[N];
double f[N][10010],a[N][N],b[N],c[N][N],ans=0;
struct edge{
    int to,next;
}data[10010];
inline void inv(){
    for(int i=1;i<=n;++i) c[i][i]=1;
    for(int i=1;i<=n;++i){
        int r=i;
        for(int j=i+1;j<=n;++j) if(fabs(a[j][i])>fabs(a[r][i])) r=j;
        if(r!=i) for(int j=1;j<=n;++j) swap(a[i][j],a[r][j]),swap(c[i][j],c[r][j]);
        double t=a[i][i];
        for(int j=1;j<=n;++j) a[i][j]/=t,c[i][j]/=t;
        for(int j=1;j<=n;++j){
            if(j==i) continue;t=a[j][i];
            for(int k=1;k<=n;++k) a[j][k]-=a[i][k]*t,c[j][k]-=c[i][k]*t;
        }
    }
}
int main(){
//  freopen("a.in","r",stdin);
    n=read();m=read();K=read();
    for(int i=1;i<=n;++i) w[i]=read();
    while(m--){
        int x=read(),y=read();
        data[++num].to=y;data[num].next=h[x];h[x]=num;du[y]++;
        if(x==y) continue;
        data[++num].to=x;data[num].next=h[y];h[y]=num;du[x]++;
    }for(int x=1;x<=n;++x){
        a[x][x]=1;if(w[x]) continue;
        for(int i=h[x];i;i=data[i].next){
            int y=data[i].to;if(y==n) continue;
            a[x][y]-=1.0/du[y];
        }
    }inv();
    for(int k=K;k>=1;--k){
        memset(b,0,sizeof(b));if(k==K) b[1]=1;
        for(int x=1;x<=n;++x){
            if(!w[x]||w[x]+k>K) continue;
            for(int i=h[x];i;i=data[i].next){
                int y=data[i].to;if(y==n) continue;
                b[x]+=f[y][w[x]+k]/du[y];
            }
        }for(int i=1;i<=n;++i)
            for(int j=1;j<=n;++j)
                f[i][k]+=c[i][j]*b[j];ans+=f[n][k];
    }printf("%.8lf\n",ans);
    return 0;
}