容斥和反演 1 / 简单反演

通过一些简单的反演公式，可以把统计“恰好”类型的限制转化为“钦定”类型的限制。

子集反演

设全集 \(U=\{1,2,\cdots N\}\)，函数 \(f(S)\) 定义在 \(U\) 的所有子集 \(S\) 上。那么：

形式 \(1\)形式 \(2\)

对于另一函数 \(g(S)\) 满足

\[ g(S)=\sum_{S\subseteq T}f(T) \tag 1 \]

有

\[ f(S)=\sum_{S\subseteq T}(-1)^{|T|-|S|}g(T) \tag 2 \]

对于另一函数 \(g(S)\) 满足

\[ g(S)=\sum_{T\subseteq S}f(T) \tag 1 \]

有

\[ f(S)=\sum_{T\subseteq S}(-1)^{|S|-|T|}g(T) \tag 2 \]

证明

注意到形式 \(1\) 和形式 \(2\) 等价，只需 \(f(S)\gets f(\overline S),\ g(S)\gets g(\overline S)\)。

现在证明形式 \(2\)。将 \((1)\) 代入 \((2)\)，即证：

\[ f(S)=\sum_{T\subseteq S}(-1)^{|S|-|T|}\sum_{V\subseteq T}f(V) \]

考虑 \(f(V)\) 对 \(f(S)\) 的贡献系数

\[ \begin{align*} \sum_{V\subseteq T\subseteq S}(-1)^{|S|-|T|}&=(-1)^{|S|}\sum_{i=|V|}^{|S|}\binom{|S|-|V|}{i-|V|}(-1)^{i}\\ &=(-1)^{|S|-|V|}\sum_{i=0}^{|S|-|V|}\binom{|S|-|V|}{i}(-1)^i\\ &=(-1)^{|S|-|V|}\big[V=S\big]\\ &=[V=S] \end{align*} \]

因此等式恒成立，得证。

二项式反演

设函数 \(f(n)\) 是定义在 \([0,N]\) 的全体整数上，那么：

形式 \(1\)形式 \(2\)

对于另一函数 \(g(n)\) 满足

\[ g(n)=\sum_{k=n}^{N}\binom{k}{n}f(k) \]

有

\[ f(n)=\sum_{k=n}^{N}\binom{k}{n}(-1)^{k-n}g(k) \]

对于另一函数 \(g(n)\) 满足

\[ g(n)=\sum_{k=0}^{n}\binom{n}{k}f(k) \]

有

\[ f(n)=\sum_{k=0}^{n}\binom{n}{k}(-1)^{n-k}g(k) \]

证明：构造 \(f(n)=\sum_{|S|=n}F(S)\)，然后套用子集反演的证明即可。

如果钦定一个集合，会把这个集合的子集或者超集也统计一遍，就可以使用子集反演或者二项式反演求得原问题的答案。可以结合下面的一些例题理解。

错排问题

求有多少个长度为 \(n\) 的排列，满足 \(\forall i,\ p_i\ne i\)。

有 \(n\) 个条件，第 \(i\) 个条件为 \(p_i=i\)，问恰好满足 \(0\) 个的方案数。

记钦定满足 \(k\) 个条件的方案数为 \(s_k\)，容易得到 \(s_k=\binom nk(n-k)!\)。

向后（形式 \(1\)）二项式反演即可。

代码

#include<iostream>
#define int long long

using namespace std;
const int N = 50;

int n, ans;
int a[N];

int c(int a, int b) {
    int res = 1;
    for(int i = a; i >= a - b + 1; i--) res *= i;
    for(int i = b; i >= 1; i--) res /= i;
    return res;
}

int fact(int x) {
    int res = 1;
    for(int i = 1; i <= x; i++) res = res * i;
    return res;
}

signed main() {

    cin >> n;
    for(int i = 0; i <= n; i++) a[i] = c(n, i) * fact(n - i);

    for(int i = 0; i <= n; i += 2) ans += a[i];
    for(int i = 1; i <= n; i += 2) ans -= a[i];

    cout << ans << '\n';

    return 0;
}

P10596 BZOJ2839 集合计数

