Kirchhoff's Matrix Tree 定理与 BEST 定理随记
Kirchhoff 矩阵树定理适用于解决图的生成树相关计数问题。BEST 定理则适用于解决有向欧拉图的欧拉路径计数问题。
本文中的图允许重边，不允许自环（其实可以有自环，甚至不会对定理内容产生任何影响，但是为方便叙述不允许自环）。下文记矩阵 \(A_{[S], [T]}\) 为选取 \(A_{i, j} (i \in S, j \in T)\) 组成的子矩阵，\([n] = \{1, 2, \dots, n\}\)。
Matrix Tree 定理（无向图形式）
记 \(A\) 为邻接矩阵，\(D\) 为度数矩阵，其中 \(D_{i, i} = \deg_i\)，当 \(i \neq j\) 时 \(D_{i, j} = 0\)。记 \(\mathcal{T}\) 表示图的生成树集合，定义一棵生成树 \(T\) 的权值为 \(\omega(T) = \prod \limits_{e \in T} \omega(e)\)，其中 \(\omega(e)\) 为边 \(e\) 的边权。
定义 Laplace 矩阵 \(L = D - A\)，则 \(\det L_{[n] \backslash \{k\}, [n] \backslash \{k\}} = \sum \limits_{T \in \mathcal{T}} \omega(T)\)，其中 \(k\) 任意。用自然语言描述即，生成树的权值之和等于 \(L\) 任一 \(n - 1\) 阶主子式的 \(\det\)。
下面对此定理展开证明。
首先有一个重要引理，Cauchy-Binet 引理：对于一个 \(n \times m\) 的矩阵 \(A\) 与一个 \(m \times n\) 的矩阵 \(B\)，有
当 \(m < n\) 时显然 \(AB\) 不满秩，等式右侧也可被认为是 \(0\)，引理自然成立。
这里主要探讨 \(m \geq n\) 的情况。作为线性代数中的一条重要公式，学界已经给出了各式各样的代数证明。此处参考 OI wiki，给出一份简短而巧妙的证明。前置知识：[NOI2021] 路径交点。
考虑“路径交点”的模型，构造一个层数为 \(3\) 的 DAG：左部点与右部点均有 \(n\) 个，中部点 \(m\) 个。左部点 \(i\) 向中部点 \(j\) 连权为 \(a_{i, j}\) 的边，中部点 \(i\) 向右部点 \(j\) 连权为 \(b_{i, j}\) 的边。
容易发现 \(AB\) 即为 LGV 引理中的 \(\mathbf{M}\)，而引理中等式右侧计算的值等于“路径交点”模型中交点数为偶数的不交路径组的权值和，减去交点数为奇数的不交路径组的权值和。根据前置经验，我们知道 \(\det \mathbf{M}\) 与之相等，故等式成立。\(\square\)
回到 Matrix Tree 定理的证明上。一个事实是，\(L\) 矩阵有一个很好的关于“关联矩阵”的性质。其中关联矩阵 \(P\) 为一大小为 \(|V| \times |E|\) 的矩阵，对于边 \(e_i = (u, v, w), u < v\)，\(P_{u, i} = 1, P_{v, i} = -1\)，其余元素均为 \(0\)。定义另一个关联矩阵 \(Q\) 为一大小为 \(|E| \times |V|\) 的矩阵，对于边 \(e_i = (u, v, w), u < v\)，\(Q_{i, u} = w, Q_{i, v} = -w\)，其余元素均为 \(0\)。则 \(L\) 矩阵具有这样的性质：\(L = PQ\)。
上述转化是证明的核心步骤之一。记 \(n = |V|, m = |E|\)，记 \(L' = L_{[n] \backslash \{1\}, [n] \backslash \{1\}}\) 为 \(L\) 去掉第一行第一列的子矩阵，\(P', Q'\) 分别为 \(P, Q\) 去掉第一行的子矩阵，则 \(L' = P'Q'\)。对 \(L' = P'Q'\) 应用 Cauchy-Binet 引理，得到 \(\det L' = \sum \limits_{S \subseteq [m], |S| = n - 1} (\det P'_{[n - 1], S}) (\det Q'_{S, [n - 1]})\)。
考虑 \(\det P'_{[n - 1], S}\) 的组合意义：对标号为 \(2, 3, \dots, n\) 的点，分别选择一条 \(S\) 中的边（若点 \(i\) 选择边 \(j\)，则认为 \(i\) 是边 \(j\) 的起点），所有方案的带符号权值和。考虑边集 \(S\) 生成了至少一个环的情况。由于我们的选择是有向的，因此任意一个无向环都有两种定向方式，对应行列式上两种不同贡献。将点编号最小值最小的环反向，记环长为 \(k\)，由于所有边的起点终点交换，相当于 \(P\) 上对应列的 \(1\) 与 \(-1\) 交换，因此贡献乘以了 \((-1) ^ k\)。另一方面是行列式逆序对数的变化，交换对应列的 \(1\) 与 \(-1\) 后形成的新排列可以通过原排列执行 \(k - 1\) 次 swap 操作形成，因此贡献还需乘以 \((-1) ^ {k - 1}\)。无论 \(k\) 是奇数还是偶数，将最小环反转后贡献都会乘以 \(-1\)，而我们的反转操作显然构成对合映射，因此若 \(S\) 生成至少一个环，则 \(\det P'_{[n - 1], S}\) 为 \(0\)。
否则边集 \(S\) 成功生成了一棵树，而我们对这张基图也有唯一的定向方式：以 \(1\) 为根的根向树。记在计算 \(\det P'_{[n - 1], S}\) 时这种定向方式下产生的逆序对数为 \(t\)，每一列选择的权值的积为 \((-1) ^ x\)，则 \(\det P'_{[n - 1], S} = (-1) ^ {t + x}\)。显然 \(Q'\) 和 \(P'\) 的情况是对称的，它的 \(\det\) 值是 \((-1) ^ {t + x} \times \prod \limits_{e \in S} \omega(e)\)。故最后算出的值等于 \(S\) 边集对应的 \(\omega(T)\)。
求和后自然得到 \(\sum \limits_{T \in \mathcal{T}} \omega(T)\)。由于无向情况下答案与钦定的有根树根节点并无关系，因此任一 \(n - 1\) 阶主子式的 \(\det\) 值均相等，至此完成了无向图部分的证明。\(\square\)
其实矩阵树定理还有另一种形式：所求即为 Laplace 矩阵所有非零特征值（共 \(n - 1\) 个）的乘积除以 \(n\)，证明咕。一些情况下这种形式可以更便于推导问题，下面是一个例子，不过这种形式应该不怎么常见。
给定一张 \(n\) 部完全图的 \(n\) 个部分的顶点数量 \(a_1, a_2, \dots, a_n\)，记 \(N = \sum a_i\)，那么该图的生成树数量为 \(N^{n - 2}\prod\limits_{i = 1} ^ n (N - a_i) ^ {a_i - 1}\)。
Matrix Tree 定理（有向图形式）
较无向图下的矩阵树定理而言，有向图下的定理乃至证明均只有一些细微的区别。
记 \(A\) 为邻接矩阵，\(D^{\text{in}}\) 为入度矩阵，\(D^{\text{out}}\) 为出度矩阵。定义 Laplace 矩阵 \(L^{\text{in}} = D^{\text{in}} - A\)，\(L^{\text{out}} = D^{\text{out}} - A\)。
则以 \(k\) 为根的根向树数量为 \(\det {L^{\text{out}}}_{[n] \backslash \{k\}, [n] \backslash \{k\}}\)，叶向树数量为 \(\det {L^{\text{in}}}_{[n] \backslash \{k\}, [n] \backslash \{k\}}\)。
考虑证明。仍然类似地定义关联矩阵 \(P, Q\)，对于每条边 \(e_i = (u, v, w)\)，
-
若所求为叶向树，\(P_{u, i} = -1, P_{v, i} = 1, Q_{i, u} = 0, Q_{i, v} = w\)。此时有 \(L ^ {\text{in}} = PQ\)。
-
否则对于根向树，\(P_{u, i} = w, P_{v, i} = 0, Q_{i, u} = 1, Q_{i, v} = -1\)。此时有 \(L ^ {\text{out}} = PQ\)。
接下来的证明流程完全同前，如果要求树以 \(k\) 为根，仍然令 \(P', Q'\) 分别为 \(P, Q\) 去掉第 \(k\) 行形成的子矩阵，那么 \(L' = P'Q'\)。应用 Cauchy-Binet 引理得到 \(\det L' = \sum \limits_{S \subseteq [m], |S| = n - 1} (\det P'_{[n - 1], S}) (\det Q'_{S, [n - 1]})\)。
以叶向树为例，\(\det P'_{[n - 1], S}\) 的组合意义与无向情况完全一致，而由于 \(\det Q'_{S, [n - 1]}\) 的 \(n - 1\) 行分别只有一项有值，因此这一项非 \(0\) 当且仅当每列都有恰好一项有值，对应到图上意味着有向边集 \(S\) 生成的子图中每个点的入度恰好为 \(1\)。如果满足，则说明边集 \(S\) 构成以 \(k\) 为根的叶向树，此时定向方式唯一，\(\det P'_{[n - 1], S}\) 在每一列上也一定会选择正值项，因此最后求出的即为 \(\omega(T)\)。求和后自然得到 \(\sum \limits_{T \in \mathcal{T}} \omega(T)\)。\(\square\)
BEST 定理
记 \(f ^ {\text{root}}(G, k)\) 表示有向图 \(G\) 的以 \(k\) 为根的根向生成树数量，由 Matrix Tree 定理我们知道它等于 \(\det {L ^ {\text{out}}}_{[n] \backslash \{k\}, [n] \backslash \{k\}}\)。定义有向欧拉图 \(G\) 的不同欧拉回路数（这里认为重边算作互不相同的边）为 \(\mathrm{ec}(G)\)。那么 BEST 定理表明
这里的 \(k\) 可以任取，在证明中它代表了欧拉回路的起点。有向欧拉图上 \(\deg^{\text{in}}_u = \deg^{\text{out}}_u\)，统称之为 \(\deg_u\)。
注意这是考虑循环同构的情况。如果认为欧拉回路从 \(k\) 点出发，只考虑边的经过顺序是否相同而不考虑是否循环同构时，此式需要乘以 \(\deg_k\)。下文始终考虑循环同构。
证明需要考虑这样一种刻画欧拉回路的方法：对于所有的以 \(k\) 为根节点的根向生成树，我们对所有结点的非树边出边钦定顺序，对 \(k\) 结点而言有 \((\deg_k - 1)!\) 种（减一是为了去除循环同构），其余为 \((\deg_u - 1)!\) 种。从点 \(k\) 出发，如果树边是唯一一条未经过的出边，那么走树边。否则，按照钦定的顺序走下一条未经过的出边。
如果我们能证明这样做一定能唯一地得到欧拉回路，且每一条欧拉回路都能被唯一对应，那么双射自然成立。
一个显然的事实是，在有向欧拉图上的任意不经过重复边的游走过程中，如果路径起点为 \(k\)，游走过程中行至 \(u\) 点，若 \(u = k\)，则此时每个点的入度与出度相等；若 \(u \neq k\)，则 \(k\) 点的入度为出度减一，\(u\) 点的出度为入度减一，其余点均满足入度出度相等。游走将在没有未经过的出边时终止。
那么第一个命题的证明就简单了。首先游走的过程只可能在 \(k\) 点终止，否则当前点一定满足出度为入度减一，即还有至少一条出边可行。而会不会有边没有被经过呢？由于树边一定是每个非 \(k\) 点的最后一条出边，若存在一条树边没有被经过，它会给父结点带来一个入度，因此父结点的树边也一定未被经过。这样归纳下去，\(k\) 结点将有一条入边未被经过，而由欧拉图可知，此时对应地有一条出边未被经过，矛盾，故所有边都将被经过，这种方法一定得到了欧拉回路。容易发现树边即为每个点最后一条经过的出边，所以这样得到的欧拉回路是互不相同的。
第二个命题的证明仍然考虑提出所有非 \(k\) 点的最后一条出边，这样形成的图一定是一棵根向树，其余边按照访问顺序对应即可。故每条欧拉回路都能被对应且有唯一的对应方式。
综合以上两个命题，BEST 定理成立。\(\square\)
扩展：如果认为重边是本质相同的，而不考虑循环同构时，只需将所求式除以 \(\prod A_{i, j}!\) 即可，其中 \(A_{i, j}\) 为 \(i \to j\) 的有向边数量。
既要去除重边，又要考虑循环同构的情况下，由于重边的缘故，从某个点出发的回路可能本质相同。假如我们将回路的异同关系用编号表示，则这个编号的可重集可以是多种的。因此我们很难确定每条欧拉回路被统计的次数，在这种情况下可能不具有通用的公式。
但事实上这种问题已经有了很成熟的解决方法：Polya 定理。鉴于这不是本文的主体因此只说明一下做法。
记 \(A_{i, j}\) 表示 \(i \to j\) 的有向边数量，\(t = \gcd(A_{i, j})\)。若 \(d \nmid t\)，定义 \(g(d) = 0\)，否则 \(g(d)\) 表示将所有重边的重复次数除以 \(d\) 形成的新图 \(G'\) 中，从任意点出发，不考虑循环同构，重边本质相同时的欧拉回路数：\(m' \times f ^ {\text{root}}(G', k) \times \prod (\deg'_u - 1)! \times \prod \dfrac{1}{A'_{i, j}!}\)。则所求等价于下式：
可以以 \(\Theta(t ^ {\frac{1}{3}} n ^ 3)\) 的时间复杂度计算。
接下来是一些题目。
P6178 【模板】Matrix-Tree 定理
模板题，行变换求行列式即可。时间复杂度 \(\Theta(n ^ 3)\)。
Code
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
constexpr int mod = 1e9 + 7, N = 300 + 5;
namespace basic {
inline int add(int x, int y) {return (x + y >= mod ? x + y - mod : x + y);}
inline int dec(int x, int y) {return (x - y < 0 ? x - y + mod : x - y);}
inline void ad(int &x, int y) {x = add(x, y);}
inline void de(int &x, int y) {x = dec(x, y);}
inline int qpow(int a, int b) {
int r = 1;
while(b) {
if(b & 1) r = 1ll * r * a % mod;
a = 1ll * a * a % mod; b >>= 1;
}
return r;
}
inline int inv(int x) {return qpow(x, mod - 2);}
int fac[N], ifac[N];
inline void fac_init(int n = N - 1) {
fac[0] = 1;
for(int i = 1; i <= n; i++)
fac[i] = 1ll * fac[i - 1] * i % mod;
ifac[n] = inv(fac[n]);
for(int i = n - 1; i >= 0; i--)
ifac[i] = 1ll * ifac[i + 1] * (i + 1) % mod;
}
}
using namespace basic;
struct Determinant {
int a[N][N], n;
Determinant() {memset(a, 0, sizeof(a));}
inline int* operator [] (int x) {return a[x];}
inline int det() {
int coef = 1;
for(int i = 1; i <= n; i++) {
int it = -1;
for(int j = i; j <= n; j++) {
if(a[j][i] != 0) {
it = j;
break;
}
}
if(it == -1) {
return 0;
}
if(it != i) {
swap(a[it], a[i]);
coef = mod - coef;
}
for(int j = i + 1; j <= n; j++) {
if(a[j][i] == 0) {
continue;
}
int C = 1ll * a[j][i] * inv(a[i][i]) % mod;
for(int k = i; k <= n; k++) {
de(a[j][k], 1ll * C * a[i][k] % mod);
}
}
}
int ret = coef;
for(int i = 1; i <= n; i++) {
ret = 1ll * ret * a[i][i] % mod;
}
return ret;
}
};
int n, m, type;
namespace Undirected {
int D[N][N], A[N][N];
void Main() {
for(int i = 1; i <= m; i++) {
int u, v, w; cin >> u >> v >> w;
if(u == v) {
continue;
}
ad(A[u][v], w), ad(A[v][u], w);
ad(D[u][u], w), ad(D[v][v], w);
}
Determinant L;
for(int i = 1; i <= n; i++) {
for(int j = 1; j <= n; j++) {
L[i][j] = dec(D[i][j], A[i][j]);
}
}
L.n = n - 1;
cout << L.det() << "\n";
}
}
namespace Directed {
int D_in[N][N], D_out[N][N], A[N][N];
void Main() {
for(int i = 1; i <= m; i++) {
int u, v, w; cin >> u >> v >> w;
if(u == v) {
continue;
}
ad(A[u][v], w);
ad(D_out[u][u], w), ad(D_in[v][v], w);
}
Determinant L;
for(int i = 2; i <= n; i++) {
for(int j = 2; j <= n; j++) {
L[i - 1][j - 1] = dec(D_in[i][j], A[i][j]);
}
}
L.n = n - 1;
cout << L.det() << "\n";
}
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr), cout.tie(nullptr);
cin >> n >> m >> type;
if(type == 0) {
Undirected::Main();
} else {
Directed::Main();
}
}
P5807 【模板】BEST 定理 | Which Dreamed It
模板题。注意是求从 \(1\) 出发的经过顺序不同的欧拉回路数，因此要用循环同构意义下不同的欧拉回路数乘以 \(\deg_1\)。本题还需判断图是否是欧拉图，等价于所有结点出入度相等，且含有 \(1\) 的极大强连通分量是否包含了所有非零度点。后者也即非零度点是否强连通，可以用 \(\Theta(n ^ 3)\) 的暴力实现。注意去除孤立点。单组时间复杂度 \(\Theta(n ^ 3)\)。
Code
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
constexpr int mod = 1e6 + 3, N = 100 + 5, M = 2e5 + 10;
namespace basic {
// ...
}
using namespace basic;
struct Determinant {
// ...
};
int D_in[N][N], D_out[N][N], A[N][N]; bool f[N][N];
void Main() {
memset(f, 0, sizeof(f));
memset(D_in, 0, sizeof(D_in)), memset(D_out, 0, sizeof(D_out));
memset(A, 0, sizeof(A));
int n; cin >> n;
for(int i = 1; i <= n; i++) {
f[i][i] = 1;
int E; cin >> E;
for(int j = 1; j <= E; j++) {
int v; cin >> v;
A[i][v]++, D_out[i][i]++, D_in[v][v]++;
f[i][v] = 1;
}
}
for(int k = 1; k <= n; k++) {
for(int i = 1; i <= n; i++) {
for(int j = 1; j <= n; j++) {
f[i][j] |= f[i][k] && f[k][j];
}
}
}
bool fl = 1;
for(int i = 1; i <= n; i++) {
fl &= D_in[i][i] == D_out[i][i];
for(int j = 1; j <= n; j++) {
if(D_in[i][i] && D_in[j][j]) {
fl &= f[i][j];
}
}
}
if(!fl) {
return cout << 0 << "\n", void();
}
vector<int> p;
for(int i = 1; i <= n; i++) {
if(D_out[i][i]) {
p.push_back(i);
}
}
if(p.size() == 0) {
return cout << 1 << "\n", void();
}
if(p[0] != 1) {
return cout << 0 << "\n", void();
}
int res = fac[D_out[p[0]][p[0]]];
for(int i = 1; i < p.size(); i++) {
res = 1ll * res * fac[D_out[p[i]][p[i]] - 1] % mod;
}
Determinant L; L.n = p.size() - 1;
for(int i = 1; i <= L.n; i++) {
for(int j = 1; j <= L.n; j++) {
L[i][j] = dec(D_out[p[i]][p[j]], A[p[i]][p[j]]);
}
}
res = 1ll * res * L.det() % mod;
cout << res << "\n";
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr), cout.tie(nullptr);
fac_init();
int T; cin >> T;
while(T--) {
Main();
}
}
*P6624 [省选联考 2020 A 卷] 作业题
这个关于 \(\gcd\) 的式子可以直接莫比乌斯反演掉，这里直接给出反演后的结果：记 \(f(k)\) 表示取所有 \(k \mid w\) 的边，其所能形成的所有生成树 \(T\) 的 \(\sum \limits_{e \in T} \omega(e)\) 之和。则所求内容可以表示为下式：
问题来到了如何计算所有的 \(f(k)\)。涉及到生成树相关计数，考虑 Matrix Tree 定理。但该定理要求 \(\omega(T) = \prod \limits_{e \in T} \omega(e)\)，此处我们遇到的是 \(\omega(T) = \sum \limits_{e \in T} \omega(e)\)，故尝试用乘法来刻画加法。
通常加乘转化可以用到 \(\exp, \ln\)，但在这里没什么用。为了用乘法刻画加法，一种更强大的工具是多项式。考虑两个多项式 \(1 + px, 1 + qx\)，它们的乘积为 \((1 + px)(1 + qx) \equiv 1 + (p + q)x \pmod{x ^ 2}\)。因此我们可以令一条边的权值多项式为 \(W_e = 1 + \omega(e)x\)，那么 \([x](\prod \limits_{e \in T} W_e)\) 即为 \(\sum \limits_{e \in T} \omega(e)\)。依此建立行列式并应用矩阵树定理即可。
如果对所有的 \(f(k), k \in [1, V]\) 计算答案，复杂度是 \(\Theta(n ^ 3 V)\) 的，较高。如果将边集大小小于 \(n - 1\) 的 \(k\) 跳过，合法的 \(k\) 将只有 \(\Theta(V ^ {\frac{1}{3}})\) 种，复杂度降低为 \(\Theta(n ^ 3 V ^ {\frac{1}{3}})\)，绰绰有余。
当然本题还有一个很大的注意点，这在官方数据里没有体现，但存在于 uoj 的 hack 数据中。我们知道多项式存在乘法逆当且仅当其常数项存在逆元，因此我们会选择当前列常数项非 \(0\) 的行，可以通过它完成行变换。如果不存在常数项非 \(0\) 的元素，也并不说明行列式的结果为 \(\mathbf{0}\)。正确的做法是，若该列存在至少一个元素的一次项系数非 \(0\)，则取之完成行变换，否则返回 \(\mathbf{0}\)。
Code
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
constexpr int mod = 998244353, N = 30 + 5, M = 2e5;
namespace basic {
// ...
}
using namespace basic;
int phi[M], vis[M];
vector<int> pr;
void sieve(int n = M - 1) {
// ...
}
struct Poly_2 {
int a, b;
inline friend Poly_2 operator + (Poly_2 x, Poly_2 y) {
return {add(x.a, y.a), add(x.b, y.b)};
}
inline friend Poly_2 operator - (Poly_2 x, Poly_2 y) {
return {dec(x.a, y.a), dec(x.b, y.b)};
}
inline friend Poly_2 operator * (Poly_2 x, Poly_2 y) {
return {1ll * x.a * y.a % mod, add(1ll * x.a * y.b % mod, 1ll * x.b * y.a % mod)};
}
inline Poly_2 Inv() {
assert(a);
int inva = inv(a);
return {inva, 1ll * dec(0, b) * inva % mod * inva % mod};
}
inline friend Poly_2 operator / (Poly_2 x, Poly_2 y) {
return x * y.Inv();
}
};
struct Determinant {
Poly_2 a[N][N]; int n;
Determinant() {
for(int i = 0; i < N; i++) {
for(int j = 0; j < N; j++) {
a[i][j] = {0, 0};
}
}
}
inline Poly_2* operator [] (int x) {return a[x];}
inline Poly_2 det() {
int coef = 1;
for(int i = 1; i <= n; i++) {
int it = -1;
for(int j = i; j <= n; j++) {
if(a[j][i].a) {
it = j;
break;
}
}
if(it == -1) {
for(int j = i; j <= n; j++) {
if(a[j][i].b) {
it = j;
break;
}
}
if(it == -1) {
return {0, 0};
}
if(it != i) {
swap(a[it], a[i]);
coef = mod - coef;
}
for(int j = i + 1; j <= n; j++) {
int C = 1ll * a[j][i].b * inv(a[i][i].b) % mod;
for(int k = i; k <= n; k++) {
de(a[j][k].b, 1ll * C * a[i][k].b % mod);
}
}
} else {
if(it != i) {
swap(a[it], a[i]);
coef = mod - coef;
}
for(int j = i + 1; j <= n; j++) {
Poly_2 C = a[j][i] / a[i][i];
for(int k = i; k <= n; k++) {
a[j][k] = a[j][k] - C * a[i][k];
}
}
}
}
Poly_2 ret = {coef, 0};
for(int i = 1; i <= n; i++) {
ret = ret * a[i][i];
}
return ret;
}
};
int n, m;
struct Edge {int u, v, w;} e[N * N];
Poly_2 A[N][N], D[N][N];
int f(vector<int> E) {
if(E.size() < n - 1) {
return 0;
}
for(int i = 1; i <= n; i++) {
for(int j = 1; j <= n; j++) {
A[i][j] = D[i][j] = {0, 0};
}
}
for(auto id : E) {
auto [u, v, w] = e[id];
Poly_2 W = {1, w};
A[u][v] = A[u][v] + W, A[v][u] = A[v][u] + W;
D[u][u] = D[u][u] + W, D[v][v] = D[v][v] + W;
}
Determinant L;
for(int i = 1; i <= n; i++) {
for(int j = 1; j <= n; j++) {
L[i][j] = D[i][j] - A[i][j];
}
}
L.n = n - 1;
return L.det().b;
}
vector<int> buc[M], rec[M];
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr), cout.tie(nullptr);
cin >> n >> m;
int V = 0;
for(int i = 1; i <= m; i++) {
auto &[u, v, w] = e[i];
cin >> u >> v >> w;
V = max(V, w);
}
sieve(V);
for(int i = 1; i <= V; i++) {
for(int j = i; j <= V; j += i) {
buc[j].push_back(i);
}
}
for(int i = 1; i <= m; i++) {
int w = e[i].w;
for(auto d : buc[w]) {
rec[d].push_back(i);
}
}
int ans = 0;
for(int i = 1; i <= V; i++) {
ad(ans, 1ll * phi[i] * f(rec[i]) % mod);
}
cout << ans << "\n";
}
P4111 [HEOI2015] 小 Z 的房间
这个就是练手题了。对题目建立图论模型，发现所求即为生成树数量，直接应用 Matrix Tree 定理。当然行列式求值的部分需要注意，因为这里 \(p\) 不是质数，因此要规避逆元，用类似辗转相除法的方式消元，可以先做做【模板】行列式求值。
这里给出行列式求值部分的实现，时间复杂度是 \(\Theta(n ^ 3 + n ^ 2\log p)\)：
Code
struct Determinant {
int a[M][M], n;
Determinant() {}
Determinant(int _) {n = _;}
inline int* operator [] (int x) {return a[x];}
inline int det() {
int coef = 1;
for(int i = 1; i <= n; i++) {
int it = -1;
for(int j = i; j <= n; j++) {
if(a[j][i] != 0) {
it = j;
break;
}
}
if(it == -1) {
return 0;
}
if(it != i) {
coef = mod - coef;
swap(a[it], a[i]);
}
for(int j = i + 1; j <= n; j++) {
while(1) {
if(a[j][i] > a[i][i]) {
swap(a[j], a[i]);
coef = mod - coef;
}
if(a[j][i] == 0) {
break;
}
int d = a[i][i] / a[j][i];
for(int k = i; k <= n; k++) {
de(a[i][k], 1ll * a[j][k] * d % mod);
}
}
}
}
int ret = coef;
for(int i = 1; i <= n; i++) {
ret = 1ll * ret * a[i][i] % mod;
}
return ret;
}
};
P2144 [FJOI2007] 轮状病毒
还是练手题。可以直接运用 Matrix Tree 定理建行列式，然后用高精度进行行变换，推荐使用上一题的辗转相消法，时间复杂度 \(\tilde {O}(n ^ 3)\)。当然这个行列式比较特殊，可以作更进一步的分析得到递推式，此处略。
Code
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
struct BigInt {
int sign;
std::vector<int> v;
BigInt() : sign(1) {}
BigInt(const std::string &s) { *this = s; }
BigInt(int v) {
char buf[21];
sprintf(buf, "%d", v);
*this = buf;
}
void zip(int unzip) {
if (unzip == 0) {
for (int i = 0; i < (int)v.size(); i++)
v[i] = get_pos(i * 4) + get_pos(i * 4 + 1) * 10 + get_pos(i * 4 + 2) * 100 + get_pos(i * 4 + 3) * 1000;
} else
for (int i = (v.resize(v.size() * 4), (int)v.size() - 1), a; i >= 0; i--)
a = (i % 4 >= 2) ? v[i / 4] / 100 : v[i / 4] % 100, v[i] = (i & 1) ? a / 10 : a % 10;
setsign(1, 1);
}
int get_pos(unsigned pos) const { return pos >= v.size() ? 0 : v[pos]; }
BigInt &setsign(int newsign, int rev) {
for (int i = (int)v.size() - 1; i > 0 && v[i] == 0; i--)
v.erase(v.begin() + i);
sign = (v.size() == 0 || (v.size() == 1 && v[0] == 0)) ? 1 : (rev ? newsign * sign : newsign);
return *this;
}
std::string val() const {
BigInt b = *this;
std::string s;
for (int i = (b.zip(1), 0); i < (int)b.v.size(); ++i)
s += char(*(b.v.rbegin() + i) + '0');
return (sign < 0 ? "-" : "") + (s.empty() ? std::string("0") : s);
}
bool absless(const BigInt &b) const {
if (v.size() != b.v.size()) return v.size() < b.v.size();
for (int i = (int)v.size() - 1; i >= 0; i--)
if (v[i] != b.v[i]) return v[i] < b.v[i];
return false;
}
BigInt operator-() const {
BigInt c = *this;
c.sign = (v.size() > 1 || v[0]) ? -c.sign : 1;
return c;
}
BigInt &operator=(const std::string &s) {
if (s[0] == '-')
*this = s.substr(1);
else {
for (int i = (v.clear(), 0); i < (int)s.size(); ++i)
v.push_back(*(s.rbegin() + i) - '0');
zip(0);
}
return setsign(s[0] == '-' ? -1 : 1, sign = 1);
}
bool operator<(const BigInt &b) const {
return sign != b.sign ? sign < b.sign : (sign == 1 ? absless(b) : b.absless(*this));
}
bool operator==(const BigInt &b) const { return v == b.v && sign == b.sign; }
BigInt &operator+=(const BigInt &b) {
if (sign != b.sign) return *this = (*this) - -b;
v.resize(std::max(v.size(), b.v.size()) + 1);
for (int i = 0, carry = 0; i < (int)b.v.size() || carry; i++) {
carry += v[i] + b.get_pos(i);
v[i] = carry % 10000, carry /= 10000;
}
return setsign(sign, 0);
}
BigInt operator+(const BigInt &b) const {
BigInt c = *this;
return c += b;
}
void add_mul(const BigInt &b, int mul) {
v.resize(std::max(v.size(), b.v.size()) + 2);
for (int i = 0, carry = 0; i < (int)b.v.size() || carry; i++) {
carry += v[i] + b.get_pos(i) * mul;
v[i] = carry % 10000, carry /= 10000;
}
}
BigInt operator-(const BigInt &b) const {
if (b.v.empty() || b.v.size() == 1 && b.v[0] == 0) return *this;
if (sign != b.sign) return (*this) + -b;
if (absless(b)) return -(b - *this);
BigInt c;
for (int i = 0, borrow = 0; i < (int)v.size(); i++) {
borrow += v[i] - b.get_pos(i);
c.v.push_back(borrow);
c.v.back() -= 10000 * (borrow >>= 31);
}
return c.setsign(sign, 0);
}
BigInt operator*(const BigInt &b) const {
if (b < *this) return b * *this;
BigInt c, d = b;
for (int i = 0; i < (int)v.size(); i++, d.v.insert(d.v.begin(), 0))
c.add_mul(d, v[i]);
return c.setsign(sign * b.sign, 0);
}
BigInt operator/(const BigInt &b) const {
BigInt c, d;
BigInt e = b;
e.sign = 1;
d.v.resize(v.size());
double db = 1.0 / (b.v.back() + (b.get_pos((unsigned)b.v.size() - 2) / 1e4) +
(b.get_pos((unsigned)b.v.size() - 3) + 1) / 1e8);
for (int i = (int)v.size() - 1; i >= 0; i--) {
c.v.insert(c.v.begin(), v[i]);
int m = (int)((c.get_pos((int)e.v.size()) * 10000 + c.get_pos((int)e.v.size() - 1)) * db);
c = c - e * m, c.setsign(c.sign, 0), d.v[i] += m;
while (!(c < e))
c = c - e, d.v[i] += 1;
}
return d.setsign(sign * b.sign, 0);
}
BigInt operator%(const BigInt &b) const { return *this - *this / b * b; }
bool operator>(const BigInt &b) const { return b < *this; }
bool operator<=(const BigInt &b) const { return !(b < *this); }
bool operator>=(const BigInt &b) const { return !(*this < b); }
bool operator!=(const BigInt &b) const { return !(*this == b); }
};
constexpr int N = 100 + 5;
struct Determinant {
BigInt a[N][N]; int n;
Determinant() {}
Determinant(int _) {n = _;}
inline BigInt* operator [] (int x) {return a[x];}
inline BigInt det() {
BigInt coef = 1;
for(int i = 1; i <= n; i++) {
int it = -1;
for(int j = i; j <= n; j++) {
if(a[j][i] != 0) {
it = j;
break;
}
}
if(it == -1) {
return 0;
}
if(it != i) {
coef = -coef;
swap(a[it], a[i]);
}
for(int j = i + 1; j <= n; j++) {
while(1) {
BigInt x = a[j][i], y = a[i][i];
if(x < 0) {
x = -x;
}
if(y < 0) {
y = -y;
}
if(x > y) {
swap(a[j], a[i]);
coef = -coef;
}
if(a[j][i] == 0) {
break;
}
BigInt d = a[i][i] / a[j][i];
for(int k = i; k <= n; k++) {
a[i][k] = a[i][k] - a[j][k] * d;
}
}
}
}
BigInt ret = coef;
for(int i = 1; i <= n; i++) {
ret = ret * a[i][i];
}
return ret;
}
};
int n, D[N][N], A[N][N];
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr), cout.tie(nullptr);
cin >> n;
if(n == 1) {
cout << 1 << "\n";
return 0;
}
for(int i = 1; i <= n; i++) {
D[i][i] = 3;
int l = i == 1 ? n : i - 1;
A[l][i]++, A[i][l]++;
A[n + 1][i]++, A[i][n + 1]++;
}
Determinant L(n);
for(int i = 1; i <= n; i++) {
for(int j = 1; j <= n; j++) {
L[i][j] = D[i][j] - A[i][j];
}
}
cout << L.det().val() << "\n";
}
*[ABC336G] 16 Integers
首先进行图论建模：对于 \((i, j, k, l)\)，从 \(v_{i, j, k}\) 向 \(v_{j, k, l}\) 建 \(X_{i, j, k, l}\) 条有向边，则此图的欧拉路径（回路）与原问题中的合法序列成双射。
首先判掉此图不是欧拉图或半欧拉图的情况，考虑此图是欧拉图的情况。欧拉回路计数，自然运用 BEST 定理。这里的欧拉回路显然可以从任一顶点开始，且重边是本质相同的，因此所求为 \(m \times f ^ {\text{root}}(G, k) \times \prod (\deg_u - 1)! \times \prod \frac{1}{A_{i, j}!}\)。
若此图是半欧拉图，我们要对欧拉路径进行计数。考虑从终点向起点连一条边，并钦定这条边一定是最后走的，则问题转化为对欧拉回路进行计数。容易发现这等价于 \(f ^ {\text{root}}(G, k) \times \prod (\deg_u - 1)! \times \prod \frac{1}{A_{i, j}!}\)，注意处理重边时不要将新加的那条边算在内，因为我们已经钦定它被最后经过。
至此我们以 \(\Theta(n ^ 3 + m)\) 的复杂度解决了本题，在本题中 \(n = 8\)。
一个注意点是自环的问题，其实我们不用刻意管它，可以发现它无论在 BEST 定理中还是 Matrix Tree 定理中都完全没有影响。
Code
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
constexpr int mod = 998244353, N = 1e6 + 10;
namespace basic {
// ...
}
using namespace basic;
constexpr int M = 8;
struct Determinant {
// ...
};
int e[2][2][2][2], m;
int D_in[M][M], D_out[M][M], A[M][M];
namespace Cycle {
void Main() {
int ans = m;
vector<int> p;
for(int i = 0; i < 8; i++) {
if(D_out[i][i]) {
p.push_back(i);
}
}
Determinant L(p.size() - 1);
for(int i = 0; i < L.n; i++) {
for(int j = 0; j < L.n; j++) {
L[i][j] = dec(D_out[p[i]][p[j]], A[p[i]][p[j]]);
}
}
ans = 1ll * ans * L.det() % mod;
for(int i = 0; i < p.size(); i++) {
ans = 1ll * ans * fac[D_out[p[i]][p[i]] - 1] % mod;
}
for(int i = 0; i < 8; i++) {
for(int j = 0; j < 8; j++) {
ans = 1ll * ans * ifac[A[i][j]] % mod;
}
}
cout << ans << "\n";
}
}
namespace Path {
void Main() {
int S = -1, T = -1;
for(int i = 0; i < 8; i++) {
if(D_out[i][i] == D_in[i][i] + 1) {
S = i;
}
if(D_in[i][i] == D_out[i][i] + 1) {
T = i;
}
}
int ans = 1;
for(int i = 0; i < 8; i++) {
for(int j = 0; j < 8; j++) {
ans = 1ll * ans * ifac[A[i][j]] % mod;
}
}
A[T][S]++, D_in[S][S]++, D_out[T][T]++;
vector<int> p;
for(int i = 0; i < 8; i++) {
if(D_out[i][i]) {
p.push_back(i);
}
}
Determinant L(p.size() - 1);
for(int i = 0; i < L.n; i++) {
for(int j = 0; j < L.n; j++) {
L[i][j] = dec(D_out[p[i]][p[j]], A[p[i]][p[j]]);
}
}
ans = 1ll * ans * L.det() % mod;
for(int i = 0; i < p.size(); i++) {
ans = 1ll * ans * fac[D_out[p[i]][p[i]] - 1] % mod;
}
cout << ans << "\n";
}
}
bool C1[M][M], C2[M][M];
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr), cout.tie(nullptr);
fac_init();
for(int i = 0; i < 2; i++) {
for(int j = 0; j < 2; j++) {
for(int k = 0; k < 2; k++) {
for(int l = 0; l < 2; l++) {
int u = (i << 2) | (j << 1) | k, v = (j << 2) | (k << 1) | l;
cin >> e[i][j][k][l];
if(e[i][j][k][l]) {
C1[u][v] = 1, C2[u][v] = C2[v][u] = 1;
}
D_in[v][v] += e[i][j][k][l], D_out[u][u] += e[i][j][k][l];
A[u][v] += e[i][j][k][l];
m += e[i][j][k][l];
}
}
}
}
for(int i = 0; i < 8; i++) {
C1[i][i] = C2[i][i] = 1;
}
for(int k = 0; k < 8; k++) {
for(int i = 0; i < 8; i++) {
for(int j = 0; j < 8; j++) {
C1[i][j] |= C1[i][k] && C1[k][j];
C2[i][j] |= C2[i][k] && C2[k][j];
}
}
}
bool fl = 1;
for(int i = 0; i < 8; i++) {
fl &= D_in[i][i] == D_out[i][i];
for(int j = 0; j < 8; j++) {
if(D_in[i][i] && D_in[j][j]) {
fl &= C1[i][j];
}
}
}
if(fl) {
Cycle::Main();
return 0;
}
fl = 1; int S = 0, T = 0;
for(int i = 0; i < 8; i++) {
if(abs(D_in[i][i] - D_out[i][i]) > 1) {
fl = 0;
}
S += D_out[i][i] == D_in[i][i] + 1;
T += D_in[i][i] == D_out[i][i] + 1;
for(int j = 0; j < 8; j++) {
if((D_in[i][i] + D_out[i][i]) && (D_in[j][j] + D_out[j][j])) {
fl &= C2[i][j];
}
}
}
if(fl && S == 1 && T == 1) {
Path::Main();
return 0;
}
cout << 0 << "\n";
}
*P5296 [北京省选集训2019] 生成树计数
前言：我做这个题时先入为主地想到了斯特林拆幂，其实二项式定理化成 EGF 的形式好像更简单！
发现是生成树相关内容的计数，直接考虑 Matrix Tree 定理。容易发现难点在于如何表示边权，使得边权之积包含了 \((\sum \omega(e)) ^ k\) 这个信息。
根据经典思路，先试试斯特林拆幂：\(n ^ k = \sum \limits_{i = 0} ^ k {k \brace i} \binom{n}{i} i!\)。则我们可以把问题进一步转化为，找到一种表示边权的方式，使得边权之积等于 \(\binom{\sum \omega(e)}{k}\)。
这等价于，我们需要通过关于 \(\binom{p}{k}\) 与 \(\binom{q}{k}\) 的信息，求积得到关于 \(\binom{p + q}{k}\) 的信息。这个 \(\binom{p + q}{k}\) 的形式容易让我们联想到范德蒙德卷积：\(\binom{p + q}{k} = \sum \limits_{i = 0} ^ {k} \binom{p}{i} \binom{q}{k - i}\)，至此思路已经一目了然：将边权定为 \(W(e) = \sum \limits_{i = 0} ^ k \binom{\omega(e)}{i} x ^ i\)，则 \([x ^ k]\prod \limits_{e \in T} W(e) = \binom{\sum \omega(e)}{k}\)。
运用 Matrix Tree 定理建立行列式后，通过行变换对行列式求值即可。将结果带回至斯特林拆幂的形式中即可得到答案。由于数据范围很小，暴力 \(\Theta(k ^ 2)\) 的卷积和求逆是比大常数的 \(\Theta(k \log k)\) 快的，因此采用暴力。时间复杂度 \(\Theta(k ^ 2 n ^ 3)\)。
行变换时需要用到多项式乘法逆，注意常数项为 \(0\) 的问题。
Code
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
constexpr int mod = 998244353, N = 30 + 5;
namespace basic {
// ...
}
using namespace basic;
struct TinyPoly {
vector<int> a;
TinyPoly() {}
inline void resize(int n) {
a.resize(n, 0);
}
TinyPoly(int n) {resize(n);}
TinyPoly(vector<int> b) {a = b;}
inline int & operator [] (int x) {
// assert(x < a.size());
return a[x];
}
inline int size() {
return a.size();
}
inline void Rev() {
reverse(a.begin(), a.end());
}
inline int fir() {
for(int i = 0; i < a.size(); i++) {
if(a[i]) {
return i;
}
}
return -1;
}
inline void popfront(int k) {
assert(k <= a.size());
reverse(a.begin(), a.end());
while(k--) {
a.pop_back();
}
reverse(a.begin(), a.end());
}
inline friend TinyPoly operator + (TinyPoly x, TinyPoly y) {
TinyPoly ret(max(x.size(), y.size()));
for(int i = 0; i < x.size() || i < y.size(); i++) {
if(i < x.size()) {
ad(ret[i], x[i]);
}
if(i < y.size()) {
ad(ret[i], y[i]);
}
}
return ret;
}
inline friend TinyPoly operator - (TinyPoly x, TinyPoly y) {
TinyPoly ret(max(x.size(), y.size()));
for(int i = 0; i < x.size() || i < y.size(); i++) {
if(i < x.size()) {
ad(ret[i], x[i]);
}
if(i < y.size()) {
de(ret[i], y[i]);
}
}
return ret;
}
inline friend TinyPoly operator * (TinyPoly x, TinyPoly y) {
if(x.size() == 0 || y.size() == 0) {
return {};
}
TinyPoly ret(x.size() + y.size() - 1);
for(int i = 0; i < x.size(); i++) {
for(int j = 0; j < y.size(); j++) {
ad(ret[i + j], 1ll * x[i] * y[j] % mod);
}
}
return ret;
}
inline TinyPoly Inv(int n) {
assert(a[0]);
TinyPoly ret(n);
ret[0] = inv(a[0]);
for(int i = 1; i < n; i++) {
int C = 0;
for(int j = 0; j < i; j++) {
de(C, 1ll * a[i - j] * ret[j] % mod);
}
ret[i] = 1ll * C * ret[0] % mod;
}
return ret;
}
inline TinyPoly operator += (TinyPoly y) {
return (*this) = (*this) + y;
}
inline TinyPoly operator -= (TinyPoly y) {
return (*this) = (*this) - y;
}
inline TinyPoly operator *= (TinyPoly y) {
return (*this) = (*this) * y;
}
};
int k;
struct Determinant {
int n; TinyPoly a[N][N];
Determinant() {};
Determinant(int _) {n = _;}
inline TinyPoly* operator [] (int x) {
return a[x];
}
inline TinyPoly det() {
TinyPoly ret(k + 1); ret[0] = 1;
for(int i = 1; i <= k; i++) {
ret[i] = 0;
}
for(int i = 1; i <= n; i++) {
int it = -1, fir = k + 1;
for(int j = i; j <= n; j++) {
int p = a[j][i].fir();
if(p != -1 && p < fir) {
fir = p, it = j;
}
}
if(it == -1) {
return ret[0] = 0, ret;
}
if(it != i) {
swap(a[i], a[it]);
ret[0] = mod - ret[0];
}
TinyPoly Inv = a[i][i]; Inv.popfront(fir), Inv.resize(k + 1);
Inv = Inv.Inv(k + 1);
for(int j = i + 1; j <= n; j++) {
TinyPoly p = a[j][i]; p.popfront(fir);
TinyPoly D = p * Inv; D.resize(k + 1);
for(int t = i; t <= n; t++) {
a[j][t] -= D * a[i][t];
a[j][t].resize(k + 1);
}
}
}
for(int i = 1; i <= n; i++) {
ret *= a[i][i];
ret.resize(k + 1);
}
return ret;
}
};
int n; TinyPoly D[N][N], A[N][N];
inline TinyPoly generate(int w) {
TinyPoly ret(k + 1); ret[0] = 1;
for(int i = 0; i < k; i++) {
ret[i + 1] = 1ll * ret[i] * dec(w, i) % mod * inv(i + 1) % mod;
}
return ret;
}
int S2[N][N];
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr), cout.tie(nullptr);
fac_init();
cin >> n >> k;
for(int i = 1; i <= n; i++) {
for(int j = 1; j <= n; j++) {
int w; cin >> w;
if(i < j) {
TinyPoly W = generate(w);
D[i][i] += W, D[j][j] += W;
A[i][j] += W, A[j][i] += W;
}
}
}
Determinant L(n - 1);
for(int i = 1; i <= n; i++) {
for(int j = 1; j <= n; j++) {
L[i][j] = D[i][j] - A[i][j];
}
}
TinyPoly res = L.det();
// for(int i = 0; i <= k; i++) {
// cout << res[i] << " \n"[i == k];
// }
S2[0][0] = 1;
for(int i = 1; i <= k; i++) {
for(int j = 1; j <= i; j++) {
S2[i][j] = add(S2[i - 1][j - 1], 1ll * S2[i - 1][j] * j % mod);
}
}
int ans = 0;
for(int i = 0; i <= k; i++) {
ad(ans, 1ll * S2[k][i] * res[i] % mod * fac[i] % mod);
}
cout << ans << "\n";
}
浙公网安备 33010602011771号