跳转至

圖論計數

在組合數學中,圖論計數(Graph Enumeration)是研究滿足特定性質的圖的計數問題的分支。生成函數波利亞計數定理符號化方法OEIS 是解決這類問題時最重要的數學工具。圖論計數可分為有標號和無標號兩大類問題,大多數情況下1有標號版本的問題都比其對應的無標號問題更加簡單,因此我們將先考察有標號問題的計數。

有標號樹

即 Cayley 公式,參見 Prüfer 序列 一文,我們也可以使用 Kirchhoff 矩陣樹定理生成函數拉格朗日定理 得到這一結果。

習題

有標號連通圖

例題「POJ 1737」Connected Graph

例題 「POJ 1737」Connected Graph

題目大意:求有 \(n\) 個結點的有標號連通圖的方案數(\(n \leq 50\))。

這類問題最早出現於樓教主的男人八題系列中,我們設 \(g_n\)\(n\) 個點有標號圖的方案數,\(c_n\) 為待求序列。\(n\) 個點的圖至多有 \(\binom{n}{2}\) 條邊,每條邊根據其出現與否有兩種狀態,每種狀態之間獨立,因而有 \(g_n = 2^{\binom{n}{2}}\)。我們固定其中一個節點,枚舉其所在連通塊的大小,那麼還需要從剩下的 \(n-1\) 個節點中選擇 \(i-1\) 個節點組成一個連通塊。連通塊之外的節點可以任意連邊,因而有如下遞推關係:

\[ \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} \]

移項得到 \(c_n\) 序列的 \(O(n^2)\) 遞推公式,可以通過此題。

例題「集訓隊作業 2013」城市規劃

例題 「集訓隊作業 2013」城市規劃

題目大意:求有 \(n\) 個結點的有標號連通圖的方案數(\(n \leq 130000\))。

對於數據範圍更大的序列問題,往往我們需要構造這些序列的生成函數,以使用高效的多項式算法。

方法一:分治 FFT

上述的遞推式可以看作一種自卷積形式,因而可以使用分治 FFT 進行計算,複雜度 \(O(n\log^2n)\)

方法二:多項式求逆

我們將上述遞推式中的組合數展開,並進行變形:

\[ \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} \]

構造多項式:

\[ \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 = H\),使用 多項式求逆 後再卷積解出 \(C(x)\) 即可。

方法三:多項式 exp

另一種做法是使用 EGF 中多項式 exp 的組合意義,我們設有標號連通圖和簡單圖序列的 EGF 分別為 \(C(x)\)\(G(x)\),那麼它們將有下列關係:

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

使用 多項式 ln 解出 \(C(x)\) 即可。

有標號歐拉圖、二分圖

例題「SPOJ KPGRAPHS」Counting Graphs

例題 「SPOJ KPGRAPHS」Counting Graphs

題目大意:求有 \(n\) 個結點的分別滿足下列性質的有標號圖的方案數(\(n \leq 1000\))。

本題限制代碼長度,因而無法直接使用多項式模板,但生成函數依然可以幫助我們進行分析。

連通圖問題在之前的例題中已被解決,考慮歐拉圖。注意到上述對連通圖計數的幾種方法,均可以在滿足任意性質的有標號連通圖進行推廣。例如我們可以將連通圖遞推公式中的 \(g_n\),從任意圖替換成滿足頂點度數均為偶數的圖,此時得到的 \(c_n\) 即為歐拉圖。

我們將 POJ 1737 的遞推過程封裝成連通化函數,

1
2
3
4
5
6
7
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];
  }
}

前兩問即可輕鬆解決:

1
2
3
4
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。

1
2
3
4
5
6
7
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];
  }
}

下面討論有標號二分圖計數,

我們設 \(b_n\) 表示 n 個結點的二分圖方案數,\(g_n\) 表示 \(n\) 個結點對結點進行 2 染色,滿足相同顏色的結點之間不存在邊的圖的方案數。枚舉其中一種顏色節點的數量,有2

\[ g_n = \sum_{i=0}^{n} \binom{n}{i}2^{i(n-i)} \]

接下來我們用兩種不同的方法建立 \(g_n\)\(b_n\) 之間的關係。

方法一:算兩次

我們設 \(c_{n, k}\) 表示有 k 個連通分量的二分圖方案數,那麼不難得到如下關係:

\[ \begin{align} b_n &= \sum_{i=1}^{n} c_{n, i} \\ g_n &= \sum_{i=1}^{n} c_{n, i} 2^i \end{align} \]

比較兩種 \(g_n\) 的表達式,展開得:

\[ \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} \]

不難得到 \(b_n\) 的遞推關係,複雜度 \(O(n^3)\),進一步使用容斥原理,可以優化到 \(O(n^2)\) 通過本題。

方法二:連通化遞推

方法二和方法三均使用連通二分圖 \(b1_n\) A001832 來建立 \(g_n\)\(b_n\) 之間的橋樑。

注意到對於每個連通二分圖,我們恰好有兩種不同的染色方法,對應到兩組不同的連通 2 染色圖, 因而對 \(g_n\) 進行連通化,得到的序列恰好是 \(b1_n\) 的兩倍,而 \(b_n\) 則由 \(b1_n\) 進行逆連通化得到。

因此:

1
2
3
4
5
6
7
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(n^2)\),可以通過本題。

方法三:多項式 exp

我們注意到也可以使用 EGF 理解上面的遞推過程。

\(G(x)\)\(g_n\) 的 EGF,\(B1(x)\)\(b1_n\) 的 EGF,\(B(x)\)\(b_n\) 的 EGF,應用做法二的方法,我們有:

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

