专栏文章

快速 Walsh 变换 / Mobius 变换简介

铃铃悬
算法·理论参与者 40已保存评论 42

文章操作

快速查看文章及其快照的属性，并进行相关操作。

当前评论
42 条
当前快照
1 份
快照标识符
@mhz5rv0v
此快照首次捕获于
2025/11/15 01:55
4 个月前
此快照最后确认于
2025/11/29 05:24
3 个月前
查看原文
\gdef\and{\mathbin{\&}} \gdef\or{\mathbin{|}} \gdef\xor{\mathbin{\mathrm{xor}}} \gdef\abs#1{\lvert #1 \rvert} \gdef\cF{\mathcal{F}} \gdef\dd{\mathord{.}}
在本篇文章中, 我将介绍两种算法: 快速 Walsh 变换与快速 Mobius 变换.

引子

我们时不时会遇到一些位运算或者集合运算问题. 这两者在很大意义上有共同之处: 如果我们把一个二进制数 m=∑jaj2j−1m = \sum_j a_j 2^{j - 1} 同集合 {j∣aj=1}\{ j \mid a_j = 1 \} 对应起来, 例如把 19=24+21+2019 = 2^4 + 2^1 + 2^0 对应到 {1,2,5}\{ 1, 2, 5 \}. 那么二进制数的按位与, 或, 非, 异或, 恰好对应集合的交, 并, 补, 对称差.
两集合 A,BA, B 的对称差定义为恰好属于两者中的某一个集合的元素所构成的集合. 即 {x∣x∈A,x∉B 或者 x∉A,x∈B}\{ x \mid x \in A, x \notin B \text{ 或者 } x \notin A, x \in B \}.
本文中, 我们对这两种表示不加区分. 但是为了不混淆, 我们会用 i,ji, j 表示位数, x,y,zx, y, z 表示二进制数. 其中位数是从 00 开始的, 也就是说第 ii 位的 11 就是 2i2^i. 例如说, 我们会说 x⊆yx \subseteq y 表示 xx 对应的集合是 yy 对应的集合的子集, 或者说 x\andy=yx \and y = y. 又比如说, 我们会说 xx 包含 ii 这一位, 或者 i∈xi \in x, 就是说 xx 的二进制表示里第 ii 位为 11. 再比如说, 我们会用 \absx\abs{x} 表示 xx 二进制中 11 的个数 (即所谓的 popcount⁡(x)\operatorname{popcount}(x)).

问题

快速 Walsh 变换和快速 Mobius 变换处理的是某种"位运算卷积". 设 nn 是正整数, 我们接下来只考虑 nn 位二进制数. 记 N=2nN = 2^n. A=(Ax)0≤a<NA = (A_x)_{0 \leq a < N} 是一个长为 NN 的数组, B=(Bx)B = (B_x) 同理.
我们定义
Cz=∑x⊕y=zAxBy=∑x=0N−1∑y=0N−1[x⊕y=z]AxBy, C_z = \sum_{x \oplus y = z} A_x B_y = \sum_{x = 0}^{N - 1} \sum_{y = 0}^{N - 1} [x \oplus y = z] A_x B_y,
将 C=(Cz)C = (C_z) 称为 AA 和 BB 的 ⊕\oplus-卷积. 其中 ⊕\oplus 可能是 \and,\or,\xor\and, \or, \xor 之一. 在不引起混淆的时候 (⊕\oplus 明确的时候), 我们记 C=A∗BC = A \ast B.
此外, 我们还记 A⋅BA \cdot B 和 A+BA + B 分别表示逐项相加和逐项相乘.
为了解决这个问题, 我们有多种想法.

容斥原理与高维前缀和

