高斯消元

高斯消元法通常用于求解如下形式的 $n$ 元线性方程组：

{\begin{cases} a_{1, 1} x_{1} + a_{2, 2} x_{2} + \dots + a_{1, n} x_{n} = b_{1} \\ a_{2, 1} x_{1} + a_{2, 2} x_{2} + \dots + a_{2, n} x_{n} = b_{2} \\ \dots \\ a_{n, 1} x_{1} + a_{n, 2} x_{2} + \dots + a_{n, n} x_{n} = b_{n} \end{cases}

$\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} + \dots + a_{1, n} x_{n} = b_{1} \\ a_{2, 1} x_{1} + a_{2, 2} x_{2} + \dots + a_{2, n} x_{n} = b_{2} \\ ⋮ \\ a_{m, 1} x_{1} + a_{m, 2} x_{2} + \dots + a_{m, 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}$

可以表示为下面的矩阵：

[\begin{array}{ccccc} a_{1, 1} & a_{1, 2} & \dots & a_{1, n} & b_{1} \\ a_{2, 1} & a_{2, 2} & \dots & a_{2, n} & b_{2} \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ \\ a_{m, 1} & a_{m, 2} & \dots & a_{m, n} & b_{m} \end{array}]

$\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]$

整个矩阵称为增广矩阵，而没有右边系数的矩阵称为系数矩阵。

考虑把一个线性方程组的 $n$ 个元的解写在一个列向量（仅含一列的矩阵）中，如：

X = [\begin{matrix} x_{1} \\ x_{2} \\ ⋮ \\ x_{n} \end{matrix}]

$X = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix}$

于是可以得到向量方程 $AX = B$ 。

解线性方程组

P3389 【模板】高斯消元法

考虑求解如下方程组：

{\begin{cases} x_{1} + 3 x_{2} + 4 x_{3} = 5 \\ x_{1} + 4 x_{2} + 7 x_{3} = 3 \\ 9 x_{1} + 3 x_{2} + 2 x_{3} = 2 \end{cases}

$\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{matrix} 1 & 3 & 4 & 5 \\ 1 & 4 & 7 & 3 \\ 9 & 3 & 2 & 2 \end{matrix}]

$\begin{bmatrix} 1 & 3 & 4 & 5 \\ 1 & 4 & 7 & 3 \\ 9 & 3 & 2 & 2 \\ \end{bmatrix}$

把未知项的系数单独拿出来考虑：

[\begin{matrix} 1 & 3 & 4 \\ 1 & 4 & 7 \\ 9 & 3 & 2 \end{matrix}]

$\begin{bmatrix} 1 & 3 & 4 \\ 1 & 4 & 7 \\ 9 & 3 & 2 \\ \end{bmatrix}$

如果能够将其化为：

[\begin{matrix} ? & ? & ? \\ 0 & ? & ? \\ 0 & 0 & ? \end{matrix}]

$\begin{bmatrix} ? & ? & ? \\ 0 & ? & ? \\ 0 & 0 & ? \\ \end{bmatrix}$

就可以从最后一行逆推上去。设第一行第一列为主元，利用加减消元法配系数消掉后两行的第一列：

[\begin{matrix} 1 & 3 & 4 \\ 0 & 1 & 3 \\ 0 & - 24 & - 34 \end{matrix}]

$\begin{bmatrix} 1 & 3 & 4 \\ 0 & 1 & 3 \\ 0 & -24 & -34 \\ \end{bmatrix}$

再以第二行第二列为主元，加减消元得到：

[\begin{matrix} 1 & 3 & 4 \\ 0 & 1 & 3 \\ 0 & 0 & 38 \end{matrix}]

$\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 mxp = i;
 
        for (int j = i + 1; j <= n; ++j)
            if (fabs(g[j][i]) > fabs(g[mxp][i]))
                mxp = j;
 
        if (fabs(g[mxp][i]) < eps)
            return false;
 
        if (mxp != i)
            swap(g[mxp], g[i]);
 
        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] -= div * g[i][k];
        }
    }
 
    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;
}
 
signed 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{matrix} ? & 0 & 0 & ? \\ 0 & ? & 0 & ? \\ 0 & 0 & ? & ? \end{matrix}]

