说在前面
很久没有写过高斯消元了,看着别人的代码yy了一会。一遍敲出来感觉成就感++;
而且gauss这个单词感觉特别帅有没有!!
题目
BZOJ2337传送门
(突然发现一张图就把所有题目信息包含完了,不用手打还有点不习惯…)
解法
因为原题是要求异或期望,位与位之间没有影响,所以拆开处理。
( 以下的分析均是针对某一二进制位上的值 )
定义f[u]表示从u点走到N点时,该位为1的期望
对于点u和u能到的点v来说(deg[u]表示u的度数):
如果边上的权值为0,那么f[u] += f[v] / deg[u]
如果边上的权值为1,那么f[u] += ( 1 - f[v] ) / deg[u]
边界:f[N] = 0 【显然,因为N就是终点了】
这张图是无向的,不能直接靠拓扑关系求解,由是采用高斯消元
实现的时候为了方便,并没有采用小数的形式,而是等式左右两边同时乘上了deg[u]。
需要注意的地方有两个
一是在方程矩阵里对于每个F[u][u],它的初值是deg[u],因此要么在枚举完v之后F[u][u]+=deg[u],要么是在枚举值前赋值。
二是边界f[N]为0,无论F[N][N]是什么值,最后解出来N的答案都是0。在我的程序中,F[N][N]必须赋值,不然会被gauss判定无解(虽然这题根本就不用判断)。
自带大常数的代码
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std ;
int N , M , deg[105] , head[105] , tp ;
double F[105][105] , ans ;
struct Path{
int pre , to , val ;
}p[20005] ;
void In( int t1 , int t2 , int t3 ){
p[++tp].pre = head[t1] ;
p[ head[t1] = tp ].to = t2 ;
p[tp].val = t3 ;
}
template <typename T>
T abs( T x ){
return x > 0 ? x : -x ;
}
void gauss_Ele(){
for( int i = 1 , s ; i <= N ; i ++ ){
s = i ;
for( int j = i + 1 ; j <= N ; j ++ )
if( abs( F[s][i] ) < abs ( F[j][i] ) ) s = j ;
if( s != i )
for( int j = i ; j <= N + 1 ; j ++ )
swap( F[i][j] , F[s][j] ) ;
for( int j = i + 1 ; j <= N ; j ++ ){
double rate = F[j][i] / F[i][i] ;
for( int k = i ; k <= N + 1 ; k ++ )
F[j][k] -= rate * F[i][k] ;
}
}
for( int i = N ; i ; i -- ){
double delta = 0 ;
for( int j = i + 1 ; j <= N ; j ++ )
delta += F[i][j] * F[j][N+1] ;
F[i][N+1] = ( F[i][N+1] - delta ) / F[i][i] ;
}
}
int main(){
scanf( "%d%d" , &N , &M ) ;
for( int i = 1 , u , v , x ; i <= M ; i ++ ){
scanf( "%d%d%d" , &u , &v , &x ) ;
In( u , v , x ) ; deg[u] ++ ;
if( u == v ) continue ;
In( v , u , x ) ; deg[v] ++ ;
}
for( int w = 0 ; w <= 30 ; w ++ ){
for( int i = 1 ; i <= N ; i ++ )
for( int j = 1 ; j <= N + 1 ; j ++ )
F[i][j] = 0 ;
for( int u = 1 ; u < N ; u ++ ){
F[u][u] = deg[u] ;
for( int i = head[u] ; i ; i = p[i].pre ){
int v = p[i].to , val = p[i].val ;
if( val & ( 1 << w ) ) F[u][v] += 1 , F[u][N+1] += 1 ;
else F[u][v] -= 1 ;
}
}
F[N][N] = 1234321.0;
gauss_Ele() ;
ans += F[1][N+1] * ( 1 << w ) ;
}
printf( "%.3f" , ans ) ;
}