二项式反演学习笔记

アズマ比起丘比特更像魅魔吧（笑

最近算法溜大了，来点纯数学去去火。

---

概述

二项式反演常用来将“至少/至多 $k$ 个”转为“恰好 $k$ 个”，可以解决一类计数问题。

常用形式有以下两种：

形式一

$f(i)$ 表示恰好 $i$ 个的数量，$g(i)$ 表示至多 $i$ 个的数量，则二者满足以下关系：
$$
g(k)=\sum\limits_{i=0}^{k}\binom{k}{i}f(i) \Longleftrightarrow f(k)=\sum\limits_{i=0}^{k}(-1)^{k-i}\binom{k}{i}g(i)
$$

形式二

$f(i)$ 表示恰好 $i$ 个的数量，$g(i)$ 表示至少 $i$ 个的数量，则二者满足以下关系：
$$
g(k)=\sum\limits_{i=k}^{n}\binom{i}{k}f(i) \Longleftrightarrow f(k)=\sum\limits_{i=k}^{n}(-1)^{i-k}\binom{i}{k}g(i)
$$
这启示我们，如果知道其中一个序列，就可以在 $O(n)$ 的时间内计算另一个序列的某一项。

通常情况下，“至少/至多”比较好算，我们可以用这个信息来计算“恰好”的值。

读者可能会发现二项式反演形式上与广义容斥原理十分相似。本质上，广义容斥原理就是二项式反演在集合计数领域的一个具体应用，而二项式反演是一个更通用的数学工具，在生成函数等领域也有其他应用。

证明

既然是纯数学，有完整的证明才舒服。

两种形式的证明方法很类似，以下仅给出形式一的证明。
$$
g(k)=\sum\limits_{i=0}^{k}\binom{k}{i}f(i) \Longleftrightarrow f(k)=\sum\limits_{i=0}^{k}(-1)^{k-i}\binom{k}{i}g(i)
$$

我们假设右边公式成立，将 $g(i)$ 的定义代入：
$$
f(k)=\sum\limits_{i=0}^{k}(-1)^{k-i}\binom{k}{i}\sum\limits_{j=0}^{i}\binom{i}{j}f(j)
$$
交换求和顺序：
$$
f(k)=\sum\limits_{j=0}^{k}f(j)\sum\limits_{i=j}^{k}(-1)^{k-i}\binom{k}{i}\binom{i}{j}
$$
由恒等式 $\binom{k}{i}\binom{i}{j}=\binom{k}{j}\binom{k-j}{k-i}$ 得：
$$
f(k)=\sum\limits_{j=0}^{k}\binom{k}{j}f(j)\sum\limits_{i=j}^{k}(-1)^{k-i}\binom{k-j}{k-i}
$$
由二项式定理，注意到
$$
\sum\limits_{i=j}^{k}(-1)^{k-i}\binom{k-j}{k-i}=(1-1)^{k-j}=[k=j]
$$
则原式化为
$$
f(k)=\binom{k}{k}f(k)=f(k)
$$
恒成立，即形式一得证。

例题

P10596 BZOJ2839 集合计数

定义 $g(i)$ 为选取集合交集元素个数至少为 $i$ 的方案数，此处钦定元素不同则认为方案不同。

则 $g(i)=\binom{n}{i}(2^{2^{n-i}}-1)$ ，代入形式二计算即可。

#include <bits/stdc++.h>
#define FIO cin.tie(0); ios::sync_with_stdio(false)
#define all(x) (x).begin(), (x).end()
#define fi first
#define se second
#define TEST
#define TESTS int t = 1; cin >> t; while (t--)

#if 0
#define int i64
#define inf 0x3f3f3f3f3f3f3f3fLL
#else
#define inf 0x3f3f3f3f
#endif

using namespace std;
using i64 = long long;
using u32 = unsigned;
using u64 = unsigned long long;
using pii = std::pair<int, int>;

constexpr int N = 1e6 + 10;
constexpr int MOD = 1e9 + 7;

i64 f[N], g[N];

i64 ksm(i64 a, i64 b) {
    i64 res = 1;
    while (b) {
        if (b & 1) res = res * a % MOD;
        a = a * a % MOD;
        b >>= 1;
    }
    return res;
}

i64 C(int n, int m) {
    return f[n] * g[m] % MOD * g[n - m] % MOD;
}

void init() {
    f[0] = 1;
    for (int i = 1; i < N; ++i) f[i] = f[i - 1] * i % MOD;
    g[N - 1] = ksm(f[N - 1], MOD - 2);
    for (int i = N - 2; i >= 0; --i) g[i] = g[i + 1] * (i + 1) % MOD;
}

void solve() {
    int n, k;
    cin >> n >> k;
    vector<i64> g(n + 1);
    i64 tmp = 2;
    for (int i = n; i >= 0; --i) {
        g[i] = tmp;
        tmp = tmp * tmp % MOD;
    }
    for (int i = 0; i <= n; ++i) {
        g[i] = (g[i] + MOD - 1) * C(n, i) % MOD;
    }
    i64 ans = 0, o = 1;
    for (int i = k; i <= n; ++i) {
        ans = (ans + o * C(i, k) % MOD * g[i]) % MOD;
        o = (o * (MOD - 1)) % MOD;
    }
    cout << ans << "\n";
}

