[HNOI2011]XOR和路径——高斯消元、随机游走模型

最新推荐文章于 2023-11-05 07:10:33 发布

偶耶XJX

最新推荐文章于 2023-11-05 07:10:33 发布

阅读量286

点赞数

分类专栏：信息竞赛解题文章标签：算法线性代数二进制拆位随机游走高斯消元

本文链接：https://blog.csdn.net/weixin_43960287/article/details/118999461

版权

信息竞赛解题专栏收录该内容

195 篇文章 5 订阅

订阅专栏

[HNOI2011]XOR和路径

思路

拆位来看，边的权值就变为了只有0和1，而且只有经过权值为1的边的次数为奇数次时会对答案有贡献，所以我们只需算到达节点 n 时经过奇数次权值为1的边的概率即可。

先设 $P_i$ 表示从1走到 n 中经过点 i 的概率，那么应该有转移关系式
$P_i=\sum_{(i,v)\in E,v\neq n}\frac{P_v}{d_v}$
其中 $d_v$ 表示节点 v 的度数，这里的边不考虑边权。

在推概率关系的时候我们会发现一个问题，即节点 n 必须从相邻的点过来，节点 2~n-1 必须从除了n以外的相邻点过来，而节点1既可以在初始位置，又可以从不为n的相邻点过来。这样一来，应该满足 $P_1=1$ ，即 $P_1$ 不满足上面的关系式。如果不满足转移关系式，不就相当于从1走出后不能走回来了吗？

所以为了解决这个问题，我们假设1不为起点，另设一个起点0。起点0可以走到与1相邻的节点但不能走回来，节点1只能从相邻节点走过来。这样拆开后只有 $P_0=1$ 固定，其它的都可以转移了。

然而我们发现这个转移有环，不能直接枚举点转移，但是无论有无环都满足关系式。所以我们把关系式列出来，就得到了一个完整的 n 元一次方程组，用高斯消元求解即可。

把这个 $P_i$ 搞到手后就好办了。设 $p_i$ 表示走到节点 i 时经过奇数次权值为1的边的概率，那么经过偶数次的概率就为 $P_i-p_i$ 。如果我们不先求出经过每个点的概率，而直接把经过偶数次的概率也设出来求，那么复杂度就是 $O(8n^3\log w)≈O(240000000)$ ，再加上高斯消元的大常数肯定过不了。先求出 $P_i$ 可以把未知数个数缩为一半，时间缩为 $\frac{1}{8}$ 。

显然有 $p_0=0$ ，然后当 $i > 0$ 时， $p_i$ 满足以下关系
$p_i=\sum_{(i,v,1)\in E,v\neq n}\frac{P_v-p_v}{d_v} +\sum_{(i,v,0)\in E,v\neq n}\frac{p_v}{d_v}$
同样用高斯消元求解即可。因为 $P_i$ 无关边权，所以只用在最开始求一次。总复杂度为 $O(n^3\log w)$ 。

事后把数据输出发现， $P_i$ 居然很多都大于1。

P1: 1.800
P2: 2.400
P3: 2.000
P4: 1.800
P5: 0.800
...

~~(看起来真不错)~~
这是怎么回事呢？反正我A了，管它呢，不知道。

有可能我定义的根本就不是概率，
有可能高斯消元产生了精度问题（？上千倍的精度误差？），
有可能我两边求概率的方法都错了、凑一起就对了，
还有可能是由于我经常虔诚地膜拜 JZM ，老天放我一马…

代码

47ms rank1代码（还可优化）

#include<cstdio>//JZM yyds!!
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<stack>
#include<ctime>
#include<set>
#define ll long long
#define MAXN 105
#define uns unsigned
#define INF 1e17
#define MOD 998244353ll
#define lowbit(x) ((x)&(-(x)))
using namespace std;
inline ll read(){
	ll x=0;bool f=1;char s=getchar();
	while((s<'0'||s>'9')&&s>0){if(s=='-')f^=1;s=getchar();}
	while(s>='0'&&s<='9')x=(x<<1)+(x<<3)+s-'0',s=getchar();
	return f?x:-x;
}
const double MN=1e-13;
int n,m;
struct edge{
	int v,w;edge(){}
	edge(int V,int W){v=V,w=W;}
};
vector<edge>G[MAXN];
double p[MAXN];
double c[MAXN][MAXN],x[MAXN];
inline int Gauss(int n,int m,int&tm){
	for(int i=0;i<=m;i++)x[i]=0;
	int col,row;
	for(row=0,col=0;row<n&&col<m;row++,col++){
		int mxr=row;
		for(int i=row+1;i<n;i++)
			if(fabs(c[i][col])>fabs(c[mxr][col]))mxr=i;
		if(mxr!=row&&(++tm))
			for(int j=col;j<=m;j++)swap(c[row][j],c[mxr][j]);
		if(fabs(c[row][col])<MN){row--;continue;}
		for(int i=row+1;i<n;i++){
			if(fabs(c[i][col])<MN)continue;
			double g=c[i][col]/c[row][col];
			c[i][col]=0;
			for(int j=col+1;j<=m;j++)c[i][j]-=c[row][j]*g;
		}
	}
	for(int i=row;i<n;i++)//无解 
		if(fabs(c[i][col])>=MN)return -1;
	if(row<m)return m-row;//无数解 
	for(int i=n-1;i>=0;i--){//唯一解 
		double b=c[i][m];
		for(int j=i+1;j<m;j++)b-=c[i][j]*x[j];
		x[i]=b/c[i][i];
	}
	return 0;
}
signed main()
{
	n=read(),m=read();
	int mx=0;
	for(int i=1;i<=m;i++){
		int u=read(),v=read(),w=read();
		G[u].push_back(edge(v,w)),mx=max(mx,w);
		if(v!=u)G[v].push_back(edge(u,w));
	}
	c[0][0]=c[0][n+1]=1;
	for(int i=1;i<=n;i++){
		c[i][i]=1;
		for(uns j=0;j<G[i].size();j++){
			int v=G[i][j].v;
			if(v==n)continue;
			c[i][v]-=1.0/G[v].size();
			if(v==1)c[i][0]-=1.0/G[v].size();
		}
	}int tm=0;
	Gauss(n+1,n+1,tm);
	for(int i=0;i<=n;i++)p[i]=x[i];
	double ans=0;
	for(int D=0;D<30&&(1<<D)<=mx;D++){
		for(int i=n+1;i>=0;i--)
			for(int j=n+1;j>=0;j--)c[i][j]=0;
		for(int i=n+1;i>=0;i--)x[i]=0;
		c[0][0]=1,c[0][n+1]=0;
		for(int u=1;u<=n;u++){
			c[u][u]=1;
			for(uns i=0;i<G[u].size();i++){
				int v=G[u][i].v,w=G[u][i].w;
				if(v==n)continue;
				if((w>>D)&1)c[u][v]+=1.0/G[v].size(),
				c[u][n+1]+=p[v]/G[v].size();
				else c[u][v]-=1.0/G[v].size();
				if(v==1){
					if((w>>D)&1)c[u][0]+=1.0/G[v].size(),
					c[u][n+1]+=p[0]/G[v].size();
					else c[u][0]-=1.0/G[v].size();
				}
			}
		}
		Gauss(n+1,n+1,tm=0);
		ans+=x[n]*(1<<D);
	}
	printf("%.3f\n",ans);
	return 0;
}