$\begin{bmatrix} ? & 0 & 0 & ? \\ 0 & ? & 0 & ? \\ 0 & 0 & ? & ? \\ \end{bmatrix}$

首先，我们依照朴素的高斯消元不难得到：

[\begin{matrix} ? & ? & ? & ? \\ 0 & ? & ? & ? \\ 0 & 0 & ? & ? \end{matrix}]

$\begin{bmatrix} ? & ? & ? & ? \\ 0 & ? & ? & ? \\ 0 & 0 & ? & ? \\ \end{bmatrix}$

观察上下两个矩阵，不难得到消元时同时消去主元上下方的方程即可。

注意一下特殊情况：

无解：当左边系数为 $0$ 而右边系数不为 $0$ 时即为无解。
多解：当左边系数为 $0$ 而右边系数也为 $0$ 时即为多解。

 inline bool Gauss() {
	for (int i = 1; i <= n; ++i) {
		int mxp = i;
		
		for (int j = i + 1; j <= n; ++j)
			if (fabs(g[j][i]) > fabs(g[mxp][i]))
				mxpos = j;
		
		if (fabs(g[mxp][i]) < eps)
			return false;
       	
        if (i != mxp)
            swap(g[i], g[mxp]);
		
		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 \dots \oplus a_{1, n} x_{n} = b_{1} \\ a_{2, 1} x_{1} \oplus a_{2, 2} x_{2} \oplus \dots \oplus a_{2, n} x_{n} = b_{2} \\ \dots \\ a_{n, 1} x_{1} \oplus a_{n, 2} x_{2} \oplus \dots \oplus a_{n, n} x_{n} = b_{n} \end{cases}

$\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 \in \{ 0, 1 \}$ 。

由于 $\oplus$ 符合交换律与结合律，故可以按照高斯消元法逐步消元求解。值得注意的是，在消元的时候应使用异或消元而非加减消元，且不需要进行乘除改变系数（因为系数均为 $0$ 和 $1$ ）。

可以用 bitset 优化到 $O(\frac{n^3}{\omega})$ 。

 inline bool Gauss() {
	for (int i = 1; i <= n; ++i) {
		int mxp = i;
 
		while (mxpos <= n && !g[mxpos].test(i))
			++mxp;
 
		if (mxpos > n)
			return false;
		
        if (i != mxp)
			swap(g[i], g[mxp]);
 
		for (int j = 1; j <= n; ++j)
			if (j != i && g[j].test(i))
				g[j] ^= g[i];
	}
 
	return true;
}

P3429 [POI 2005] DWA-Two Parties

给出一张无向图，需要给每个点染上黑色或白色。定义一个点合法当且仅当其所有同色邻居数量为偶数，最大化合法点的数量，并构造方案。

$n \leq 200$ ， $m \leq \frac{n(n - 1)}{2}$

可以证明答案始终为 $n$ 。构造就是求解如下异或方程组：

⨁_{(u, v) \in G} (x_{u} \oplus x_{v}) = \deg_{u} mod 2

$\bigoplus_{(u, v) \in G} (x_u \oplus x_v) = \deg_u \bmod 2$

这个方程始终有解，否则若出现能将其消为左边没有元而右边为 $1$ 的情况，其说明存在奇数个奇度数的点。而每个奇度数的点要选奇数个邻居在集合中，偶度数点要选偶数个邻居在集合中，从而其导出子图的度数为奇数，矛盾。

 #include <bits/stdc++.h>
using namespace std;
const int N = 2e2 + 7;
 
bitset<N> g[N];
 
int n;
 
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 void Gauss() {
    for (int i = 1; i <= n; ++i) {
        if (!g[i].test(i)) {
            for (int j = i + 1; j <= n; ++j)
                if (g[j].test(i)) {
                    swap(g[i], g[j]);
                    break;
                }
        }
 
        for (int j = 1; j <= n; ++j)
            if (j != i && g[j].test(i))
                g[j] ^= g[i];
    }
}
 
signed main() {
    n = read();
 
    for (int i = 1; i <= n; ++i) {
        int d = read();
 
        if (d & 1)
            g[i].set(i), g[i].set(n + 1);
 
        while (d--)
            g[i].set(read());
    }
 
    Gauss();
    int cnt = 0;
 
    for (int i = 1; i <= n; ++i)
        if (g[i].test(n + 1))
            ++cnt;
 
    printf("%d\n", cnt);
 
    for (int i = 1; i <= n; ++i)
        if (g[i].test(n + 1))
            printf("%d ", i);
    
    return 0;
}

矩阵求逆

P4783 【模板】矩阵求逆

对于矩阵 $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 mxp = i;
 
        while (mxp <= n && !g[mxp][i])
            ++mxp;
 
        if (mxp > n)
            return false;
 
        if (mxp != i)
            swap(g[i], g[mxp]);
 
        int inv = mi(g[i][i], Mod - 2);
 
        for (int j = 1; j <= n; ++j)
            if (j != i) {
                int div = 1ll * g[j][i] * inv % Mod;
 
                for (int k = i; k <= n * 2; ++k)
                    g[j][k] = dec(g[j][k], 1ll * div * g[i][k] % Mod);
            }
    }
 
    for (int i = 1; i <= n; ++i) {
        int inv = mi(g[i][i], Mod - 2);
 
        for (int j = n + 1; j <= n * 2; ++j)
            g[i][j] = 1ll * g[i][j] * inv % Mod;
    }
 
    return true;
}

