【笔记】集合幂级数

这篇笔记可以说从今年年初就开始动笔，然后烂尾了几次又断断续续写，然后还经历了几次完全看不懂之前写的内容的情况（雾

上班上多了感觉脑子都不灵光了。

这篇笔记几乎没什么前置知识，如果有可能就是二维前缀和吧。

高维前缀和 & 高维差分

考虑朴素的不使用容斥的二维前缀和，可以写成

for (int i = 1; i <= n; i++) {
    for (int j = 1; j <= n; j++) {
        a[i][j] += a[i - 1][j];
    }
}
for (int i = 1; i <= n; i++) {
    for (int j = 1; j <= n; j++) {
        a[i][j] += a[i][j - 1];
    }
}

进一步，下面是三维前缀和

for (int i = 1; i <= n; i++) {
    for (int j = 1; j <= n; j++) {
        for (int k = 1; k <= n; k++) { 
            a[i][j][k] += a[i - 1][j][k];
        }
    }
}
for (int i = 1; i <= n; i++) {
    for (int j = 1; j <= n; j++) {
        for (int k = 1; k <= n; k++) {
            a[i][j][k] += a[i][j - 1][k];
        }
    }
}
for (int i = 1; i <= n; i++) {
    for (int j = 1; j <= n; j++) {
        for (int k = 1; k <= n; k++) {
            a[i][j][k] += a[i][j][k - 1];
        }
    }
}

假如按这种方式写 \(m\) 维前缀和，那就有 \(m\) 个循环体，每个循环体 \(m\) 层循环，复杂度 \(O(m\cdot n^m)\) 。

这就是高维前缀和了，只不过，一般我们遇到的是 “每一维长度均为 \(2\)” 的情况，在上述例子中也就是 \(n=2\)。此时，可以状态压缩写成下面这样，二进制下每一位都代表一个维度。

for (int j = 0; j < m; j++) { // 遍历每一维
    for (int i = 0; i < (1 << m); i++) {
        if (i >> j & 1) {
            a[i] += a[i ^ (1 << j)];
        }
    }
}

接下来，考虑高维前缀和后 \(\{a\}\) 的意义。以下将高维前缀和后的 \(\{a\}\) 称作 \(\{f\}\)。

• 一维前缀和中，\(f[x]=\sum_{i\le x} a[i]\)。
• 二维前缀和中，\(f[x][y]=\sum_{i\le x,j\le y}a[i][j]\)。

因此，\(m\) 维前缀和后，\(f[x]\)（这里 \(x\) 是 \(m\) 位二进制数）就表示满足 \(\forall i\in [0,m)\)，\(y_i\le x_i\)（\(y\) 二进制下第 \(i\) 位 $$ \(x\) 二进制下第 \(i\) 位）的所有 \(a[y]\) 的和，即 \[ f[S]=\sum_{T\subseteq S}a[T] \] 既然高维前缀和能通过 \(a\) 算 \(f\)，那有办法能通过 \(f\) 算 \(a\) 吗？

有的，与之相对的就是高维差分。这就好比一维前缀和是 \(f[i]=f[i-1]+a[i]\)，想反过来算出 \(a\) 就是 \(a[i]=f[i]-f[i-1]\)。

代码写出来只用改一个正负号。

for (int j = 0; j < m; j++) {
    for (int i = 0; i < (1 << m); i++) {
        if (i >> j & 1) {
            a[i] -= a[i ^ (1 << j)];
        }
    }
}

这段代码实际上隐含了容斥原理。手玩一下 \(m=2\) 就会发现，所有的 a[i] -= a[i ^ (1 << j)] 合并在一起形如 \(a[11] = f[11] - f[10] -f[01] + f[00]\)，每个元素的正负号在反复相减中变化。

所以高维差分的式子写出来就是容斥的样子。 \[ a[S]=\sum_{T\subseteq S}(-1)^{|S|-|T|}f[T] \] 上述从 \(\{a\}\) 到 \(\{f\}\) 的操作也被称为 “子集和” 或 “快速莫比乌斯变换” 或 “Zeta 变换”，从 \(\{f\}\) 到 \(\{a\}\) 的过程也被称作 “子集反演” 或 “快速莫比乌斯反演”。

超集和 & 超集差分

