图论计数

在组合数学中，图论计数（Graph Enumeration）是研究满足特定性质的图的计数问题的分支．生成函数、波利亚计数定理 与 符号化方法 和 OEIS 是解决这类问题时最重要的数学工具．图论计数可分为有标号和无标号两大类问题，大多数情况下1有标号版本的问题都比其对应的无标号问题更加简单，因此我们将先考察有标号问题的计数．

有标号树

即 Cayley 公式，参见 Prüfer 序列 一文，我们也可以使用 Kirchhoff 矩阵树定理 或 生成函数 和 拉格朗日定理 得到这一结果．

习题

有标号连通图

例题「POJ 1737」Connected Graph

例题 「POJ 1737」Connected Graph

题目大意：求有 nn 个结点的有标号连通图的方案数（n≤50n \leq 50）．

这类问题最早出现于楼教主的男人八题系列中，我们设 gng_n 为 nn 个点有标号图的方案数，cnc_n 为待求序列．nn 个点的图至多有 (n2)\binom{n}{2} 条边，每条边根据其出现与否有两种状态，每种状态之间独立，因而有 gn=2(n2)g_n = 2^{\binom{n}{2}}．我们固定其中一个节点，枚举其所在连通块的大小，那么还需要从剩下的 n−1n-1 个节点中选择 i−1i-1 个节点组成一个连通块．连通块之外的节点可以任意连边，因而有如下递推关系：

∑i=1n(n−1i−1)cign−i=gncn=gn−∑i=1n−1(n−1i−1)cign−i\begin{align} \sum_{i=1}^{n} \binom{n-1}{i-1} c_i g_{n-i} &= g_n \\ c_n &= g_n - \sum_{i=1}^{n-1} \binom{n-1}{i-1} c_i g_{n-i} \end{align}

移项得到 cnc_n 序列的 O(n2)O(n^2) 递推公式，可以通过此题．

例题「集训队作业 2013」城市规划

例题 「集训队作业 2013」城市规划

题目大意：求有 nn 个结点的有标号连通图的方案数（n≤130000n \leq 130000）．

对于数据范围更大的序列问题，往往我们需要构造这些序列的生成函数，以使用高效的多项式算法．

方法一：分治 FFT

上述的递推式可以看作一种自卷积形式，因而可以使用分治 FFT 进行计算，复杂度 O(nlog⁡2n)O(n\log^2n)．

方法二：多项式求逆

我们将上述递推式中的组合数展开，并进行变形：

∑i=1n(n−1i−1)cign−i=gn∑i=1nci(i−1)!gn−i(n−i)!=gn(n−1)!\begin{align} \sum_{i=1}^{n} \binom{n-1}{i-1} c_i g_{n-i} &= g_n \\ \sum_{i=1}^{n} \frac{c_i}{(i-1)!} \frac{g_{n-i}}{(n-i)!} &= \frac{g_n}{(n-1)!} \end{align}

构造多项式：

C(x)=∑n=1cn(n−1)!xnG(x)=∑n=0gnn!xnH(x)=∑n=1gn(n−1)!xn\begin{align} C(x) &= \sum_{n=1} \frac{c_n}{(n-1)!} x^n \\ G(x) &= \sum_{n=0} \frac{g_n}{n!} x^n \\ H(x) &= \sum_{n=1} \frac{g_n}{(n-1)!} x^n \end{align}

代换进上式得到 CG=HCG = H，使用 多项式求逆 后再卷积解出 C(x)C(x) 即可．

方法三：多项式 exp

另一种做法是使用 EGF 中多项式 exp 的组合意义，我们设有标号连通图和简单图序列的 EGF 分别为 C(x)C(x) 和 G(x)G(x)，那么它们将有下列关系：

exp⁡(C(x))=G(x)C(x)=ln⁡(G(x))\begin{align} \exp(C(x)) &= G(x) \\ C(x) &= \ln(G(x)) \end{align}

使用 多项式 ln 解出 C(x)C(x) 即可．

有标号欧拉图、二分图

