快速沃尔什变换与子集卷积

前置知识：FFT（快速傅里叶变换）。

快速沃尔什变换

Luogu P4717 【模板】快速莫比乌斯/沃尔什变换 (FMT/FWT)

快速沃尔什变换（Fast Walsh–Hadamard transform）解决二进制运算下的卷积。

给定序列 \(f,g\)，求以下三个序列 \(A,B,C\)：

\[A_i=\sum_{j \operatorname{or} k = i} f_j \times g_k \nonumber \]

\[A_i=\sum_{j \operatorname{and} k = i} f_j \times g_k \nonumber \]

\[A_i=\sum_{j \operatorname{xor} k = i} f_j \times g_k \nonumber \]

其中 \(\operatorname{or}\)、\(\operatorname{and}\)、\(\operatorname{xor}\) 为二进制的或、与、异或运算。

显然可以 \(\mathcal{O}(3^n)\) 做平凡的子集枚举，但这不是我们想要的，我们的目标是 \(\mathcal{O}(n \log n)\) 求，其中 \(n\) 是序列长度。

FWT运用了和FFT相同的思想，将序列变换后做运算再逆变换回去来加速卷积。

或卷积

考虑构造一个运算 \(\mathscr{F}\) 使得若

\[A_i=\sum_{j \operatorname{or} k = i} f_j \times g_k \nonumber \]

则

\[(\mathscr{F}A)_i=(\mathscr{F}f)_i \times (\mathscr{F}g)_i \nonumber \]

并且给定序列 \(f\)，我们可以快速地计算出 \(\mathscr{F}A\) 以及 \(\mathscr{F}^{-1}A\)，其中 \(\mathscr{F}^{-1}\) 表示 \(\mathscr{F}\) 的逆变换，即 \(\mathscr{F}^{-1} (\mathscr{F}A)=A\)。

构造算子 \(\mathscr{F}\) 满足：

\[(\mathscr{F}A)_i=\sum_{j \operatorname{or} i = i}A_j \nonumber \]

也就是 \((\mathscr{F}A)_i\) 是所有 \(A_j\) 的和，其中下标 \(j\) 是 \(i\) 表示的集合的子集。

容易验证 \((\mathscr{F}A)_i=(\mathscr{F}f)_i \times (\mathscr{F}g)_i\)：

\[\begin{aligned} (\mathscr{F}f)_i \times (\mathscr{F}g)_i=\sum_{j \operatorname{or} i = i}\sum_{k \operatorname{or} i = i} f_jg_k=\sum_{(j \operatorname{or} k) \operatorname{or}i = i} f_jg_k=(\mathscr{F}A)_i \end{aligned} \nonumber \]

于是我们只需要快速对一个序列进行变换及逆变换。

对一个序列 \(A\) 进行变换是容易的，先将序列补齐为 \(2^k\) 项，考虑分治：设 \(A_0\) 表示 \(A\) 的前 \(2^{k-1}\) 项（即下标首位为 \(0\) 的位置）构成的序列，\(A_1\) 表示 \(A\) 的后 \(2^{k-1}\) 项即下标首位为 \(1\) 的位置构成的序列，根据 \(\mathscr{F}A_0\) 和 \(\mathscr{F}A_1\) 求出 \(\mathscr{F}A\)。

对于 \(\mathscr{F}A\) 的前 \(2^{k-1}\) 项，它只能由前 \(2^{k-1}\) 项贡献过来，因为后 \(2^{k-1}\) 项首位都是 \(1\)，不可能或出 \(0\) 来；对于后 \(2^{k-1}\) 项，整个序列都可以贡献，因此得到

\[\mathscr{F}A=\operatorname{merge}\Big( \mathscr{F}A_0,\mathscr{F}A_0+\mathscr{F}A_1 \Big) \nonumber \]

其中 \(\operatorname{merge}(P,Q)\) 就是将两个序列像字符串那样拼起来，加法就是对两个序列对应位置做加法得到的新序列。

逆变换也就呼之欲出：

\[\mathscr{F}^{-1}A=\operatorname{merge}\Big( \mathscr{F}^{-1}A_0,\mathscr{F}A_1-\mathscr{F}^{-1}A_0 \Big) \nonumber \]

时间复杂度 \(\mathcal{O}(n \log n)\)，代码如下：

inline void OR_FWT(int A[], int n, int type) {
// type = -1 表示逆变换，否则type = 1
// n 是一个 2^k
    for(int h = 2, H = 1; h <= n; h <<= 1, H <<= 1) {
        for(int i = 0; i < n; i += h) 
            for(int j = 0; j < H; j ++) 
                A[i + j + H] = ((A[i + j + H] + type * A[i + j]) % MOD + MOD) % MOD;
    }
}

与卷积

构造和或卷积类似的算子 \(\mathscr{F}\) 使得

\[(\mathscr{F}A)_i=\sum_{j \operatorname{and} i = i} A_j \nonumber \]

容易验证若 \(A_i=\sum\limits_{j \operatorname{and} k = i} f_jg_k\)，那么 \((\mathscr{F}A)_i=(\mathscr{F}f)_i \times (\mathscr{F}g)_i\)。