和子集和类似，定义超集和是 \[ f[S]=\sum_{S\subseteq T}a[T] \] 即所有包含 \(S\) 的状态 \(T\) 的 \(a[T]\) 之和。

对比子集和，代码上唯一的区别在于：当我们在处理第 \(j\) 维时，要把 “坐标为 \(1\)” 的值加到 “坐标为 \(0\)” 的位置上。也就是

for (int j = 0; j < m; j++) {
    for (int i = 0; i < (1 << m); i++) {
        if (!(i >> j & 1)) {
            a[i] += a[i ^ (1 << j)];
        }
    }
}

同理如果想要通过 \(f\) 还原 \(a\)，就是

if (!(i >> j & 1)) {
    a[i] -= a[i ^ (1 << j)];
}

OR 卷积

定义子集和变换 \(\hat f=\text{FMT}(f)\)，为 \[ \hat f_i=\sum\limits_{j\subseteq i}f_j \]

其逆变换为 \[ f_i=\text{IFMT}(\hat f)_i=\sum\limits_{j\subseteq i}(-1)^{|i|-|j|}\hat f_j \] 也就是子集反演。

所谓 OR 卷积，就是 \[ h_k=\sum_{i|j=k}f_i\cdot g_j \] 要做 \(f\) 和 \(g\) 的 OR 卷积，我们只要把 \(f\) 和 \(g\) 分别做 \(\text{FMT}\)，点乘起来，再做一遍 \(\text{IFMT}\) 即可。 \[ h=\text{IFMT}(\text{FMT}(f)\cdot \text{FMT}(g)) \] 这个不难理解，因为 \[ \sum_{i \subseteq S}\sum_{j \subseteq S} f_i g_j = \sum_{i|j \subseteq S} f_i g_j = \sum_{k\subseteq S}\sum_{i|j=k} f_i g_j = \sum_{k\subseteq S}h_k \] 所以变换后点乘得到的是 \(h\) 的子集和，再做 IFMT 就能得到 OR 卷积。

代码在下方 AND 卷积一同给出。

AND 卷积

同样的，注意到 \[ \sum_{S \subseteq i}\sum_{S \subseteq j} f_i g_j = \sum_{S \subseteq (i\& j)} f_i g_j = \sum_{S\subseteq k}\sum_{i\&j=k} f_i g_j = \sum_{S\subseteq k}h_k \] 所以变换后点乘得到的是 \(h\) 的超集和。换句话说，点乘后的数组在状态 \(S\) 上存的是 \(\sum_{S\subseteq k}h_k\)，也就是 \(h\) 的超集和。因此我们把子集变换改成超集变换就行。

定义 \(\text{FMT}_{\text{sup}}\) 和 \(\text{IFMT}_{\text{sup}}\) 分别表示超集和以及超集反演，那么 \[ h=\text{IFMT}_{\text{sup}}(\text{FMT}_{\text{sup}}(f)\cdot \text{FMT}_{\text{sup}}(g)) \]

// op = 0 表示做 OR 卷积, op = 1 表示做 AND 卷积
vector<int> conv(vector<int> a, vector<int> b, int op) {
    auto fmt = [&] (vector<int> &c, int inv) {
        int N = c.size();
        for (int j = 1; j < N; j <<= 1) {
            for (int i = 0; i < N; i++) {
                if ((op == 0 && i & j) || (op == 1 && !(i & j))) {
                    if (inv) {
                        c[i] = (c[i] - c[i ^ j] + mod) % mod;
                    } else {
                        c[i] = (c[i] + c[i ^ j]) % mod;
                    }
                }
            }
        }
    };
    fmt(a, 0);
    fmt(b, 0);
    for (int i = 0; i < a.size(); i++) {
        a[i] = a[i] * b[i] % mod;
    }
    fmt(a, 1);
    return a;
}

分治写法

除了上面基于高维前缀和的迭代写法，分治写法也比较常见。

以 OR 卷积为例，设当前数组 \(\{a\}\) 长度为 \(2^n\)，我们可以平分成两半：

• \(a_0\)：前一半元素，下标最高位为 \(0\)。
• \(a_1\)：后一半元素，下标最高位为 \(1\)。