例题「SPOJ KPGRAPHS」Counting Graphs

例题 「SPOJ KPGRAPHS」Counting Graphs

题目大意：求有 nn 个结点的分别满足下列性质的有标号图的方案数（n≤1000n \leq 1000）．

本题限制代码长度，因而无法直接使用多项式模板，但生成函数依然可以帮助我们进行分析．

连通图问题在之前的例题中已被解决，考虑欧拉图．注意到上述对连通图计数的几种方法，均可以在满足任意性质的有标号连通图进行推广．例如我们可以将连通图递推公式中的 gng_n，从任意图替换成满足顶点度数均为偶数的图，此时得到的 cnc_n 即为欧拉图．

我们将 POJ 1737 的递推过程封装成连通化函数，

void ln(Int C[], Int G[]) {
  for (int i = 1; i <= n; ++i) {
    C[i] = G[i];
    for (int j = 1; j <= i - 1; ++j)
      C[i] -= binom[i - 1][j - 1] * C[j] * G[i - j];
  }
}

前两问即可轻松解决：

for (int i = 1; i <= n; ++i) G[i] = pow(2, binom[i][2]);
ln(C, G);
for (int i = 1; i <= n; ++i) G[i] = pow(2, binom[i - 1][2]);
ln(E, G);

注意到这里的连通化递推过程其实等价于对其 EGF 求多项式 ln，同理我们也可以写出逆连通化函数，它等价于对其 EGF 求多项式 exp．

void exp(Int G[], Int C[]) {
  for (int i = 1; i <= n; ++i) {
    G[i] = C[i];
    for (int j = 1; j <= i - 1; ++j)
      G[i] += binom[i - 1][j - 1] * C[j] * G[i - j];
  }
}

下面讨论有标号二分图计数，

我们设 bnb_n 表示 n 个结点的二分图方案数，gng_n 表示 nn 个结点对结点进行 2 染色，满足相同颜色的结点之间不存在边的图的方案数．枚举其中一种颜色节点的数量，有2：

gn=∑i=0n(ni)2i(n−i)g_n = \sum_{i=0}^{n} \binom{n}{i}2^{i(n-i)}

接下来我们用两种不同的方法建立 gng_n 与 bnb_n 之间的关系．

方法一：算两次

我们设 cn,kc_{n, k} 表示有 k 个连通分量的二分图方案数，那么不难得到如下关系：

bn=∑i=1ncn,ign=∑i=1ncn,i2i\begin{align} b_n &= \sum_{i=1}^{n} c_{n, i} \\ g_n &= \sum_{i=1}^{n} c_{n, i} 2^i \end{align}

比较两种 gng_n 的表达式，展开得：

∑i=0n(ni)2i(n−i)=∑i=1ncn,i2icn,i=∑i=0n−1(n−1i−1)cn,1cn−i,k−1\begin{align} \sum_{i=0}^{n} \binom{n}{i}2^{i(n-i)} &= \sum_{i=1}^{n} c_{n, i} 2^i \\ c_{n, i} &= \sum_{i=0}{n-1} \binom{n-1}{i-1} c_{n, 1}c_{n-i,k-1} \end{align}

不难得到 bnb_n 的递推关系，复杂度 O(n3)O(n^3)，进一步使用容斥原理，可以优化到 O(n2)O(n^2) 通过本题．

方法二：连通化递推

方法二和方法三均使用连通二分图 b1nb1_n A001832 来建立 gng_n 与 bnb_n 之间的桥梁．

注意到对于每个连通二分图，我们恰好有两种不同的染色方法，对应到两组不同的连通 2 染色图， 因而对 gng_n 进行连通化，得到的序列恰好是 b1nb1_n 的两倍，而 bnb_n 则由 b1nb1_n 进行逆连通化得到．

因此：

for (int i = 1; i <= n; ++i) {
  G[i] = 0;
  for (int j = 0; j < i + 1; ++j) G[i] += binom[i][j] * pow(2, j * (i - j));
}
ln(B1, G);
for (int i = 1; i <= n; ++i) B1[i] /= 2;
exp(B, B1);