设全集 \(U={1,2,\cdots n}\)，它有 \(2^n\) 个子集。现在要从这些子集中选取若干个（至少一个），使得它们交集的大小为 \(k\)。问选取的方案数，对 \(10^9+7\) 取模。

\(n\le 10^6,\ 0\le k\le n\)

注意到恰好 \(k\) 个不好做，但是钦定 \(k\) 个作为交集中的元素是显然容易的，方案数为 \(s_k=\binom{n}{k}(2^{2^{n-k}}-1)\)。直接向后二项式反演即可。

注意：指数中的模数为原模数减 \(1\)（费马小定理），即 \(10^9+6\)。

代码

#include<iostream>
#define int long long
using namespace std;
const int N = 1e6 + 10;
const int MOD = 1e9 + 7;

int n, k, ans;
int inv[N], fact[N], ifact[N];
int pw2[N];

#define add(a, b) a = (a + b) % MOD;
#define sub(a, b) a = (a + MOD - b) % MOD;

inline int c(int a, int b) {
    return fact[a] * ifact[a - b] % MOD * ifact[b] % MOD;
}

inline int qpow(int a, int b) {
    int res = 1;
    while(b) {
        if(b & 1) res = res * a % MOD;
        a = a * a % MOD;
        b >>= 1;
    }
    return res;
}

int a[N];

signed main() {

    cin >> n >> k;

    inv[1] = 1;
    fact[0] = ifact[0] = 1;
    pw2[0] = 1;
    for(int i = 2; i <= n; i++) inv[i] = MOD - inv[MOD % i] * (MOD / i) % MOD;
    for(int i = 1; i <= n; i++) fact[i] = fact[i - 1] * i % MOD;
    for(int i = 1; i <= n; i++) ifact[i] = ifact[i - 1] * inv[i] % MOD;
    for(int i = 1; i <= n; i++) pw2[i] = pw2[i - 1] * 2 % (MOD - 1);

    for(int i = k; i <= n; i++) a[i] = c(n, i) * (qpow(2, pw2[n - i]) + MOD - 1) % MOD;

    for(int i = k; i <= n; i += 2) add(ans, c(i, k) * a[i] % MOD);
    for(int i = k + 1; i <= n; i += 2) sub(ans, c(i, k) * a[i] % MOD);

    cout << ans << '\n';

    return 0;
}

P5505 [JSOI2011] 分特产

有 \(m\) 种物品，第 \(i\) 种物品有 \(a_i\) 个，同种物品之间完全相同。你需要把 \(m\) 种物品分给 \(n\) 个人，每个人都必须至少分到一个物品，问分配方案数。

\(n,m,a_i\le 1000\)

由于同种物品之间完全相同，因此我们必须分开单独考虑每种物品，用插板法计算方案数。总方案数应是 \(\prod_{i=1}^{m}{\binom{a_i+n}{n}}\)。然而在分配单种物品的过程中，可以分给一个人 \(0\) 件物品，因此最终有人为空的情况也会被统计到。

考虑二项式反演，钦定 \(k\) 个人为空的方案数 \(s_k=\prod_{i=1}^{m}{\binom{a_i+n-k}{n-k}}\) 是容易计算的。直接向后二项式反演计算恰好 \(0\) 人时的方案数即可。

时间复杂度 \(O(nm)\)。

代码

#include<iostream>
#define int long long
using namespace std;
const int N = 2010;
const int MOD = 1e9 + 7;

int n, m, ans;
int a[N];
int inv[N], fact[N], ifact[N];
int s[N];

int c(int x, int y) {
    return fact[x] * ifact[x - y] % MOD * ifact[y] % MOD;
}

signed main() {

    inv[1] = fact[0] = ifact[0] = 1;
    for(int i = 2; i < N; i++) inv[i] = MOD - inv[MOD % i] * (MOD / i) % MOD;
    for(int i = 1; i < N; i++) fact[i] = fact[i - 1] * i % MOD;
    for(int i = 1; i < N; i++) ifact[i] = ifact[i - 1] * inv[i] % MOD;

    cin >> n >> m;
    for(int i = 1; i <= m; i++) cin >> a[i];

    for(int k = 0; k < n; k++) {
        s[k] = c(n, k);
        for(int i = 1; i <= m; i++) s[k] = s[k] * c(a[i] + n - k - 1, n - k - 1) % MOD;
    }
    for(int i = 0; i <= n; i += 2) ans = (ans + s[i]) % MOD;
    for(int i = 1; i <= n; i += 2) ans = (ans + MOD - s[i]) % MOD;

    cout << ans << '\n';

    return 0;
}

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

