bzoj 3270: 博物馆【dp+高斯消元】

好像是高斯消元解互相推(?)的dp的例子
首先考虑dp，设f[i][j]为一人在i一人在j的概率，点i答案显然就是f[i][i]；
然后根据题意，得到转移是

\[f[i][j]=f[i][j]*p_i*p_j+\sum_{edge(x,i)\in E}f[x][j]*p_j*\frac{1-p[x]}{d[x]}+\sum_{edge(y,j)\in E}f[i][y]*p_i*\frac{1-p[y]}{d[y]}++\sum_{edge(x,i)\in E,edge(y,j)\in E}f[x][y]*\frac{1-p[x]}{d[x]}*\frac{1-p[y]}{d[y]} \]

\[f[i][j]*(p_i*p_j-1)+\sum_{edge(x,i)\in E}f[x][j]*p_j*\frac{1-p[x]}{d[x]}+\sum_{edge(y,j)\in E}f[i][y]*p_i*\frac{1-p[y]}{d[y]}++\sum_{edge(x,i)\in E,edge(y,j)\in E}f[x][y]*\frac{1-p[x]}{d[x]}*\frac{1-p[y]}{d[y]}=0 \]

这样他就可以用高斯消元解了
注意，当f[i][j]中i==j时，就不能去转移其他情况了，因为这样已经是结束状态了

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int N=405;
int n,m,sx,sy,h[N],cnt,id[25][25],tot;
double p[N],a[N][N],d[N];
struct qwe
{
	int ne,to;
}e[N<<1];
void add(int u,int v)
{
	cnt++;
	e[cnt].ne=h[u];
	e[cnt].to=v;
	h[u]=cnt;
}
void gaosi(int n)
{
	for(int i=1;i<=n;i++)
	{
		int nw=i;
		for(int j=i+1;j<=n;j++)
			if(fabs(a[nw][i])<fabs(a[j][i]))
				nw=j;
		for(int j=i;j<=n+1;j++)
			swap(a[nw][j],a[i][j]);
		for(int j=i+1;j<=n+1;j++)
			a[i][j]/=a[i][i];
		a[i][i]=1;
		for(int j=i+1;j<=n;j++)
		{
			for(int k=i+1;k<=n+1;k++)
				a[j][k]-=a[j][i]*a[i][k];
			a[j][i]=0;
		}
	}
	for(int i=n-1;i>=1;i--)
		for(int j=i+1;j<=n;j++)
			a[i][n+1]-=a[j][n+1]*a[i][j];
}
int main()
{
	scanf("%d%d%d%d",&n,&m,&sx,&sy);
	for(int i=1,x,y;i<=m;i++)
	{
		scanf("%d%d",&x,&y);
		d[x]++,d[y]++;
		add(x,y),add(y,x);
	}
	for(int i=1;i<=n;i++)
		scanf("%lf",&p[i]);
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			id[i][j]=++tot;
	for(int x=1;x<=n;x++)
		for(int y=1;y<=n;y++)
		{
			a[id[x][y]][id[x][y]]--;
			if(x!=y)
				a[id[x][y]][id[x][y]]+=p[x]*p[y];
			for(int i=h[x];i;i=e[i].ne)
				for(int j=h[y];j;j=e[j].ne)
					if(e[i].to!=e[j].to)
						a[id[x][y]][id[e[i].to][e[j].to]]+=(1-p[e[i].to])*(1-p[e[j].to])/d[e[i].to]/d[e[j].to];
			for(int i=h[x];i;i=e[i].ne)
				if(e[i].to!=y)
					a[id[x][y]][id[e[i].to][y]]+=p[y]*(1-p[e[i].to])/d[e[i].to];
			for(int i=h[y];i;i=e[i].ne)
				if(e[i].to!=x)
					a[id[x][y]][id[x][e[i].to]]+=p[x]*(1-p[e[i].to])/d[e[i].to];
		}
	a[id[sx][sy]][n*n+1]--;
	// for(int i=1;i<=n*n;i++)
	// {
		// for(int j=1;j<=n*n+1;j++)
			// printf("%.6f ",a[i][j]);
		// puts("");
	// }
	gaosi(n*n);
	for(int i=1;i<=n;i++)
		printf("%.6f ",a[id[i][i]][n*n+1]/a[id[i][i]][id[i][i]]);
	return 0;
}