左半部分的状态最高位为 \(0\)，它的子集仍在左半部分；右半部分的状态最高位为 \(1\)，它的子集分成两类：不含最高位的子集在左半部分，含最高位的子集在右半部分。 \[ \text{FMT}(a)=\text{merge}(\text{FMT}(a_0),\text{FMT}(a_0)+\text{FMT}(a_1)) \] 其中 \(\text{merge}\) 表示数组首尾相连。于是我们只用把右半部分加上左半部分就行。

void FMT_OR(vector<int> &f) {
    int N = f.size();
    for (int o = 2, k = 1; o <= N; o <<= 1, k <<= 1) {
        for (int i = 0; i < N; i += o) {
            for (int j = 0; j < k; j++) {
                f[i + j + k] = (f[i + j + k] + f[i + j]) % mod;
            }
        }
    }
}

IFMT 同理，改成 f[i + j + k] -= f[i + j]; 即可。写在一起就是

void conv_OR(vector<int> &f, int v) {
    int N = f.size();
    for (int o = 2, k = 1; o <= N; o <<= 1, k <<= 1) {
        for (int i = 0; i < N; i += o) {
            for (int j = 0; j < k; j++) {
                f[i + j + k] = (f[i + j + k] + f[i + j] * v + mod) % mod;
            }
        }
    }
}
vector<int> conv_OR(vector<int> f, vector<int> g) {
    conv_OR(f, 1);
    conv_OR(g, 1);
    for (int i = 0; i < f.size(); i++) {
        f[i] = f[i] * g[i] % mod;
    }
    conv_OR(f, -1);
    return f;
}

AND 卷积同理，子集变超集，也就是把 f[i + j + k] 和 f[i + j] 换个位置。

void conv_AND(vector<int> &f, int v) {
    int N = f.size();
    for (int o = 2, k = 1; o <= N; o <<= 1, k <<= 1) {
        for (int i = 0; i < N; i += o) {
            for (int j = 0; j < k; j++) {
                f[i + j] = (f[i + j] + f[i + j + k] * v + mod) % mod;
            }
        }
    }
}
vector<int> conv_AND(vector<int> f, vector<int> g) {
    conv_AND(f, 1);
    conv_AND(g, 1);
    for (int i = 0; i < f.size(); i++) {
        f[i] = f[i] * g[i] % mod;
    }
    conv_AND(f, -1);
    return f;
}

快速沃尔什变换 & XOR 卷积

快速沃尔什变换（FWT），通常指基于 Hadamard 矩阵的变换，这个之后再提。

总之最常见的用法是 XOR 卷积，即求 \[ h_k=\sum_{i\oplus j=k}f_i\cdot g_j \] 基于和 FMT 类似的思想，我们想要构造变换 \(\text{FWT}(f)\) 使得 \[ h=\text{IFWT}(\text{FWT}(f)\cdot \text{FWT}(g)) \]

还是考虑分治。设当前数组长度为 \(2^n\)，把最高位为 \(0\) 和最高位为 \(1\) 的部分分开，记作 \(f_0,f_1\)。

记结果的两半为 \(h_0,h_1\)。如果 \(k\) 的最高位为 \(0\)，那么 \(i,j\) 的最高位相同；如果 \(k\) 的最高位为 \(1\)，那么 \(i,j\) 的最高位不同。

固定除最高位外的低位部分为 \(t\)，有 \[ h_0[t]=\sum_{i\oplus j=t}f_0[i]g_0[j]+\sum_{i\oplus j=t}f_1[i]g_1[j] \] \[ h_1[t]=\sum_{i\oplus j=t}f_0[i]g_1[j]+\sum_{i\oplus j=t}f_1[i]g_0[j] \]

仔细观察这两个式子，如果我们把 \(f_0\) 和 \(f_1\) 先变成 \(f_0+f_1\)，把 \(g_0\) 和 \(g_1\) 先变成 \(g_0+g_1\)，那么它们低一维 XOR 卷积的结果就是 \(h_0+h_1\)；同理，把两边分别变成 \(f_0-f_1\) 和 \(g_0-g_1\)，低一维 XOR 卷积的结果就是 \(h_0-h_1\)。

综上，可以定义一层 FWT 变换为 \[ (x,y)\rightarrow (x+y,x-y) \] 反变换则为 \[ (x,y)\rightarrow \left(\frac{x+y}{2},\frac{x-y}{2}\right) \]