给定两个长为 \(n\) 的序列 \(a,b\)（保证不出现重复的数字），你需要将 \(b\) 重新排列，使得 \(a_i> b_i\) 的位置 \(i\) 的数量比 \(a_i<b_i\) 的位置数量恰好多 \(k\) 个。问重排（匹配）方案数，答案对 \(10^9+9\) 取模。

\(n\le 2000, 0\le k\le n\)

先 \(k\gets \frac 12(n+k)\)，表示 \(a_i>b_i\) 的位置恰好有 \(k\) 个。

我们只需要求出钦定 \(x\) 个位置满足 \(a_i>b_i\) 的方案数 \(s_x\)，然后向后二项式反演即可。

注意到 \(s_x\) 不易用组合数直接计算。我们先将 \(a,b\) 分别排好序，考虑 dp，设 \(f_{i,j}\) 表示 \(a\) 的前 \(i\) 项种钦定了 \(j\) 项匹配了一个 \(b_i<a_i\)，对应的方案数。由于序列已经排好序，每个 \(a_i\) 能匹配的位置是一个前缀 \(p\)，使用 lower_bound 求出即可。容易写出转移：

\[ f_{i,j}=f_{i-1,j}+(p-j+1)f_{i-1,j-1} \]

最后 \((n-x)!f_{n,x}\) 就是 \(s_x\)。

代码

#include<iostream>
#include<algorithm>
#define int long long
using namespace std;
const int N = 2010;
const int MOD = 1e9 + 9;

#define add(a, b) a = (a + b) % MOD

int fact[N], ifact[N], inv[N];

int c(int a, int b) {
    return fact[a] * ifact[a - b] % MOD * ifact[b] % MOD;
}

int n, k, ans;
int a[N], b[N];
int f[N];

signed main() {

    inv[1] = fact[0] = ifact[0] = 1;
    for(int i = 2; i < N; i++) inv[i] = MOD - inv[MOD % i] * (MOD / i) % MOD;
    for(int i = 1; i < N; i++) fact[i] = fact[i - 1] * i % MOD;
    for(int i = 1; i < N; i++) ifact[i] = ifact[i - 1] * inv[i] % MOD;

    cin >> n >> k;
    for(int i = 1; i <= n; i++) cin >> a[i];
    for(int i = 1; i <= n; i++) cin >> b[i];

    if((n & 1) != (k & 1)) { cout << "0\n"; return 0; }
    k = (n + k) >> 1;

    sort(a + 1, a + 1 + n);
    sort(b + 1, b + 1 + n);

    for(int i = 1; i <= n; i++) {
        int p = lower_bound(b + 1, b + 1 + n, a[i]) - b - 1;
        f[0] = 1;
        for(int j = min(i, p); j >= 1; j--) {
            add(f[j], f[j - 1] * (p - j + 1) % MOD);
        }
    }

    for(int i = 0; i <= n; i++) f[i] = f[i] * fact[n - i] % MOD;

    for(int i = k; i <= n; i += 2) add(ans, c(i, k) * f[i] % MOD);
    for(int i = k + 1; i <= n; i += 2) add(ans, MOD - c(i, k) * f[i] % MOD);

    cout << ans << '\n';

    return 0;
}

P1450 [HAOI2008] 硬币购物

有四种硬币，面值分别为 \(c_0,c_1,c_2,c_3\)。

有 \(q\) 次询问，每次询问给定 \(d_0,d_1,d_2,d_3,s\)，表示四种硬币各有 \(d_0,d_1,d_2,d_3\) 个，询问拼出面额 \(s\) 的方案数。

\(c_i,d_i,s\le 10^5,\ q\le 1000\)

直接背包过不了。注意到如果没有 \(d_i\) 的限制，那么可以先跑一遍完全背包，然后 \(O(1)\) 回答。

考虑容斥，枚举钦定违反哪些限制，然后直接先拿出 \(d_i+1\) 个硬币。剩下的子问题就是完全背包预处理。