signed main() {
    
    FIO;
    init();
    solve();

    return 0;
}

P5505 [JSOI2011] 分特产

定义 $g(i)$ 为未分到特产的人数至少为 $i$ 的方案数。

根据隔板法， $g(i)=\binom{n}{i}\prod\limits_{x\in a}\binom{x+n-i-1}{n-i-1}$ ，代入形式二计算即可。

#include <bits/stdc++.h>
#define FIO cin.tie(0); ios::sync_with_stdio(false)
#define all(x) (x).begin(), (x).end()
#define fi first
#define se second
#define TEST
#define TESTS int t = 1; cin >> t; while (t--)

#if 0
#define int i64
#define inf 0x3f3f3f3f3f3f3f3fLL
#else
#define inf 0x3f3f3f3f
#endif

using namespace std;
using i64 = long long;
using u32 = unsigned;
using u64 = unsigned long long;
using pii = std::pair<int, int>;

constexpr int N = 2e3 + 10;
constexpr int MOD = 1e9 + 7;

i64 C[N][N];

void init() {
    for (int i = 0; i < N; ++i) {
        C[i][0] = 1;
        for (int j = 1; j <= i; ++j) {
            C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % MOD;
        }
    }
}

void solve() {
    int n, m;
    cin >> n >> m;
    vector<int> a(m);
    for (int i = 0; i < m; ++i) cin >> a[i];
    // 至少i个人未分到特产
    vector<i64> g(n + 1);
    for (int i = 0; i <= n; ++i) {
        g[i] = C[n][i];
        for (int x : a) {
            g[i] = (g[i] * C[x + n - i - 1][n - i - 1]) % MOD;
        }
    }
    i64 ans = 0, o = 1;
    for (int i = 0; i <= n; ++i) ans = (ans + o * g[i]) % MOD, o = (o * (MOD - 1)) % MOD;
    cout << ans << "\n";
}

signed main() {
    
    FIO;
    init();
    solve();

    return 0;
}

P4859 已经没有什么好害怕的了

转换题目条件，即糖果能量大的组数恰好为 $\frac{n+k}{2}$ 个。记 $t=\frac{n+k}{2}$ 。

定义 $g(i)$ 为糖果能量大的组数至少为 $i$ 的方案数。

使用动态规划，定义 $dp[i][j]$ 为考虑前 $i$ 个糖果，钦定了糖果能量大的组数为 $j$ 的方案数。

考虑如何进行状态转移。容易发现糖果和药片的顺序不影响答案，将其升序排序，使用双指针计算所有药片中值小于每个糖果 $A[i]$ 的数量 $cnt[i]$，则 $dp[i][j]=dp[i-1][j]+dp[i-1][j-1]\times(cnt[i]-j+1)$ 。

那么 $g(i)=dp[n][i] \times (n-i)!$ 。

代入形式二，可得 $f(t)=\sum\limits_{i=t}^{n}(-1)^{i-t}\binom{i}{t}(n-i)! \times dp[n][i]$ 。

#include <bits/stdc++.h>
#define FIO cin.tie(0); ios::sync_with_stdio(false)
#define all(x) (x).begin(), (x).end()
#define fi first
#define se second
#define TEST
#define TESTS int t = 1; cin >> t; while (t--)

#if 0
#define int i64
#define inf 0x3f3f3f3f3f3f3f3fLL
#else
#define inf 0x3f3f3f3f
#endif

using namespace std;
using i64 = long long;
using u32 = unsigned;
using u64 = unsigned long long;
using pii = std::pair<int, int>;

constexpr int N = 2e3 + 10;
constexpr int MOD = 1e9L + 9;

i64 C[N][N], fac[N];

void init() {
    for (int i = 0; i < N; ++i) {
        C[i][0] = 1;
        for (int j = 1; j <= i; ++j) {
            C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % MOD;
        }
    }
    fac[0] = 1;
    for (int i = 1; i < N; ++i) fac[i] = fac[i - 1] * i % MOD;
}

void solve() {
    int n, k;
    cin >> n >> k;
    vector<int> a(n), b(n), s(n);
    for (int i = 0; i < n; ++i) cin >> a[i];
    for (int i = 0; i < n; ++i) cin >> b[i];
    if ((n + k) & 1) {
        cout << 0 << "\n";
        return;
    }
    k = (n + k) / 2;
    ranges::sort(a);
    ranges::sort(b);
    for (int i = 0, p = -1; i < n; ++i) {
        while (p + 1 < n && b[p + 1] < a[i]) p++;
        s[i] = p + 1;
    }
    vector<i64> g(n + 1);
    vector<vector<i64>> dp(n + 1, vector<i64>(n + 1));
    for (int i = 0; i <= n; ++i) {
        dp[i][0] = 1;
        for (int j = 1; j <= i; ++j) {
            dp[i][j] = (dp[i - 1][j] + dp[i - 1][j - 1] * (s[i - 1] - j + 1)) % MOD;
        }
    }
    for (int i = 0; i <= n; ++i) g[i] = dp[n][i] * fac[n - i] % MOD;
    i64 ans = 0, o = 1;
    for (int i = k; i <= n; ++i) {
        ans = (ans + o * C[i][k] % MOD * g[i]) % MOD;
        o = (o * (MOD - 1)) % MOD;
    }
    cout << ans << "\n";
}