同样的，将 \(A\) 补齐到 \(2^k\) 项，设 \(A_0\) 表示前 \(2^{k-1}\) 项构成的序列，\(A_1\) 表示后 \(2^{k-1}\) 项构成的序列，不难得到

\[\begin{aligned} &\mathscr{F}A=\operatorname{merge}\Big( \mathscr{F}A_0+\mathscr{F}A_1,\mathscr{F}A_1 \Big)\newline &\mathscr{F}^{-1}A=\operatorname{merge}\Big( \mathscr{F}^{-1}A_0-\mathscr{F}^{-1}A_1,\mathscr{F}^{-1}A_1 \Big) \end{aligned} \nonumber \]

时间复杂度同样是 \(\mathcal{O}(n \log n)\)，代码如下：

inline void AND_FWT(int A[], int n, int type) {
// type = -1 表示逆变换，否则type = 1
// n 是一个 2^k
    for(int h = 2, H = 1; h <= n; h <<= 1, H <<= 1) {
        for(int i = 0; i < n; i += h)
            for(int j = 0; j < H; j ++)
                A[i + j] = ((A[i + j] + type * A[i + j + H]) % MOD + MOD) % MOD;
    }
}

异或卷积

异或的情况就比较不同了，我们不能照搬之前与和或的形式来做变换（\(j \operatorname{xor} i = i \Rightarrow j = 0\)，这样的变换肯定是不可行的）。

接下来给出异或的正确变换，为了方便，对于数 \(x\)，记 \(|x|\) 表示 \(x\) 的二进制表示中 \(1\) 的个数，即 \(|x|=\operatorname{popcount}(x)\)。

基于奇偶性，由于 \(\big((a \operatorname{and} b) \bmod 2\big) \operatorname{xor} \big((b \operatorname{and} c)]\bmod 2 \big)=\big( (a \operatorname{and} c) \bmod 2 \big)\)，得到变换 \(\mathscr{F}\)：

\[(\mathscr{F}A)_i=\sum_{|j \operatorname{and} i| \bmod 2 = 0}A_j-\sum_{|j \operatorname{and} i|\bmod 2=1}A_j \nonumber \]

为了方便，记 \(S_i\) 表示满足 \(|j \operatorname{and} i| \bmod 2 = 0\) 的 \(j\) 构成的集合，\(T_i\) 表示满足 \(|j \operatorname{and} i| \bmod 2 = 1\) 的 \(j\) 构成的集合，那么上式可以就写成

\[(\mathscr{F}A)_i=\sum_{j \in S_i}A_j - \sum_{j \in T_i}A_j \nonumber \]

下证若 \(A_i=\sum\limits_{j \operatorname{xor} k=i}f_jg_k\)，那么 \((\mathscr{F}A)_i=(\mathscr{F}f)_i \times (\mathscr{F}g)_i\)。

证明：

\[\begin{aligned} (\mathscr{F}f)_i \times (\mathscr{F}g)_i&=\left(\sum_{j \in S_i}f_j-\sum_{j \in T_i}f_j\right) \times \left(\sum_{j \in S_i}g_j-\sum_{j \in T_i}g_j\right)\newline &=\sum_{j,k \in S_i}f_jg_k + \sum_{j,k \in T_i}f_jg_k-\sum_{j \in S_i, k \in T_i}f_jg_k - \sum_{j \in T_i, k \in S_i}f_jg_k\newline &=\sum_{(j \operatorname{xor} k) \in S_i}f_jg_k-\sum_{(j \operatorname{xor} k) \in T_i}f_jg_k\newline &=(\mathscr{F}A)_i \end{aligned} \nonumber \]

然后还是分治求 \(\mathscr{F}A\) 和 \(\mathscr{F}^{-1}A\)，将 \(A\) 补齐为 \(2^k\) 项，设 \(A_0\) 为前 \(2^{k-1}\) 项构成的序列，\(A_1\) 为后 \(2^{k-1}\) 项构成的序列，那么

\[\begin{aligned} &\mathscr{F}A=\operatorname{merge}(\mathscr{F}A_0+\mathscr{F}A_1,\mathscr{F}A_0-\mathscr{F}A_1)\newline &\mathscr{F}^{-1}A=\operatorname{merge}(\frac{\mathscr{F}^{-1}A_0+\mathscr{F}^{-1}A_1}{2},\frac{\mathscr{F}^{-1}A_0-\mathscr{F}^{-1}A_1}{2}) \end{aligned} \nonumber \]

时间复杂度 \(\mathcal{O}(n \log n)\)，代码如下：

inline void XOR_FWT(int A[], int n, int type) {
// type = 2的逆元 表示逆变换，否则type = 1
// n 是一个 2^k
    for(int h = 2, H = 1, x, y; h <= n; h <<= 1, H <<= 1) {
        for(int i = 0; i < n; i += h)
            for(int j = 0; j < H; j ++) {
                x = A[i + j], y = A[i + j + H]; 
                A[i + j] = 1ll * (x + y) % MOD * type % MOD, A[i + j + H] = 1ll * ((x - y) % MOD + MOD) % MOD * type % MOD;
            }
    }
}

