bzoj3640 JC的小苹果

题目描述：

bz

题解：

有重边有自环

（我邻接矩阵死了）

高斯消元矩阵求逆。

设状态$dp[i][j]$指走到点$i$时正好有$j$点血的概率。

首先有个大暴力，把$n*hp$个状态放在一起跑高消。

当然会死。

然后发现，方程是：$$dp[v][h]=\sum_{(u,v)}{dp[u][h+w[v]]/ind[u]}$$

我们可以：$$dp[v][h]-\sum_{(u,v)}{dp[u][h]/ind[u]}=0,w[v]==0$$
$$dp[v][h]=\sum_{(u,v)}{dp[u][h+w[v]]/ind[u]},w[v]!=0$$

左边系数矩阵是确定的，右边一直在变。

所以我们可以求左边系数矩阵的逆矩阵，倒序枚举$hp$然后$O(n^2)$相乘，这样时间复杂度就是$O(hp*n^2+n^3)$的了。

代码：

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int N = 200;
template<typename T>
inline void read(T&x)
{
    T f = 1,c = 0;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){c=c*10+ch-'0';ch=getchar();}
    x = f*c;
}
int n,m,hp,w[N],ind[N];
int hed[N],cnt;
struct EG
{
    int to,nxt;
}e[10050];
void ae(int f,int t)
{
    e[++cnt].to = t;
    e[cnt].nxt = hed[f];
    hed[f] = cnt;
}
double a[N][N],b[N][N];
void gs()
{
    for(int i=1;i<=n;i++)
    {
        int tmp = i;
        for(int j=i+1;j<=n;j++)
            if(fabs(a[j][i])>fabs(a[tmp][i]))tmp=j;
        if(tmp!=i)
        {
            for(int j=i;j<=n;j++)swap(a[i][j],a[tmp][j]);
            for(int j=1;j<=n;j++)swap(b[i][j],b[tmp][j]);
        }
        double now = a[i][i];
        for(int j=i;j<=n;j++)a[i][j]/=now;
        for(int j=1;j<=n;j++)b[i][j]/=now;
        for(int j=1;j<=n;j++)if(j!=i)
        {
            now = a[j][i];
            for(int k=i;k<=n;k++)a[j][k]-=now*a[i][k];
            for(int k=1;k<=n;k++)b[j][k]-=now*b[i][k];
        }
    }
}
double dp[N][10050],D[N];
int main()
{
//    freopen("tt.in","r",stdin);
    read(n),read(m),read(hp);
    for(int i=1;i<=n;i++)
        read(w[i]);
    for(int f,t,i=1;i<=m;i++)
    {
        read(f),read(t);
        ind[f]++,ae(f,t);
        if(f!=t)ind[t]++,ae(t,f);
    }
    for(int i=1;i<=n;i++)
    {
        a[i][i]=1.0;
        if(!w[i])for(int j=hed[i];j;j=e[j].nxt)
        {
            int to = e[j].to;
            if(to!=n)a[i][to]-=1.0/ind[to];
        }
    }
    for(int i=1;i<=n;i++)
        b[i][i]=1;
    gs();
    for(int h=hp;h>=1;h--)
    {
        for(int j=1;j<=n;j++)D[j]=0;
        if(h==hp)D[1]=1;
        for(int i=1;i<n;i++)if(w[i]&&h+w[i]<=hp)
            for(int j=hed[i];j;j=e[j].nxt)
            {
                int to = e[j].to;
                if(to==n)continue;
                D[i]+=dp[to][h+w[i]]/ind[to];
            }
        for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)dp[i][h]+=D[j]*b[i][j];
    }
    double ans = 0;
    for(int i=1;i<=hp;i++)
        ans+=dp[n][i];
    printf("%.8f\n",ans);
    return 0;
}