行列式求值

P7112 【模板】行列式求值

考虑一个情况，当一个矩阵任意一个位置出现 $0$ ，其对行列式的影响非常大（直接没有贡献了）。

我们利用上面的一些性质显然是可以让矩阵不断变化出现 $0$ 的。考虑将矩阵一行（列）消成只有最后一个元素非 $0$ 该怎么做。也就是说：

[\begin{matrix} a_{1, 1} & a_{1, 2} & \dots & a_{1, n} \\ a_{2, 1} & a_{2, 2} & \dots & a_{2, n} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ a_{n, 1} & a_{n, 2} & \dots & a_{n, n} \end{matrix}] \Rightarrow [\begin{matrix} a_{1, 1} & a_{1, 2} & \dots & a_{1, n} \\ a_{2, 1} & a_{2, 2} & \dots & a_{2, n} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & a_{n, n} \end{matrix}]

$\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;
    }
}

CF24D Broken robot

$n$ 行 $m$ 列的矩阵，现在在 $(x,y)$ ，每次等概率向左、右、下走或原地不动，但不能走出去，求走到最后一行期望的步数。

$n, m \leq 10^3$

记 $f_{i, j}$ 表示机器人在 $(i, j)$ 时走到最后一行的期望步数，则：

$m = 1$ 时有（省略第二维）：

f_{i} = 1 + \frac{1}{2} (f_{i + 1} + f_{i})

$f_i = 1 + \dfrac{1}{2} (f_{i + 1} + f_i)$

即：

f_{i} = f_{i + 1} + 2

$f_i = f_{i + 1} + 2$

$m > 1$ 时有：

f_{i, j} = 1 + {\begin{cases} \frac{1}{3} (f_{i, j} + f_{i + 1, j} + f_{i, j + 1}), & j = 1 \\ \frac{1}{4} (f_{i, j} + f_{i + 1, j} + f_{i, j - 1} + f_{i, j + 1}), & 1 < j < m \\ \frac{1}{3} (f_{i, j} + f_{i + 1, j} + f_{i, j - 1}), & j = m \end{cases}

$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;
}

CF963E Circles of Waiting

平面直角坐标系上有一个点，一开始在 $(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;
}

P4457 [BJOI2018] 治疗之雨

你有 $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)} (\frac{m}{m + 1} f_{i - j} + \frac{1}{m + 1} f_{i - j + 1}) \times (\binom{k}{j}) (\frac{1}{m + 1})^{j} (\frac{m}{m + 1})^{k - j}

$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;
}

P6899 [ICPC2014 WF] Pachinko

有一个宽度为高度为的方格纸， $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)$ 。

P5643 [PKUWC2018] 随机游走

给定一棵有根树， $q$ 次询问从根走到 $S$ 中的所有点至少一次的期望步数。

$n \leq 18$ ， $q \leq 5000$

求经过 $S$ 里所有元素的期望时间，即到达 $S$ 中最后一个点的期望步数（ $\max$ ），那么可以转化为枚举 $S$ 的子集 $T$ ，求到达 $T$ 中第一个元素的期望时间（ $\min$ ）。

