题解 BZOJ3812 / P11714【[清华集训 2014] 主旋律】
题目描述
给定一个 \(n\) 个点的有向图，问这个图有多少个强连通子图。保证 \(n\leq 15\)。
前置问题
一个有向图有 \(n\) 个点，求它的 DAG 子图数量。\(n\leq 15\)。
解：令 \(f_S\) 表示点集 \(S\) 有多少个 DAG 子图。枚举这个 DAG 中入度为 0 的部分，然后递归下去。由于不能一次枚举完所有入度为 0 的点，我们对着这些点做容斥。对于一张合法的 DAG 来说，假如它的零度点集合为 \(R\)，则我们会计算它 \(\sum_{T\subseteq R,T\neq\varnothing}(-1)^{|T|+1}\) 次，由二项式定理我们知道这个系数恰好为 \([|R|>0]\)，即它是 DAG 时系数为 \(1\)，不是 DAG 时系数为 \(0\)。
solution
一个非强连通的子图，缩点后是一个至少两个点的 DAG（这个 DAG 可以不连通）。
对于点集 \(S\)，考虑强连通子图 \(=\) 所有图 \(-\) 非强连通子图。
令只考虑点集 \(S\) 的强连通子图的数量为 \(f_S\)，枚举 \(T\subset S\)，使 \(T\) 是一个强连通分量，\(S-T\) 就随意，而且强制规定 \(T\) 在缩点图的入度为零。
- \(h_S\) 表示点集为 \(S\) 导出的所有图的数量。
- \(\text{ways}(S,T)\) 表示 \(\sum_{u\in S,v\in T}[(u, v)\in E]\)。注意到 \(S\cap T=\varnothing\)，它的复杂度顶满是 \(O(3^nn^2)\)。
这是错误的。考虑图为 \(\{(1, 3), (2, 3)\}\) 时会算至少两次。
考虑问题在哪里呢？枚举 \(T\) 的时候，我们没有保证 \(T\) 一定是所有的入度为 \(0\) 的点。我们对着 \(T\) 容斥：
- 将 \(T\) 划分为若干个强连通分量，奇数个的贡献是 \(1\)，偶数个的贡献是 \(-1\)，这个系数和前置问题中的系数是一样的。
- 令 \(g_T\) 表示所有划分方案的和。
将 \(f_T\) 换成 \(g_T\) 就行了。
怎么算 \(g_S\)？
- 枚举一个 \(\varnothing\neq T\subset S\)，然后划分成 \(-f_Tg_{S\setminus T}\)。注意 \(T\) 不能和 \(S\) 取等，因为整个 \(S\) 都是强连通分量不是我们需要减去的东西（非强连通子图）。
- 为了避免判重等问题，我们钦定一个 \(x\in S\)，然后枚举 \(x\in T\land T\subset S\) 的 \(T\)。
- 转移顺序为：先计算 \(g_S\) 但是不用加上 \(f_S\)，然后计算 \(f_S\)，最后再将 \(f_S\) 加到 \(g_S\)，保证 \(f_S\) 的转移合法。
复杂度 \(O(3^nn^2)\)。
以上复杂度无法通过。事实上只需要优化 \(\text{ways}(S,T)\) 这个瓶颈的计算。提前计算 \(g_u\) 为一个 \(2^n\) 大小的数，第 \(v\) 位上是 \(1\) 表示存在 \((u, v)\in E\)。计算 \(\text{ways}(S,T)\) 枚举 \(u\in S\)，然后计算 \(g_u\) 与 \(T\) 的按位与，再求这个结果的 __builtin_popcount，即可将复杂度削减一个 \(n\)。
以上复杂度无法通过另一道依赖本问题的题目，这里就不点明是什么题目了。事实上只需要优化 \(\text{ways}(S,T)\) 这个瓶颈的计算。预先处理两个数组 \(o_{S}=\sum_{i\in S}3^i\)，然后将 \((S,T)\) 映射为 \(o_S+2o_T\)，形成一个 \(3^n\) 的数组，计算时只需要随意选定 \(x\in S\)，拿出 \(\text{ways}(S\setminus\{x\},T)\) 的结果和 \(\text{ways}(\{x\},T)\) 相加即可，复杂度 \(O(n2^n+3^n)\)。
code
代码中 b 对应 \(g\) 数组，pw2[h[S]] 对应 \(h_S\)。
#include <bits/stdc++.h>
using namespace std;
#ifdef LOCAL
#define debug(...) fprintf(stderr, ##__VA_ARGS__)
#else
#define endl "\n"
#define debug(...) void(0)
#endif
using LL = long long;
template <unsigned umod>
struct modint {/*{{{*/
static constexpr int mod = umod;
unsigned v;
modint() = default;
template <class T, enable_if_t<is_integral<T>::value, int> = 0>
modint(const T& y) : v((unsigned)(y % mod + (is_signed<T>() && y < 0 ? mod : 0))) {}
modint& operator+=(const modint& rhs) { v += rhs.v; if (v >= umod) v -= umod; return *this; }
modint& operator-=(const modint& rhs) { v -= rhs.v; if (v >= umod) v += umod; return *this; }
modint& operator*=(const modint& rhs) { v = (unsigned)(1ull * v * rhs.v % umod); return *this; }
modint& operator/=(const modint& rhs) { assert(rhs.v); return *this *= qpow(rhs, mod - 2); }
friend modint operator+(modint lhs, const modint& rhs) { return lhs += rhs; }
friend modint operator-(modint lhs, const modint& rhs) { return lhs -= rhs; }
friend modint operator*(modint lhs, const modint& rhs) { return lhs *= rhs; }
friend modint operator/(modint lhs, const modint& rhs) { return lhs /= rhs; }
template <class T> friend modint qpow(modint a, T b) {
modint r = 1;
for (assert(b >= 0); b; b >>= 1, a *= a) if (b & 1) r *= a;
return r;
}
friend int raw(const modint& self) { return self.v; }
friend ostream& operator<<(ostream& os, const modint& self) { return os << raw(self); }
explicit operator bool() const { return v != 0; }
modint operator-() const { return modint(0) - *this; }
bool operator==(const modint& rhs) const { return v == rhs.v; }
bool operator!=(const modint& rhs) const { return v != rhs.v; }
};/*}}}*/
using mint = modint<(int)1e9 + 7>;
template <class T>
constexpr auto lowbit(T x) { return x & -x; }
constexpr int pw3[] = {1, 3, 9, 27, 81, 243, 729, 2187, 6561, 19683, 59049, 177147, 531441, 1594323, 4782969, 14348907, 43046721, 129140163, 387420489, 1162261467};
constexpr int inv2 = (mint::mod + 1) / 2;
struct subset {/*{{{*/
int s;
struct iterator {
int s, t;
constexpr bool operator!=(const iterator& rhs) const { return t != rhs.t; }
constexpr iterator& operator++() { t = t ? (t - 1) & s : -1; return *this; }
constexpr int operator*() const { return t; }
};
constexpr auto begin() const { return iterator{s, s}; }
constexpr auto end() const { return iterator{s, -1}; }
};/*}}}*/
struct eachbit {/*{{{*/
int s;
struct iterator {
int s;
constexpr bool operator!=(const iterator& rhs) const { return s != rhs.s; }
constexpr iterator& operator++() { s -= lowbit(s); return *this; }
constexpr int operator*() const { return __builtin_ctz(s); }
};
constexpr auto begin() const { return iterator{s}; }
constexpr auto end() const { return iterator{0}; }
};/*}}}*/
mint pw2[100010], ipw2[100010];
vector<mint> solve(int n, int g[]) {
vector<mint> f(1 << n), b(1 << n);
vector<int> h(1 << n), w(pw3[n]), ot1(1 << n), ot2(1 << n);
for (int i = 0; i < n; i++) {
for (int j : eachbit{g[i]}) h[1 << i | 1 << j] += 1;
}
for (int i = 0; i < n; i++) {
for (int S = 0; S < 1 << n; S++) if (S >> i & 1) h[S] += h[S ^ (1 << i)];
}
for (int S = 1; S < 1 << n; S++) {
int i = __builtin_ctz(S);
ot1[S] = ot1[S ^ lowbit(S)] + pw3[i];
ot2[S] = ot2[S ^ lowbit(S)] + 2 * pw3[i];
}
for (int S = 1; S < 1 << n; S++) {
for (int T : subset{S}) if (T) {
int i = __builtin_ctz(T);
w[ot1[T] + ot2[S ^ T]] = w[ot1[T ^ lowbit(T)] + ot2[S ^ T]] + __builtin_popcount(g[i] & (S ^ T));
}
}
for (int S = 1; S < 1 << n; S++) {
f[S] = pw2[h[S]];
for (int T : subset{S ^ lowbit(S)}) if (T) b[S] -= b[T] * f[S ^ T];
for (int T : subset{S}) if (T) f[S] -= b[T] * pw2[w[ot1[T] + ot2[S ^ T]] + h[S ^ T]];
b[S] += f[S];
}
return f;
}
int main() {
#ifndef LOCAL
cin.tie(nullptr)->sync_with_stdio(false);
#endif
pw2[0] = 1;
for (int i = 1; i <= 1e5; i++) pw2[i] = pw2[i - 1] + pw2[i - 1];
ipw2[0] = 1;
for (int i = 1; i <= 1e5; i++) ipw2[i] = ipw2[i - 1] * inv2;
int n, g[15]{}, m;
cin >> n >> m;
for (int i = 1, u, v; i <= m; i++) cin >> u >> v, g[--u] |= 1 << --v;
cout << solve(n, g)[(1 << n) - 1] << endl;
return 0;
}
本文来自博客园，作者：caijianhong，转载请注明原文链接：https://www.cnblogs.com/caijianhong/p/solution-bzoj3812.html