代码如下。

void FWT_XOR(vector<int> &f, int v) {
	int N = f.size();
	for (int o = 2, k = 1; o <= N; o <<= 1, k <<= 1) {
		for (int i = 0; i < N; i += o) {
			for (int j = 0; j < k; j++) {
				int x = f[i + j], y = f[i + j + k];
				f[i + j] = (x + y) % mod;
				f[i + j + k] = (x - y + mod) % mod;
			}
		}
	}
	if (v) {
		int invN = qpow(N, mod - 2);
		for (int &x : f) {
			x = x * invN % mod;
		}
	}
}
vector<int> conv_XOR(vector<int> f, vector<int> g) {
    FWT_XOR(f, 0);
    FWT_XOR(g, 0);
    for (int i = 0; i < f.size(); i++) {
        f[i] = f[i] * g[i] % mod;
    }
    FWT_XOR(f, 1);
    return f;
}

从另一个角度看，FWT 实际上是在用下面这个容斥函数做变换： \[ \hat f[S]=\sum_T(-1)^{|S\cap T|}f[T] \] 它满足 \[ (-1)^{|S\cap A|}(-1)^{|S\cap B|}=(-1)^{|S\cap (A\oplus B)|} \] 所以 XOR 卷积在变换后就变成了点乘。

至此，OR / AND / XOR 三类按位卷积都可以写成 “变换、点乘、反变换” 的形式。

接下来考虑另一种集合上的乘法。前面的 OR 卷积要求 \(i|j=S\)，允许 \(i\) 和 \(j\) 有交；但很多问题里，我们希望把一个集合 \(S\) 拆成两个不相交的部分 \(T\) 和 \(S\setminus T\)，分别从 \(f\) 和 \(g\) 中取贡献。为了刻画这种 “不交合并” 的乘法，可以引入集合幂级数。

集合幂级数

给全集中的每个元素 \(i\) 对应一个变量 \(x_i\)。对于集合 \(S\)，记 \[ x^S=\prod_{i\in S}x_i \] 也就是说，集合 \(S\) 对应的单项式就是把 \(S\) 中所有元素对应的变量乘起来。因此集合幂级数可以写成 \[ f=\sum_S f[S]x^S \] 其中乘法规定为 \[ x^A x^B=\begin{cases} x^{A\cup B} & A\cap B=\varnothing\\ 0 & A\cap B\ne \varnothing \end{cases} \] 由乘法定义可知，一旦两个单项式包含相同变量，乘积就是 \(0\)。特别地，\(x_i^2=0\)。因此每个非零单项式中每个变量最多出现一次，所有级数都是有限的，最高次数不超过全集的元素个数 \(n\)。

子集卷积

由上面的乘法定义，两个集合幂级数 \(f,g\) 的乘积 \(h=fg\) 满足 \[ h[S]=\sum_{T\subseteq S}f[T]g[S\setminus T] \] 这也叫子集卷积。

注意它和 OR 卷积很像，但并不一样。OR 卷积枚举的是 \(i|j=S\)，允许 \(i\) 和 \(j\) 有交；而子集卷积要求两边不交，且并起来正好是 \(S\)。

一个朴素写法是枚举子集，复杂度 \(O(3^n)\)。

vector<int> conv_subset_slow(vector<int> f, vector<int> g) {
    int N = f.size();
    vector<int> h(N);
    for (int s = 0; s < N; s++) {
        for (int t = s; ; t = (t - 1) & s) {
            h[s] = (h[s] + f[t] * g[s ^ t]) % mod;
            if (t == 0) break;
        }
    }
    return h;
}

想优化它，需要额外记录集合大小。令 \[ F_i[S]=\begin{cases} f[S] & |S|=i\\ 0 & |S|\ne i \end{cases} \] 对每个 \(i\) 分别做子集和。这样变换后，固定一个 \(S\)，所有 \(T\subseteq S\) 都被加进来了。

如果 \(T\) 和 \(S\setminus T\) 不交，那么大小满足 \[ |T|+|S\setminus T|=|S| \] 所以我们可以在大小这一维上做普通卷积。

具体地：

