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

119 阅读 0 评论 79 点赞

我是靠谱客的博主激动乌龟，这篇文章主要介绍[HNOI2011]XOR和路径——高斯消元、随机游走模型[HNOI2011]XOR和路径思路代码，现在分享给大家，希望可以做个参考。

[HNOI2011]XOR和路径

思路

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

先设 $P_i$ 表示从1走到 n 中经过点 i 的概率，那么应该有转移关系式
$P_i=sum_{(i,v)in E,vneq 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^3log w)≈O(240000000)$ ，再加上高斯消元的大常数肯定过不了。先求出 $P_i$ 可以把未知数个数缩为一半，时间缩为 $1 8 frac{1}{8}$ 。

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

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

复制代码

1
2
3
4
5
6
7
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("%.3fn",ans);
return 0;
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
#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("%.3fn",ans);
return 0;
}