考虑 Min-Max 容斥，设 $f_{u, S}$ 表示 $u$ 第一次走到 $S$ 中的点的期望步数， $d_u$ 表示 $u$ 的度数，则：

f_{u, S} = \frac{f_{f a_{u}, S} + \sum_{v \in s o n (u)} f_{v, S}}{d_{u}} + 1 (x \notin S) f_{u, S} = 0 (x \in S)

$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 s o n (u)} k_{v}, B_{u} = \sum_{v \in s o n (u)} b_{v}

$K_u = \sum_{v \in son(u)} k_v, B_u = \sum_{v \in son(u)} b_v$

则得到：

f_{u, S} = \frac{1}{d_{u} - K_{u}} \times f_{f a_{u}, S} + \frac{d_{u} + B_{u}}{d_{u} - K_{u}}

$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} = \frac{1}{d_{u} - K_{u}}, b_{u} = \frac{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;
}

P11736 [集训队互测 2015] 胡策的小树

有一棵树，点权 $a_{1 \sim n}$ 构成了一个 $0 \sim n - 1$ 的排列，并且满足 $a_1 = 0$ 。

初始时，每个节点各有一只猴子，每一秒第 $i$ 个点上的猴子会有 $p(i)$ 的概率跳到父亲， $\frac{1 - p(i)}{siz_i}$ 的概率跳到子树内的任意点，其中 $p(i) = \begin{cases} 0 & i = 1 \\ \frac{a_i}{n} & 2 \leq i \leq n \end{cases}$ 。

初始时将所有点的权值变为 $(a_i + x) \bmod n$ ，其中 $x$ 为自行选定任意非负整数。

记第 $i$ 秒成功跳到父亲的猴子数量为 $g_i$ ，幸福指数被定义为 $g_{0 \sim +\infty}$ 的平均数。选取适当的 $x$ ，求幸福指数期望的最大值。

$n \leq 5 \times 10^5$ ，节点 $i$ 的父亲从 $1 \sim i - 1$ 中等概率选取

先考虑初始操作 $x = 0$ 的情况。可以发现每只猴子对答案的贡献是独立的，进一步发现每只猴子的初始位置是不重要的。因为时间无穷，因此一只猴子一定会跳到根，并可以忽略跳到根之前的贡献。下面讨论从根开始跳的答案。

记 $p_i = p(i), q_i = \frac{1 - p_i}{siz_i}$ ，设猴子处于 $u$ 的概率为 $f_u$ ，则：

f_{u} = (\sum_{w \in a n c (u) \cup {u}} q_{w} f_{w}) + (\sum_{v \in s o n (u)} p_{v} f_{v}) \sum f_{u} = 1

$f_u = \left( \sum_{w \in anc(u) \cup \{ u \}} q_w f_w \right) + \left( \sum_{v \in son(u)} p_v f_v \right) \\ \sum f_u = 1$

暴力高消可以做到 $O(n^3)$ ，但是没有利用树上高消的性质。

枚举祖先并不好树上高消，记 $F_u = \sum_{w \in anc(u)} q_w f_w$ ，则 $f_u = \frac{1}{q_u} (F_u - F_{fa_u})$ ，带入得到：

\frac{1}{q_{u}} (F_{u} - F_{f a_{u}}) = F_{u} + \sum_{v \in s o n (u)} \frac{p_{v}}{q_{v}} (F_{v} - F_{u}) (\frac{1}{q_{u}} - 1 + \sum_{v \in s o n (u)} \frac{p_{v}}{q_{v}}) F_{u} = \frac{1}{q_{u}} F_{f a_{u}} + \sum_{v \in s o n (u)} \frac{p_{v}}{q_{v}} F_{v}

$\frac{1}{q_u} (F_u - F_{fa_u}) = F_u + \sum_{v \in son(u)} \frac{p_v}{q_v} (F_v - F_u) \\ (\frac{1}{q_u} - 1 + \sum_{v \in son(u)} \frac{p_v}{q_v}) F_u = \frac{1}{q_u} F_{fa_u} + \sum_{v \in son(u)} \frac{p_v}{q_v} F_v$

记：

