高斯消元
高斯消元
高斯消元法通常用于求解如下的 \(n\) 元线性方程组：
\[\begin{cases}
a_{1, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{1, n} x_n = b_1
\\
a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n = b_2
\\
\cdots
\\
a_{n, 1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n = b_n
\end{cases}
\]
线性方程组与矩阵的联系
对于一个线性方程组：
\[\begin{cases}
a_{1, 1} x_1 + a_{1, 2} x_2 + \cdots + a_{1, n} x_n = b_1 \\
a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n = b_2 \\
\vdots \\
a_{m, 1} x_1 + a_{m, 2} x_2 + \cdots + a_{m, n} x_n = b_n
\end{cases}
\]
就可以表示为下面的矩阵：
\[\left[
\begin{array}{cccc|c}
a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} & b_1 \\
a_{2, 1} & a_{2, 2} & \cdots & a_{2, n} & b_2 \\
\vdots & \vdots & \ddots & \vdots & \vdots \\
a_{m, 1} & a_{m, 2} & \cdots & a_{m, n} & b_m \\
\end{array}
\right]
\]
矩阵里每个对应的系数整理出来就是这个矩阵。
整个矩阵称为增广矩阵，而没有右边系数的矩阵称为系数矩阵。
因为解均形如 \(x_i = a_i\) ，考虑把一个线性方程组的 \(n\) 个元的解写在一个列向量（仅含一列的矩阵）中，如：
\[X = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix}
\]
引入一个线性组合的概念：对于 \(n\) 个同型向量 \(\vec{v_i}\) 与标量 \(c_i\) ，向量 \(\vec{x} = \sum^{n}_{i = 1} c_i \vec{v_i}\) 称为向量 \(\vec{v_i}\) 以 \(c_i\) 为权的线性组合。
于是可以得到向量方程：
\[A \vec{x} = \begin{bmatrix} \vec{a_1}, \vec{a_2}, \cdots, \vec{a_n} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} = x_1 \begin{bmatrix} a_{1, 1} \\ a_{1, 2} \\ \vdots \\ a_{1, n} \end{bmatrix} + x_2 \begin{bmatrix} a_{2, 1} \\ a_{2, 2} \\ \vdots \\ a_{2, n} \end{bmatrix} + \cdots + x_n \begin{bmatrix} a_{n, 1} \\ a_{n, 2} \\ \vdots \\ a_{n, n} \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{bmatrix}
\]
即 \(AX = B\) 。
高斯消元解线性方程组
考虑求解如下方程组：
\[\begin{cases}
x_1 + 3x_2 + 4x_3 = 5
\\
x_1 + 4x_2 + 7x_3 = 3
\\
9x_1 + 3x_2 + 2x_3 = 2
\end{cases}
\]
把系数写到矩阵里：
\[\begin{bmatrix}
1 & 3 & 4 & 5 \\
1 & 4 & 7 & 3 \\
9 & 3 & 2 & 2 \\
\end{bmatrix}
\]
把未知项的系数单独拿出来考虑：
\[\begin{bmatrix}
1 & 3 & 4 \\
1 & 4 & 7 \\
9 & 3 & 2 \\
\end{bmatrix}
\]
如果能够将其化为：
\[\begin{bmatrix}
? & ? & ? \\
0 & ? & ? \\
0 & 0 & ? \\
\end{bmatrix}
\]
就可以从最后一行逆推上去。设第一行第一列为主元，利用加减消元法配系数消掉后两行的第一列：
\[\begin{bmatrix}
1 & 3 & 4 \\
0 & 1 & 3 \\
0 & -24 & -34 \\
\end{bmatrix}
\]
再以第二行第二列为主元，加减消元得到：
\[\begin{bmatrix}
1 & 3 & 4 \\
0 & 1 & 3 \\
0 & 0 & 38 \\
\end{bmatrix}
\]
即得目标状态，实现时把 \(b\) 的那一列一起带入消元即可。
每次选取系数最大的为主元可以有效降低精度误差。
本质就是用初等变换消为一个下三角矩阵。
#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-9;
const int N = 1e2 + 7;
double g[N][N], ans[N];
int n;
inline bool Gauss() {
for (int i = 1; i <= n; ++i) {
int mxpos = i;
for (int j = i + 1; j <= n; ++j)
if (fabs(g[j][i]) > fabs(g[mxpos][i]))
mxpos = j;
if (fabs(g[mxpos][i]) < eps)
return false;
if (mxpos != i)
swap(g[i], g[mxpos]);
for (int j = i + 1; j <= n; ++j) {
double div = g[j][i] / g[i][i];
for (int k = i; k <= n + 1; ++k)
g[j][k] -= g[i][k] * div;
}
}
for (int i = n; i; --i) {
ans[i] = g[i][n + 1];
for (int j = i + 1; j <= n; ++j)
ans[i] -= g[i][j] * ans[j];
ans[i] /= g[i][i];
}
return true;
}
int main() {
scanf("%d", &n);
for (int i = 1; i <= n; ++i)
for (int j = 1; j <= n + 1; ++j)
scanf("%lf", g[i] + j);
if (!Gauss())
return puts("No Solution"), 0;
for (int i = 1; i <= n; ++i)
printf("%.2lf\n", ans[i]);
return 0;
}
约旦消元法解线性方程组
其为高斯消元的进阶版，能避免回带来求出答案，也就是要把矩阵变为：
\[\begin{bmatrix}
? & 0 & 0 & ? \\
0 & ? & 0 & ? \\
0 & 0 & ? & ? \\
\end{bmatrix}
\]
首先，我们依照朴素的高斯消元不难得到：
\[\begin{bmatrix}
? & ? & ? & ? \\
0 & ? & ? & ? \\
0 & 0 & ? & ? \\
\end{bmatrix}
\]
观察上下两个矩阵，不难得到，我们在消元时不仅仅消去主元所在式子下方的式子，而对于上方的式子也应当予以处理。
注意一下特殊情况：
-
无解：当左边系数为 \(0\) 而右边系数不为 \(0\) 时即为无解。
-
多解：当左边系数为 \(0\) 而右边系数也为 \(0\) 时即为多解。
inline bool Gauss() {
for (int i = 1; i <= n; ++i) {
int mxpos = i;
for (int j = i + 1; j <= n; ++j)
if (fabs(g[j][i]) > fabs(g[mxpos][i]))
mxpos = j;
if (i != mxpos)
swap(g[i], g[mxpos]);
if (fabs(g[i][i]) < eps)
return false;
for (int j = 1; j <= n; ++j)
if (j != i) {
double div = g[j][i] / g[i][i];
for (int k = i + 1; k <= n + 1; ++k)
g[j][k] -= g[i][k] * div;
}
}
for (int i = 1; i <= n; ++i)
ans[i] = g[i][n + 1] / g[i][i];
return true;
}
约旦消元法解异或方程组
异或方程组是指形如：
\[\begin{cases}
a_{1, 1} x_1 \oplus a_{2, 2} x_2 \oplus \cdots \oplus a_{1, n} x_n = b_1
\\
a_{2, 1} x_1 \oplus a_{2, 2} x_2 \oplus \cdots \oplus a_{2, n} x_n = b_2
\\
\cdots
\\
a_{n, 1} x_1 \oplus a_{n, 2} x_2 \oplus \cdots \oplus a_{n, n} x_n = b_n
\end{cases}
\]
的方程组，其中 \(a_{i, j}, b_i\) 的值为 \(0\) 或 \(1\) 。
由于 \(\oplus\) 符合交换律与结合律，故可以按照高斯消元法逐步消元求解。值得注意的是，我们在消元的时候应使用异或消元而非加减消元，且不需要进行乘除改变系数（因为系数均为 \(0\) 和 \(1\) ）。
注意到异或方程组的增广矩阵是 \(01\) 矩阵，可以用 bitset 优化到 \(O(\dfrac{n^2 m}{\omega})\) 。
inline bool Gauss() {
for (int i = 1; i <= n; ++i) {
int mxpos = i;
while (mxpos <= n && !g[mxpos].test(i))
++mxpos;
if (mxpos > n)
return false;
if (i != mxpos)
swap(g[i], g[mxpos]);
for (int j = 1; j <= m; ++j)
if (j != i && g[j].test(i))
g[j] ^= g[i];
}
return true;
}
矩阵求逆
对于矩阵 \(A\) ，若存在矩阵 \(A^{-1}\) 使得 \(A \times A^{-1} = A^{-1} \times A = I\) ，则称矩阵 \(A\) 可逆，\(A^{-1}\) 为其逆矩阵。
给出 \(n\) 阶方阵 \(A\) ，求解其逆矩阵的方法如下：
- 构造一个 \(n \times 2n\) 的矩阵 \((A, I_n)\)
- 用高斯消元法将其化简为最简形 \((I_n, A^{-1})\) ，即可得到 \(A\) 的逆矩阵 \(A^{-1}\) 。
- 如果最终最简形的左半部分不是单位矩阵 \(I_n\) ，则矩阵 \(A\) 不可逆。
inline bool Gauss() {
for (int i = 1; i <= n; ++i)
g[i][i + n] = 1;
for (int i = 1; i <= n; ++i) {
int mxpos = i;
while (mxpos <= n && !g[mxpos][i])
++mxpos;
if (mxpos > n)
return false;
if (i != mxpos)
swap(g[i], g[mxpos]);
int inv = mi(g[i][i], Mod - 2);
for (int j = 1; j <= n; ++j) {
if (j == i)
continue;
int div = 1ll * g[j][i] * inv % Mod;
for (int k = i; k <= n * 2; ++k)
g[j][k] = dec(g[j][k], 1ll * g[i][k] * div % Mod);
}
}
for (int i = 1; i <= n; ++i) {
int inv = mi(g[i][i], Mod - 2);
for (int j = 1; j <= 2 * n; ++j)
g[i][j] = 1ll * g[i][j] * inv % Mod;
}
return true;
}
行列式求值
考虑一个情况，当一个矩阵任意一个位置出现 \(0\) ，其对行列式的影响非常大（直接没有贡献了）。
我们利用上面的一些性质显然是可以让矩阵不断变化出现 \(0\) 的。考虑将矩阵一行（列）消成只有最后一个元素非 \(0\) 该怎么做。也就是说：
\[\begin{bmatrix}
a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\
a_{2, 1} & a_{2, 2} & \cdots & a_{2, n} \\
\vdots & \vdots & \ddots & \vdots \\
a_{n, 1} & a_{n, 2} & \cdots & a_{n, n} \\
\end{bmatrix}
\Rightarrow
\begin{bmatrix}
a_{1, 1} & a_{1, 2} & \cdots & a_{1, n} \\
a_{2, 1} & a_{2, 2} & \cdots & a_{2, n} \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & a_{n, n} \\
\end{bmatrix}
\]
对于 \(1 \sim n - 1\) 列中的第 \(i\) 列，我们只要让第 \(i\) 列整列加上第 \(n\) 列的 \(-\dfrac{a_{n, i}}{a_{n, n}}\) 倍即可在 \(\det A\) 不变的情况下消掉点 \(n - 1\) 个元素
运用行列式展开发现 \(\det A = a_{n, n} \times A_{n, n}\)
去掉 \(a_{n, n}\) 所在的行和列，继续消元，则发现 \(\det A = a_{n - 1, n - 1} \times A_{n - 1, n - 1}\)
以此类推，如果我们一直对当前行列式消元，取最末（右下角）位的指和其余子式，余子式作为新行列式重新做。
一直递归下去做，我们发现最后我们得到一个下三角行列式。
我们就会发现这样一个矩阵的行列式是其对角线所有元素的乘积，即 \(\prod_{i = 1}^n a_{i, i}\) 。
inline int Gauss(int n) {
int res = 1;
for (int i = 1; i <= n; ++i)
for (int j = i + 1; j <= n; ++j) {
while (g[i][i]) {
int div = g[j][i] / g[i][i];
for (int k = i; k <= n; ++k)
g[j][k] = dec(g[j][k], 1ll * g[i][k] * div % Mod);
swap(g[i], g[j]), res = Mod - res;
}
swap(g[i], g[j]), res = Mod - res;
}
for (int i = 1; i <= n; ++i)
res = 1ll * res * g[i][i] % Mod;
return res;
}
band-matrix
长这个样子：
空白部分都为 \(0\) ，橙色部分可以为任何数，这样中间就形成了一个宽度大约为 \(d\) 的带。
可以发现任意一个 \(i\) 满足从 \((i, i)\) 向右或向下拓展都有不超过 \(d - 1\) 个非零数字，即很多位置根本不需要消。
具体地，假设现在要消第 \(i\) 列，那么从第 \(i\) 行开始往下枚举 \(d - 1\) 行，每行往右消 \(d\) 个数字即可，最后仍能得到一个上三角矩阵。
与普通高斯消元有点不一样的地方在于当主元为 \(0\) 的时候的处理方法。在 band-matrix 中，若直接交换行会破坏 band-matrix 。注意到每次交换完后交换的行右边最多多出 \(d\) 个数，于是每次往右消元 \(2d\) 个数即可。
时间复杂度 \(O(nd^2)\) 。
inline bool Gauss(int n, int band) {
for (int i = 1; i <= n; ++i) {
int mxpos = i;
while (mxpos <= min(i + band, n) && fabs(g[mxpos][i]) < eps)
++mxpos;
if (mxpos > min(i + band, n))
return false;
if (i != mxpos)
swap(g[i], g[mxpos]);
for (int j = i + 1; j <= min(i + band, n); ++j) {
double div = g[j][i] / g[i][i];
for (int k = i; k <= min(i + 2 * band, n); ++k)
g[j][k] -= div * g[i][k];
g[j][n + 1] -= div * g[i][n + 1];
}
}
for (int i = n; i; --i) {
ans[i] = g[i][n + 1];
for (int j = i + 1; j <= min(i + 2 * band, n); ++j)
ans[i] -= g[i][j] * ans[j];
ans[i] /= g[i][i];
}
return true;
}
还有一种解决主元为 \(0\) 的方法，普通高消是交换行，这里只要交换列就可以保持 band-matrix 的性质了，时间复杂度也是 \(O(nd^2)\) 。
inline void Gauss(int n, int band) {
iota(id + 1, id + 1 + n, 1);
for (int i = 1; i <= n; ++i) {
if (!g[i][i]) {
for (int j = i + 1; j <= min(n, i + band); ++j)
if (g[i][j]) {
for (int k = 1; k <= n; ++k)
swap(g[k][i], g[k][j]);
swap(id[i], id[j]);
break;
}
}
int Inv = mi(g[i][i], Mod - 2);
for (int j = i + 1; j <= n; ++j) {
int div = 1ll * g[j][i] * Inv % Mod;
for (int k = i; k <= min(n, i + band); ++k)
g[j][k] = dec(g[j][k], 1ll * g[i][k] * div % Mod);
g[j][n + 1] = dec(g[j][n + 1], 1ll * g[i][n + 1] * div % Mod);
}
}
for (int i = n; i; --i) {
ans[id[i]] = g[i][n + 1];
for (int j = i + 1; j <= min(n, i + band); ++j)
ans[id[i]] = dec(ans[i], 1ll * ans[id[j]] * g[i][j] % Mod);
ans[id[i]] = 1ll * ans[id[i]] * mi(g[i][i], Mod - 2) % Mod;
}
}
应用
\(n\) 行 \(m\) 列的矩阵，现在在 \((x,y)\)，每次等概率向左、右、下走或原地不动，但不能走出去，求走到最后一行期望的步数。
\(n, m \leq 10^3\)
记 \(f_{i, j}\) 表示机器人在 \((i, j)\) 时走到最后一行的期望步数，则：
\(m = 1\) 时有（省略第二维）：
\[f_i = 1 + \dfrac{1}{2} (f_{i + 1} + f_i)
\]
即：
\[f_i = f_{i + 1} + 2
\]
\(m > 1\) 时有：
\[f_{i, j} = 1 +
\begin{cases}
\dfrac{1}{3} (f_{i, j} + f_{i + 1, j} + f_{i, j + 1}), & j = 1 \\
\dfrac{1}{4} (f_{i, j} + f_{i + 1, j} + f_{i, j - 1} + f_{i, j + 1}), & 1 < j < m \\
\dfrac{1}{3} (f_{i, j} + f_{i + 1, j} + f_{i, j - 1}), & j = m
\end{cases}
\]
注意到这是一个 \(d = 2\) 的 band-matrix ，直接使用 band-matrix 消元即可，时间复杂度 \(O(nmd^2)\) 。
#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-12;
const int N = 1e3 + 7;
double g[N][N], ans[N][N];
int n, m, x, y;
inline void Gauss(int n, int band, double *ans) {
for (int i = 1; i <= n; ++i) {
int mxpos = i;
while (mxpos <= min(i + band, n) && fabs(g[mxpos][i]) < eps)
++mxpos;
if (mxpos > min(i + band, n))
continue;
if (i != mxpos)
swap(g[i], g[mxpos]);
for (int j = i + 1; j <= min(i + band, n); ++j) {
double div = g[j][i] / g[i][i];
for (int k = i; k <= min(i + 2 * band, n); ++k)
g[j][k] -= div * g[i][k];
g[j][n + 1] -= div * g[i][n + 1];
}
}
for (int i = n; i; --i) {
ans[i] = g[i][n + 1];
for (int j = i + 1; j <= min(i + 2 * band, n); ++j)
ans[i] -= g[i][j] * ans[j];
ans[i] /= g[i][i];
}
}
signed main() {
scanf("%d%d%d%d", &n, &m, &x, &y);
if (m == 1)
return printf("%.10lf\n", 2.0 * (n - x)), 0;
for (int i = n - 1; i; --i) {
g[1][1] = 2, g[1][2] = -1, g[1][m + 1] = 3 + ans[i + 1][1];
for (int j = 2; j < m; ++j)
g[j][j] = 3, g[j][j - 1] = g[j][j + 1] = -1, g[j][m + 1] = 4 + ans[i + 1][j];
g[m][m] = 2, g[m][m - 1] = -1, g[m][m + 1] = 3 + ans[i + 1][m];
Gauss(m, 1, ans[i]);
}
printf("%.10lf\n", ans[x][y]);
return 0;
}
平面直角坐标系上有一个点，一开始在 \((0, 0)\) ，每秒钟这个点都会随机移动，如果它在 \((x, y)\) ，下一秒：
- 在 \((x - 1, y)\) 的概率是 \(p_1\)
- 在 \((x, y - 1)\) 的概率是 \(p_2\)
- 在 \((x + 1, y)\) 的概率是 \(p_3\)
- 在 \((x, y + 1)\) 的概率是 \(p_4\)
保证 \(p_1 + p_2 + p_3 + p_4 = 1\) ，求该点移动至距离原点距离为大于 \(R\) 的点的期望步数
\(0 \leq R \leq 50\)
把所有满足 \(i^2 + j^2 \leq R^2\) 的点依次编号，显然有 \(O(R^2)\) 个点。
设 \(f_{id(i, j)}\) 表示 \((i, j)\) 走出圆的期望步数，\((i, j)\) 只要转移到 \((i, j - 1), (i - 1, j), (i, j + 1), (i + 1, j)\) 。因为是依次编号，所以建出来的矩阵带宽 \(\leq 2R + 1\) 。套用 band-matrix ，可以做到 \(O(R^2 d^2) = O(R^4)\) 。
#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int Mod = 1e9 + 7;
const int R = 5e1 + 7, SIZE = 8e3 + 7;
vector<pair<int, int> > Pos;
int g[SIZE][SIZE];
int id[R << 1][R << 1];
int ans[SIZE];
int r, p1, p2, p3, p4;
template <class T = int>
inline T read() {
char c = getchar();
bool sign = c == '-';
while (c < '0' || c > '9')
c = getchar(), sign |= c == '-';
T x = 0;
while ('0' <= c && c <= '9')
x = (x << 1) + (x << 3) + (c & 15), c = getchar();
return sign ? (~x + 1) : x;
}
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline int mi(int g, int b) {
int res = 1;
for (; b; b >>= 1, g = 1ll * g * g % Mod)
if (b & 1)
res = 1ll * res * g % Mod;
return res;
}
inline int getid(int x, int y) {
return id[x + R][y + R];
}
inline void Gauss(int n, int band) {
for (int i = 1; i <= n; ++i) {
int mxpos = i;
while (mxpos <= min(i + band, n) && !g[mxpos][i])
++mxpos;
if (mxpos > min(i + band, n))
continue;
if (i != mxpos)
swap(g[i], g[mxpos]);
int inv = mi(g[i][i], Mod - 2);
for (int j = i + 1; j <= min(i + band, n); ++j) {
int div = 1ll * g[j][i] * inv % Mod;
for (int k = i; k <= min(i + 2 * band, n); ++k)
g[j][k] = dec(g[j][k], 1ll * div * g[i][k] % Mod);
g[j][n + 1] = dec(g[j][n + 1], 1ll * div * g[i][n + 1] % Mod);
}
}
for (int i = n; i; --i) {
ans[i] = g[i][n + 1];
for (int j = i + 1; j <= min(i + 2 * band, n); ++j)
ans[i] = dec(ans[i], 1ll * g[i][j] * ans[j] % Mod);
ans[i] = 1ll * ans[i] * mi(g[i][i], Mod - 2) % Mod;
}
}
signed main() {
r = read(), p1 = read(), p2 = read(), p3 = read(), p4 = read();
int invsum = mi(add(add(p1, p2), add(p3, p4)), Mod - 2);
p1 = 1ll * p1 * invsum % Mod, p2 = 1ll * p2 * invsum % Mod;
p3 = 1ll * p3 * invsum % Mod, p4 = 1ll * p4 * invsum % Mod;
for (int i = -r; i <= r; ++i)
for (int j = -r; j <= r; ++j)
if (i * i + j * j <= r * r)
Pos.emplace_back(make_pair(i, j)), id[i + R][j + R] = Pos.size();
int n = Pos.size(), band = 0;
for (auto it : Pos) {
int x = it.first, y = it.second;
g[getid(x, y)][getid(x, y)] = g[getid(x, y)][n + 1] = 1;
if (getid(x - 1, y)) {
g[getid(x, y)][getid(x - 1, y)] = Mod - p1;
band = max(band, abs(getid(x, y) - getid(x - 1, y)));
}
if (getid(x, y - 1)) {
g[getid(x, y)][getid(x, y - 1)] = Mod - p2;
band = max(band, abs(getid(x, y) - getid(x, y - 1)));
}
if (getid(x + 1, y)) {
g[getid(x, y)][getid(x + 1, y)] = Mod - p3;
band = max(band, abs(getid(x, y) - getid(x + 1, y)));
}
if (getid(x, y + 1)) {
g[getid(x, y)][getid(x, y + 1)] = Mod - p4;
band = max(band, abs(getid(x, y) - getid(x, y + 1)));
}
}
Gauss(n, band);
printf("%d", ans[getid(0, 0)]);
return 0;
}
你有 \(p\) 滴血量，血量上限为 \(n\) 。每轮操作如下：
- 先以 \(\dfrac{1}{m + 1}\) 的概率增加 \(1\) 滴血，满血时则概率为 \(0\) 。
- \(k\) 次判定，每次以 \(\dfrac{1}{m + 1}\) 的概率减少一滴血，死了则概率为 \(0\) 。
求期望几轮死亡。
\(n \leq 1500\) ，\(m, k \leq 10^9\)
设 \(f_i\) 表示血量为 \(i\) 时期望多少轮结束，则：
\[f_i = 1 + \sum_{j = 0}^{\min(i, k)} (\dfrac{m}{m + 1} f_{i - j} + \dfrac{1}{m + 1} f_{i - j + 1}) \times \dbinom{k}{j} (\dfrac{1}{m + 1})^j (\dfrac{m}{m + 1})^{k - j}
\]
注意一下 \(i = n\) 时的情况即可。
观察到 \(f_i\) 只与 \(f_{0 \sim i + 1}\) 有关，所以矩阵应该是一个类似下三角矩阵的东西。因为我们高斯消元的时候是拿自己这行去减下面的，所以每一行中只有 \(2\) 个系数要去和下面的相减，时间复杂度就可以做到 \(O(n^2)\) 了。
#include <bits/stdc++.h>
using namespace std;
const int Mod = 1e9 + 7;
const int N = 1.5e3 + 7;
int g[N][N];
int inv[N], d[N], id[N], ans[N];
int n, p, m, k;
template <class T = int>
inline T read() {
char c = getchar();
bool sign = (c == '-');
while (c < '0' || c > '9')
c = getchar(), sign |= (c == '-');
T x = 0;
while ('0' <= c && c <= '9')
x = (x << 1) + (x << 3) + (c & 15), c = getchar();
return sign ? (~x + 1) : x;
}
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline int mi(int a, int b) {
int res = 1;
for (; b; b >>= 1, a = 1ll * a * a % Mod)
if (b & 1)
res = 1ll * res * a % Mod;
return res;
}
inline void prework() {
inv[0] = inv[1] = 1;
for (int i = 2; i < N; ++i)
inv[i] = 1ll * (Mod - Mod / i) * inv[Mod % i] % Mod;
}
inline void Gauss(int n, int band) {
iota(id + 1, id + 1 + n, 1);
for (int i = 1; i <= n; ++i) {
if (!g[i][i]) {
for (int j = i + 1; j <= min(n, i + band); ++j)
if (g[i][j]) {
for (int k = 1; k <= n; ++k)
swap(g[k][i], g[k][j]);
swap(id[i], id[j]);
break;
}
}
int Inv = mi(g[i][i], Mod - 2);
for (int j = i + 1; j <= n; ++j) {
int div = 1ll * g[j][i] * Inv % Mod;
for (int k = i; k <= min(n, i + band); ++k)
g[j][k] = dec(g[j][k], 1ll * g[i][k] * div % Mod);
g[j][n + 1] = dec(g[j][n + 1], 1ll * g[i][n + 1] * div % Mod);
}
}
for (int i = n; i; --i) {
ans[id[i]] = g[i][n + 1];
for (int j = i + 1; j <= min(n, i + band); ++j)
ans[id[i]] = dec(ans[i], 1ll * ans[id[j]] * g[i][j] % Mod);
ans[id[i]] = 1ll * ans[id[i]] * mi(g[i][i], Mod - 2) % Mod;
}
}
signed main() {
prework();
int T = read();
while (T--) {
n = read(), p = read(), m = read(), k = read();
if (!k) {
puts("-1");
continue;
} else if (!m) {
if (k == 1)
puts("-1");
else
printf("%d\n", (p + k - 2 - (p == n)) / (k - 1));
continue;
}
d[0] = 1ll * mi(m, k) * mi(mi(m + 1, k), Mod - 2) % Mod;
for (int i = 1, invm = mi(m, Mod - 2); i <= min(n, k); ++i)
d[i] = 1ll * d[i - 1] * invm % Mod * (k - i + 1) % Mod * inv[i] % Mod;
memset(g, 0, sizeof(g));
for (int i = 1, invm1 = mi(m + 1, Mod - 2); i < n; ++i) {
g[i][n + 1] = g[i][i] = 1;
for (int j = 0; j <= min(i, k); ++j) {
g[i][i - j] = dec(g[i][i - j], 1ll * dec(1, invm1) * d[j] % Mod);
g[i][i - j + 1] = dec(g[i][i - j + 1], 1ll * invm1 * d[j] % Mod);
}
}
g[n][n + 1] = g[n][n] = 1;
for (int i = 0; i <= min(n, k); ++i)
g[n][n - i] = dec(g[n][n - i], d[i]);
Gauss(n, 1);
printf("%d\n", ans[p]);
}
return 0;
}
有一个宽度为 \(w\) 高度为 \(h\) 的方格纸， $ w \times h$ 的格子中，有一些是空的，有一些是洞，有一些是障碍物。从第一行的空的格子中随机选一个放置一个球，向上下左右移动的概率比为 \(p_u : p_d : p_l : p_r\) ，不能移动到有障碍物的格子上。对于每个洞，输出落入该洞的概率。
\(2 \leq 20, h \leq 10^4, p_u + p_d + p_l + p_r = 100\) 。
和上面比较类似。
#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-12;
const int dx[] = {-1, 1, 0, 0};
const int dy[] = {0, 0, -1, 1};
const int N = 1e4 + 7, M = 2e1 + 7;
double g[N * M][M << 1], ans[N * M];
int id[N][M];
char str[N][M];
int n, m, p[4], tot;
inline bool check(int x, int y) {
return x && x <= n && y && y <= m;
}
inline void Gauss(int n, int band) {
for (int i = 1; i <= n; ++i) {
if (fabs(g[i][i]) < eps)
continue;
for (int j = i + 1; j <= min(i + band, n); ++j)
g[i][j] /= g[i][i];
ans[i] /= g[i][i], g[i][i] = 1;
for (int j = i + 1; j <= min(i + band, n); ++j) {
double div = g[j][i] / g[i][i];
for (int k = i; k <= min(i + band, n); ++k)
g[j][k] -= g[i][k] * div;
ans[j] -= ans[i] * div;
}
}
for (int i = n; i; --i)
for (int j = i + 1; j <= min(i + band, n); ++j)
ans[i] -= g[i][j] * ans[j];
}
signed main() {
scanf("%d%d", &m, &n);
for (int i = 0; i < 4; ++i)
scanf("%d", p + i);
for (int i = 1; i <= n; ++i) {
scanf("%s", str[i] + 1);
for (int j = 1; j <= m; ++j)
if (str[i][j] != 'X')
id[i][j] = ++tot;
}
int sum = m - count(id[1] + 1, id[1] + 1 + m, 0);
for (int i = 1; i <= n; ++i)
for (int j = 1; j <= m; ++j) {
if (!id[i][j])
continue;
if (i == 1)
ans[id[i][j]] = 1.0 / sum;
if (str[i][j] == 'T')
continue;
g[id[i][j]][id[i][j]] = 1;
int base = 100;
for (int k = 0; k < 4; ++k) {
int x = i + dx[k], y = j + dy[k];
if (!check(x, y) || !id[x][y])
base -= p[k];
}
for (int k = 0; k < 4; ++k) {
int x = i + dx[k], y = j + dy[k];
if (check(x, y) && id[x][y])
g[id[x][y]][id[i][j]] = -1.0 * p[k] / base;
}
}
Gauss(tot, m);
for (int i = 1; i <= n; ++i)
for (int j = 1; j <= m; ++j)
if (str[i][j] == 'T')
printf("%.9lf\n", ans[id[i][j]]);
return 0;
}
树上高消
若在树上进行高斯消元，往往一个叶子和其父亲的值呈一次函数关系。
那么类似数学归纳法地向上递推，可以发现所有点的值和根的值都是一次函数关系。
此时会发现所有值之和为定值，或者某些点（如根、叶子）的值容易求，那么就可以递推优化到 \(O(n)\) 。
给定一棵有根树，\(q\) 次询问从根走到 \(S\) 中的所有点至少一次的期望步数。
\(n \leq 18\) ，\(q \leq 5000\)
记 \(S'\) 表示根走到 \(S\) 中每个点分别的期望步数，考虑 Min-Max 容斥：
\[E(\max(S')) = \sum_{T' \subseteq S'} (-1)^{|T'| + 1} \min(T')
\]
考虑最值的意义，\(\max(S')\) 表示从根走完 \(S\) 中的所有点的期望步数，\(\min(S')\) 表示第一次走到 \(S\) 中的点的期望步数。设 \(f_{u, S}\) 表示 \(u\) 第一次走到 \(S\) 中的点的期望步数，\(d_u\) 表示 \(u\) 的度数则：
\[f_{u, S} = \frac{f_{fa_u, S} + \sum_{v \in son(u)} f_{v, S}}{d_u} + 1 \quad (x \not \in S) \\
f_{u, S} = 0 \quad (x \in S)
\]
考虑将每个点的值都写作 \(f_{u, S} = k_u \times f_{fa_u, S} + b_u\) 的形式，记：
\[K_u = \sum_{v \in son(u)} k_v, B_u = \sum_{v \in son(u)} b_v
\]
则得到：
\[f_{u, S} = \dfrac{1}{d_u - K_u} \times f_{fa_u, S} + \dfrac{d_u + B_u}{d_u - K_u}
\]
即：
\[k_u = \dfrac{1}{d_u - K_u}, b_u = \dfrac{d_u + B_u}{d_u - K_u}
\]
答案即为 \(\sum_{T \subseteq S} (-1)^{|T| + 1} f_{r, T}\) ，不难用高维前缀和预处理后 \(O(1)\) 查询。时间复杂度 \(O(n 2^n + q)\) 。
#include <bits/stdc++.h>
using namespace std;
const int Mod = 998244353;
const int N = 19;
struct Graph {
vector<int> e[N];
inline void insert(int u, int v) {
e[u].emplace_back(v);
}
} G;
int f[1 << N], g[1 << N];
int deg[N];
int n, q, r;
template <class T = int>
inline T read() {
char c = getchar();
bool sign = (c == '-');
while (c < '0' || c > '9')
c = getchar(), sign |= (c == '-');
T x = 0;
while ('0' <= c && c <= '9')
x = (x << 1) + (x << 3) + (c & 15), c = getchar();
return sign ? (~x + 1) : x;
}
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline int mi(int a, int b) {
int res = 1;
for (; b; b >>= 1, a = 1ll * a * a % Mod)
if (b & 1)
res = 1ll * res * a % Mod;
return res;
}
inline int sgn(int n) {
return n & 1 ? Mod - 1 : 1;
}
struct Node {
int k, b;
inline Node(const int _b = 0) : k(0), b(_b) {}
inline Node(const int _k, const int _b) : k(_k), b(_b) {}
inline Node operator + (const Node &rhs) const {
return Node(add(k, rhs.k), add(b, rhs.b));
}
} nd[N];
void dfs(int u, int f, int state) {
Node res = 0;
for (int v : G.e[u])
if (v != f)
dfs(v, u, state), res = res + nd[v];
if (~state >> u & 1) {
nd[u].k = mi(dec(deg[u], res.k), Mod - 2);
nd[u].b = 1ll * add(deg[u], res.b) * nd[u].k % Mod;
} else
nd[u] = 0;
}
signed main() {
n = read(), q = read(), r = read() - 1;
for (int i = 1; i < n; ++i) {
int u = read() - 1, v = read() - 1;
G.insert(u, v), G.insert(v, u);
++deg[u], ++deg[v];
}
for (int i = 1; i < (1 << n); ++i)
dfs(r, -1, i), f[i] = 1ll * nd[r].b * sgn(__builtin_popcount(i) + 1) % Mod;
for (int i = 0; i < n; ++i)
for (int j = 0; j < (1 << n); ++j)
if (j >> i & 1)
f[j] = add(f[j], f[j ^ (1 << i)]);
while (q--) {
int k = read(), state = 0;
while (k--)
state |= 1 << (read() - 1);
printf("%d\n", f[state]);
}
return 0;
}
