跳转至

容斥和反演

容斥和反演能够将一些不易直接统计的数字,转化为一些能够直接计算的数字乘以一些容斥系数。和二分一样,容斥反演也是反人类直觉的一种首先应当尝试的思路。

简单反演

有些问题可以直接套用公式,实现“钦定”和“恰好”之间的转化。

子集反演

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

对于另一函数 \(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(S)\) 可以表示恰好满足 \(S\) 集合对应的所有条件,对应的方案数。由于“恰好”的限制过于严格,我们很难保证当前方案在满足 \(S\) 的情况下是否会满足其他的若干条件。这时候,我们就可以使用形式 \(1\) 中的 \(g(S)\),它表示钦定至少)满足 \(S\) 中的所有条件,对应的方案数。一般 \(g\) 都是容易求出的。然后,我们再根据反演的公式,就可以得到对应的 \(f\)。形式 \(2\) 的应用同理。

二项式反演

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

对于另一函数 \(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(S)\) 只和 \(|S|\) 有关时的特殊情况。

需要注意的是,在二项式反演(形式 \(1\))中 \(g(k)\) 不是“至少满足 \(k\) 个条件的方案数”(不然直接差分完事了),而是先随意钦定 \(k\) 个条件满足,剩下的随意;也就是说,对于每个 \(f(k)\) 都应该对 \(g(n)\) 产生 \(\binom{k}{n}\) 倍的贡献。

错排问题

求有多少个长度为 \(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;
}

loj575 不等关系

NTT

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;
}

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;
}

复杂容斥系数

容斥原理的核心是,通过设置一系列容斥系数,使得每种合法方案被恰好统计 \(1\) 次,不合法方案被恰好统计 \(0\) 次。

如果我们计算一个集合 \(T\),就会把所有 \(T_0\supseteq T\)\(T_0\subseteq T\))的 \(T_0\) 统计一遍,那么我们可以列一个表,手动递推容斥系数。

