高斯消元，就是 $O(n^3)$ 解方程组
bzoj3270 博物馆
一个无向图，两个人一个在 A 一个在 B，这两个人开始随机走，求这两个人在每个点相遇的概率（在边上不会相遇），每个点有一个自环，每次有 $P_i$ 的概率走自环，剩下 $1 - P_i$ 的概率等概率选一个相邻点走过去
$n leq 20$
sol：
我是不是...开错题了？
原题等价于求这两个人在每个点相遇的期望步数
设 $f_{(i,j)}$ 为第一个人在 $i$ ，第二个人在 $j$ 的期望步数，枚举这两个人走的情况转移
注意：
1.$i==j$ 的点不转移，因为相遇了就停了
2.因为一开始就在起点，所以起点的概率要 $+1$ ，概率 $+1$ 相当于方程解出来右边是 $-1$
#include<bits/stdc++.h> #define LL long long using namespace std; inline int read() { int x = 0,f = 1;char ch = getchar(); for(;!isdigit(ch);ch = getchar())if(ch == '-') f = -f; for(;isdigit(ch);ch = getchar())x = 10 * x + ch - '0'; return x * f; } const int maxn = 410; int n,m,A,B; int d[maxn];double pr[maxn],a[maxn][maxn]; vector<int> G[maxn]; int mem[25][25],dfn; inline int f(int x,int y){return mem[x][y] ? mem[x][y] : (mem[x][y] = ++dfn);} void gauss() { int now=1,tot = n * n; for(int i=1;i<=tot;i++) { int j; for(j=now;!a[j][now]&&j<=tot;j++); for(int k=1;k<=tot+1;k++)swap(a[now][k],a[j][k]); for(int j=1;j<=tot;j++) if(j!=now) { double t=a[j][now]/a[now][now]; for(int k=1;k<=tot+1;k++) a[j][k]-=t*a[now][k]; } now++; } } int main() { n = read(),m = read(),A = read(),B = read(); for(int i=1;i<=n;i++)G[i].push_back(i);for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)f(i,j); a[f(A,B)][n * n + 1] = -1; for(int i=1;i<=m;i++) { int u = read(),v = read(); d[u]++;G[u].push_back(v); d[v]++;G[v].push_back(u); } for(int i=1;i<=n;i++)cin >> pr[i]; for(int x=1;x<=n;x++) for(int y=1;y<=n;y++) { a[f(x,y)][f(x,y)]--; for(int i=0;i<G[x].size();i++) for(int j=0;j<G[y].size();j++) { int tx = G[x][i],ty = G[y][j]; int p = f(x,y),np = f(tx,ty); if(tx == ty)continue; if(tx == x && ty == y)a[p][np] += pr[x] * pr[y]; else if(tx == x && ty != y)a[p][np] += pr[x] * (1 - pr[ty]) / d[ty]; else if(tx != x && ty == y)a[p][np] += (1 - pr[tx]) / d[tx] * pr[y]; else if(tx != x && ty != y)a[p][np] += (1 - pr[tx]) / d[tx] * (1 - pr[ty]) / d[ty]; } } gauss(); for(int i=1;i<=n;i++) { int fd = f(i,i); printf("%.6lf ",a[fd][n * n + 1] / a[fd][fd]); } }
bzoj2337 XOR 和路径
一个无向图，从 $1$ 出发每次随机选一条出边走，走到 $n$ 停止，求经过每条边权异或和的期望
$n leq 100$
sol：
XOR 的期望...XOR 看起来是整数，期望看起来是小数，没法直接算（两个小数没法 XOR ）
可以拆开每一位，算出每一位最终结果的期望，再按位值加起来就可以了
#include<bits/stdc++.h> #define LL long long using namespace std; inline int read() { int x = 0,f = 1;char ch = getchar(); for(;!isdigit(ch);ch = getchar())if(ch == '-') f = -f; for(;isdigit(ch);ch = getchar())x = 10 * x + ch - '0'; return x * f; } const int maxn = 310; int n,m; int d[maxn];double a[maxn][maxn]; int first[maxn],to[maxn * maxn],nx[maxn * maxn],val[maxn * maxn],cnt; inline void add(int u,int v,int w) { to[++cnt] = v; nx[cnt] = first[u]; first[u] = cnt; val[cnt] = w; d[u]++; } void gauss() { int now=1,tot = n; for(int i=1;i<=tot;i++) { int j; for(j=now;!a[j][now]&&j<=tot;j++); for(int k=1;k<=tot+1;k++)swap(a[now][k],a[j][k]); for(int j=1;j<=tot;j++) if(j!=now) { double t=a[j][now]/a[now][now]; for(int k=1;k<=tot+1;k++) a[j][k]-=t*a[now][k]; } now++; } } int main() { n = read(),m = read(); for(int i=1;i<=m;i++) { int u = read(),v = read(),w = read(); add(u,v,w); if(u != v)add(v,u,w); }double ans = 0; for(int T=0;T<=30;T++) { memset(a,0,sizeof(a)); for(int i=1;i<n;i++)a[i][i] = 1.0; for(int x=1;x<n;x++) { for(int i=first[x];i;i=nx[i]) { if(val[i] & (1 << T))a[x][to[i]] += (1.0 / d[x]),a[x][n + 1] += (1.0 / d[x]); else a[x][to[i]] -= (1.0 / d[x]); } }a[n][n] = 1; gauss(); ans += (1 << T) * (a[1][n + 1] / a[1][1]); } printf("%.3lf ",ans); }
bzoj4820 硬币游戏
给出 $n$ 个长度均为 $m$ 的不同 $01$ 串,随机生成一个无限长的 $01$ 串,对 $n$ 个 $01$ 串中的每个,求出它最先在随机串中出现的概率.
$n,m leq 300$
sol：
如果直接上 AC 自动机的话，会有 $O(n^2)$ 个方程，所以要考虑优化方程个数
可以把“没有任何给定 $01$ 串出现”作为一个状态 $N$，如果一个 $01$ 串 $A$，然后某次转移从 $N$ 转移到了 $N + A$，那么 $A$ 就有可能获胜
至于“有可能”是因为有可能 $A$ 还没有出现完，别的就已经出现了
这时候就要考虑处理出这种情况发生的概率
容易知道这种情况发生是因为 $A$ 的前缀等于 $B$ 的后缀
之后我们可以对于每一对 $(i,j)$ 做一次 kmp ，求出所有 $S_i$ 和 $S_j$ 的相同前后缀
然后概率就是 $P_{(i,j)} = sumlimits_k 2^{k-m}$
设走到 $N$ 的期望步数是 $x_0$，走到第 $i$ 个人的终止节点的期望步数是 $x_i$，则可以对每一个人列一个方程
$2^{-m} space x_0 space = space x_i+sumlimits_{j=1}^n P_{(i,j)} space x_j$
最后还有一个方程 $sumlimits_{i=1}^n x_i = 1$ （因为最后总是要有人赢的）
这样就是 $n+1$ 个方程，高斯消元即可
#include <bits/stdc++.h> #define LL long long #define rep(i, s, t) for (register int i = (s), i##end = (t); i <= i##end; ++i) #define dwn(i, s, t) for (register int i = (s), i##end = (t); i >= i##end; --i) using namespace std; inline int read() { int x = 0, f = 1;char ch; for (ch = getchar(); !isdigit(ch); ch = getchar()) if (ch == '-') f = -f; for (; isdigit(ch); ch = getchar()) x = 10 * x + ch - '0'; return x * f; } const int maxn = 310; int n, m; int fail[maxn << 1]; double f[maxn][maxn]; char ch[maxn][maxn], a[maxn << 1]; void kmp() { int j = 0; int len = strlen(a + 1); rep(i, 2, len) { while(j && a[j + 1] != a[i]) j = fail[j]; if(a[j + 1] == a[i]) j++; fail[i] = j; } } void gauss() { int now = 1; for(int i=1;i<=n+1;i++) { int j; for(j=now;!f[j][now] && j<=n+1;j++); swap(f[now], f[j]); for(int j=1;j<=n+1;j++) { if(j != now) { double t = f[j][now] / f[now][now]; for(int k=1;k<=n+2;k++) f[j][k] -= t * f[now][k]; } } now++; } } int main() { n = read(), m = read(); rep(i, 1, n) scanf("%s", ch[i] + 1); rep(i, 1, n) rep(j, 1, n) { rep(k, 1, m) a[k] = ch[i][k], a[k + m] = ch[j][k]; kmp(); for(int k = fail[m << 1]; k; k = fail[k]) if(k < m) f[i][j] += pow(0.5, m - k); } rep(i, 1, n) f[i][i] += 1.0, f[i][n+1] -= 1.0; f[n + 1][n + 2] = 1.0; rep(i, 1, n) f[n + 1][i] = 1.0; gauss(); rep(i, 1, n) printf("%.10f ", f[i][n + 2] / f[i][i]); } //2019.3.11 终于填上了...
UPD：真正的高斯消元
大概就是从左到右考虑每一个变量，找到系数最大的一行当主元，换到第 $i$ 行
然后这一行从第 $i$ 个开始每个除以 $mat(i,i)$
之后从第 $i+1$ 行到第 $n$ 行，第 $i$ 个到第 $n+1$ 个开始每个减去 $mat(j,i) div mat(i,i) imes mat(i,k)$
这样能消成一个下三角矩阵，把主对角线乘起来就是行列式
求答案的话搞一个 $ans$ 数组，一开始 $ans[i] = mat(i,n+1)$
之后从 $n-1$ 到 $1$ 枚举 $i$，$ans[i] = ans[i] - sumlimits_{j=i+1}^n ans[j] imes mat(i,j)$
最后 $ans[i]$ 就是 $i$ 的值
（然而大概是要用来求行列式吧，之前的 Gauss-Jordan 解方程似乎更好，然而 Gauss-Jordan 行列式简直吔屎
#include <bits/stdc++.h> #define LL long long using namespace std; #define rep(i, s, t) for (register int i = (s), i##end = (t); i <= i##end; ++i) #define dwn(i, s, t) for (register int i = (s), i##end = (t); i >= i##end; --i) inline int read() { int x = 0, f = 1; char ch = getchar(); for (; !isdigit(ch); ch = getchar()) if (ch == '-') f = -f; for (; isdigit(ch); ch = getchar()) x = 10 * x + ch - '0'; return x * f; } int n; double a[110][110], ans[110]; int Gauss() { rep(i, 1, n) { int now = i; rep(k, i+1, n) if(fabs(a[k][i]) > fabs(a[now][i])) now = k; if(a[now][i] == 0) return 1; if(now != i) rep(k, i, n+1) swap(a[i][k], a[now][k]); double t = a[i][i]; rep(j, i, n+1) a[i][j] /= t; rep(j, i+1, n) { double t = a[j][i] / a[i][i]; rep(k, i, n+1) a[j][k] -= t * a[i][k]; } } return 0; } int main() { n = read(); rep(i, 1, n) rep(j, 1, n+1) cin >> a[i][j]; if(Gauss()) cout << "No Solution" << endl; else { rep(i, 1, n) ans[i] = a[i][n+1]; dwn(i, n-1, 1) rep(j, i+1, n) ans[i] -= a[i][j] * ans[j]; rep(i, 1, n) cout << fixed << showpoint << setprecision(2) << ans[i] << endl; } }
UPD：新板子
int Gauss() { int i, j; for(i = 1, j = 1; i <= n && j <= m; ++j) { int pos = 0; rep(k, i, n) if(a[k][j]) { pos = k; break; } if(!pos) continue; if(pos != i) rep(k, 1, m) swap(a[i][k], a[pos][k]); int _inv = ksm(a[i][j], mod - 2); rep(k, i+1, n) if(a[k][j]) { int tmp = 1LL * a[k][j] * _inv % mod; rep(l, j, m) (a[k][l] -= (1LL * tmp * a[i][l] % mod)) %= mod; } i++; } return i - 1; }