两种递推的过程复杂度均为 O(n2)O(n^2)，可以通过本题．

方法三：多项式 exp

我们注意到也可以使用 EGF 理解上面的递推过程．

设 G(x)G(x) 为 gng_n 的 EGF，B1(x)B1(x) 为 b1nb1_n 的 EGF，B(x)B(x) 为 bnb_n 的 EGF，应用做法二的方法，我们有：

G(x)=exp⁡(2B1(x))B(x)=exp⁡(B1(x))=exp⁡(ln⁡G(x)2)=G\begin{align} G(x) &= \exp(2B1(x)) \\ B(x) &= \exp(B1(x)) \\ &= \exp(\frac{\ln{G(x)}}{2}) \\ &= \sqrt{G} \end{align}

我们可以对等式两边分别进行求导并比较两边系数，以得到易于编码的递推公式，通过此题． 注意到做法二与做法三本质相同，且一般情况下做法三可以得到更优的时间复杂度．

Bn2=G2BnBn′=G′\begin{align} B_n^2 &= G \\ 2B_nB_n' &= G' \end{align}
参考代码 
#include <iostream>
using namespace std;

using LL = long long;
constexpr int MOD = int(1e9) + 7;

// <<= '2. Number Theory .,//{
namespace NT {
void INC(int& a, int b) {
  a += b;
  if (a >= MOD) a -= MOD;
}

int sum(int a, int b) {
  a += b;
  if (a >= MOD) a -= MOD;
  return a;
}

void DEC(int& a, int b) {
  a -= b;
  if (a < 0) a += MOD;
}

int dff(int a, int b) {
  a -= b;
  if (a < 0) a += MOD;
  return a;
}

void MUL(int& a, int b) { a = (LL)a * b % MOD; }

int pdt(int a, int b) { return (LL)a * b % MOD; }

int _I(int b) {
  int a = MOD, x1 = 0, x2 = 1, q;
  while (1) {
    q = a / b, a %= b;
    if (!a) return x2;
    DEC(x1, pdt(q, x2));

    q = b / a, b %= a;
    if (!b) return x1;
    DEC(x2, pdt(q, x1));
  }
}

void DIV(int& a, int b) { MUL(a, _I(b)); }

int qtt(int a, int b) { return pdt(a, _I(b)); }

int pow(int a, LL b) {
  int c(1);
  while (b) {
    if (b & 1) MUL(c, a);
    MUL(a, a), b >>= 1;
  }
  return c;
}

template <class T>
T pow(T a, LL b) {
  T c(1);
  while (b) {
    if (b & 1) c *= a;
    a *= a, b >>= 1;
  }
  return c;
}

template <class T>
T pow(T a, int b) {
  return pow(a, (LL)b);
}

struct Int {
  int val;

  operator int() const { return val; }

  Int(int _val = 0) : val(_val) {
    val %= MOD;
    if (val < 0) val += MOD;
  }

  Int(LL _val) : val(_val) {
    _val %= MOD;
    if (_val < 0) _val += MOD;
    val = _val;
  }

  Int& operator+=(const int& rhs) {
    INC(val, rhs);
    return *this;
  }

  Int operator+(const int& rhs) const { return sum(val, rhs); }

  Int& operator-=(const int& rhs) {
    DEC(val, rhs);
    return *this;
  }

  Int operator-(const int& rhs) const { return dff(val, rhs); }

  Int& operator*=(const int& rhs) {
    MUL(val, rhs);
    return *this;
  }

  Int operator*(const int& rhs) const { return pdt(val, rhs); }

  Int& operator/=(const int& rhs) {
    DIV(val, rhs);
    return *this;
  }

  Int operator/(const int& rhs) const { return qtt(val, rhs); }

  Int operator-() const { return MOD - *this; }
};

}  // namespace NT

using namespace NT;

constexpr int N = int(1e3) + 9;
Int binom[N][N], C[N], E[N], B[N], B1[N], G[N];
int n;

