BZOJ3640 JC的小苹果 概率DP+高斯消元灵活运用 题解

题目描述

让我们继续JC和DZY的故事。

“你是我的小丫小苹果，怎么爱你都不嫌多！”

“点亮我生命的火，火火火火火！”

话说JC历经艰辛来到了城市B，但是由于他的疏忽DZY偷走了他的小苹果！没有小苹果怎么听歌！他发现邪恶的DZY把他的小苹果藏在了一个迷宫里。JC在经历了之前的战斗后他还剩下hp点血。开始JC在1号点，他的小苹果在N号点。DZY在一些点里放了怪兽。当JC每次遇到位置在i的怪兽时他会损失Ai点血。当JC的血小于等于0时他就会被自动弹出迷宫并且再也无法进入。

但是JC迷路了，他每次只能从当前所在点出发等概率的选择一条道路走。所有道路都是双向的，一共有m条，怪兽无法被杀死。现在JC想知道他找到他的小苹果的概率。

P.S.大家都知道这个系列是提高组模拟赛，所以这是一道送分题balabala

输入格式

第一行三个整数表示n，m，hp。
接下来一行整数，第i个表示jc到第i个点要损失的血量。保证第1个和n个数为0。
接下来m行每行两个整数a，b表示ab间有一条无向边。

输出格式

仅一行，表示JC找到他的小苹果的期望概率，保留八位小数。

样例

样例输入

 3 3 2
 0 1 0
 1 2
 1 3
 2 3

样例输出

0.87500000

数据范围与提示

对于100%的数据 2<=n<=150，hp<=10000，m<=5000，保证图联通。

　　首先看到概率期望，我们要分辨它是真是假，好吧它是真的 。。。在看一眼无向图，完了消元。。看n的范围，好像有谱。。。先把方程搞出来：$f[w,x]$表示到x点还剩w点血量的概率，那转移也是很套路的，$f[w,x]=\sum \limits_{x to y}\frac{f[w+a[x]]}{du[y]}$，这个应该很好理解，然后。。。。。我*,这**不是$O(hpn^{3})$的吗，然后就没然后了，接着想复杂度里$O(hp)$一定要有，就是那个$n^{3}$太难看了。。。突然想到是不是高斯消元可以乱搞一下，因为不管血量多少，图都是长一样的，那高斯消元的系数矩阵都是一样的。理论是应该是一些模板灵活运用的题目，然后就不会了。koala学长讲过以后，突然明白了，同时对高斯消元有了一定的理解。

　　$AX=B$首先这是线性方程组的矩阵写法，$A$表示$n*n$系数矩阵，$X$表示$n*1$的未知量矩阵，$B$表示右侧方程的值。然后就可以yy一下消元的过程，不断变化A使之成为单位矩阵，右边B剩啥就是啥了。其实我们可以把右边看做一个$n*n$的单位矩阵和$B$的乘积，那么对左边$A$进行一系列操作，右边对单位矩阵做同样的操作，然后将右边剩下的矩阵和$B$乘起来，这是一样的。然而这样两种思路已经有了区别，第二种可以形象的写为$X=B*A^{-1}$，就是当操作完后，右边剩下的矩阵其实是A的逆矩阵，这个东西在题里是一成不变的，那么我们在循环外面就可以处理出来，然后循环里面直接将$B$与$A$向乘，由于$B$是一列，复杂度只有$n^{2}$，那么总复杂度就能变成$O(n^{3}+hp*n^{2})$，算一下的话是1e9级别，可以接受。

　　附上代码，没有看懂实现的同学想一想再看

　　

#include<iostream>
#include<cmath>
#include<cstdlib>
#include<algorithm>
#include<vector>
#include<cstdio>
#include<cstring>
using namespace std;
const int N=200,HP=10020,M=5020;
int du[N],t[N],c[N],fr[M*2],tt;
double f[HP][N],a[N][N],b[N][N],p[N],ans[N];
struct node{int fr,to,pr;}mo[M*2];
int rd()
{
    char cc=getchar();
    int s=0,w=1;
    while(cc<'0'||cc>'9') {if(cc=='-') w=-1;cc=getchar();}
    while(cc>='0'&&cc<='9') s=(s<<3)+(s<<1)+(cc^48),cc=getchar();
    return s*w;
}
void add(int x,int y)
{
    mo[++tt].fr=x;mo[tt].to=y;
    mo[tt].pr=fr[x];fr[x]=tt;
}
void solve(int n)
{
    for(int i=1;i<=n;i++)
    {
        int k=i;
        for(int j=i+1;j<=n;j++)if(fabs(a[j][i])>fabs(a[k][i]))k=j;
        if(k!=i)for(int j=1;j<=n;j++)swap(a[i][j],a[k][j]),swap(b[i][j],b[k][j]);
        for(int j=i+1;j<=n;j++)
        {
            double t=a[j][i]/a[i][i];
            for(int k=1;k<=n;k++) a[j][k]-=a[i][k]*t,b[j][k]-=b[i][k]*t;
        }
    }
    for(int i=n;i;i--)
    {
        for(int j=n;j>i;j--)
        {
            double t=a[i][j];
            for(int k=1;k<=n;k++)a[i][k]-=a[j][k]*t,b[i][k]-=b[j][k]*t;
        }
        double t=a[i][i];
        for(int j=1;j<=n;j++)a[i][j]/=t,b[i][j]/=t;
    }
}
int main()
{
    int n=rd(),m=rd(),hp=rd(),cnt=0;
    double sum=0;
    for(int i=1;i<=n;i++)c[i]=rd();
    for(int i=1,x,y;i<=m;i++)
    {
        x=rd(),y=rd();
        add(x,y);du[x]++;
        if(x!=y)add(y,x),du[y]++;
    }
    for(int i=1;i<=n;i++)
    {
        a[i][i]++,b[i][i]=1;
        if(i==n) break; 
        for(int j=fr[i];j;j=mo[j].pr) if(!c[mo[j].to]) a[mo[j].to][i]-=1.0/du[i];
    }
    solve(n);
    f[hp][1]=1;
    for(int w=hp;w>0;w--)
    {
        memset(ans,0,sizeof(ans));
        for(int j=1;j<=n;j++)
            for(int k=1;k<=n;k++)
                ans[j]+=b[j][k]*f[w][k];
        for(int j=1;j<=n;j++)f[w][j]=ans[j];
        for(int j=1;j<n;j++)
            for(int k=fr[j];k;k=mo[k].pr)
                if(c[mo[k].to]&&w>c[mo[k].to])
                    f[w-c[mo[k].to]][mo[k].to]+=f[w][j]/du[j];
        sum+=f[w][n];
    }
    printf("%.8lf\n",sum);
}
/*
g++ 1.cpp -o 1
./1
3 3 2
0 1 0
1 2
1 3
2 3
*/