bzoj 2337: [HNOI2011]XOR和路径 （概率与期望DP+高斯消元）

2337: [HNOI2011]XOR和路径

Time Limit: 10 Sec   Memory Limit: 128 MB
Submit: 1117   Solved: 653
[ Submit][ Status][ Discuss]

Description

Input

Output

Sample Input

Sample Output

HINT

Source

Day2

[ Submit][ Status][ Discuss]

题解：概率与期望DP+高斯消元

这道题刚拿到手的时候毫无思路，然后想到之前一道跟边权异或有关的网络流，于是就想对于每个二进制位分开考虑试试。

那么因为是异或，所以如果路径上经过了奇数个1，那么2^i*路径上经过奇数个1的概率就是这一位的贡献。

于是就对于每个点维护两个概率，分别表示当前位到当前点经过奇数个1的概率，经过偶数个1的概率，然后建立方程组，用高斯消元求解即可。

注意子环只需要加一次。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#define N 20003
#define M 203
using namespace std;
int tot,point[N],v[N],c[N],w[32][N],nxt[N],n,m,du[N];
double a[M][M],b[M],ans[M];
void add(int x,int y,int z)
{
	tot++; nxt[tot]=point[x]; point[x]=tot; v[tot]=y; c[tot]=z;
	if (x==y) return;
	tot++; nxt[tot]=point[y]; point[y]=tot; v[tot]=x; c[tot]=z;
}
void guass(int n)
{
	//for (int i=1;i<=n;i++,cout<<b[i-1]<<endl)
	// for (int j=1;j<=n;j++) cout<<a[i][j]<<" ";
    for (int i=1;i<=n;i++) {
    	int num=i;
    	for (int j=i+1;j<=n;j++)
    	 if (fabs(a[j][i])>fabs(a[num][i])) num=j;
    	for (int j=1;j<=n;j++) swap(a[num][j],a[i][j]);
    	swap(b[num],b[i]);
    	for (int j=i+1;j<=n;j++)
    	 if(a[j][i]) {
    	 	double t=a[j][i]/a[i][i];
    	 	for (int k=1;k<=n;k++)
    	 	 a[j][k]-=a[i][k]*t;
    	 	b[j]-=b[i]*t;
		 }
	}	
	for (int i=n;i>=1;i--) {
	    double t=1.0*b[i];
		for (int j=i+1;j<=n;j++)
		 if (a[i][j]) t-=ans[j]*a[i][j];
		ans[i]=t/a[i][i];
	}
}
int main()
{
	freopen("xorpath.in","r",stdin);
	freopen("xorpath.out","w",stdout);
	scanf("%d%d",&n,&m);
	for (int i=1;i<=m;i++) {
		int x,y,z; scanf("%d%d%d",&x,&y,&z);
		add(x,y,i); du[x]++; 
		if (x!=y) du[y]++;
		for (int j=0;j<=30;j++)
		 w[j][i]=(z>>j)&1;
	}
	double E=0; double mi=1;
	for (int pos=0;pos<=30;pos++) {
		memset(a,0,sizeof(a));
		memset(b,0,sizeof(b));
		a[n*2][n*2]=1; b[n*2]=1;
		a[n*2-1][n*2-1]=1; b[n*2-1]=0;
		for (int i=1;i<n;i++) {
         double t=1.0/du[i];
         a[i*2][i*2]=1;
         a[i*2-1][i*2-1]=1;
		 for (int j=point[i];j;j=nxt[j]) 
		  if (w[pos][c[j]])  {
		  	 int x=v[j]*2-1;
		     a[i*2][x]-=t;
	         a[i*2-1][x+1]-=t;	 
		  }
		  else {
		  	int x=v[j]*2-1;
		  	a[i*2-1][x]-=t;
		  	a[i*2][x+1]-=t;
		  }
	    }
	    guass(n*2);
	    E+=ans[1]*mi;
		mi*=2.0;   
	}
	printf("%.3lf\n",E);
}