void ln(Int C[], Int G[]) {
  for (int i = 1; i <= n; ++i) {
    C[i] = G[i];
    for (int j = 1; j <= i - 1; ++j)
      C[i] -= binom[i - 1][j - 1] * C[j] * G[i - j];
  }
}

void exp(Int G[], Int C[]) {
  for (int i = 1; i <= n; ++i) {
    G[i] = C[i];
    for (int j = 1; j <= i - 1; ++j)
      G[i] += binom[i - 1][j - 1] * C[j] * G[i - j];
  }
}

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);

  n = 1000;
  for (int i = 0; i < n + 1; ++i) {
    binom[i][0] = 1;
    for (int j = 0; j < i; ++j)
      binom[i][j + 1] = binom[i - 1][j] + binom[i - 1][j + 1];
  }

  for (int i = 1; i <= n; ++i) G[i] = pow(2, binom[i][2]);
  ln(C, G);
  for (int i = 1; i <= n; ++i) G[i] = pow(2, binom[i - 1][2]);
  ln(E, G);
  for (int i = 1; i <= n; ++i) {
    G[i] = 0;
    for (int j = 0; j < i + 1; ++j) G[i] += binom[i][j] * pow(2, j * (i - j));
  }
  ln(B1, G);
  for (int i = 1; i <= n; ++i) B1[i] /= 2;
  exp(B, B1);

  int T;
  cin >> T;
  while (T--) {
    cin >> n;
    cout << "Connected: " << C[n] << '\n'
         << "Eulerian: " << E[n] << '\n'
         << "Bipartite: " << B[n] << "\n\n";
  }
}

习题

Riddell's Formula

上述关于 EGF 的 exp 的用法，有时又被称作 Riddell's formula for labeled graphs，生成函数的 欧拉变换，有时也被称为 Riddell's formula for unlabeled graphs，后者最早出现在欧拉对分拆数的研究中，除了解决图论计数问题之外，也在完全背包问题中出现．

对于给定序列 aia_i，和对应的 OGF A(x)A(x)，定义 A(x)A(x) 的欧拉变换为：

E(A(x))=∏i(1−xi)−ai=exp⁡(∑iA(xi)i)\begin{align} \mathcal{E}(A(x)) &= \prod_{i} (1-x^i)^{-a_i} \\ &= \exp (\sum_{i} \frac{A(x^i)}{i}) \end{align}

设 E(A(x))\mathcal{E}(A(x)) 的各项系数为 bib_i，定义辅助数组 ci=∑d∣ndadc_i = \sum_{d|n} d a_d，则有递推公式

nbn=cn+∑i=1n−1cibn−in b_n = c_n + \sum_{i=1}^{n-1} c_i b_{n-i}

无标号树

例题「SPOJ PT07D」Let us count 1 2 3

例题 「SPOJ PT07D」Let us count 1 2 3

题目大意：求有 n 个结点的分别满足下列性质的树的方案数．

有根树

有标号情况以在前文中解决，下面考察无标号有根树，设其 OGF 为 F(x)F(x)，应用欧拉变换，可得：

F(x)=xE(F(x))F(x) = x\mathcal{E}(F(x))

取出系数即可．

无根树

考虑容斥，我们用有根树的方案中减去根不是重心的方案，并对 nn 的奇偶性进行讨论．

当 nn 是奇数时：

必然存在一棵子树大小 ≥⌈n2⌉\geq \left\lceil \frac{n}{2}\right\rceil，枚举这棵子树的大小有．

gn=fn−∑i=⌈n2⌉n−1fifn−ig_n = f_n - \sum_{i=\left\lceil\frac{n}{2}\right\rceil}^{n-1} f_i f_{n-i}

当 nn 是偶数时：

注意到当有两个重心的情况时，上面的过程只会减去一次，因此还需要减去

gn=fn−∑i=n2+1n−1fifn−i−(fn22)g_n = f_n - \sum_{i=\frac{n}{2}+1}^{n-1} f_i f_{n-i} - \binom{f_{\frac{n}{2}}}{2}