\begin{aligned} (1) & A_{u} & = \frac{1}{q_{u}} - 1 + \sum_{v \in s o n (u)} \frac{p_{v}}{q_{v}} \\ (2) & B_{u} & = \frac{1}{q_{u}} \\ (3) & C_{u} & = \frac{p_{u}}{q_{u}} \end{aligned}

$\begin{align} A_u &= \frac{1}{q_u} - 1 + \sum_{v \in son(u)} \frac{p_v}{q_v} \\ B_u &= \frac{1}{q_u} \\ C_u &= \frac{p_u}{q_u} \end{align}$

则：

A_{u} F_{u} = B_{u} F_{f a_{u}} + \sum_{v \in s o n (u)} C_{v} F_{v}

$A_u F_u = B_u F_{fa_u} + \sum_{v \in son(u)} C_v F_v$

考虑设 $F_u = k_u F_{fa_u} + b_u$ ，则：

\begin{aligned} (4) & A_{u} (k_{u} F_{f a_{u}} + b_{u}) & = B_{u} F_{f a_{u}} + \sum_{v \in s o n (u)} C_{v} (k_{v} F_{u} + b_{v}) \\ (5) & (A_{u} - \sum_{v \in s o n (u)} C_{v} k_{v}) F_{u} & = B_{u} F_{f a_{u}} + \sum_{v \in s o n (u)} C_{v} b_{v} \end{aligned}

$\begin{align} A_u (k_u F_{fa_u} + b_u) &= B_u F_{fa_u} + \sum_{v \in son(u)} C_v (k_v F_u + b_v) \\ (A_u - \sum_{v \in son(u)} C_v k_v) F_u &= B_u F_{fa_u} + \sum_{v \in son(u)} C_v b_v \end{align}$

于是可以 $O(n)$ 高消求解 $x = 0$ 时幸福指数的期望 $\sum p_i f_i$ 。

考虑 $x \neq 0$ 的情况，此时一定存在一个点 $i$ 满足 $p(i) = 0$ ，也就是说跳到 $i$ 子树内后就跳不出来了。于是对于每个点的子树做一遍树上高消即可。

时间复杂度 $O(\sum siz_i)$ ，因为树的形态随机，因此时间复杂度 $O(n \log n)$ 。

 #include <bits/stdc++.h>
typedef long double ld;
using namespace std;
const int N = 5e5 + 7;
 
struct Graph {
    vector<int> e[N];
    
    inline void insert(int u, int v) {
        e[u].emplace_back(v);
    }
} G;
 
struct Node {
    ld k, b;
 
    inline Node() {}
 
    inline Node(ld _b) : b(_b) {}
 
    inline Node(ld _k, ld _b) : k(_k), b(_b) {}
 
    inline Node operator + (Node rhs) {
        return Node(k + rhs.k, b + rhs.b);
    }
 
    inline Node operator - (Node rhs) {
        return Node(k - rhs.k, b - rhs.b);
    }
 
    inline Node operator * (ld x) {
        return Node(k * x, b * x);
    }
 
    inline Node operator / (ld x) {
        return Node(k / x, b / x);
    }
} nd[N], F[N];
 
ld p[N], q[N], A[N], B[N], C[N], f[N];
int fa[N], a[N], siz[N], in[N], out[N], id[N];
 
int n, dfstime;
 
void dfs(int u) {
    siz[u] = 1, id[in[u] = ++dfstime] = u;
 
    for (int v : G.e[u])
        dfs(v), siz[u] += siz[v];
 
    out[u] = dfstime;
}
 