我们从或卷积开始. 那么就有 Cz=∑x\ory=zAxBy.C_z = \sum_{x \or y = z} A_x B_y.
想象如何计算 Cz=∑x\ory=zAxByC_z = \sum_{x \or y = z} A_x B_y.
可以发现, 如果 x\ory=zx \or y = z, 那么一定有 x,y⊆zx, y \subseteq z. 因此我们可以先把所有 x,y⊆zx, y \subseteq z 这样的 AxByA_x B_y 都加起来 (也就是 (∑x⊆zAx)(∑y⊆zBy)(\sum_{x \subseteq z} A_x)(\sum_{y \subseteq z} B_y)), 然后对每个 i∈zi \in z, 减掉那些 xx 和 yy 都不包含第 ii 位时的 AxByA_x B_y.
但是这样的话, 如果 xx 和 yy 都不包含的位有两位, 我们就重复加了, 所以还需要再减掉. 然后还需要加上都不包含某三位的; 依此类推.
如果把这个式子写开 (不难发现它其实就是容斥原理), 就是
Cz=∑w⊆z(−1)\absz−w(∑x⊆wAx)(∑y⊆wBy). C_z = \sum_{w \subseteq z} (-1)^{\abs{z - w}} \bigl(\sum_{x \subseteq w} A_x\bigr) \bigl(\sum_{y \subseteq w} B_y \bigr).
如果我们不做容斥, 那么上式等价于说
∑w⊆zCw=(∑x⊆wAx)(∑y⊆wBy). \sum_{w \subseteq z} C_w = \bigl(\sum_{x \subseteq w} A_x\bigr) \bigl(\sum_{y \subseteq w} B_y \bigr).
这件事情是容易的: x\oryx \or y 是 zz 的子集当且仅当 xx 和 yy 都是 zz 的子集.
类似地, 如果 CC 是 AA 和 BB 的与卷积, 就有
∑z⊆wCw=(∑w⊆xAx)(∑w⊆yBy). \sum_{z \subseteq w} C_w = \bigl(\sum_{w \subseteq x} A_x\bigr) \bigl(\sum_{w \subseteq y} B_y \bigr).
因此
Cz=∑z⊆w(−1)\absw−z(∑w⊆xAx)(∑w⊆yBy). C_z = \sum_{z \subseteq w} (-1)^{\abs{w - z}}\bigl(\sum_{w \subseteq x} A_x\bigr) \bigl(\sum_{w \subseteq y} B_y \bigr).
如果读者熟悉数论函数的话, 可以把上面的式子和 h(d)=∑gcd⁡(i,j)=df(i)g(j)h(d) = \sum_{\gcd(i, j) = d} f(i) g(j) 时的下述等式做比较: ∑d∣kh(d)=(∑k∣if(i))(∑k∣jg(j)).\sum_{d \mid k} h(d) = \bigl(\sum_{k \mid i} f(i)\bigr)\bigl(\sum_{k \mid j} g(j)\bigr). h(d)=∑d∣kμ(k/d)(∑k∣if(i))(∑k∣jg(j)).h(d) = \sum_{d \mid k} \mu(k / d) \bigl(\sum_{k \mid i} f(i)\bigr)\bigl(\sum_{k \mid j} g(j)\bigr). 如果我们把 ∣\mid 和 ⊆\subseteq 对应起来, μ\mu 与 (−1)\abs⋅(-1)^{\abs{\cdot}} 对应起来, 那这两个式子基本上是一样的!
可能这就是 "Mobius 变换" 名称的由来? 笔者并未查询到这个名字的来源.
不过这种办法似乎不能用来给出 \xor\xor 卷积的一种形式.
那么, 剩下的问题就是已知 AxA_x 之后如何对每个 ww 算出 ∑w⊆xAx\sum_{w \subseteq x} A_x, 以及对每个 ww 都知道 ∑w⊆xCx\sum_{w \subseteq x} C_x 如何还原出 CxC_x. (⊇\supseteq 的情况是都对称的).
把 AA 想象成一个 nn 维, 每一维都是 0/10/1 的数组, 即 A0/1,0/1,…,0/1A_{0/1, 0/1, \dots, 0/1}. 那么不难看出, ∑w⊇xAx\sum_{w \supseteq x} A_x 就是 AA 的 (nn 维) 前缀和, 而 ∑w⊆xAx\sum_{w \subseteq x} A_x 就是 AA 的后缀和. 而反过来的过程, 就是 nn 维差分.
我们知道, 高维 (比如说二维) 前缀和可以每一维分别做. 比如说要求 Ai,jA_{i, j} 的前缀和, 就可以
CPP
for (int i = 0; i < n; ++i)
  for (int j = 1; j < m; ++j)
    A[i][j] += A[i][j - 1];

for (int j = 0; j < m; ++i)
  for (int i = 1; i < n; ++j)
    A[i][j] += A[i - 1][j];