例题「Luogu P5900」无标号无根树计数

例题 「Luogu P5900」无标号无根树计数

题目大意：求有 n 个结点的无标号无根树的方案数（n≤200000n \leq 200000）．

对于数据范围更大的情况，做法同理，欧拉变换后使用多项式模板即可．

无标号简单图

例题「SGU 282. Isomorphism」Isomorphism

例题 「SGU 282. Isomorphism」Isomorphism

题目大意：求有 n 个结点的无标号完全图的边进行 m 染色的方案数．

注意到当 m = 2 时，所求对象就是无标号简单图 A000088，考察波利亚计数定理，

1∣G∣∑g∈Gmc(g)\frac{1}{|G|}\sum_{g\in G} m^{c(g)}

本题中置换群 GG 为顶点的 nn 阶对称群生成的边集置换群，但暴力做法的枚举量为 O(n!)O(n!)，无法通过此题．

考虑根据按照置换的循环结构进行分类，每种循环结构对应一种数的分拆，我们用 dfs() 生成分拆，那么问题即转化为求每一种分拆 pp 所对应的置换数目 w(p)w(p) 和每一类置换中的循环个数 c(p)c(p)，答案为

1∣G∣∑p∈Pw(p)mc(p)\frac{1}{|G|} \sum_{p \in P} w(p) m^{c(p)}

考虑 w(p)w(p)，每一个分拆对应一个循环排列，同时同一种大小的分拆之间的顺序无关，因而我们有：

w(p)=n!∏i(pi)∏i(qi!)w(p) = \frac{n!}{\prod_{i}(p_i)\prod_{i}(q_i!)}

这里 qiq_i 表示大小为 ii 的分拆在 pp 中出现的次数．

考虑 c(p)c(p)，pp 所影响的点集的循环即为 ∣p∣|p|，但题目考察的是边染色，所以还需要考察点置换所生成的边置换，

如果一条边关联的顶点处在同一个循环内，设该循环大小为 pip_i，那么边所生成的循环数恰好为 ⌊pi2⌋\left\lfloor \frac{p_i}{2} \right\rfloor．

如果一条边关联的顶点处在两个不同的循环中，设分别为 pip_i,pjp_j，每个循环节的长度均为 lcm⁡(pi,pj)\operatorname{lcm}(p_i,p_j)，因而边所生成的循环数恰好为 pipjlcm⁡(pi,pj)=gcd⁡(pi,pj)\frac{p_i p_j}{\operatorname{lcm}(p_i,p_j)} = \gcd(p_i, p_j)．

参考代码 
#include <iostream>
#include <vector>
using namespace std;

using LL = long long;
int MOD = int(1e9) + 7;

namespace NT {
void INC(int& a, int b) {
  a += b;
  if (a >= MOD) a -= MOD;
}

int sum(int a, int b) {
  a += b;
  if (a >= MOD) a -= MOD;
  return a;
}

void DEC(int& a, int b) {
  a -= b;
  if (a < 0) a += MOD;
}

int dff(int a, int b) {
  a -= b;
  if (a < 0) a += MOD;
  return a;
}

void MUL(int& a, int b) { a = (LL)a * b % MOD; }

int pdt(int a, int b) { return (LL)a * b % MOD; }

int _I(int b) {
  int a = MOD, x1 = 0, x2 = 1, q;
  while (1) {
    q = a / b, a %= b;
    if (!a) return x2;
    DEC(x1, pdt(q, x2));

    q = b / a, b %= a;
    if (!b) return x1;
    DEC(x2, pdt(q, x1));
  }
}

void DIV(int& a, int b) { MUL(a, _I(b)); }

int qtt(int a, int b) { return pdt(a, _I(b)); }

int pow(int a, LL b) {
  int c(1);
  while (b) {
    if (b & 1) MUL(c, a);
    MUL(a, a), b >>= 1;
  }
  return c;
}

template <class T>
T pow(T a, LL b) {
  T c(1);
  while (b) {
    if (b & 1) c *= a;
    a *= a, b >>= 1;
  }
  return c;
}

template <class T>
T pow(T a, int b) {
  return pow(a, (LL)b);
}

struct Int {
  int val;