inline ld solve(int r) {
    for (int i = in[r]; i <= out[r]; ++i) {
        int u = id[i];
        p[u] = (ld)((a[u] + n - a[r]) % n) / n;
        q[u] = (1 - p[u]) / siz[u];
        A[u] = 1 / q[u] - 1, B[u] = 1 / q[u], C[u] = p[u] / q[u];
 
        if (u != r)
            A[fa[u]] += p[u] / q[u];
    }
 
    for (int i = out[r]; i > in[r]; --i) {
        int u = id[i];
        nd[u] = Node(B[u], 0);
        ld div = A[u];
 
        for (int v : G.e[u])
            nd[u].b += C[v] * nd[v].b, div -= C[v] * nd[v].k;
 
        nd[u] = nd[u] / div;
    }
 
    F[r] = Node(1, 0);
    Node all = F[r] / q[r];
 
    for (int i = in[r] + 1; i <= out[r]; ++i) {
        int u = id[i];
        F[u] = F[fa[u]] * nd[u].k + nd[u].b;
        all = all + (F[u] - F[fa[u]]) / q[u];
    }
 
    f[r] = (1 - all.b) / all.k / q[r];
 
    for (int i = in[r] + 1; i <= out[r]; ++i) {
        int u = id[i];
        Node now = (F[u] - F[fa[u]]) / q[u];
        f[u] = now.k * f[r] * q[r] + now.b;
    }
 
    ld res = 0;
 
    for (int i = in[r]; i <= out[r]; ++i) {
        int u = id[i];
        res += f[u] * p[u];
    }
 
    return res;
}
 
signed main() {
    scanf("%d", &n);
 
    for (int i = 1; i <= n; ++i)
        scanf("%d", fa + i), G.insert(fa[i], i);
 
    for (int i = 1; i <= n; ++i)
        scanf("%d", a + i);
 
    dfs(1);
    ld ans = 0;
 
    for (int i = 1; i <= n; ++i)
        ans = max(ans, solve(i));
 
    printf("%.9LF", ans);
    return 0;
}

posted @ 2024-07-22 19:48 wshcl 阅读(20) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

相关博文：

· 点分治相关

· 连通性相关

· 高斯消元法求线性方程组

· 高斯消元学习笔记

· 「Note」高斯消元（未完）

阅读排行：
· 地球OL攻略 —— 某应届生求职总结
· 周边上新：园子的第一款马克杯温暖上架
· Open-Sora 2.0 重磅开源！
· 提示词工程——AI应用必不可少的技术
· .NET周刊【3月第1期 2025-03-02】

公告

昵称： wshcl
园龄： 4年7个月
粉丝： 12
关注： 13

+加关注

2025年3月

日

一

二

三

四

五

六

wshcl

awa

高斯消元

高斯消元

线性方程组与矩阵的联系

解线性方程组

约旦消元法

解异或方程组

P3429 [POI 2005] DWA-Two Parties

矩阵求逆

行列式求值

band-matrix

CF24D Broken robot

CF963E Circles of Waiting

P4457 [BJOI2018] 治疗之雨

P6899 [ICPC2014 WF] Pachinko

树上高消

P5643 [PKUWC2018] 随机游走

P11736 [集训队互测 2015] 胡策的小树

公告

搜索

常用链接

最新随笔

积分与排名

随笔档案 (28)

阅读排行榜