1. 按照 popcount 分层得到 \(F_i,G_i\)。
2. 对每一层做 FMT。
3. 对每个集合 \(S\)，计算 \[ H_k[S]=\sum_{i=0}^k F_i[S]G_{k-i}[S] \]
4. 对每一层做 IFMT。
5. 答案是 \(H_{|S|}[S]\)。

代码如下，复杂度 \(O(n^2 2^n)\)。

void FMT(vector<int> &a, int inv) {
    int N = a.size();
    for (int j = 1; j < N; j <<= 1) {
        for (int i = 0; i < N; i++) {
            if (i & j) {
                if (inv) {
                    a[i] = (a[i] - a[i ^ j] + mod) % mod;
                } else {
                    a[i] = (a[i] + a[i ^ j]) % mod;
                }
            }
        }
    }
}

vector<int> conv_subset(vector<int> f, vector<int> g) {
    int N = f.size(); // N 是 2 的幂
    int n = __builtin_ctz(N);
    vector F(n + 1, vector<int>(N));
    vector G(n + 1, vector<int>(N));
    vector H(n + 1, vector<int>(N));

    for (int s = 0; s < N; s++) {
        int c = __builtin_popcountll(s);
        F[c][s] = f[s];
        G[c][s] = g[s];
    }
    for (int i = 0; i <= n; i++) {
        FMT(F[i], 0);
        FMT(G[i], 0);
    }
    for (int s = 0; s < N; s++) {
        for (int i = 0; i <= n; i++) {
            for (int j = 0; i + j <= n; j++) {
                H[i + j][s] = (H[i + j][s] + F[i][s] * G[j][s]) % mod;
            }
        }
    }
    for (int i = 0; i <= n; i++) {
        FMT(H[i], 1);
    }

    vector<int> h(N);
    for (int s = 0; s < N; s++) {
        h[s] = H[__builtin_popcountll(s)][s];
    }
    return h;
}

为什么这样能去掉 “不交” 的限制？