三个FWT合并起来的代码如下：

code for Luogu P4717

#include <bits/stdc++.h>

using namespace std;

const int maxn = (1 << 18) + 10, MOD = 998244353, inv2 = 499122177;
int n, A[maxn], B[maxn], C[maxn];

inline void print(int A[], int n) {
    for(int i = 0; i < n; i ++) cout << A[i] << ' ';
    cout << '\n';
}

inline void OR_FWT(int A[], int n, int type) {
// type = -1 表示逆变换，否则type = 1
// n 是一个 2^k
    for(int h = 2, H = 1; h <= n; h <<= 1, H <<= 1) {
        for(int i = 0; i < n; i += h) 
            for(int j = 0; j < H; j ++) 
                A[i + j + H] = ((A[i + j + H] + type * A[i + j]) % MOD + MOD) % MOD;
    }
}
inline void AND_FWT(int A[], int n, int type) {
// type = -1 表示逆变换，否则type = 1
// n 是一个 2^k
    for(int h = 2, H = 1; h <= n; h <<= 1, H <<= 1) {
        for(int i = 0; i < n; i += h)
            for(int j = 0; j < H; j ++)
                A[i + j] = ((A[i + j] + type * A[i + j + H]) % MOD + MOD) % MOD;
    }
}
inline void XOR_FWT(int A[], int n, int type) {
// type = 2的逆元 表示逆变换，否则type = 1
// n 是一个 2^k
    for(int h = 2, H = 1, x, y; h <= n; h <<= 1, H <<= 1) {
        for(int i = 0; i < n; i += h)
            for(int j = 0; j < H; j ++) {
                x = A[i + j], y = A[i + j + H]; 
                A[i + j] = 1ll * (x + y) % MOD * type % MOD, A[i + j + H] = 1ll * ((x - y) % MOD + MOD) % MOD * type % MOD;
            }
    }
}

int main() {
    ios::sync_with_stdio(false), cin.tie(0);

    cin >> n; n = (1 << n);
    for(int i = 0; i < n; i ++) cin >> A[i];
    for(int i = 0; i < n; i ++) cin >> B[i];

    OR_FWT(A, n, 1), OR_FWT(B, n, 1);
    for(int i = 0; i < n; i ++) C[i] = 1ll * A[i] * B[i] % MOD;
    OR_FWT(C, n, -1); print(C, n); OR_FWT(A, n, -1), OR_FWT(B, n, -1);

    AND_FWT(A, n, 1), AND_FWT(B, n, 1);
    for(int i = 0; i < n; i ++) C[i] = 1ll * A[i] * B[i] % MOD;
    AND_FWT(C, n, -1); print(C, n); AND_FWT(A, n, -1), AND_FWT(B, n, -1);

    XOR_FWT(A, n, 1), XOR_FWT(B, n, 1);
    for(int i = 0; i < n; i ++) C[i] = 1ll * A[i] * B[i] % MOD;
    XOR_FWT(C, n, inv2); print(C, n); XOR_FWT(A, n, inv2), XOR_FWT(B, n, inv2);
    
    return 0;
}

集合幂级数

就是将集合放到指数位置上，例如 \(f(x)=x^{\emptyset}+2x^{\{1\}}+4x^{\{2\}}+3x^{\{1,2\}}\)。

设全集为 \(U\)，那么集合幂级数可以表示为

\[f(x)=\sum_{S \subseteq U} f_Sx^S \nonumber \]

定义集合幂级数的加、减法为

\[f(x) \pm g(x)=\sum_{S \subseteq U} (f_S \pm g_S)x^S \nonumber \]

子集卷积

Luogu P6097 【模板】子集卷积

基本形式

给定两个长度为 \(2^n\) 的序列 \(a_0,a_1,\cdots,a_{2^n-1}\) 和 \(b_0,b_1,\cdots,b_{2^n-1}\)，你需要求出一个序列 \(c_0,c_1,\cdots,c_{2^n-1}\)，其中 \(c_k\) 满足：

\[c_k=\sum_{\substack{ {i \operatorname{and} j=0}\newline{i \operatorname{or} j=k} } } a_i b_j \nonumber \]

其中 \(\operatorname{and}\) 和 \(\operatorname{or}\) 为位运算。

尝试将它转变成我们熟悉的与、或卷积的形式：注意到若 \(i \operatorname{or} j=k\) 并且 \(|i|+|j|=|k|\)（其中 \(|x|\) 表示 \(\operatorname{popcount}(x)\)，即 \(x\) 的二进制中 \(1\) 的个数），那么就自然满足 \(i \operatorname{and} j=k\)，于是设

\[A_{i,j}=\begin{cases} a_j, &|j|=i\newline 0, &|j|\ne i \end{cases} \nonumber \]

同理定义 \(B_{i,j},C_{i,j}\)，上式就可以改写为