推荐排行榜

	#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 mxp = i;

	for (int j = i + 1; j <= n; ++j)
	if (fabs(g[j][i]) > fabs(g[mxp][i]))
	mxp = j;

	if (fabs(g[mxp][i]) < eps)
	return false;

	if (mxp != i)
	swap(g[mxp], g[i]);

	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] -= div * g[i][k];
	}
	}

	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;
	}

	signed 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;
	}

	#include <bits/stdc++.h>
	using namespace std;
	const int N = 2e2 + 7;

	bitset<N> g[N];

	int n;

	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 void Gauss() {
	for (int i = 1; i <= n; ++i) {
	if (!g[i].test(i)) {
	for (int j = i + 1; j <= n; ++j)
	if (g[j].test(i)) {
	swap(g[i], g[j]);
	break;
	}
	}

	for (int j = 1; j <= n; ++j)
	if (j != i && g[j].test(i))
	g[j] ^= g[i];
	}
	}

	signed main() {
	n = read();

	for (int i = 1; i <= n; ++i) {
	int d = read();

	if (d & 1)
	g[i].set(i), g[i].set(n + 1);

	while (d--)
	g[i].set(read());
	}

	Gauss();
	int cnt = 0;

	for (int i = 1; i <= n; ++i)
	if (g[i].test(n + 1))
	++cnt;

	printf("%d\n", cnt);

	for (int i = 1; i <= n; ++i)
	if (g[i].test(n + 1))
	printf("%d ", i);

	return 0;
	}

	inline bool Gauss() {
	for (int i = 1; i <= n; ++i)
	g[i][i + n] = 1;

	for (int i = 1; i <= n; ++i) {
	int mxp = i;

	while (mxp <= n && !g[mxp][i])
	++mxp;

	if (mxp > n)
	return false;

	if (mxp != i)
	swap(g[i], g[mxp]);

	int inv = mi(g[i][i], Mod - 2);

	for (int j = 1; j <= n; ++j)
	if (j != i) {
	int div = 1ll * g[j][i] * inv % Mod;

	for (int k = i; k <= n * 2; ++k)
	g[j][k] = dec(g[j][k], 1ll * div * g[i][k] % Mod);
	}
	}

	for (int i = 1; i <= n; ++i) {
	int inv = mi(g[i][i], Mod - 2);

	for (int j = n + 1; j <= n * 2; ++j)
	g[i][j] = 1ll * g[i][j] * inv % Mod;
	}

	return true;
	}

	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;
	}

	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;
	}

	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;
	}
	}

	#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;
	}

	#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;
	}

	#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;
	}

	#include <bits/stdc++.h>
	typedef long double ld;
	using namespace std;
	const int N = 5e5 + 7;

	struct Graph {
	vector<int> e[N];

	inline void insert(int u, int v) {
	e[u].emplace_back(v);
	}
	} G;

	struct Node {
	ld k, b;

	inline Node() {}

	inline Node(ld _b) : b(_b) {}

	inline Node(ld _k, ld _b) : k(_k), b(_b) {}

	inline Node operator + (Node rhs) {
	return Node(k + rhs.k, b + rhs.b);
	}

	inline Node operator - (Node rhs) {
	return Node(k - rhs.k, b - rhs.b);
	}

	inline Node operator * (ld x) {
	return Node(k * x, b * x);
	}

	inline Node operator / (ld x) {
	return Node(k / x, b / x);
	}
	} nd[N], F[N];

	ld p[N], q[N], A[N], B[N], C[N], f[N];
	int fa[N], a[N], siz[N], in[N], out[N], id[N];

	int n, dfstime;

	void dfs(int u) {
	siz[u] = 1, id[in[u] = ++dfstime] = u;

	for (int v : G.e[u])
	dfs(v), siz[u] += siz[v];

	out[u] = dfstime;
	}

	inline ld solve(int r) {
	for (int i = in[r]; i <= out[r]; ++i) {
	int u = id[i];
	p[u] = (ld)((a[u] + n - a[r]) % n) / n;
	q[u] = (1 - p[u]) / siz[u];
	A[u] = 1 / q[u] - 1, B[u] = 1 / q[u], C[u] = p[u] / q[u];

	if (u != r)
	A[fa[u]] += p[u] / q[u];
	}

	for (int i = out[r]; i > in[r]; --i) {
	int u = id[i];
	nd[u] = Node(B[u], 0);
	ld div = A[u];

	for (int v : G.e[u])
	nd[u].b += C[v] * nd[v].b, div -= C[v] * nd[v].k;

	nd[u] = nd[u] / div;
	}

	F[r] = Node(1, 0);
	Node all = F[r] / q[r];

	for (int i = in[r] + 1; i <= out[r]; ++i) {
	int u = id[i];
	F[u] = F[fa[u]] * nd[u].k + nd[u].b;
	all = all + (F[u] - F[fa[u]]) / q[u];
	}

	f[r] = (1 - all.b) / all.k / q[r];

	for (int i = in[r] + 1; i <= out[r]; ++i) {
	int u = id[i];
	Node now = (F[u] - F[fa[u]]) / q[u];
	f[u] = now.k * f[r] * q[r] + now.b;
	}

	ld res = 0;

	for (int i = in[r]; i <= out[r]; ++i) {
	int u = id[i];
	res += f[u] * p[u];
	}

	return res;
	}

	signed main() {
	scanf("%d", &n);

	for (int i = 1; i <= n; ++i)
	scanf("%d", fa + i), G.insert(fa[i], i);

	for (int i = 1; i <= n; ++i)
	scanf("%d", a + i);

	dfs(1);
	ld ans = 0;

	for (int i = 1; i <= n; ++i)
	ans = max(ans, solve(i));

	printf("%.9LF", ans);
	return 0;
	}