FMT 后的 \(F_i[S]\) 表示 \(S\) 的所有大小为 \(i\) 的子集 \(T\) 的 \(f[T]\) 之和。因此点乘得到的 \(H'_k[S]\) 会统计所有 \(A,B\subseteq S\) 且 \(|A|+|B|=k\) 的贡献。

这些贡献可以按 \(A\cup B\) 分组，所以 \(H'_k[S]\) 实际上是某个数组 \(H_k\) 的子集和： \[ H'_k[S]=\sum_{U\subseteq S}H_k[U] \] 其中 \(H_k[U]\) 表示所有满足 \(A\cup B=U\) 且 \(|A|+|B|=k\) 的贡献。因此对 \(H'_k\) 做 IFMT 后，\(H_k[S]\) 只保留 \(A\cup B=S\) 的贡献。

此时再取 \(k=|S|\)，由于 \(|A|+|B|=|A\cup B|\) 当且仅当 \(A\cap B=\varnothing\)，所以 \(H_{|S|}[S]\) 正好是子集卷积的答案。

集合幂级数运算

逆元

若 \(f[\varnothing]\ne 0\)，则 \(f\) 有乘法逆元 \(g=f^{-1}\)，满足 \[ f g=1 \] 其中 \(1\) 表示 \(a[\varnothing]=1\)，其余位置为 \(0\)。

对非空集合 \(S\)，有 \[ \sum_{T\subseteq S}f[T]g[S\setminus T]=0 \] 把 \(T=\varnothing\) 的项单独拿出来，就得到 \[ g[S]=-\frac{1}{f[\varnothing]}\sum_{\varnothing\ne T\subseteq S}f[T]g[S\setminus T] \]

朴素递推复杂度也是 \(O(3^n)\)。

// f[0] != 0
vector<int> inv_set(vector<int> f) {
    int N = f.size();
    vector<int> g(N);
    int inv_f0 = qpow(f[0], mod - 2);
    g[0] = inv_f0;
    for (int s = 1; s < N; s++) {
        int sum = 0;
        for (int t = s; t; t = (t - 1) & s) {
            sum = (sum + f[t] * g[s ^ t]) % mod;
        }
        g[s] = (mod - sum) * inv_f0 % mod;
    }
    return g;
}

exp

在普通的形式幂级数中，\(\text{exp}\) 的定义直接源于微积分里的泰勒展开 \[ e^x=\sum_{k=0}^{\infty} \frac{x^k}{k!} \] 我们如果将自变量 \(x\) 替换为集合幂级数 \(f\)，就得到了集合幂级数 \(\text{exp}\) 的定义 \[ \exp(f)=\sum_{k\ge 0}\frac{f^k}{k!} \] 注意这里必须要求 \(f[\varnothing]=0\)。

由于集合幂级数最高只有 \(n\) 次，所以这个和实际上到 \(n\) 就结束了。

设 \(g=\exp(f)\)。为了方便推导，我们定义一个类似求导的 \('\) 操作，表示 “每一项乘上它自己包含的元素个数（集合大小）”：若 \[ F=\sum_S F[S]x^S \] 则 \[ F'=\sum_S |S|F[S]x^S \] 类比链式求导法则，有 \(g'=f'g\)，据此可以得到递推式 \[ |S|g[S]=\sum_{\varnothing\ne T\subseteq S}|T|f[T]g[S\setminus T] \] 所以 \[ g[S]=\frac{1}{|S|}\sum_{\varnothing\ne T\subseteq S}|T|f[T]g[S\setminus T] \]

// f[0] == 0
vector<int> exp_set(vector<int> f) {
    int N = f.size();
    vector<int> g(N);
    g[0] = 1;
    for (int s = 1; s < N; s++) {
        int sum = 0;
        for (int t = s; t; t = (t - 1) & s) {
            sum = (sum + __builtin_popcountll(t) * f[t] % mod * g[s ^ t]) % mod;
        }
        g[s] = sum * qpow(__builtin_popcountll(s), mod - 2) % mod;
    }
    return g;
}

ln

同样地，微积分中 \(\ln(1+x)\) 在 \(x=0\) 处的泰勒展开为 \[ \ln(1+x)=\sum_{k\ge 1}\frac{(-1)^{k-1}x^k}{k} \] 令 \(1+x=f\)，带入上式就得到了集合幂级数 \(\ln(f)\) 的定义 \[ \ln(f)=\sum_{k\ge 1}\frac{(-1)^{k-1}(f-1)^k}{k} \]

同样这里要求 \(f[\varnothing]=1\)。

设 \(g=\ln(f)\)。沿用上面对 \('\) 的定义，两边求导得到 \(g'=\frac{f'}{f}\)，即 \[ f'=g'f \] 展开提取系数 \[ |S|f[S]=\sum_{\varnothing\ne T\subseteq S}|T|g[T]f[S\setminus T] \] 我们要算的是 \(g[S]\)，所以把 \(T=S\) 的项单独拿出来 \[ |S|f[S]=|S|g[S]f[\varnothing] + \sum_{\varnothing\ne T\subsetneq S}|T|g[T]f[S\setminus T] \]

因为前提条件是 \(f[\varnothing]=1\)，所以再整理一下，得到 \[ g[S]=f[S]-\frac{1}{|S|}\sum_{\varnothing\ne T\subsetneq S}|T|g[T]f[S\setminus T] \]

// f[0] == 1
vector<int> ln_set(vector<int> f) {
    int N = f.size();
    vector<int> g(N);
    for (int s = 1; s < N; s++) {
        int sum = 0;
        for (int t = (s - 1) & s; t; t = (t - 1) & s) {
            sum = (sum + __builtin_popcountll(t) * g[t] % mod * f[s ^ t]) % mod;
        }
        int c = __builtin_popcountll(s);
        g[s] = (f[s] - sum * qpow(c, mod - 2)) % mod;
        if (g[s] < 0) g[s] += mod;
    }
    return g;
}

pow

有了 \(\exp\) 和 \(\ln\)，算集合幂级数的幂 \(f^k\)，就很简单了。 \[ f^k=\exp(k\ln f) \]

// f[0] == 1
vector<int> pow_set(vector<int> f, int k) {
    f = ln_set(f);
    for (int i = 0; i < f.size(); i++) {
        f[i] = f[i] * k % mod;
    }
    return exp_set(f);
}

注：上面给出的集合幂级数 exp, ln 和 pow 的代码，都是 \(O(3^n)\)。

但实际上，利用子集卷积的思想，都能优化到 \(O(n^2 2^n)\)。

但我脑子已经不够用了，就这样吧.jpg