\[C_{l,k}=\sum_{i=0}^{l}\left( \sum_{x \operatorname{or} y = k} A_{i,x} \times B_{l-i,y} \right) \nonumber \]

后面括号里面是或卷积的形式，可以 \(\mathcal{O}(n2^n)\) 算，对于每个 \(i \in [0,|k|]\) 都算一次，时间复杂度 \(\mathcal{O}(n^2 2^n)\)。

集合幂级数下的解释

在集合幂级数的意义下，相当与加一个变量 \(z^{|S|}\)，两个集合幂级数的乘积为

\[\left( \sum_{S \subseteq U} f_Sx^Sz^{|S|} \right) \times \left( \sum_{S \subseteq U} g_Sx^Sz^{|S|} \right) \nonumber \]

最终提取出 \(x^{T}z^{|T|}\) 的系数就是答案。

code for Luogu P6097

#include <bits/stdc++.h>

using namespace std;

const int N = (1 << 20) + 10, MOD = 1e9 + 9;
inline int Plus(int a, int b) {return a + b >= MOD ? a + b - MOD : a + b; }
inline int Minus(int a, int b) {return a - b < 0 ? a - b + MOD : a - b; }
int n, a[N], b[N], c[N];
int A[21][N], B[21][N], C[21][N];

inline void FWT(int A[], int n, int type) {
    for(int h = 2; h <= n; h <<= 1) 
        for(int i = 0; i < n; i += h)
            for(int j = i; j < i + (h >> 1); j ++) {
                if(type == 1) A[j + (h >> 1)] = Plus(A[j + (h >> 1)], A[j]);
                else A[j + (h >> 1)] = Minus(A[j + (h >> 1)], A[j]);
            }
}

int main() {
    ios::sync_with_stdio(false), cin.tie(0);

    cin >> n;
    for(int i = 0; i < (1 << n); i ++)
        cin >> a[i], A[__builtin_popcount(i)][i] = a[i];
    for(int i = 0; i < (1 << n); i ++)
        cin >> b[i], B[__builtin_popcount(i)][i] = b[i];
    for(int i = 0; i <= n; i ++)
        FWT(A[i], 1 << n, 1), FWT(B[i], 1 << n, 1);
    
    for(int i = 0; i <= n; i ++) {
        for(int k = 0; k <= i; k ++)
            for(int j = 0; j < (1 << n); j ++)
                C[i][j] = Plus(C[i][j], 1ll * A[k][j] * B[i - k][j] % MOD);
        FWT(C[i], 1 << n, -1);
    }
    for(int i = 0; i < (1 << n); i ++)
        c[i] = C[__builtin_popcount(i)][i];
    for(int i = 0; i < (1 << n); i ++)
        cout << c[i] << ' ';
    cout << '\n';
    
    return 0;
}

注意看这个式子

\[C_{l,k}=\sum_{i=0}^{l}\left( \sum_{x \operatorname{or} y = k} A_{i,x} \times B_{l-i,y} \right) \nonumber \]

改写成

\[C_i=\sum_{i_1+i_2=i}A_{i_1} *_{\operatorname{or}} B_{i_2} \nonumber \]

一个卷积？是不是也可以做求逆、\(\exp\)、\(\ln\)！

定义集合幂级数 \(\exp f\) 为满足 \(\exp f =\sum_{i=0}^{+\infty}\frac{f^i}{i!}\) 的集合幂级数，\(\ln f\) 为满足 \(\exp (\ln f) = f\) 的集合幂级数。

实际上是当然可以的，对每个 \(A_i,B_i\) 都FWT一下，现在对每个 \(x\)，\(C_{i, x}\) 是互相独立的，因此可以对每个 \(x \in [0,2^n-1]\) 求出这一维的答案，然后最终合并，问题转化为长度为 \(n\) 的序列的求逆、\(\ln\)、\(\exp\)，当然可以上多项式全家桶做到 \(\mathcal{O}(2^n n \log n)\)，但是这是完全没必要的，\(n\) 很小，也就 \(20\) 左右，多项式全家桶的大常数是得不偿失的，可以直接 \(\mathcal{O}(n^2)\) 暴力做这些操作，时间复杂度是 \(\mathcal{O}(n^2 2^n)\)。

再细说一下的话，求逆用递推 \(g_n=\frac{1-\sum_{j=0}^{i-1}g_jf_{i-j}}{f_0}\)，\(\exp\) 用 \(g'=gf'\)，写成系数形式就是 \(g_n=\frac{1}{n}\sum_{i=1}^{n}if_ig_{n-i}\)。\(\ln\) 可以利用它是 \(\exp\) 的逆变换得到 \(g_n=f_n-\frac{1}{n}\sum_{i=1}^{n-1}ig_if_{n-i}\)。

求逆

for(int i = 0; i <= n; i ++)
    FWT(f[i], 1 << n, 1);
for(int mask = 0; mask < (1 << n); mask ++) {
    g[0][mask] = 1;
    for(int i = 1; i <= n; i ++)
        for(int j = 1; j <= i; j ++)
            g[i][mask] = Minus(g[i][mask], 1ll * f[j][mask] * g[i - j][mask] % MOD);
}
for(int i = 0; i <= n; i ++)
    FWT(g[i], 1 << n, -1);