P4515 [COCI 2009/2010 #6] XOR

给定 \(n\) 个平面上的等腰直角三角形,三点分别为 \((x_i,y_i),(x_i+r_i,y_i),(x_i,y_i+r_i),~(r_i>0)\)。询问恰好被奇数个三角形覆盖的部分的面积。

\(n\le 10\)

暴力容斥不难做到 \(O(3^n)\)。这里讲如何手动递推容斥系数并做到 \(O(n2^n)\)

\(S\) 为一个三角形的集合,\(f_S\)\(S\) 内三角形交集的面积。设 \(w_i\)\(f_S,~|S|=i\) 的容斥系数。那么一块恰好被集合 \(S_0\) 覆盖的三角形会被统计:

\[ \sum_{\varnothing\subset S\subseteq S_0}w_{|S|}=\sum_{i=1}^{|S_0|}\binom{|S|}{i}w_i \]

而我们希望上式满足

\[ \sum_{\varnothing\subset S\subseteq S_0}w_{|S|}=\sum_{i=1}^{|S_0|}\binom{|S|}{i}w_i=tar_{|S|}=|S|\bmod 2 \]

注意到:

  • \(w_1\) 满足 \(w_1=tar_1=1\)
  • \(w_2\) 满足 \(2w_1+w_2=tar_2=0\),解得 \(w_2=-2\)
  • \(w_3\) 满足 \(3w_1+3w_2+w_3=tar_3=1\),解得 \(w_3=4\)
  • \(w_4\) 满足 \(4w_1+6w_2+4w_3+w_4=tar_4=0\),解得 \(w_4=-8\)
  • ...

对于这种容斥,手动计算前 \(4\) 项大部分可以发现规律,然后再尝试证明即可。我们可以列出一个表格,然后手动递推即可:

\(\vert S\vert\) \(1\) \(2\) \(3\) \(4\) \(5\) ...
\(tar\) \(1\) \(0\) \(1\) \(0\) \(1\) ...
\(w\) \(1\) \(-2\) \(4\) \(-8\) \(16\) ...

时间复杂度 \(O(n2^n)\)

代码
#include<iostream>
#include<iomanip>
#define int long long
using namespace std;
const int N = 12;

int n, ans;
int px[N], py[N], pr[N];
int cnt[1 << N], w[N];

int calc(int s) {
    int x = 0, y = 0, r = 3e6 + 10;
    for(int i = 0; i < n; i++) {
        if(s >> i & 1) {
            int nx = max(x, px[i]), ny = max(y, py[i]);
            int nr = min(x + y + r - nx - ny, px[i] + py[i] + pr[i] - nx - ny);
            x = nx, y = ny, r = nr;
            if(r < 0) return 0;
        }
    }
    return r * r;
}

signed main() {

    cin >> n;
    for(int i = 0; i < n; i++) cin >> px[i] >> py[i] >> pr[i];

    w[1] = 1;
    for(int i = 2; i <= n; i++) w[i] = w[i - 1] * (-2);

    for(int i = 1; i < (1 << n); i++) {
        cnt[i] = cnt[i >> 1] + (i & 1);
        ans += w[cnt[i]] * calc(i);
    }

    cout << fixed << setprecision(1) << (double)ans / 2 << '\n';

    return 0;
}

P11363 [NOIP2024] 树的遍历

给定一棵 \(n\) 个节点的无根树。称两条边是相邻的,当且仅当它们有一个公共端点。

考虑一个 dfs 的过程:初始选择一条边 \(s\) 作为起点,每次找到与当前边相邻且未被访问过的边,递归进入访问。如果不存在,则回溯。每次递归儿子时,都在两条边之间连一条“树边”,最后得到一棵以边为节点的搜索树。称两棵搜索树不同当且仅当一条边 \((u,v)\) 包含于一棵树但不包含于另一棵树。

给定 \(k\) 条关键边,问:可以从某一条关键边出发得到的搜索树 有多少个。答案对 \(10^9+7\) 取模。

\(n\le 10^5\)

首先发现一个点的所有邻边在搜索树上形成一条链,并且顺序(几乎)任意排列。先考虑从一条关键边出发会得到什么样的搜索树。不难发现每个节点的邻边中,必须有一个待在链头,因此答案为 \(\prod (deg_i-1)!\)

再考虑多条关键边的情况。发现从不同的关键边出发,可能得到同一棵搜索树。而给定一棵搜索树,合法的起始边会形成一条链(也就是一个固定的集合 \(S_0\),任选一个子集都能统计到它)。因此如果我们直接枚举关键边的集合 \(S\),再计数能同时被它们得到的搜索树数量,则会算重。因此考虑容斥,对于关键边集合 \(S\),设 \(g_S\) 表示能被关键边集合 \(S\) 中的每条边同时得到的搜索树数量。给 \(g_S\) 赋予 \((-1)^{|S|+1}\) 的容斥系数,即可使每棵合法的搜索树恰好被统计 \(1\) 次:

\(\vert S\vert\) \(1\) \(2\) \(3\) \(4\) \(5\) ...
\(tar\) \(1\) \(1\) \(1\) \(1\) \(1\) ...
\(w\) \(1\) \(-1\) \(1\) \(-1\) \(1\) ...

因为它满足

\[ \sum_{i=1}^{k}\binom ki w_i=tar_k=1 \]

至于如何求 \(g_S\),我们发现链上的点的贡献变为了 \((deg_i-2)!\)。因此对于链上的点都给答案乘以 \(\frac{1}{deg_i-1}\) 即可。这里可以用树上 dp 简单维护。

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

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

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

int inv[N], fact[N];

int T;
int n, k, ans;
int deg[N], dep[N], key[N];
int f[N]; 

void dfs1(int u, int fa) {
    dep[u] = dep[fa] + 1;
    for(int e = head[u]; e; e = pool[e].next) {
        int v = pool[e].v;
        if(v == fa) continue;
        dfs1(v, u);
    }
}

void dfs2(int u, int fa) {
    f[u] = 0;
    for(int e = head[u]; e; e = pool[e].next) {
        int v = pool[e].v;
        if(v == fa) continue;
        dfs2(v, u);
        if(!key[v]) {
            ans = (ans + f[u] * f[v] % MOD) % MOD;
            f[u] = (f[u] + f[v] * inv[deg[u] - 1] % MOD) % MOD;
        } else {
            ans = (ans + MOD - f[v]) % MOD;
            ans = (ans + MOD - 1) % MOD;
            ans = (ans + MOD - f[u]) % MOD;
            f[u] = (f[u] + MOD - inv[deg[u] - 1]) % MOD;
        }
    }
}

void clear() {
    for(int i = 0; i <= n + 5; i++) head[i] = key[i] = deg[i] = 0;
    ne = 1; ans = 0;
}

signed main() {

    inv[1] = fact[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;

    cin >> T >> T;
    while(T--) {
        cin >> n >> k;
        clear();
        for(int i = 1; i <= n - 1; i++) {
            int u, v;
            cin >> u >> v;
            ++deg[u], ++deg[v];
            addEdge(u, v);
            addEdge(v, u);
        }
        dfs1(1, 0);
        for(int i = 1; i <= k; i++) {
            int e, u, v;
            cin >> e;
            u = pool[2 * e].v, v = pool[2 * e + 1].v;
            if(dep[v] < dep[u]) swap(u, v);
            key[v] = 1;
        }
        dfs2(1, 0);
        ans = MOD - ans;
        for(int i = 1; i <= n; i++) ans = ans * fact[deg[i] - 1] % MOD;
        cout << ans << '\n';
    }

    return 0;
}

DAG 子图计数

给定一个有向图,请问有多少个边集是无环的。

\(n\le 17\)

考虑状压,设 \(f_S\) 表示点集 \(S\) 的答案。对于 DAG,我们可以枚举它的顶层节点(入度为 \(0\))组成的集合 \(T\),然后从 \(f_{S-T}\) 转移。然而这样做会将一个 DAG 算重多次:考虑拆贡献,对于一个 DAG,设它的顶层节点集合为 \(T_0\);那么只要我们枚举的 \(T\) 是真实的 \(T_0\) 的子集,就会将这个 DAG 计算一遍。

可以和上一题一样,构造出容斥系数 \((-1)^{|T|+1}\),使得每个 \(|T_0|\) 恰好被统计一次:

\(\vert T\vert\) \(1\) \(2\) \(3\) \(4\) \(5\) ...
\(tar\) \(1\) \(1\) \(1\) \(1\) \(1\) ...
\(w\) \(1\) \(-1\) \(1\) \(-1\) \(1\) ...

或者形式化的:

\[ \sum_{\varnothing\subset T\subseteq S}(-1)^{|T|+1}=\sum_{i=1}^{|S|}(-1)^{i+1}\binom{|S|}{i}=1 \]

因此转移为:

\[ f_S=\sum_{\varnothing\subset T\subseteq S}(-1)^{|T|+1}f_{S-T}2^{E(T,S-T)} \]

至于如何在 \(O(3^n)\) 内计算出两个点集之间的边数,考虑刷表法,从当前集合 \(S\) 转移到 \(S+T\)。这样我们只需在 \(S\) 固定的情况下,动态维护 \(E(S,T)\) 即可。

P11714 [清华集训 2014] 主旋律

给定一个 \(n\) 个节点的有向图,询问它有多少个边集是强联通的(所有节点之间都相互可达)。

\(n\le 15\)

\(ans_S\) 表示点集 \(S\) 构成恰好一个 SCC 的方案数。考虑正难则反,统计图缩点之后的 DAG 包含大于 \(1\) 个节点的情况。这里可以使用经典的 DAG 容斥方法,枚举 DAG 的顶层节点集合 \(T\)

\[ ans_S=2^{E(S,S)}-\sum_{\varnothing\subset T\subseteq S}{(-1)^{?}f(T)2^{E(T,S-T)+E(S-T,S-T)}} \]

其中 \(f_{T}\) 是将点集 \(T\) 划分为若干 SCC 且它们之间没有边的方案数。但是此时我们没法确定上式中的容斥系数,因为 \(-1\) 的指数应该等于 \(T\) 中 SCC 的数量加一。因此我们将容斥系数放到 \(f\) 中计算。考虑如何计算 \(f\),枚举它的划分方案:

\[ f_{S}=-\sum_{\{x\}\subseteq T\subseteq S}ans_Tf_{S-T} \]

其中钦定一个元素 \(x\)\(T\) 中是为了防止算重。那么 \(ans\) 的转移可以写成:

\[ ans_S=2^{E(S,S)}+\sum_{\varnothing\subset T\subseteq S}{f(T)2^{E(T,S-T)+E(S-T,S-T)}} \]

然而还会出现一个问题:\(f_S\)\(ans_S\) 之间似乎存在循环依赖的问题。但马上我们就会注意到,\(ans_S\) 其实不应该统计 \(T=S\)\(T\) 恰好形成一个 SCC 的情况。因此计算 \(f_S\) 时先忽略 \(T=S\) 的项,然后再计算 \(ans_S\),最后再让 \(f_S\) 加上 \(ans_S\) 的贡献,去转移别人。通过刷表法可以动态维护 \(E(T,S-T)\),时间复杂度 \(O(3^n)\)

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

int n, m;
int adj[1 << N];

int pw[300], cnt[1 << N];
int f[1 << N], ans[1 << N];
int ess[1 << N], est[1 << N];

signed main() {

    pw[0] = 1;
    for(int i = 1; i < 300; i++) pw[i] = pw[i - 1] * 2 % MOD;
    for(int i = 1; i < (1 << N); i++) cnt[i] = cnt[i >> 1] + (i & 1);

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

    for(int s = 1; s < (1 << n); s++) {
        for(int i = 0; i < n; i++) {
            if(!(s >> i & 1)) continue;
            for(int j = 0; j < n; j++) {
                if((s >> j & 1) && (adj[1 << i] >> j & 1)) ++ess[s];
            }
        }
    }

    f[0] = ans[0] = 1;
    for(int s = 1; s < (1 << n); s++) {
        int x = s & -s, is = (~s) & ((1 << n) - 1);
        for(int t = (s - 1) & s; t > 0; t = (t - 1) & s) {
            if(~t & x) continue;
            f[s] = (f[s] + MOD - ans[t] * f[s - t] % MOD) % MOD;
        }
        ans[s] = (ans[s] + f[s]) % MOD;
        ans[s] = (pw[ess[s]] + ans[s]) % MOD;
        f[s] = (f[s] + MOD - ans[s]) % MOD;
        for(int tt = (is - 1) & is; tt > 0; tt = (tt - 1) & is) {
            int t = is ^ tt;
            est[t] = est[t - (t & -t)] + cnt[adj[t & -t] & s];
        }
        est[is] = est[is - (is & -is)] + cnt[adj[is & -is] & s];
        for(int t = is; t > 0; t = (t - 1) & is) {
            ans[s + t] = (ans[s + t] + f[s] * pw[est[t]] % MOD * pw[ess[t]] % MOD);
        }
    }

    cout << ans[(1 << n) - 1] << '\n';

    return 0;
}