[BZOJ2337][HNOI2011]XOR和路径（高斯消元解期望方程）

先来介绍一下期望方程：
考虑问题：一个$N$个点$M$条边的无向图，边有权值，到每个节点都会等概率地选择与这个点关联的下一个点走。求从$1$走到$N$的期望路径权值和。
考虑DP。设$f[u]$为从$u$到$N$的期望，$d_u$为$u$的度，$v$为与$u$有边相连的点集，$w(u,v)$表示边$(u,v)$的权值，可以发现：
$f[u]=\frac{1}{d_u}\sum(w(u,v)+f[v]),u<>N$
$f[N]=0$
但是由于存在环，所以这个DP的转移是有后效性的。所以把$f[1],f[2],...,f[N]$看成$N$个未知数，把上面的方程移项后用高斯消元求解。
回到原题。可以发现，异或不容易直接求解，所以这里按位统计：
首先枚举位$k$，设$f[i]$表示从$i$到$N$，第$k$位的期望异或和，也就是从$i$到$N$，经过的边中有奇数条边权第$k$位为1的概率。
同样可以列出：
$f[N]=0$
当$u<>N$时，对于任何一个与$u$相邻的点$v$，
如果$w(u,v)$的第$k$位为$0$，则$f[u]+=\frac{f[v]}{d_u}$，
否则$f[u]+=\frac{1-f[v]}{d_u}$。
同样移项后高斯消元。
代码：

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
inline int read() {
    int res = 0; bool bo = 0; char c;
    while (((c = getchar()) < '0' || c > '9') && c != '-');
    if (c == '-') bo = 1; else res = c - 48;
    while ((c = getchar()) >= '0' && c <= '9')
        res = (res << 3) + (res << 1) + (c - 48);
    return bo ? ~res + 1 : res;
}
const int N = 105, M = 2e4 + 5;
int n, m, ecnt, nxt[M], adj[N], go[M], val[M], cnt[N];
double a[N][N];
void add_edge(int u, int v, int w) {
    nxt[++ecnt] = adj[u]; adj[u] = ecnt; go[ecnt] = v; val[ecnt] = w;
    if (u != v) nxt[++ecnt] = adj[v], adj[v] = ecnt, go[ecnt] = u, val[ecnt] = w;
    cnt[u]++; if (u != v) cnt[v]++;
}
void Gauss() {
    int i, j, k;
    for (i = 1; i <= n; i++) {
        int u = i;
        for (j = i + 1; j <= n; j++)
            if (fabs(a[j][i]) > fabs(a[u][i])) u = j;
        if (u != i) for (j = i; j <= n + 1; j++)
            swap(a[i][j], a[u][j]);
        double tmp = a[i][i];
        for (j = i; j <= n + 1; j++) a[i][j] /= tmp;
        for (j = 1; j <= n; j++) if (j != i) {
            tmp = a[j][i];
            for (k = i; k <= n + 1; k++)
                a[j][k] -= a[i][k] * tmp;
        }
    }
}
double solve(int x) {
    memset(a, 0, sizeof(a));
    int i; for (i = 1; i < n; i++) {
        a[i][i] = 1; double d = 1.0 / cnt[i];
        for (int e = adj[i], v; e; e = nxt[e]) {
            int op = (val[e] >> x) & 1; v = go[e];
            if (op) a[i][n + 1] += d, a[i][v] += d;
            else a[i][v] -= d;
        }
    }
    return a[n][n] = 1, Gauss(), a[1][n + 1] * pow(2, x);
}
int main() {
    int i, x, y, z; n = read(); m = read();
    for (i = 1; i <= m; i++) {
        x = read(); y = read(); z = read();
        add_edge(x, y, z);
    }
    double res = 0; for (i = 0; i < 31; i++)
        res += solve(i); printf("%.3lf\n", res);
    return 0;
}