$\exp$

inline void exp(int f[], int g[], int n) {
    for(int i = 0; i <= n; i ++) g[i] = 0;
    g[0] = 1;
    for(int i = 1; i <= n; i ++) {
        for(int j = 1; j <= i; j ++)
            g[i] = Plus(g[i], 1ll * f[j] * j % MOD * g[i - j] % MOD);
        g[i] = 1ll * g[i] * inv[i] % MOD;   // inv[i] 是 i 的逆元
    }
}

// 主函数中
for(int i = 0; i <= n; i ++)
    FWT(f[i], 1 << n, 1);
for(int mask = 0; mask < (1 << n); mask ++) {
    static int F[N + 1], G[N + 1];
    for(int i = 0; i <= n; i ++) 
        F[i] = f[i][mask];
    exp(F, G, n);
    for(int i = 0; i <= n; i ++)
        g[i][mask] = G[i];
}
for(int i = 0; i <= n; i ++)
    FWT(g[i], 1 << n, -1);

$\ln$

inline void ln(int f[], int g[], int n) {
    for(int i = 0; i <= n; i ++) g[i] = 0;
    for(int i = 1; i <= n; i ++) {
        for(int j = 1; j < i; j ++)
            g[i] = Plus(g[i], 1ll * g[j] * j % MOD * f[i - j] % MOD);
        g[i] = Minus(f[i], 1ll * g[i] * inv[i] % MOD); // inv[i] 是 i 的逆元
    }
}

// 主函数中
for(int i = 0; i <= n; i ++)
    FWT(f[i], 1 << n, 1);
for(int mask = 0; mask < (1 << n); mask ++) {
    static int F[N + 1], G[N + 1];
    for(int i = 0; i <= n; i ++)
        F[i] = f[i][mask];
    ln(F, G, n);
    for(int i = 0; i <= n; i ++)
        g[i][mask] = G[i];
}
for(int i = 0; i <= n; i ++)
    FWT(g[i], 1 << n, -1);

例题

「WC2018」 州区划分

题目描述

小 S 现在拥有 \(n\) 座城市，第 \(i\) 座城市的人口为 \(w_i\)，城市与城市之间可能有双向道路相连。

现在小 S 要将这 \(n\) 座城市划分成若干个州，每个州由至少一个城市组成，每个城市在恰好一个州内。

假设小 S 将这些城市划分成了 \(k\) 个州，设 \(V_i\) 是第 \(i\) 个州包含的所有城市组成的集合。定义一条道路是一个州的内部道路，当且仅当这条道路的两个端点城市都在这个州内。如果一个州内部存在一条起点终点相同，不经过任何不属于这个州的城市，且经过这个州的所有内部道路都恰好一次并且经过这个州的所有城市至少一次的路径（路径长度可以为 \(0\)），则称这个州是不合法的。

定义第 \(i\) 个州的满意度为：第 \(i\) 个州的人口在前 \(i\) 个州的人口中所占比例的 \(p\) 次幂，即：

\[\left(\dfrac{\sum _ {x \in V _ i} w _ x}{\sum _ {j = 1} ^ i \sum _ {x \in V _ j} w _ x}\right) ^ p \]

定义一个划分的满意度为所有州的满意度的乘积。

求所有合法的划分方案的满意度之和。

答案对 \(998244353\) 取模。
两个划分 \(\{V_1, V _ 2, \cdots, V_k\}\) 和 \(\{C_1, C _ 2, \cdots, C_s\}\) 是不同的，当且仅当 \(k \neq s\)，或存在某个 \(1 \leq i \leq k\)，使得 \(V_i \neq C_i\)。

保证对于所有数据有：\(0 \leq n \leq 21\)， \(0 \leq m \leq \dfrac{n\times (n-1)}{2}\) ， \(0 \leq p \leq 2\)， \(1 \leq w_i \leq 100\)。

题解

注意到一个州不合法当且仅当它连通且存在欧拉回路，状压dp，设 \(f_S\) 表示集合 \(S\) 的答案，\(g_S\) 表示集合 \(S\) 内所有点的 \(w\) 和，就有转移

\[f_S=\sum_{T \subset S}f_{S-T}\left(\frac{g_T}{g_S}\right)^p=\frac{1}{g_S^p}\sum_{T \subset S} f_{S-T}g_T^p \nonumber \]

这是一个子集卷积的形式，按照 \(|S|\) 升序，依次算出 \(f\) 即可。

时间复杂度 \(\mathcal{O}(n^22^n)\)。

code

#include <bits/stdc++.h>

using namespace std;

const int N = 21, M = (1 << N), MOD = 998244353;
inline int Plus(int a, int b) {return a + b >= MOD ? a + b - MOD : a + b; }
inline int Minus(int a, int b) {return a - b < 0 ? a - b + MOD : a - b; }
inline int ksm(long long a, int b) {
    long long r = 1;
    for(; b; b >>= 1, a = a * a % MOD)
        if(b & 1) r = r * a % MOD;
    return r;
}