时间复杂度 \(O(4n+2^4q)\)。

代码

#include<iostream>
#define int long long
using namespace std;
const int N = 1e5 + 10;

int c0, c1, c2, c3;
int d0, d1, d2, d3;
int q, s, ans;

int f[N];

int calc(int x) {
    int tmp = s;
    if(x & 1) tmp -= c0 * d0;
    if(x & 2) tmp -= c1 * d1;
    if(x & 4) tmp -= c2 * d2;
    if(x & 8) tmp -= c3 * d3;
    if(tmp < 0) return 0;
    return (__builtin_popcount(x) & 1) ? -f[tmp] : f[tmp];
}

signed main() {

    cin >> c0 >> c1 >> c2 >> c3;

    f[0] = 1;
    for(int i = 0; i <= 100000 - c0; i++) f[i + c0] += f[i];
    for(int i = 0; i <= 100000 - c1; i++) f[i + c1] += f[i];
    for(int i = 0; i <= 100000 - c2; i++) f[i + c2] += f[i];
    for(int i = 0; i <= 100000 - c3; i++) f[i + c3] += f[i];

    cin >> q;
    while(q--) {
        cin >> d0 >> d1 >> d2 >> d3 >> s;
        ++d0; ++d1; ++d2; ++d3;
        ans = 0;
        for(int i = 0; i < 16; i++) ans += calc(i);
        cout << ans << '\n';
    }

    return 0;
}

P3349 [ZJOI2016] 小星星

给你一棵 \(n\) 个节点的树和 \(n\) 个节点的无向连通图。问存在多少种给树重标号的方案，使得重标号之后的树边都存在于图上。

\(n\le 17\)

因为我们要维护两个集合的匹配，考虑在树上进行 dp，将一个集合钦定为子树。设 \(f_{i,S,j}\) 表示，\(i\) 子树匹配了图上的 \(S\) 集合，钦定树上节点 \(i\) 匹配图上节点 \(j\) 的方案数。容易做到 \(O(n^33^n)\) 的时间复杂度，无法通过。

考虑将 dp 中一个集合的状态容斥掉。如果我们不记录子树内的点匹配了哪个集合，可能会出现树上多个点匹配图上一个点的情况。考虑容斥，我们可以先钦定树上的点都是匹配图的一个子集 \(S\) 之内的点，然后直接子集反演即可求出恰好匹配图上的点集 \(\{1,2,\cdots,n\}\) 的方案数。

去掉集合的状态后，一次 dp 的复杂度为 \(O(n^3)\)。总时间复杂度 \(O(n^32^n)\)，可以通过。

代码

#include<iostream>
#define int long long
using namespace std;
const int N = 17;

struct Edge {
    int v, next;
} pool[2 * N + 10];
int ne, head[N + 5];

void addEdge(int u, int v) {
    pool[++ne] = {v, head[u]};
    head[u] = ne;
}

int n, m, ans;
int adj[N + 5][N + 5];

int f[N + 5][N + 5], vld[N + 5];
void dfs(int u, int fa) {
    for(int i = 0; i <= n; i++) f[u][i] = vld[i];
    for(int e = head[u]; e; e = pool[e].next) {
        int v = pool[e].v;
        if(v == fa) continue;
        dfs(v, u);
        for(int i = 1; i <= n; i++) {
            int tmp = 0;
            for(int j = 1; j <= n; j++) {
                if(adj[i][j]) tmp += f[v][j];
            }
            f[u][i] = f[u][i] * tmp;
        }
    }
}

signed main() {

    cin >> n >> m;
    for(int i = 1; i <= m; i++) {
        int u, v;
        cin >> u >> v;
        adj[u][v] = adj[v][u] = 1;
    }

    for(int i = 1; i <= n - 1; i++) {
        int u, v;
        cin >> u >> v;
        addEdge(u, v);
        addEdge(v, u);
    }

    for(int i = 1; i < (1 << n); i++) {
        for(int j = 1; j <= n; j++) vld[j] = (i >> (j - 1)) & 1;
        dfs(1, 0);
        int res = 0;
        for(int i = 1; i <= n; i++) res += f[1][i];
        if(__builtin_popcountll(((1 << n) - 1) ^ i) & 1) ans -= res;
        else ans += res;
    }

    cout << ans << '\n';

    return 0;
}

