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
",sum);
}
/*
g++ 1.cpp -o 1
./1
3 3 2
0 1 0
1 2
1 3
2 3
*/
View Code
Zeit und Raum trennen dich und mich.时空将你我分开。
原文地址:https://www.cnblogs.com/starsing/p/11241878.html