int n, m, p, w[N];
vector<int> G[N];
int f[N + 1][M], g[N + 1][M], ig[M];

bool vis[N];
void dfs(int u, int mask, int c) {
    vis[u] = c;
    for(auto v : G[u]) 
        if((mask >> v & 1) && (vis[v] != c))
            dfs(v, mask, c);
}
inline int check(int mask) {
    static int deg[N];
    bool flag = false;
    for(int i = 0; i < n; i ++) if(mask >> i & 1) {
        for(auto j : G[i]) if(mask >> j & 1) deg[i] ++;
        if(deg[i] & 1) flag = true;
    }
    for(int i = 0; i < n; i ++) if(mask >> i & 1) {
        for(auto j : G[i]) if(mask >> j & 1) deg[i] --;
    }
    for(int i = 0; i < n; i ++) 
        if(mask >> i & 1) {dfs(i, mask, 1); break; }
    int siz = 0;
    for(int i = 0; i < n; i ++) 
        if(mask >> i & 1) siz += vis[i];
    if(siz != __builtin_popcount(mask)) flag = true;
    for(int i = 0; i < n; i ++)
        if(mask >> i & 1) {dfs(i, mask, 0); break; }
    if(!flag) return 0;
    int sum = 0;
    for(int i = 0; i < n; i ++)
        if(mask >> i & 1) sum = Plus(sum, w[i]);
    return sum;
}

inline void FWT(int A[], int n, int type) {
    for(int h = 2; h <= n; h <<= 1) 
        for(int i = 0; i < n; i += h)
            for(int j = i; j < i + (h >> 1); j ++) {
                if(type == 1) A[j + (h >> 1)] = Plus(A[j + (h >> 1)], A[j]);
                else A[j + (h >> 1)] = Minus(A[j + (h >> 1)], A[j]);
            }
}

int main() {
    ios::sync_with_stdio(false), cin.tie(0);

    cin >> n >> m >> p;
    for(int i = 0; i < m; i ++) {
        int a, b; cin >> a >> b;
        a --, b --; G[a].emplace_back(b), G[b].emplace_back(a);
    }
    for(int i = 0; i < n; i ++) cin >> w[i];

    for(int mask = 1; mask < (1 << n); mask ++) {
        int x = check(mask);
        if(x != 0) {
            if(p == 0) g[__builtin_popcount(mask)][mask] = 1;
            else if(p == 1) g[__builtin_popcount(mask)][mask] = x;
            else g[__builtin_popcount(mask)][mask] = 1ll * x * x % MOD;
        }
        int sum = 0;
        for(int i = 0; i < n; i ++)
            if(mask >> i & 1) sum = Plus(sum, w[i]);
        ig[mask] = ksm(ksm(sum, MOD - 2), p);
    }
    for(int i = 0; i <= n; i ++)
        FWT(g[i], 1 << n, 1);

    f[0][0] = 1; FWT(f[0], 1 << n, 1);
    for(int i = 1; i <= n; i ++) {
        for(int j = 0; j < i; j ++)
            for(int k = 0; k < (1 << n); k ++)
                f[i][k] = Plus(f[i][k], 1ll * f[j][k] * g[i - j][k] % MOD);
        FWT(f[i], 1 << n, -1);
        for(int k = 0; k < (1 << n); k ++)
            f[i][k] = __builtin_popcount(k) == i ? 1ll * f[i][k] * ig[k] % MOD : 0;
        FWT(f[i], 1 << n, 1);
    }
    FWT(f[n], 1 << n, -1);
    cout << f[n][(1 << n) - 1] << '\n';

    return 0;
}

「CEOI2019」 Amusement Park

题目描述

有一个含 \(n\) 个点，\(m\) 条边的有向图，图无重边，无自环，两点之间不成环。

现在我们想改变一些边的方向，使得该有向图无环。

您需要求出，每一种改变方向后使得该有向图无环的方案的需改变边的数量之和 \(\bmod\ 998244353\) 之后的答案。

对于 \(100\%\) 的数据，保证 \(1\le n\le 18\)，\(0\le m\le \frac{n\times (n-1)}{2}\)，\(1\le a_i,b_i\le n\)，\(a_i\not=b_i\)，对于 \(i\not=j\)，均有 \(a_i\not=a_j\) 或者 \(b_i\not=b_j\)，无序数对 \(\{a_i,b_i\}\) 互不相同。

题解

主要到一个DAG的所有边反向后仍然是DAG，因此答案就是方案数乘 \(\frac{m}{2}\)。

设 \(f_S\) 表示集合 \(S\) 中点的方案数，转移时钦定一些度数为 \(0\) 的点后容斥出结果，具体的有

\[f_S=\sum_{\substack{ {T \subseteq S}\newline\newline{T\text{是独立集} } } }(-1)^{|T|-1}f_{S-T} \nonumber \]

于是就可以枚举子集做到 \(\mathcal{O}(3^n)\) 了，应该可以通过，但接下来我们将它优化到 \(\mathcal{O}(n^22^n)\)。