  operator int() const { return val; }

  Int(int _val = 0) : val(_val) {
    val %= MOD;
    if (val < 0) val += MOD;
  }

  Int(LL _val) : val(_val) {
    _val %= MOD;
    if (_val < 0) _val += MOD;
    val = _val;
  }

  Int& operator+=(const int& rhs) {
    INC(val, rhs);
    return *this;
  }

  Int operator+(const int& rhs) const { return sum(val, rhs); }

  Int& operator-=(const int& rhs) {
    DEC(val, rhs);
    return *this;
  }

  Int operator-(const int& rhs) const { return dff(val, rhs); }

  Int& operator*=(const int& rhs) {
    MUL(val, rhs);
    return *this;
  }

  Int operator*(const int& rhs) const { return pdt(val, rhs); }

  Int& operator/=(const int& rhs) {
    DIV(val, rhs);
    return *this;
  }

  Int operator/(const int& rhs) const { return qtt(val, rhs); }

  Int operator-() const { return MOD - *this; }
};

}  // namespace NT

using namespace NT;

constexpr int N = int(5e1) + 9;
Int Fact[N];
vector<vector<int>> Partition;
vector<int> cur;
int n, m;

void gen(int n = ::n, int s = 1) {
  if (!n) {
    Partition.push_back(cur);
  } else if (n >= s) {
    cur.push_back(s);
    gen(n - s, s);
    cur.pop_back();
    gen(n, s + 1);
  }
}

Int w(const vector<int> P) {
  Int z = Fact[n];
  int c = 0, l = P.front();

  for (auto p : P) {
    z /= p;
    if (p != l) {
      z /= Fact[c];
      l = p;
      c = 1;
    } else {
      ++c;
    }
  }

  z /= Fact[c];
  return z;
}

int gcd(int x, int y) { return y ? gcd(y, x % y) : x; }

int c(const vector<int> P) {
  int z = 0;
  for (int i = 0; i < P.size(); ++i) {
    z += P[i] / 2;
    for (int j = 0; j < i; ++j) z += gcd(P[i], P[j]);
  }
  return z;
}

int main() {
  cin >> n >> m >> MOD;
  Fact[0] = 1;
  for (int i = 1; i <= n; ++i) Fact[i] = Fact[i - 1] * i;

  gen();

  Int res = 0;
  for (auto P : Partition) {
    res += w(P) * pow(m, c(P));
  }
  res /= Fact[n];
  cout << res << endl;
}

习题

参考资料与注释

  1. WC2015, 顾昱洲营员交流资料 Graphical Enumeration
  2. WC2019, 生成函数，多项式算法与图的计数
  3. Counting labeled graphs - Algorithms for Competitive Programming
  4. Graphical Enumeration Paperback, Frank Harary, Edgar M. Palmer
  5. The encyclopedia of integer sequences, N. J. A. Sloane, Simon Plouffe
  6. Combinatorial Problems and Exercises, László Lovász
  7. Graph Theory and Additive Combinatorics

Footnotes

  1. 也许无标号二叉树是一个反例，在结构简单的情况下，对应的置换群是恒等群（Identity Group），此时有标号版本可以直接通过乘以 n!n! 得到． ↩

  2. 粉兔的 blog 告诉我们，这个序列也可以使用 Chirp Z-Transform 优化． ↩

上一节
236 Pólya 计数
下一节
导读 线性代数简介

本文内容来自 OI Wiki《图论计数》；作者与修订历史见原文。AIOJ 对排版、站内链接和资源路径做了格式适配。

内容依 CC BY-SA 4.0 与 SATA 许可共享；附加条款以原文声明为准。感谢 OI Wiki 作者与贡献者；请按 SATA 约定前往 OI Wiki 仓库 Star 支持。