loj575 不等关系

见 NTT。

P5405 [CTS2019] 氪金手游

有一个抽卡游戏，包含 \(n\) 种卡。第 \(i\) 种卡有权值 \(w_i\)，\(w_i\) 以 \(p_{i,1}\) 的概率为 \(1\)，\(p_{i,2}\) 的概率为 \(2\)，\(p_{i,3}\) 的概率为 \(3\)。\(p_{i,j}\) 是给定的。

现在你将不断进行抽卡，直到集齐全部 \(n\) 种卡。每次抽卡，你抽到第 \(i\) 种卡的概率为 \(\frac{w_i}{\sum_{j=1}^{n}w_j}\)。记第一次抽到第 \(i\) 种卡是在第 \(t_i\) 次抽卡。

给定 \(n-1\) 个有序二元组 \((u_i,v_i)\)，保证它们对应的无向边组成一棵树。问：对于所有 \(i\in [1,n-1]\) 都有 \(t_{u_i}<t_{v_i}\) 的概率是多少。

\(n\le 1000\)

发现这种乱七八糟的树不易处理，我们希望这是一棵内向树或者外向树。对于两种完全相反的限制，考虑容斥。选择树的 \(1\) 号节点为根节点，然后将所有儿子指向父亲的边都容斥掉。这样，原树就变成了外向树森林。

显然森林中的每棵树是相互独立的。我们考虑如何对一棵外向树求答案。在第一次抽出这棵树内的卡时，我们要求这张卡必须是根节点对应的卡，其概率为：在抽出的的卡在树内的前提下抽出根节点对应的卡的概率，应为 \(\frac{w_i}{s_i}\)（\(s_i\) 表示子树内 \(w_i\) 的和）。那么一棵树和森林的答案显然都是 \(\prod\frac{w_i}{s_i}\)。在 dp 时记一下子树和即可。

由于我们不能暴力枚举内向边的子集，因此考虑在 dp 的过程中记录容斥系数。对于一条内向边，我们分讨是否将它容斥，如果钦定它变为外向边，那么执行卷积并添加一个负号；否则求和后执行点积。

时间复杂度 \(O(n^2)\)。

代码

#include<iostream>
#define int long long
using namespace std;
const int N = 1010;
const int MOD = 998244353;

inline int qpow(int a, int b) {
    int res = 1;
    while(b) {
        if(b & 1) res = res * a % MOD;
        a = a * a % MOD;
        b >>= 1;
    }
    return res;
}

struct Edge {
    int v, next;
} pool[2 * N];
int ne, head[N];

void addEdge(int u, int v) {
    pool[++ne] = {v, head[u]};
    head[u] = ne;
}

int n;
int p[N][4];
int sz[N], f[N][3 * N];
int inv[3 * N];

void dfs(int u, int fa) {
    f[u][0] = 1;
    sz[u] = 1;
    for(int e = head[u]; e; e = pool[e].next) {
        int v = pool[e].v;
        if(v == fa) continue;
        dfs(v, u);
        if(e & 1) {
            for(int i = 3 * (sz[u] + sz[v]); i >= 0; i--) {
                int sum = 0;
                for(int j = max(0ll, i - 3 * sz[u]); j <= 3 * sz[v] && j <= i; j++) {
                    sum = (sum + f[u][i - j] * f[v][j] % MOD) % MOD;
                }
                f[u][i] = sum;
            }
        } else {
            int tmp = 0;
            for(int j = 0; j <= 3 * sz[v]; j++) tmp = (tmp + f[v][j]) % MOD;
            for(int i = 3 * (sz[u] + sz[v]); i >= 0; i--) {
                int sum = tmp * f[u][i] % MOD;
                for(int j = max(0ll, i - 3 * sz[u]); j <= 3 * sz[v] && j <= i; j++) {
                    sum = (sum + MOD - f[u][i - j] * f[v][j] % MOD) % MOD;
                }
                f[u][i] = sum;
            }
        }
        sz[u] += sz[v];
    }
    for(int i = 3 * sz[u]; i >= 0; i--) {
        int sum = 0;
        for(int j = 1; j <= 3 && j <= i; j++) {
            sum = (sum + f[u][i - j] * p[u][j] % MOD * j % MOD * inv[i] % MOD) % MOD;
        }
        f[u][i] = sum;
    }
}

