• 【BZOJ3640】JC的小苹果(高斯消元)


    点此看题面

    • 一张(n)个点(m)条边的无向图,要从(1)号点走到(n)号点,初始体力为(hp)
    • 每当你走到编号为(i)的点时,体力都会失去(a_i),然后等概率选择当前点的一条边走出去。
    • 当体力值小于等于(0)的时候就失败了,求走到(n)号点的概率。
    • (nle150,mle5 imes10^3,hple10^4,0le a_ile hp)

    成环的概率(DP)

    (f_{k,i})表示在体力值为(k)时走到(i)的概率,显然有转移方程:

    [f_{k,i}=sum_{j=1}^{n-1}f_{k+a_i,j} imes frac{w_{i,j}}{deg_j} ]

    其中(w_{i,j})表示(i,j)之间的边数,注意此题有重边有自环。

    看起来这样就完事了,但(a_i)可能等于(0),也就是说我们不能简单地分层(DP),因为在(k)相同的状态之间也可能存在转移。

    这种时候就要套路地想到高斯消元,把这个转移式看作一个方程。

    然而,如果直接暴力这么去做显然会(T)飞,因此要考虑优化。

    高斯消元的优化

    考虑我们总共要做(hp)次高斯消元,但实际上每次高斯消元的系数都是相同的,区别只在于等号右边的值。

    因此我们可以在一开始先做一遍高斯消元,预处理出(p_{i,j})表示第(j)个式子等号右边的值对第(i)个式子等号右边的值的贡献系数。

    那么接下来每次求解就变成(O(n^2))了。

    代码:(O(n^3+n^2hp))

    #include<bits/stdc++.h>
    #define Tp template<typename Ty>
    #define Ts template<typename Ty,typename... Ar>
    #define Reg register
    #define RI Reg int
    #define Con const
    #define CI Con int&
    #define I inline
    #define W while
    #define N 150
    #define M 5000
    #define HP 10000
    #define DB double
    #define eps 1e-12
    #define add(x,y) (e[++ee].nxt=lnk[x],e[lnk[x]=ee].to=y)
    using namespace std;
    int n,m,hp,a[N+5],d[N+5],w[N+5][N+5];DB f[HP+5][N+5];
    namespace Gauss//高斯消元
    {
    	DB a[N+5][N+5],p[N+5][N+5],v[N+5],res[N+5];
    	I void Add(CI i,CI j)//用第i行去消第j行
    	{
    		DB t=-a[j][i]/a[i][i];for(RI k=1;k<=n;++k) a[j][k]+=t*a[i][k],p[j][k]+=t*p[i][k];
    	}
    	I void Init()//初始化
    	{
    		RI i,j,k;DB t;for(i=1;i<=n;++i) for(p[i][i]=1,j=i+1;j<=n;++j) Add(i,j);//从上往下消成三角形
    		for(i=n;i;--i) {for(j=1;j<=n;++j) p[i][j]/=a[i][i];for(a[i][i]=1,j=i-1;j;--j) Add(i,j);}//从下往上消得只剩对角线
    	}
    	I void Solve()//快速求解
    	{
    		RI i,j;for(i=1;i<=n;++i) for(j=1;j<=n;++j) res[i]+=p[i][j]*v[j];//根据预处理出的系数计算
    	}
    }
    int main()
    {
    	RI i,j,k,x,y;for(scanf("%d%d%d",&n,&m,&hp),i=1;i<=n;++i) scanf("%d",a+i);
    	for(i=1;i<=m;++i) scanf("%d%d",&x,&y),++w[x][y],++d[x],x^y&&(++w[y][x],++d[y]);//注意自环
    	for(i=1;i<=n;++Gauss::a[i][i],++i) if(!a[i]) for(j=1;j^n;++j) Gauss::a[i][j]=-1.0*w[i][j]/d[j];//求出系数矩阵
    	DB t=0;for(Gauss::Init(),Gauss::v[1]=1,k=hp;k;--k)
    	{
    		for(i=1;i<=n;++i) if(a[i]&&k+a[i]<=hp) for(j=1;j^n;++j) Gauss::v[i]+=f[k+a[i]][j]*w[i][j]/d[j];//求出等号右边的值
    		for(Gauss::Solve(),i=1;i<=n;++i) f[k][i]=Gauss::res[i],Gauss::v[i]=Gauss::res[i]=0;t+=f[k][n];//把值移到DP数组里
    	}return printf("%.8lf
    ",t),0;//输出答案
    }
    
  • 相关阅读:
    error while loading shared libraries: xxx.so.x"错误的原因和解决办法
    玩转html5(三)---智能表单(form),使排版更加方便
    玩转html5(二)----用canvas结合脚本在画布上画简单的图(html5又一强大功能)
    设计模式 命令模式 之 管理智能家电
    Android Volley 之自定义Request
    玩转html5(一)-----盘点html5新增的那些酷酷的input类型和属性
    玩转Web之servlet(五)---- 怎样解决servlet的线程安全问题
    Android 完美实现图片圆角和圆形(对实现进行分析)
    Android SwipeRefreshLayout 官方下拉刷新控件介绍
    scala akka基础编程
  • 原文地址:https://www.cnblogs.com/chenxiaoran666/p/BZOJ3640.html
Copyright © 2020-2023  润新知