FWT 和集合幂级数
高维前缀和
给定 \(a_0 \cdots a_{2^k-1}\)， 求 \(b_i = \sum_{j \subseteq i} a_j\)。
其实我们可以把这个 \(a\) 数组按二进制位拆分，做一个 \(k\) 维前缀和即可。
代码 \(\text{for} \; 0 \sim 2^n - 1\) ：
for(int i = 0； i < n; i ++){
for(int s = 0; s < (1 << n); s ++) if(s & (1 << i - 1)){
a[s] += a[s ^ (1 << i - 1)];
}
}
同样地，我们还有高维差分的做法：
for(int i = 0； i < n; i ++){
for(int s = 0; s < (1 << n); s ++) if(s & (1 << i - 1)){
a[s] -= a[s ^ (1 << i - 1)]; // 其实就是把 + 改成 -
}
}
例
给定 \(S_i\)，求 \(|\{T | T \cap S \ne \phi\}|\)。
我们考虑正难则反，求 \(\{T | T \cap S = \phi\}\)。
我们可以让 \(dp_s\) 记录有多少单词在当前集合 \(s\) 中，那么就有最终此集合 \(s\) 的补集 \(\bar{s}\) 的答案就为 \(n - dp_s\)。
接下来用 sosdp 就可以求出 \(dp_s\)。
时间复杂度：\(\mathcal O(n \cdot 2^n)\)。
代码
#include<bits/stdc++.h>
using namespace std;
using LL = long long;
int n;
int dp[(1 << 24) + 5];
int main(){
scanf("%d", &n);
for(int i = 1; i <= n; i ++){
string s; cin>>s;
int t = 0;
for(int j = 0; j < 3; j ++){
if(s[j] <= 'x') t |= (1 << (s[j] - 'a'));
}
dp[t] ++;
}
for(int i = 0; i < 24; i ++){
for(int j = 0; j < (1 << 24); j ++) if(j & (1 << i)){
dp[j] += dp[j ^ (1 << i)];
}
}
LL ans = 0;
for(int i = 0; i < (1 << 24); i ++) ans ^= 1ll * (n - dp[i]) * (n - dp[i]);
printf("%lld\n", ans);
return 0;
}
应用-卷积
or 卷积
\[c_{S \cup T} += a_S \times b_T
\]
\[c_W = \sum_{S\cup T = W} a_S \times b_T
\]
且由此，我们做一个前缀和，则会有：
\[c'_W = \sum_{S \cup T \subseteq W} a_S \times b_T
\]
\[= \sum_{S \subseteq W, T \subseteq W} a_S \times b_T
\]
\[= (\sum_{S \subseteq W} a_S)(\sum_{T \subseteq W} b_T)
\]
\[= a' \cdot b'
\]
其中 \(\cdot\) 表示对位相乘，\(a', b'\) 分别为 \(a, b\) 的高维前缀和。
上述过程可以用以下一张表表示：
sosdp
a, b -----------> a', b'
or|and | 对
卷|卷 | 位
积|积 | 相
↓ sosdp ↓ 乘
c ------------> c'
and 卷积
同理。
再看高维前缀和
我们都知道，如果直接做高位前缀和，时间复杂度会来到 \(\mathcal O(4^n)\)，这是我们所万万不能接受的。
考虑使用分治的思想：
定义 \(b = S(a)\)。
设 \(a = (x, y)\)，其中，\(x, y\) 分别表示 \(a\) 的前一半和后一半。
则 \(S(x, y) = (S(x), S(x) + S(y))\)，此时复杂度为 \(\mathcal O(k \cdot 2^k)\)，\(k\) 为维数。
void s(int a[], int n){
if(n == 1) return;
s(a, n / 2), s(a + n / 2, n / 2);
for(int i = 0; i < n / 2; i ++) a[i + n / 2] += a[i];
}
或：
for(int len = 1; len < (1 << n); len <<= 1){
for(int i = 0; i < (1 << n); i += (len << 1)){
for(int j = i; j <= i + len - 1; j ++) a[j + len] += a[j];
}
}
上面那个东西就叫做快速沃尔什变换，即 FWT。
再看卷积
or 卷积
\[(a, b) \times (c, d) = (a \times c, a \times d + b \times c + b \times d)
\]
\[= (a \times c, (a + b) \times (c + d) - a \times c)
\]
时间复杂度 \(\mathcal O(n \log{n})\)。
// 分治版
void mul(int M[], int N[], int n){
// 对位乘法
if(n == 1) return M[0] *= N[0], void();
// 先做加法防止覆盖
for(int i = 0; i < (n >> 1); i ++) M[i + (n >> 1)] += M[i];
for(int i = 0; i < (n >> 1); i ++) N[i + (n >> 1)] += N[i];
// 分治
mul(M, N, n >> 1); mul(M + (n >> 1), N + (n >> 1), n >> 1);
// 还原
for(int i = 0; i < (n >> 1); i ++) M[i + (n >> 1)] -= M[i];
}
或：
inline void OR(LL *a, int tp){
for(int x = 2; x <= (1 << n); x <<= 1){
LL k = x >> 1;
for(LL i = 0; i < (1 << n); i += x) {
for(LL j = 0; j < k; j ++){
(a[i + j + k] += a[i + j] * tp % Mod + Mod) %= Mod;
}
}
}
}
// 笔者推荐这种写法，原因是常数较小
and 卷积
\[(a, b) \times (c, d) = (a \times c + a \times d + b \times c, b \times d)
\]
\[= ((a + b) \times (c + d) - b \times d, b \times d)
\]
inline void AND(LL *a, int tp){
for(int x = 2; x <= (1 << n); x <<= 1){
LL k = x >> 1;
for(LL i = 0; i < (1 << n); i += x) {
for(LL j = 0; j < k; j ++){
(a[i + j] += a[i + j + k] * tp % Mod + Mod) %= Mod;
}
}
}
}
xor 卷积
\[(a, b) \times (c, d) = (a \times c + b \times d, a \times d + b \times c)
\]
\[= (\frac{(a + b) \times (c + d) + (a - b) \times (c - d)}{2}, \frac{(a + b) \times (c + d) - (a - b) \times (c - d)}{2})
\]
inline void XOR(LL *a, LL tp){
for(int x = 2; x <= (1 << n); x <<= 1){
LL k = x >> 1;
for(LL i = 0; i < (1 << n); i += x) {
for(LL j = 0; j < k; j ++){
((a[i + j] += a[i + j + k]) %= Mod);
(a[i + j + k] = a[i + j] - (a[i + j + k] << 1) % Mod) %= Mod;
(a[i + j] *= tp) %= Mod;
(a[i + j + k] *= tp) %= Mod;
}
}
}
}
3 进制 or/max 卷积
\[(x_0, x_1, x_2) \times (y_0, y_1, y_2)
\]
\[= (x_0y_0, x_1y_1+x_1y_0+x_0y_1, x_2y_2+x_2y_1+x_2y_0+x_1y_2+x_0y_2)
\]
\[= (x_0y_0, (x_0+x_1)(y_0+y_1) - x_0y_0, (x_0+x_1+x_2)(y_0+y_1+y_2) - (x_0+x_1)(y_0+y_1)-x_0y_0)
\]
时间复杂度 \(\mathcal O(n\log_3{n}) = \mathcal O(k \cdot 3^k)\)。
xor（不进位加法）卷积
\[(x_0, x_1, x_2) \times (y_0, y_1, y_2)
\]
\[= (x_0y_0+x_1y_2+x_2y_1, \cdots, \cdots)
\]
线性代数
我们发现这是一个线性变换，可以考虑使用矩阵描述。具体地，对于二进制 FWT ，设位矩阵：
\[\begin{bmatrix}\phi(0, 0) & \phi(0, 1) \\ \phi(1, 0) & \phi(1, 1)\end{bmatrix}
\]
下面给出三种基本的矩阵和对应逆矩阵：
| 卷积 | 位矩阵 | 逆矩阵 |
|---|---|---|
| 并卷积 | \(\begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix}\) | \(\begin{bmatrix}1 & 0 \\ -1 & 1\end{bmatrix}\) |
| 交卷积 | \(\begin{bmatrix}1 & 1 \\ 0 & 1\end{bmatrix}\) | \(\begin{bmatrix}1 & -1 \\ 0 & 1\end{bmatrix}\) |
| 对称差卷积（异或卷积） | \(\begin{bmatrix}1 & 1 \\ 1 & -1\end{bmatrix}\) | \(\begin{bmatrix}0. & 0.5 \\ 0.5 & -0.5\end{bmatrix}\) |
由上述视角可知，FWT 具有线性性，即：
\[\text{FWT}(aF + bG) = a \text{FWT}(F) + b\text{FWT}(G)
\]
P4717 【模板】快速莫比乌斯 / 沃尔什变换 (FMT / FWT)
就是板子。
#include<bits/stdc++.h>
using namespace std;
using LL = long long;
const int N = 20, Mod = 998244353;
int n;
LL a[(1 << N) + 1], b[(1 << N) + 1];
LL f[(1 << N) + 1], g[(1 << N) + 1];
inline void OR(LL *a, int tp){
for(int x = 2; x <= (1 << n); x <<= 1){
LL k = x >> 1;
for(LL i = 0; i < (1 << n); i += x) {
for(LL j = 0; j < k; j ++){
(a[i + j + k] += a[i + j] * tp % Mod + Mod) %= Mod;
}
}
}
}
inline void AND(LL *a, int tp){
for(int x = 2; x <= (1 << n); x <<= 1){
LL k = x >> 1;
for(LL i = 0; i < (1 << n); i += x) {
for(LL j = 0; j < k; j ++){
(a[i + j] += a[i + j + k] * tp % Mod + Mod) %= Mod;
}
}
}
}
inline void XOR(LL *a, LL tp){
for(int x = 2; x <= (1 << n); x <<= 1){
LL k = x >> 1;
for(LL i = 0; i < (1 << n); i += x) {
for(LL j = 0; j < k; j ++){
((a[i + j] += a[i + j + k]) %= Mod);
(a[i + j + k] = a[i + j] - (a[i + j + k] << 1) % Mod) %= Mod;
(a[i + j] *= tp) %= Mod;
(a[i + j + k] *= tp) %= Mod;
}
}
}
}
int main(){
scanf("%d", &n);
for(int i = 0; i < (1 << n); i ++) scanf("%lld", &a[i]);
for(int i = 0; i < (1 << n); i ++) scanf("%lld", &b[i]);
memcpy(f, a, sizeof f); memcpy(g, b, sizeof g);
OR(f, 1), OR(g, 1);
for(int i = 0; i < (1 << n); i ++) f[i] = f[i] * g[i] % Mod;
OR(f, -1); for(int i = 0; i < (1 << n); i ++) printf("%lld ", (f[i] + Mod) % Mod);
puts("");
memcpy(f, a, sizeof f); memcpy(g, b, sizeof g);
AND(f, 1), AND(g, 1);
for(int i = 0; i < (1 << n); i ++) f[i] = f[i] * g[i] % Mod;
AND(f, -1); for(int i = 0; i < (1 << n); i ++) printf("%lld ", (f[i] + Mod) % Mod);
puts("");
memcpy(f, a, sizeof f); memcpy(g, b, sizeof g);
XOR(f, 1), XOR(g, 1);
for(int i = 0; i < (1 << n); i ++) f[i] = f[i] * g[i] % Mod;
XOR(f, 499122177); for(int i = 0; i < (1 << n); i ++) printf("%lld ", (f[i] + Mod) % Mod);
puts("");
return 0;
}
P6097 【模板】子集卷积
求 \(c_k = \sum_{\substack{i \& j = 0 \\ i | j = c_k}} a_ib_j\)。
我们发现，这里对于 or 异或的约束多了一个 \(i \& j = 0\)，发现无法直接使用 or 卷积做。
我们考虑在 fwt 数组上多开一维 popcount，我们发现，在这道题中：
\[i\& j = 0 且 i | j = k \iff popcount(i) + popcount(j) = popcount(k) 且 i | j = k
\]
由此，我们分开 popcount 处理就可以 AC。
#include<bits/stdc++.h>
using namespace std;
using LL = long long;
const int N = (1 << 20) + 5, Mod = 1e9 + 9;
int n;
LL a[N], b[N], f[21][N], g[21][N], ans[21][N];
inline void OR(LL *a, int tp){
for(int x = 2; x <= (1 << n); x <<= 1){
LL k = x >> 1;
for(LL i = 0; i < (1 << n); i += x) {
for(LL j = 0; j < k; j ++){
(a[i + j + k] += a[i + j] * tp % Mod + Mod) %= Mod;
}
}
}
}
int main(){
scanf("%d", &n);
for(int i = 0; i < (1 << n); i ++) scanf("%lld", &a[i]);
for(int i = 0; i < (1 << n); i ++) scanf("%lld", &b[i]);
for(int s = 0; s < (1 << n); s ++) f[__builtin_popcount(s)][s] = a[s], g[__builtin_popcount(s)][s] = b[s];
for(int i = 0; i <= n; i ++) OR(f[i], 1), OR(g[i], 1);
for(int s = 0; s < (1 << n); s ++){
for(int k = 0; k <= n; k ++){
for(int i = 0; i <= k; i ++){
(ans[k][s] += f[i][s] * g[k - i][s] % Mod) %= Mod;
}
}
}
for(int i = 0; i <= n; i ++) OR(ans[i], -1);
for(int s = 0; s < (1 << n); s ++) printf("%lld ", ans[__builtin_popcount(s)][s]);
return 0;
}
半在线子集卷积
求 \(H_s = C(s)\cdot \sum_{\substack{T \subsetneq U = S \\ T \cap (S \backslash T) = \phi}} F_T G_{S\backslash T}\)。
这其实就是自己卷自己。（内卷吗？）
对于这种东西，其实我们可以对于每一层单独做 or 卷积 FWT 即可。
对于空集，我们只需要特判一下就可以。
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int N = (1 << 20) + 5, Mod = 998244353;
int n;
int a[21][N], f[21][N];
int g[N], h[N];
inline void FWT(int *a, int n, int tp){
for(int i = 0; i < n; i ++) for(int j = 0; j < (1 << n); j ++) if(j & (1 << i)) (a[j] += Mod + tp * a[j ^ (1 << i)]) %= Mod;
}
signed main(){
scanf("%lld%lld", &n, &f[0][0]);
for(int i = 0; i < (1ll << n); i ++) scanf("%lld", &a[__builtin_popcount(i)][i]);
for(int i = 0; i < (1 << n); i ++) scanf("%lld", &g[i]);
for(int i = 0; i <= n; i ++) FWT(a[i], n, 1);
FWT(f[0], n, 1);
for(int i = 0; i <= n; i ++){
if(i){ // 对于第一层，需要将 g 数组加入卷积。
FWT(f[i], n, -1);
for(int j = 0; j < (1 << n); j ++) f[i][j] = 1ll * f[i][j] * g[j] % Mod;
FWT(f[i], n, 1);
}
for(int j = 1; j <= n - i; j ++){
for(int k = 0; k < (1 << n); k ++){
(f[i + j][k] += 1ll * f[i][k] * a[j][k] % Mod) %= Mod;
}
}
}
for(int i = 0; i <= n; i ++) FWT(f[i], n, -1);
for(int i = 0; i < (1 << n); i ++) printf("%lld ", f[__builtin_popcount(i)][i]);
return 0;
}
集合幂级数
Pre-knowledge——形式幂级数
定义
\[A(x) = \sum_{k = 0}^{\infty} = a_0 + a_1x^1 + a_2x^2 + a_3x^3 + \cdots
\]
在这里，\(x\) 只是形式记号，只是用来区分系数下标。
同时，\(a_k\) 是 \(x^k\) 的系数，记作 \([x^k]A(x) = a_k\)。
定义
我们先写出普通一元幂级数的泰勒公式：
\[\exp(A) = \sum\limits_{k \ge 0} \frac{A^k}{k!}, \;\; \ln(1 + A) = \sum\limits_{k \ge 1}\frac{(-1)^{k+1}A^k}{k}
\]
而在这里，乘法表示的是普通多项式乘法。
现在，我们将乘法重新定义为子集卷积，并作如下强制限制：
- \(\exp(G) = F\)：要求 \([x^\phi](G) = 0\)（即空集系数为 \(0\)）。
- \(\ln(F) = G\)：要求 \([x^\phi]F = 1\)。
乘法规则
\[x^S \cdot x^T = \begin{cases} x^{S \cup T} & S \cap T = \phi \\ 0 & S \cap T \ne \phi \end{cases}
\]
集合幂级数的形式微积分：\(\exp\) 与 \(\ln\)
$ \exp(G)$
定义：
\[\exp(G) = \sum\limits_{k \ge 0}\frac{G^{*k}}{k !}
\]
其中 \(*\) 表示子集卷积；\(G^{*0} = x^0\)（单位元，仅空集系数为 \(1\)）。
组合意义：
$ \exp(G)$ 的系数 \(f_S\) 表示将子集 \(S\) 拆分为任意多个互不相交非空子集，所有拆分方案中，各小块的权值乘积之和。
递推公式：
\[|U| \cdot f_U = \sum\limits_{\substack{A \cup B = U \\ A \cap B = \phi, A \ne \phi}} |A|\cdot g_A \cdot f_B
\]
由此，分层卷积可求。
\(\ln(F)\)
定义：
满足 \(\exp(\ln(F)) = F\)，为 \(\exp\) 的逆运算。
组合意义：
\(g_S\) 表示子集 \(S\) 不可再拆分成多个非空不交子集的权值，即连通块生成函数。
递推公式：
\[k \cdot g_U = f_U - \sum\limits_{\substack{A \cup B = U \\ A \cap B = \phi, A \ne \phi}} |A| \cdot g_A \cdot f_B
\]
例题-集合幂级数 exp
让我们对 \(\exp(A)\) 求导：
\[\exp'(A(x)) = exp(A(X))A'(x)
\]
\[B_{n+1}(n+1) = B(x)A'(x)
\]
\[B_{n+1}(n+1)=\sum_{i=0}^nB_iA_{n-i+1}(n-i+1)
B_n=\frac{1}{n}\sum_{i=0}^nA_iB_{n-i}i\]
预处理逆元后 FWT 即可，具体见代码：
#include<bits/stdc++.h>
using namespace std;
using LL = long long;
const int N = (1 << 20) + 2;
const int Mod = 998244353;
inline LL qp(LL x, int y){
LL res = 1ll;
while(y){
if(y & 1) (res *= x) %= Mod;
(x *= x) %= Mod, y >>= 1;
}
return res;
}
int n;
LL a[21][N];
LL inv[N];
inline void FWT(LL *a, int tp){
for(int len = 1; len < (1 << n); len <<= 1){
for(int i = 0; i < (1 << n); i += (len << 1)){
for(int j = i; j <= i + len - 1; j ++) a[j + len] = (a[j + len] + 1ll * tp * a[j] + Mod) % Mod;
}
}
}
inline void ln(LL *a, LL *b){
b[0] = 0;
for(int i = 1; i <= n; i ++){
LL sum = 0;
for(int j = 1; j < i; j ++) (sum += 1ll * j * b[j] % Mod * a[i - j] % Mod) %= Mod;
b[i] = (a[i] + (Mod - sum) * inv[i]) % Mod;
}
}
inline void exp(LL *a, LL *b){
b[0] = 1;
for(int i = 1; i <= n; i ++){
LL sum = 0;
for(int j = 0; j <= i; j ++) (sum += 1ll * j * a[j] % Mod * b[i - j] % Mod) %= Mod;
(sum += Mod) %= Mod;
b[i] = sum * inv[i] % Mod;
}
}
signed main(){
scanf("%d", &n);
for(int i = 1; i <= n; i ++) inv[i] = qp(1ll * i, Mod - 2);
for(int i = 0; i < (1 << n); i ++) scanf("%lld", &a[__builtin_popcount(i)][i]);
for(int i = 0; i <= n; i ++) FWT(a[i], 1);
for(int s = 0; s < (1 << n); s ++){
LL A[21], B[21];
for(int i = 0; i <= n; i ++) A[i] = a[i][s];
//ln(A, B);
exp(A, B);
for(int i = 0; i <= n; i ++) a[i][s] = B[i];
}
for(int i = 0; i <= n; i ++) FWT(a[i], -1);
for(int i = 0; i < (1 << n); i ++) printf("%lld ", a[__builtin_popcount(i)][i]);
return 0;
}
例二【模板】集合幂级数 \(\ln\) - 洛谷
依旧求导：
\[\ln'(A(x)) = \frac{1}{A(x)}A'(x)
\]
\[\ln'(A(x))A(x) = A'(x)
\]
\[\sum_{i=0}^nA_{n-i}B_{i+1}(i+1)=A_{i+1}(i+1)
B_n=A_n-\frac{1}{n}\sum_{i=0}^{n-1}B_iA_{n-i}i\]
依旧代码：
#include<bits/stdc++.h>
using namespace std;
using LL = long long;
const int N = (1 << 20) + 2;
const int Mod = 998244353;
inline LL qp(LL x, int y){
LL res = 1ll;
while(y){
if(y & 1) (res *= x) %= Mod;
(x *= x) %= Mod, y >>= 1;
}
return res;
}
int n;
LL a[21][N];
LL inv[N];
inline void FWT(LL *a, int tp){
for(int len = 1; len < (1 << n); len <<= 1){
for(int i = 0; i < (1 << n); i += (len << 1)){
for(int j = i; j <= i + len - 1; j ++) a[j + len] = (a[j + len] + 1ll * tp * a[j] + Mod) % Mod;
}
}
}
inline void ln(LL *a, LL *b){
b[0] = 0;
for(int i = 1; i <= n; i ++){
LL sum = 0;
for(int j = 1; j < i; j ++) (sum += 1ll * j * b[j] % Mod * a[i - j] % Mod) %= Mod;
b[i] = (a[i] + (Mod - sum) * inv[i]) % Mod;
}
}
inline void exp(LL *a, LL *b){
b[0] = 1;
for(int i = 1; i <= n; i ++){
LL sum = 0;
for(int j = 0; j <= i; j ++) (sum += 1ll * j * a[j] % Mod * b[i - j] % Mod) %= Mod;
(sum += Mod) %= Mod;
b[i] = sum * inv[i] % Mod;
}
}
signed main(){
scanf("%d", &n);
for(int i = 1; i <= n; i ++) inv[i] = qp(1ll * i, Mod - 2);
for(int i = 0; i < (1 << n); i ++) scanf("%lld", &a[__builtin_popcount(i)][i]);
for(int i = 0; i <= n; i ++) FWT(a[i], 1);
for(int s = 0; s < (1 << n); s ++){
LL A[21], B[21];
for(int i = 0; i <= n; i ++) A[i] = a[i][s];
ln(A, B);
// exp(A, B);
for(int i = 0; i <= n; i ++) a[i][s] = B[i];
}
for(int i = 0; i <= n; i ++) FWT(a[i], -1);
for(int i = 0; i < (1 << n); i ++) printf("%lld ", a[__builtin_popcount(i)][i]);
return 0;
}
浙公网安备 33010602011771号