设 \(g_S=(-1)^{|S|-1}[S\text{是独立集}]\)，那么 \(f=f \times g + 1\)，这里的乘法是子集卷积。解得 \(f=\frac{1}{1-g}\)，做一个子集卷积求逆即可，时间复杂度 \(\mathcal{O}(n^22^n)\)。

code

#include <bits/stdc++.h>

using namespace std;

const int N = 18, M = (1 << N), MOD = 998244353;
inline int Plus(int a, int b) {return a + b >= MOD ? a + b - MOD : a + b; }
inline int Minus(int a, int b) {return a - b < 0 ? a - b + MOD : a - b; }
inline int ksm(long long a, int b) {
    long long r = 1;
    for(; b; b >>= 1, a = a * a % MOD)
        if(b & 1) r = r * a % MOD;
    return r;
}
int n, m;
int flag[M];

inline void FWT(int A[], int n, int type) {
    for(int h = 2; h <= n; h <<= 1)
        for(int i = 0; i < n; i += h)
            for(int j = i; j < i + (h >> 1); j ++) {
                if(type == 1) A[j + (h >> 1)] = Plus(A[j + (h >> 1)], A[j]);
                else A[j + (h >> 1)] = Minus(A[j + (h >> 1)], A[j]);
            }
}

int f[N + 1][M], g[N + 1][M];

int main() {
    ios::sync_with_stdio(false), cin.tie(0);

    cin >> n >> m;
    for(int i = 0; i < m; i ++) {
        int a, b; cin >> a >> b; a--, b --;
        flag[(1 << a) | (1 << b)] = 1;
    }
    FWT(flag, 1 << n, 1);
    for(int mask = 0; mask < (1 << n); mask ++) {
        if(flag[mask]) continue;
        if(__builtin_popcount(mask) & 1) f[__builtin_popcount(mask)][mask] = 1;
        else f[__builtin_popcount(mask)][mask] = MOD - 1;
    }

    for(int i = 0; i <= n; i ++)
        FWT(f[i], 1 << n, 1);
    for(int mask = 0; mask < (1 << n); mask ++) {
        for(int i = 1; i <= n; i ++)
            f[i][mask] = Minus(0, f[i][mask]);
        g[0][mask] = f[0][mask] = 1;
        for(int i = 1; i <= n; i ++)
            for(int j = 1; j <= i; j ++)
                g[i][mask] = Minus(g[i][mask], 1ll * f[j][mask] * g[i - j][mask] % MOD);
    }
    FWT(g[n], 1 << n, -1);
    cout << 1ll * g[n][(1 << n) - 1] * m % MOD * ksm(2, MOD - 2) % MOD << '\n';

    return 0;
}

ABC236Ex Distinct Multiples

题目描述

给定两个正整数 \(N,M\) 和一个正整数序列 \(D\)，询问满足条件的序列 \(A\) 的个数：

• \(1\leq A_i\leq M(1\leq i\leq N)\)
• \(A_i\neq A_j(1\leq i<j\leq N)\)
• \(D_i|A_i\)

满足 \(2\leq N\leq 16\),\(1\leq M\leq 10^{18}\),\(1\leq D_i\leq M\)。

题解

见 AtCoder ABC Ex乱写。

loj154. 集合划分计数

题目描述

给定一个集合 \(S = \{ {x_1}, {x_2}, \dots, {x_n} \}\) 和一个 \(S\) 上的集合族 \(\mathcal F = \{ {S_0}, {S_1}, \dots, {S_{m-1}} \}\)。

一个划分 \(\mathcal P\) 是 \(\mathcal F\) 的一个子族，满足 \(\mathcal P\) 中所有集合的并为 \(S\) ，任意两个集合不相交。

求大小不大于 \(k\) 的划分的数量 \(\text{mod } 998244353\)。

两个划分 \({\mathcal P_1}\), \({\mathcal P_2}\) 不同，当且仅当存在 \(i\) 使 \(S_i \in {\mathcal P_1} \land S_i \notin {\mathcal P_2}\) 或 \(S_i \notin {\mathcal P_1} \land S_i \in {\mathcal P_2}\)。\(S_i\) 和 \(S_j\) 不同当且仅当 \(i \neq j\)。

题解

这就是一个不完整的子集卷积 \(\exp\) 。记 \(f(x)\) 为所给集族 \(\mathcal{F}\)，所求即为

\[\sum_{i=0}^{k}\frac{f^i(x)}{i!} \nonumber \]

这里的乘法是用子集卷积定义的，所以不用担心转移时集合之间有交的问题。

怎么算这个东西呢？记上面式子为 \(g\)，求导得到

\[g'=\sum_{i=0}^{k}\frac{if^{i-1}f'}{i!}=f'\sum_{i=0}^{k-1}\frac{f^i}{i!} \nonumber \]

即

\[g'=f'(g-\frac{f^k}{k!}) \nonumber \]

用先 \(\ln\) 再乘 \(k\) 再 \(\exp\) 的方法容易求得 \(h=\frac{f^k}{k!}\)，之后系数形式就可以写作