我們可以對等式兩邊分別進行求導並比較兩邊係數,以得到易於編碼的遞推公式,通過此題。 注意到做法二與做法三本質相同,且一般情況下做法三可以得到更優的時間複雜度。

\[ \begin{align} B_n^2 &= G \\ 2B_nB_n' &= G' \end{align} \]
參考代碼
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
#include <bits/stdc++.h>
using namespace std;

#define Ts *this
#define rTs return Ts
typedef long long LL;
const 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)); }

inline 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>
inline 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>
inline 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);
    rTs;
  }

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

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

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

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

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

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

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

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

}  // namespace NT

using namespace NT;

const 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() {
#ifndef ONLINE_JUDGE
  // freopen("in.txt", "r", stdin);
#endif

  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--) {
    scanf("%d", &n);
    printf("Connected: %d\n", C[n]);
    printf("Eulerian: %d\n", E[n]);
    printf("Bipartite: %d\n", B[n]);
    puts("");
  }
}

習題

Riddell's Formula

上述關於 EGF 的 exp 的用法,有時又被稱作 Riddell's formula for labeled graphs,生成函數的 歐拉變換,有時也被稱為 Riddell's formula for unlabeled graphs,後者最早出現在歐拉對分拆數的研究中,除了解決圖論計數問題之外,也在完全揹包問題中出現。

對於給定序列 \(a_i\),和對應的 OGF \(A(x)\),定義 \(A(x)\) 的歐拉變換為:

\[ \begin{align} \mathcal{E}(A(x)) &= \prod_{i} (1-x^i)^{-a_i} \\ &= \exp (\sum_{i} \frac{A(x^i)}{i}) \end{align} \]

\(\mathcal{E}(A(x))\) 的各項係數為 \(b_i\),定義輔助數組 \(c_i = \sum_{d|n} d a_d\),則有遞推公式

\[ n 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) = x\mathcal{E}(F(x)) \]

取出係數即可。

無根樹

考慮容斥,我們用有根樹的方案中減去根不是重心的方案,並對 \(n\) 的奇偶性進行討論。

\(n\) 是奇數時:

必然存在一棵子樹大小 \(\geq \left\lceil \frac{n}{2}\right\rceil\),枚舉這棵子樹的大小有。

\[ g_n = f_n - \sum_{i=\left\lceil\frac{n}{2}\right\rceil}^{n-1} f_i f_{n-i} \]

\(n\) 是偶數時:

注意到當有兩個重心的情況時,上面的過程只會減去一次,因此還需要減去

\[ g_n = f_n - \sum_{i=\left\lceil\frac{n}{2}\right\rceil}^{n-1} f_i f_{n-i} - \binom{f_{\frac{n}{2}}}{2} \]

例題「Luogu P5900」無標號無根樹計數

例題 「Luogu P5900」無標號無根樹計數

題目大意:求有 n 個結點的無標號無根樹的方案數(\(n \leq 200000\))。

對於數據範圍更大的情況,做法同理,歐拉變換後使用多項式模板即可。

無標號簡單圖

例題「SGU 282. Isomorphism」Isomorphism

例題 「SGU 282. Isomorphism」Isomorphism

題目大意:求有 n 個結點的無標號完全圖的邊進行 m 染色的方案數。

注意到當 m = 2 時,所求對象就是無標號簡單圖 A000088,考察波利亞計數定理,

\[ \frac{1}{|G|}\sum_{g\in G} m^{c(g)} \]

本題中置換羣 \(G\) 為頂點的 \(n\) 階對稱羣生成的邊集置換羣,但暴力做法的枚舉量為 \(O(n!)\),無法通過此題。

考慮根據按照置換的循環結構進行分類,每種循環結構對應一種數的分拆,我們用 dfs() 生成分拆,那麼問題即轉化為求每一種分拆 \(p\) 所對應的置換數目 \(w(p)\) 和每一類置換中的循環個數 \(c(p)\),答案為

\[ \frac{1}{|G|} \sum_{p \in P} w(p) m^{c(p)} \]

考慮 \(w(p)\),每一個分拆對應一個循環排列,同時同一種大小的分拆之間的順序無關,因而我們有:

\[ w(p) = \frac{n!}{\prod_{i}(p_i)\prod_{i}(q_i!)} \]

這裏 \(q_i\) 表示大小為 \(i\) 的分拆在 \(p\) 中出現的次數。

考慮 \(c(p)\)\(p\) 所影響的點集的循環即為 \(|p|\),但題目考察的是邊染色,所以還需要考察點置換所生成的邊置換,

如果一條邊關聯的頂點處在同一個循環內,設該循環大小為 \(p_i\),那麼邊所生成的循環數恰好為 \(\left\lfloor \frac{p_i}{2} \right\rfloor\)

如果一條邊關聯的頂點處在兩個不同的循環中,設分別為 \(p_i\),\(p_j\),每個循環節的長度均為 \(\operatorname{lcm}(p_i,p_j)\),因而邊所生成的循環數恰好為 \(\frac{p_i p_j}{\operatorname{lcm}(p_i,p_j)} = \gcd(p_i, p_j)\)

參考代碼
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
#include <bits/stdc++.h>
using namespace std;

#define Ts *this
#define rTs return Ts
typedef long long LL;
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)); }

inline 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>
inline 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>
inline 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);
    rTs;
  }

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

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

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

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

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

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

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

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

}  // namespace NT

using namespace NT;

const 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 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

  1. 也許無標號二叉樹是一個反例,在結構簡單的情況下,對應的置換羣是恆等羣(Identity Group),此時有標號版本可以直接通過乘以 \(n!\) 得到。 

  2. 粉兔的 blog 告訴我們,這個序列也可以使用 Chirp Z-Transform 優化。