signed main() {
    
    FIO;
    init();
    solve();

    return 0;
}

P6478 [NOI Online #2 提高组] 游戏

和上一题做法很类似。

定义 $g(i)$ 为非平局回合数至少为 $i$ 的方案数。

使用动态规划，定义 $dp[u][i]$ 为考虑以 $u$ 为根的子树内，钦定了 $i$ 轮非平局的方案数。

考虑如何进行状态转移。分讨根节点是否被选取参与非平局的情况。

如果根节点不被选取，那么直接做子树间的转移即可。

如果根节点被选取，那么相当于从 $dp[u][i-1]$ 的状态转移过来，设子树内和根节点不同色的点数为 $num$ ，则可能的方案数为 $dp[u][i-1]\times (num-i+1)$ 。

那么 $g(i)=dp[1][i] \times (\frac{n}{2}-i)!$ 。

代入形式二进行计算，整体复杂度 $O(n^2)$ 。

#include <bits/stdc++.h>
#define FIO cin.tie(0); ios::sync_with_stdio(false)
#define all(x) (x).begin(), (x).end()
#define fi first
#define se second
#define TEST
#define TESTS int t = 1; cin >> t; while (t--)

#if 0
#define int i64
#define inf 0x3f3f3f3f3f3f3f3fLL
#else
#define inf 0x3f3f3f3f
#endif

using namespace std;
using i64 = long long;
using u32 = unsigned;
using u64 = unsigned long long;
using pii = std::pair<int, int>;

constexpr int N = 5e3 + 10;
constexpr int MOD = 998244353;

i64 C[N][N], fac[N];

void init() {
    for (int i = 0; i < N; ++i) {
        C[i][0] = 1;
        for (int j = 1; j <= i; ++j) {
            C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % MOD;
        }
    }
    fac[0] = 1;
    for (int i = 1; i < N; ++i) fac[i] = fac[i - 1] * i % MOD;
}

void solve() {
    int n, m;
    cin >> n;
    m = n / 2;
    string s;
    cin >> s;
    vector<int> a(n + 1), sz(n + 1), bak(n + 1);
    vector<vector<int>> e(n + 1);
    vector<i64> f(m + 1), g(m + 1);
    vector<vector<i64>> dp(n + 1, vector<i64>(m + 1));
    vector<array<int, 2>> belong(n + 1);
    for (int i = 0; i < n; ++i) {
        a[i + 1] = s[i] - '0';
    }
    for (int i = 0, u, v; i < n - 1; ++i) {
        cin >> u >> v;
        e[u].push_back(v);
        e[v].push_back(u);
    }
    auto dfs = [&](auto&& dfs, int u, int f) -> void {
        sz[u] = 1;
        dp[u][0] = 1;
        belong[u][a[u]] = 1;
        for (int v : e[u]) {
            if (v == f) continue;
            dfs(dfs, v, u);
            for (int i = 0; i <= sz[u] / 2; ++i) bak[i] = dp[u][i], dp[u][i] = 0;
            for (int l = 0; l <= sz[u] / 2; ++l) {
                for (int r = 0; r <= min(sz[v] / 2, m - l); ++r) {
                    dp[u][l + r] = (dp[u][l + r] + bak[l] * dp[v][r]) % MOD;
                }
            }
            sz[u] += sz[v];
            belong[u][0] += belong[v][0];
            belong[u][1] += belong[v][1];
        }
        int num = belong[u][a[u] ^ 1];
        for (int i = min(num, sz[u] / 2); i >= 1; --i) {
            dp[u][i] = (dp[u][i] + dp[u][i - 1] * (num - i + 1)) % MOD;
        }
    };
    dfs(dfs, 1, 0);
    for (int i = 0; i <= m; ++i) g[i] = dp[1][i] * fac[m - i] % MOD;
    for (int k = 0; k <= m; ++k) {
        i64 o = 1;
        for (int i = k; i <= m; ++i) {
            f[k] = (f[k] + o * C[i][k] % MOD * g[i]) % MOD;
            o = (o * (MOD - 1)) % MOD;
        }
    }
    debug(dp[1]);
    for (int i = 0; i <= m; ++i) cout << f[i] << "\n";
}

signed main() {
    
    FIO;
    init();
    solve();

    return 0;
}

结语

二项式反演与莫比乌斯反演不同，不涉及到大量的和式变换，只需要计算 $g(i)$ ，再代入公式，非常适合反演入门。只不过计算 $g(i)$ 的过程因题而异，借助组合数学、动态规划、生成函数等，八仙过海各显神通。