\[g_n=\frac{1}{n}\sum_{i=1}^{n}if_i\left(g_{n-i}-h_{n-i}\right) \nonumber \]

然后就可以直接求了。

code

#include <bits/stdc++.h>

using namespace std;

const int N = 21, MOD = 998244353;
inline int Plus(int a, int b) {return a + b >= MOD ? a + b - MOD : a + b; }
inline int Minus(int a, int b) {return a - b < 0 ? a - b + MOD : a - b; }
inline int ksm(long long a, int b) {
    long long r = 1;
    for(; b; b >>= 1, a = a * a % MOD)
        if(b & 1) r = r * a % MOD;
    return r;
}
int n, m, k, f[N + 1][1 << N], g[N + 1][1 << N];
int inv[N + 10], fac[N + 1], ifac[N + 1];

inline void FWT(int A[], int n, int type) {
    for(int h = 2; h <= n; h <<= 1)
        for(int i = 0; i < n; i += h)
            for(int j = i; j < i + (h >> 1); j ++) {
                if(type == 1) A[j + (h >> 1)] = Plus(A[j + (h >> 1)], A[j]);
                else A[j + (h >> 1)] = Minus(A[j + (h >> 1)], A[j]);
            }
}

inline void exp(int f[], int g[], int n) {
    for(int i = 0; i <= n; i ++) g[i] = 0;
    g[0] = 1;
    for(int i = 1; i <= n; i ++) {
        for(int j = 1; j <= i; j ++)
            g[i] = Plus(g[i], 1ll * f[j] * j % MOD * g[i - j] % MOD);
        g[i] = 1ll * g[i] * inv[i] % MOD;   // inv[i] 是 i 的逆元
    }
}
inline void ln(int f[], int g[], int n) {
    for(int i = 0; i <= n; i ++) g[i] = 0;
    for(int i = 1; i <= n; i ++) {
        for(int j = 1; j < i; j ++)
            g[i] = Plus(g[i], 1ll * g[j] * j % MOD * f[i - j] % MOD);
        g[i] = Minus(f[i], 1ll * g[i] * inv[i] % MOD); // inv[i] 是 i 的逆元
    }
}

inline void power(int f[], int n, int k) {
    static int g[N];
    for(int i = 0; i <= n; i ++) g[i] = 0;
    int pos = 0;
    while(pos <= n && f[pos] == 0) pos ++;
    if(pos > n || pos * k > n) {
        for(int i = 0; i <= n; i ++) f[i] = 0;
        return;
    }
    int d = f[pos], id = ksm(d, MOD - 2);
    for(int i = 0; i <= n; i ++)
        if(i + pos <= n) f[i] = 1ll * f[i + pos] * id % MOD;
        else f[i] = 0;
    ln(f, g, n - pos);
    for(int i = 0; i <= n - pos; i ++)
        f[i] = 1ll * g[i] * k % MOD;
    exp(f, g, n - pos);
    for(int i = n, mul = ksm(d, k); i >= 0; i --)
        if(i - pos * k >= 0) f[i] = 1ll * g[i - pos * k] * mul % MOD;
        else f[i] = 0;
}

int main() {
    ios::sync_with_stdio(false), cin.tie(0);

    cin >> n >> m >> k;
    inv[1] = 1;
    for(int i = 2; i <= n; i ++)
        inv[i] = Minus(0, 1ll * (MOD / i) * inv[MOD % i] % MOD);
    fac[0] = 1; for(int i = 1; i <= n; i ++) fac[i] = 1ll * fac[i - 1] * i % MOD;
    ifac[n] = ksm(fac[n], MOD - 2); for(int i = n; i >= 1; i --) ifac[i - 1] = 1ll * ifac[i] * i % MOD;
    while(m --) {
        int x; cin >> x;
        f[__builtin_popcount(x)][x] ++;
    }

    const int Lim = (1 << n);
    for(int i = 0; i <= n; i ++)
        FWT(f[i], 1 << n, 1);
    for(int mask = 0; mask < Lim; mask ++) {
        static int F[N + 1], G[N + 1], H[N + 1];
        for(int i = 0; i <= n; i ++)
            H[i] = F[i] = f[i][mask], G[i] = 0;
        power(H, n, k);
        for(int i = 0; i <= n; i ++)
            H[i] = 1ll * H[i] * ifac[k] % MOD;
        for(int i = 0, mul = 1; i <= k; i ++) {
            G[0] = Plus(G[0], 1ll * mul * ifac[i] % MOD);
            mul = 1ll * mul * F[0] % MOD;
        }
        for(int i = 1; i <= n; i ++) {
            for(int j = 1; j <= i; j ++)
                G[i] = Plus(G[i], 1ll * j * F[j] % MOD * Minus(G[i - j], H[i - j]) % MOD);
            G[i] = 1ll * G[i] * inv[i] % MOD;
        }
        for(int i = 0; i <= n; i ++) g[i][mask] = G[i];
    }
    FWT(g[n], Lim, - 1);
    cout << g[n][Lim - 1] << '\n';

    return 0;
}