signed main() {

    ios::sync_with_stdio(0);
    cin.tie(0);

    inv[1] = 1;
    for(int i = 2; i < 3 * N; i++) inv[i] = MOD - inv[MOD % i] * (MOD / i) % MOD;

    cin >> n;
    for(int i = 1; i <= n; i++) {
        int a1, a2, a3, tmp;
        cin >> a1 >> a2 >> a3;
        tmp = qpow(a1 + a2 + a3, MOD - 2);
        p[i][1] = a1 * tmp % MOD;
        p[i][2] = a2 * tmp % MOD;
        p[i][3] = a3 * tmp % MOD;
    }
    for(int i = 1; i <= n - 1; i++) {
        int u, v;
        cin >> u >> v;
        addEdge(u, v);
        addEdge(v, u);
    }

    dfs(1, 0);

    int ans = 0;
    for(int i = 1; i <= 3 * n; i++) ans = (ans + f[1][i]) % MOD;
    cout << ans << '\n';

    return 0;
}

AT_abc213_g [ABC213G] Connectivity 2

给定一张 \(n\) 个点的无向图。请对每个点 \(k,~(k\in [2,n])\) 求出：无向图的边集有多少个子集，满足 \(1\) 号节点与 \(k\) 节点连通。

\(n\le 17\)

我们可以枚举 \(1\) 号节点所在的连通块 \(S\)，钦定 \(S\to \overline{S}\) 之间没有连边，然后任意选择 \(\overline{S}\) 内部的边。记 \(f(S),~(1\in S)\) 表示 \(1\) 所在的连通块恰好是 \(S\)，\(S\) 内部连边的方案数；\(g(S)\) 表示 \(S\) 内部任意连边的方案数。那么答案就是

\[ ans_k=\sum_{\{1\}\subseteq S\subseteq \{1,2,\cdots,n\}}[k\in S]f(S)g(\overline{S}) \]

\(g(S)\) 显然可以 \(O(m2^n)\) 求出。考虑如何求出 \(f(S)\)。考虑正难则反，从总方案数 \(g(S)\) 中减去 \(S\) 不连通的方案数。在 \(S\) 不连通时，我们不妨设 \(1\) 号节点所在的连通块是 \(T,~(T\ne S)\)。此时方案数应为

\[ f(T)g(S-T) \]

转移时枚举子集即可。算上统计答案，时间复杂度 \(O(m2^n+3^n+n2^n)\)。

代码

#include<iostream>
#define int long long
using namespace std;
const int N = 17;
const int MOD = 998244353;

struct Edge {
    int u, v;
} edg[N * N];

int n, m;
int adj[N + 5][N + 5];

int pw2[N * N];
int f[1 << N], g[1 << N];

int ans[N];

signed main() {

    pw2[0] = 1;
    for(int i = 1; i < N * N; i++) pw2[i] = pw2[i - 1] * 2 % MOD;

    cin >> n >> m;
    for(int i = 1; i <= m; i++) cin >> edg[i].u >> edg[i].v;
    for(int i = 1; i <= m; i++) --edg[i].u, --edg[i].v;

    for(int i = 0; i < (1 << n); i++) {
        int cnt = 0;
        for(int j = 1; j <= m; j++) {
            int u = edg[j].u, v = edg[j].v;
            if((i & (1 << u)) && (i & (1 << v))) ++cnt;
        }
        g[i] = pw2[cnt];
    }

    f[1] = 1;
    for(int i = 3; i < (1 << n); i += 2) {
        f[i] = g[i];
        int j = (i - 1) & i;
        while(j) {
            (f[i] += MOD - f[j] * g[i ^ j] % MOD) %= MOD;
            j = (j - 1) & i;
        }
    }

    for(int i = 3; i < (1 << n); i += 2) {
        int res = f[i] * g[((1 << n) - 1) ^ i] % MOD;
        for(int j = 0; j < n; j++) {
            if(i & (1 << j)) (ans[j + 1] += res) %= MOD;
        }
    }

    for(int i = 2; i <= n; i++) cout << ans[i] << '\n';

    return 0;
}