同理, 差分的话, 把上面的枚举顺序反过来, 并把 += 改成 -= 即可.
因此我们要做高维前缀和, 也可以类似地每一维分别做:
CPP
// 求子集和
for (int i = 0; i < n; ++i) // 枚举 i 之后对第 i 维做前缀和
  for (int x = 0; x < N; ++x) // 需要枚举除了第 i 维之外的所有下标.
                              // 为了方便, 这里枚举所有 x, 并只考虑不包含 i 的那些 x
    if ((x >> i) & 1 == 0) A[x + (1 << i)] += A[x]; // 固定其他维之后对第 i 维做前缀和, 可以理解为 "A[x][1] += A[x][0]".
如果要逆运算 (知道子集和求原本的值), 就只需要把 += 改成 -= 即可.

"待定系数法"

我们发现, 如果卷积时的运算符是 ++ (准确地说, 是改成 z=(x+y) mod Nz = (x + y) \bmod N), 那么这个卷积就是我们通常意义的卷积, 也就是多项式乘法. 对于多项式乘法, 我们有 Fourier 变换 (如果读者不熟悉 Fourier 变换也没关系, 我们所需要的只有下面这两句话说明的思想).
形式地说, Fourier 变换就是对一个数组 AA, 算出了一个新的数组 \cF(A)\cF(A) (并且可以用逆变换从 \cF(A)\cF(A) 变回 AA), 使得若 AA 和 AA 卷积得到 CC, 那么 \cF(A)⋅\cF(B)=\cF(C)\cF(A) \cdot \cF(B) = \cF(C); 其中 ⋅\cdot 表示逐项相乘, 即 \cF(A)i×\cF(B)i=\cF(C)i\cF(A)_i \times \cF(B)_i = \cF(C)_i.
并且如果 AA 按下标分成两段, 那么我们可以在递归计算这两段的 Fourier 变换之后快速地合并成整个数组的 Fourier 变换;
此外, \cF\cF 还是线性的: 也就是说两个数组逐项相加之后 Fourier 变换和 Fourier 之后逐项相加是相等的, \cF(A+B)=\cF(A)+\cF(B)\cF(A + B) = \cF(A) + \cF(B).
在这个意义下, 我们的位运算卷积 (至少看上去) 有得天独厚的优势: 因为位运算是每一位互相独立的, 所以如果真的存在这样一个变换, 那么他很可能确实是可以分治递归之后合并的.
我们先假设这样的变换是存在的, 仍然记作 \cF\cF (上面我们用这个记号表示 Fourier 变换, 不过我们接下来再也用不到它了, 所以不用担心符号混淆), 然后想象这个 \cF\cF 在分治后怎么合并.
对于 0,…,N−10, \dots, N - 1 里的每一个数 x′x', 我们都可以把 x′x' 写成两部分: 二进制第 n−1n-1 位, 和剩下的位. 假设两部分分别是 aa 和 xx, 我们就简记 x′=a\ddxx' = a \dd x. 换句话说, a\ddxa \dd x 就表示 a×2n−1+xa \times 2^{n - 1} + x.
由于位运算是逐位做的, 我们有 (a\ddx′)⊕(b\ddy′)=(a⊕b)\dd(x′⊕y′)(a \dd x') \oplus (b \dd y') = (a \oplus b) \dd (x' \oplus y'). 因此,
Cc\ddz=∑a⊕b=c∑x⊕y=zAa\ddxBb\ddy.C_{c \dd z} = \sum_{a \oplus b = c} \sum_{x \oplus y = z} A_{a \dd x} B_{b \dd y}.
如果我们把 A,B,CA, B, C 分别按下标分成两段, 记作 A(0),A(1),B(0),B(1),C(0),C(1)A^{(0)}, A^{(1)}, B^{(0)}, B^{(1)}, C^{(0)}, C^{(1)}, 其中 Ax(0)=A0\ddx=Ax,Ax(1)=A1\ddx=A2n−1+xA^{(0)}_x = A_{0 \dd x} = A_x, A^{(1)}_x = A_{1 \dd x} = A_{2^{n - 1} + x}; B,CB, C 类似.
那么把上面的式子改写一下, 就有
&&C^{(c)}_z &= \sum_{a \oplus b = c} \sum_{x \oplus y = z} A^{(a)}_x B^{(b)}_y \\ &\implies& C^{(c)} &= \sum_{a \oplus b = c} A^{(a)} \ast B^{(b)}. \end{aligned}$$ 所以如果我们前面说的那个变换真的存在, 那我们应该有 $$\begin{aligned} && \cF(C^{(c)}) &= \sum_{a \oplus b = c} \cF(A^{(a)}) \cdot \cF(B^{(b)}) \\ &\implies& \cF(C^{(c)})_i &= \sum_{a \oplus b = c} \cF(A^{(a)})_i \cF(B^{(b)})_i. \end{aligned}$$ 最后这个式子中, 如果我们想象 $i$ 是一个常数, 而 $a, b, c$ 才是下标: 换句话说, 固定一个 $i$, 然后定义 $A'_0 = \cF(A^{(0)})_i, A'_1 = \cF(A^{1})_i$; $B, C$ 类似. 那么相当于说 $C' = A' \ast B'$, 其中 $\ast$ 是 $n = 1, N = 2$ 时候的卷积. 这样的话, 如果我们对每个 $i$, 再 "横着做一次 $n = 1$ 的变换" —— 也就是说, 按上面这个定义, 固定一个 $i$ 然后定义 $A'$ 之后, 我们就直接令 $\cF(A)_{0 \dd i} = \cF(A')_0, \cF_{1 \dd i} = \cF(A')_1$, 其中等号右侧的 $\cF$ 是对 $n = 1$ 做的变换. 那么这个 $\cF$ 就是我们想要的结果. **一句话总结**: 如果我们通过 $x = a \dd x'$ 的方法把数组看作二维数组: $A_x = A_{a, x'}$, 那么上面的结论告诉我们, 只要固定 $a$ 之后对 $x'$ 做一次变换, 再固定 $x'$ 对 $a$ 做一次变换, 就可以得到最终的变换. > 如果读者比较掌握 Fourier 变换, 那么可以有一个类比. > > 假设我们要对数组 $A_{i, j}$ 做 "二维 Fourier 变换", 那么我们可以固定 $i$, 对 $A_{i, \cdot}$ 做 Fourier 变换; > 然后再固定 $j$, 对 $A_{\cdot, i}$ 做 Fourier 变换. 那我们递归使用这个方法, 可以考虑把数组看作 $n$ 维数组 $A_x = A_{0/1, \dots, 0/1}$, 然后对每一维, 固定其余下标之后做一次 $N = 2$ 的变换. 这样就可以得到整个数组的变换. (伪)代码如下: ```cpp // 求子集和 for (int i = 0; i < n; ++i) // 枚举 i 之后对第 i 维做变换 for (int x = 0; x < N; ++x) // 需要枚举除了第 i 维之外的所有下标. // 为了方便, 这里枚举所有 x, 并只考虑不包含 i 的那些 x if ((x >> i) & 1 == 0) { (*) 对 A[x] 和 A[x + (1 << i)] 做 $N = 2$ 的变换. } ``` 那么我们只需要知道 $N = 2$ 的变换该怎么做. 我们可以用各种方法 (包括但不限于猜和问别人) 得知与, 或, 异或的 $N = 2$ 的变换分别是 $$\begin{aligned} &(\text{与}) & (a_0, a_1) &\to (a_0 + a_1, a_1) \\ &(\text{或}) & (a_0, a_1) &\to (a_0, a_0 + a_1) \\ &(\text{异或}) & (a_0, a_1) &\to (a_0 + a_1, a_0 - a_1) \\ \end{aligned}$$ 因此比方说, 把上述代码中 (\*) 部分改成 `int a0 = A[x], a1 = A[x + (1 << i)]; A[x] = a0 + a1; A[x + (1 << i)] = a0 - a1;` 就可以做异或的变换了. 而选手可能会发现, 通过这种办法得到的与和或的求法是和上一节通过容斥原理的求法是完全相同的! 我们将与/或的变换叫做 "快速 Mobius 变换" (未知出处); 而异或的变换叫做 "快速 Walsh 变换" (全名是 "Walsh–Hadamard 变换"). ## 结语 原本笔者还想写一些数学上的想法. 但是写了一千字之后感觉篇幅过长了, 不太符合本文"简介"的标题. 事实上如果对熟悉抽象代数的读者(我想应该基本上没有吧), 上面其实就是在说群直积的群代数同构于群代数的张量积. 不过即使是简要介绍张量积的概念, 也需要过多篇幅. 遂放弃. 希望读完本文的读者和选手们能有所启发! 此外, 位运算相关还有子集卷积, 乃至"集合幂级数"等内容, 感兴趣的读者可以自行了解! 感谢观看.

评论

共 42 条评论，欢迎与作者交流。

正在加载评论...