BZOJ 3640: JC的小苹果

貌似之前在那次讲题的时候听到过这题，结果现在还是只能想点暴力的做法的说

首先容易设出一个DP，(f_{i,j})表示还剩(i)滴血时在(j)点的概率，显然(ans=sum_{i=1}^{hp} f_{i,n})

然后我们根据走路可以写出一个转移：

[f_{i-w_v,v}+=frac{f_{i,u}}{deg_u}; u ot =n,(u,v)in E ]

乍一看woc真TM简单，直接分层DP就好了。然后写道一半发现(w_i)会等于(0)，那这个转移就成环了

不慌不慌成环都是套路，一波大力高斯消元就好了，然后就有了一个优秀的(O(n^3 imes hp))的做法，成功GG

那么这个算法那里可以优化呢，相信明眼人都能看出来那个高斯消元的地方很浪费，每次都要重新做，因为其实转移的系数矩阵其实都是一样的，只是后面的常数项不同

然后就有了一种暴力的做法，我们刚开始只做一次高斯消元，然后把路径记录下来，每次只要把常数项按着路径走一遍即可，消元的部分复杂度就是(O(n^2))的了

这种做法详细可以参照Sengxian's Blog，因为我没写这种做法

考虑一种更加优美的做法，我们考虑把每一层看作一个列向量(X)（行向量也行，个人习惯），然后(A)是转移矩阵，(B)是答案矩阵（常数矩阵），那么有：

[A imes X=BLeftrightarrow X=A^{-1} imes B ]

显然(A)每一层都是一样的，(B)也可以由之前的推导出来，那么我们只要把(A)的逆求出来即可

EXPAND：本人从来没写过矩阵求逆，这里简单地讲一下做法

首先我们对于(n)阶矩阵(A)，构造一个增广矩阵(W=[A|E])（表示在(A)的右边放上一个(n)阶的单位矩阵，注意是把矩阵变成(n imes 2n)的而不是乘上去）

然后我们对于(W)进行行变换，把它变成([E|B])的形式，那么此时(B)就是(A)的逆

证明如下：

因为我们现在仅对增广矩阵(W)进行行变换，那么这个变化过程其实可以看作左乘了一个矩阵(P)，满足(PW=P[A|E]=[PA|P]=[E|B])，那么我们可以得出(PA=E)且(P=B)，则(AB=E,B=A^{-1})

那么实现的过程也很简单了，只要对(A)进行高斯消元，同时对于伴随矩阵(B)一同做变换即可，具体可参考代码

然后这样复杂度就是(O(n^3+n^2 imes hp))的了，足以通过此题

PS：后来翻了翻《线性代数》，发现对方程的系数矩阵求逆本来就是一种解方程的方法，因此上面的两种做法本质是一样的

PPS：TMD注意自环，度数只用加(1)就好了！！！

#include<cstdio>
#include<iostream>
#include<cmath>
#define RI register int
#define CI const int&
using namespace std;
const int N=155,M=10005;
struct edge
{
    int to,nxt;
}e[M]; int n,m,hp,cnt,head[N],x,y,w[N],deg[N]; double ans,f[M][N],A[N][N],B[N][N],C[N];
inline void addedge(CI x,CI y)
{
    e[++cnt]=(edge){y,head[x]}; head[x]=cnt; ++deg[x];
    if (x==y) return; e[++cnt]=(edge){x,head[y]}; head[y]=cnt; ++deg[y];
}
inline void Inverse_Matrix(CI n)
{
    RI i,j,k; double tp; for (i=1;i<=n;++i) B[i][i]=1;
    for (i=1;i<=n;++i)
    {
        for (k=i,j=i+1;j<=n;++j) if (fabs(A[j][i])>fabs(A[k][i])) k=j;
        for (swap(A[i],A[k]),swap(B[i],B[k]),j=i+1;j<=n;++j)
        for (tp=A[j][i]/A[i][i],k=1;k<=n;++k)
        A[j][k]-=A[i][k]*tp,B[j][k]-=B[i][k]*tp;
    }
    for (i=n;i;--i)
    {
        for (j=n;j>i;--j) for (tp=A[i][j],k=1;k<=n;++k) A[i][k]-=A[j][k]*tp,B[i][k]-=B[j][k]*tp;
        for (tp=A[i][i],j=1;j<=n;++j) A[i][j]/=tp,B[i][j]/=tp;
    }
}
int main()
{
    RI i,j,k; for (scanf("%d%d%d",&n,&m,&hp),i=1;i<=n;++i) scanf("%d",&w[i]);
    for (i=1;i<=m;++i) scanf("%d%d",&x,&y),addedge(x,y);
    for (i=1;i<=n;++i) if (A[i][i]=1,i!=n) for (j=head[i];j;j=e[j].nxt)
    if (!w[e[j].to]) A[e[j].to][i]-=1.0/deg[i];
    for (Inverse_Matrix(n),f[hp][1]=1,i=hp;i;--i)
    {
        for (j=1;j<=n;++j) for (C[j]=0,k=1;k<=n;++k) C[j]+=B[j][k]*f[i][k];
        for (ans+=C[n],j=1;j<n;++j) for (k=head[j];k;k=e[k].nxt)
        if (w[e[k].to]&&i>w[e[k].to]) f[i-w[e[k].to]][e[k].to]+=C[j]/deg[j];
    }
    return printf("%.8lf",ans),0;
}