正文

【BZOJ3640】JC的小苹果 概率DP+高斯消元

IP云
IP云V管理员/3 阅读
12/02
提示：扫一扫查出行
【扫一扫了解最新限行尾号】
复制提示

【BZOJ3640】JC的小苹果

Description

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

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

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

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

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

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

Input

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

Output

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

Sample Input
3 3 2
0 1 0
1 2
1 3
2 3

Sample Output
0.87500000

HINT

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

题解 ：如果没有Ai=0，直接DP，如果hp很小，直接高斯消元，但是都有，所以用DP+高斯消元。

发现，方程组中只有Ai=0的点有系数，Ai！=0的都可以直接拿过来变成常数项，所以每次我们的方程组只有一列是变化的，所以我们将方程组表示成Ax=b，x=A -1 b。x和b都是列向量。所以我们只需要预处理出A的逆，然后每次O（n 2 ）乘一下就行了。

矩阵求逆的方法：先将(A|I)拼一起，然后通过行变换将左边的A消成I，右边剩下的就是A -1 。 

#include <cstdio>

#include <cstring>

#include <iostream>

#include <cmath>

using namespace std;

int n,m,hp;

int dam[160],pa[5010],pb[5010],d[160];

double f[160][10000],B[160],ans;

struct M

{

	double v[160][160];

	M (){memset(v,0,sizeof(v));}

	double* operator [](int x) {return v[x];}

	void I()

	{

		for(int i=1;i<=n;i++)	for(int j=1;j<=n;j++)	v[i][j]=(i==j)?1:0;

	}

	M getinv()

	{

		int i,j,k;

		M re;

		double t;

		re.I();

		for(i=1;i<=n;i++)

		{

			for(j=i;j<=n;j++)	if(fabs(v[j][i])>fabs(v[i][i]))	for(k=1;k<=n;k++)

				swap(v[i][k],v[j][k]),swap(re[i][k],re[j][k]);

			t=v[i][i];

			for(j=1;j<=n;j++)	v[i][j]/=t,re[i][j]/=t;

			for(j=1;j<=n;j++)	if(i!=j)

			{

				t=v[j][i];

				for(k=1;k<=n;k++)	v[j][k]-=t*v[i][k],re[j][k]-=t*re[i][k];

			}

		}

		return re;

	}

	void operator * (int x)

	{

		for(int i=1;i<=n;i++)	for(int j=1;j<=n;j++)	f[i][x]+=v[i][j]*B[j];

	}

};

M A,A1;

int rd()

{

	int ret=0,f=1;	char gc=getchar();

	while(gc<'0'||gc>'9')	{if(gc=='-')f=-f;	gc=getchar();}

	while(gc>='0'&&gc<='9')	ret=ret*10+gc-'0',gc=getchar();

	return ret*f;

}

int main()

{

	n=rd(),m=rd(),hp=rd();

	int i,j;

	for(i=1;i<=n;i++)	dam[i]=rd();

	for(i=1;i<=m;i++)

	{

		pa[i]=rd(),pb[i]=rd();

		d[pa[i]]++;

		if(pa[i]!=pb[i])	d[pb[i]]++;

	}

	for(i=1;i<=m;i++)

	{

		if(!dam[pb[i]])	A[pb[i]][pa[i]]-=1.0/d[pa[i]];

		if(pa[i]==pb[i])	continue;

		if(!dam[pa[i]])	A[pa[i]][pb[i]]-=1.0/d[pb[i]];

	}

	for(i=1;i<=n;i++)	A[i][n]=0;

	for(i=1;i<=n;i++)	A[i][i]++;

	A1=A.getinv();

	for(j=hp;j;j--)

	{

		for(i=1;i<=n;i++)	B[i]=0;

		if(j==hp)	B[1]=1;

		else	for(i=1;i<=m;i++)

		{

			if(dam[pb[i]]&&dam[pb[i]]+j<=hp&&pa[i]!=n)	B[pb[i]]+=f[pa[i]][j+dam[pb[i]]]*1.0/d[pa[i]];

			if(pa[i]==pb[i])	continue;

			if(dam[pa[i]]&&dam[pa[i]]+j<=hp&&pb[i]!=n)	B[pa[i]]+=f[pb[i]][j+dam[pa[i]]]*1.0/d[pb[i]];

		}

		A1*j,ans+=f[n][j];

	}

	printf("%.8lf",ans);

	return 0;

}

相关阅读

目录[+]