最通俗易懂的FWT变换讲解（快速沃尔什变换）

前言

其实我还没有退役，又回来更新算法笔记啦~~

我们知道 FFT 可以在 $O(n\log n)$ 的时间内求解 $C_k={\sum }_{i+j=k}{A}_i\times B_j$。

而类似的，我们将求和符号中的加号换成其余的位运算符号，即得到了 FWT 。

FWT 可以在 $O(n\log n)$ 的时间内求解：

• 按位或运算 $\cup$，$C_k={\sum }_{i\cup j=k}{A}_i\times B_j$
• 按位与运算 $\cap$，$C_k={\sum }_{i\cap j=k}{A}_i\times B_j$
• 按位异或运算 $\oplus$，$C_k={\sum }_{i\oplus j=k}{A}_i\times B_j$

且总体流程跟 FWT 一模一样：

• 先对 $A,B$ 做正变换，在 $O(n\log n)$ 的时间内求解得到 $A^{\prime },B^{\prime }$
• 随后在 $O(n)$ 的时间内求解 $C_i^{\prime }=A_i^{\prime }\times B_i^{\prime }$
  • 也就是两个数组按位置做点积；
  • 即两个数组中每个位置分别相乘。
• 最后对 $C^{\prime }$ 做逆变换，在 $O(n\log n)$ 的时间内求解得到 $C$。

总时间复杂度还是 $O(n\log n)$ 的。

于是我们需要关注的点就是正变换、逆变换怎么做。

• $A^{\prime }=FWT[A]$ 表示对 $A$ 做正变换。
• $UFWT[A^{\prime }]=A$ 表示对 $A^{\prime }$ 做逆变换。

而上面的第三步，其实就是：

• $UFWT[FWT[A]]=A$

现在我们需要解决的问题，就是：

• 如何构造 $FWT,UFWT$ 。
• 以及知道定义后，如何快速计算 $FWT,UFWT$。

需要注意一个点，位运算始终只对下标，元素的运算还是普通乘法。

符号约定

本文中，默认多项式长度 $n$ 均为 $2$ 的非负整数次幂。文中多项式、数组其实含义基本相同，都是对一个序列处理。

为方便陈述，约定如下符号：

多项式 $A$

• 多项式，长度为 $n$，次数为 $n-1$ 。

• 具体有多项式 $A(x)=a_0+a_{1}x+a_2{x}^2\cdots +a_{n-1}{x}^{n-1}$，

• 简写为 $A$，更像把其理解为一个数组 $A=[a_0,a_1,\cdots ,a_{n-1}]$

系数 $a_i$

• 为多项式 $A$ 中 $x^i$ 前的系数，即 $[x^i]A(x)$
• 本文中大写字母表示多项式，小写字母表示具体的元素/系数。

多项式 $A_0,A_1$

• 分别表示 $A$ 的前一半和后一半。
• $A_0=[a_0,a_1,\cdots ,a_{n/2-1}]$
• $A_1=[a_{n/2},a_{n/2+1},\cdots ,a_n]$

可以注意到，前两条都涉及到下标，为了含义不冲突。对于单字母的多项式，笔者采用大小写来区分。

多项式加法 $A+B$

• 表示对两个多项式做加法。
• $A+B=[a_0+b_0,a_1+b_1,\cdots ,a_{n-1}+b_{n-1}]$

按位卷积 $A\oplus B$

• 此处 $\oplus$ 指的是某种二元位运算，表面含义是异或，也可以换成按位与、或、同或。
• 有 $A\oplus B={\sum }_{k=0}^{n-1}{\sum }_{i\oplus j=k}{a}_i{b}_j$

FWT 变换 $FWT(A)$

• 对 $A$ 做 FWT 正变换，

FWT 逆变换 $UFWT(A)$

• 对 $A$ 做 FWT 逆变换

因为不存在歧义，所以我们直接用下标 $FWT(A{)}_i$ 表示其 $x^i$​ 前的系数

正变换后的数组 $A^{\prime }$

• 即 $A^{\prime }=FWT(A)$
• 二者等价

数组拼接 $[A,B]$

• 表示对两个数组进行拼接
• 即 $[A,B]=[a_0,a_1,\cdots ,a_{n-1},b_0,b_1,\cdots ,b_{m-1}]$​

点积 $FWT(A)\cdot FWT(B)$

• 运算符 $\cdot$​ 用在两个多项式之间时，表示点积，即相同位置乘积。

切片 $A[l:r]$

• 表示多项式 $A$ 下标在 $[l,r)$ 之间的系数。
• 即 $A[l:r]=[a_l,a_{l+1},\cdots ,a_{r-1}]$
• 类似 python 编程语言

或运算

即求解 $c_k={\sum }_{i\cup j=k}{a}_i{b}_j$ 。

FWT 变换

我们固定下标 $k=i\cup j$，去寻找满足条件的 $i,j$。

• 显然 $i$ 的二进制是 $1$ 的位，也必然是 $k$ 中二进制为 $1$ 的位。

• 考虑二进制状态压缩，即 $i$ 是 $k$ 的子集。

• 对于 $j$​ 同理。

以下构造 $FWT(A)$ 的关键规则：

• $k$ 的子集的子集，仍然是 $k$​​ 的子集。
• 子集的子集之间两两乘积最后统一求和，等价于子集的和和与另一个子集的总和的乘积，即满足交换律。

这样说可能有点抽象，我们继续往下看。

我们令 $A^{\prime }=FWT[A]$，

我们让 $a_i^{\prime }$ 等于所有是 $i$ 的子集的下标 $j$ 的系数和 ${\sum }_{j\subseteq i}{a}_j$，

即 $FWT(A{)}_i=a_i^{\prime }={\sum }_{j\subseteq i}{a}_j$

• 其中 $j\subseteq i$ 指的是，枚举 $i$ 的所有二进制子集 $j$。
• 在编程时，表示为 j or i == i

那么对于卷积 $FWT(A)\cdot FWT(B)$，考虑其下标 $k$ 有：

$$\begin{array}{rl}FWT(A{)}_k\times FWT(B{)}_k& =a_k^{\prime }\cdot b_k^{\prime }\\ & =(\sum _{i\subseteq k}{a}_i)(\sum _{j\subseteq k}{b}_j)\\ & =\sum _{i\subseteq k}\sum _{j\subseteq k}{a}_i{b}_j\\ & =\sum _{i\subseteq k}\sum _{j\subseteq k}{a}_i{b}_j\\ & =\sum _{(i\cup j)\subseteq k}{a}_i{b}_j\\ & =FWT(C{)}_k\end{array}$$
所以 $FWT(C)=FWT(A)\cdot FWT(B)$ 成立。

那么 $FWT(\cdot )={\sum }_{j\subseteq i}{a}_j$ 如何求解呢？

• 先考虑暴力，对于每个 $i$ 判断所有比它小的 $j$ 是否是它的子集。
• 复杂度为 $O(n^2)$。
• 需要继续优化。

对于多项式计算的一个常规套路，就是考虑分治。

• 当 $2$ 个区间长度为 $\frac{n}{2}$ 的都求解出来之后，
• 考虑合并，得到一个区间长度为 $n$ 的答案。

我们引入多项式 $A_0,A_1$，分别表示 $A$ 的前一半和后一半。

即下标最高位分别为 $0$ 和 $1$ 的。

• $A_0=[a_0,a_1,\cdots ,a_{n/2-1}]$
• $A_1=[a_{n/2},a_{n/2+1},\cdots ,a_n]$

假设已经求解出了 $FWT(A_0),FWT(A_1)$，那么如何合并求解 $FWT(A)$​ 呢？

因为当区间长度为 $n=1$ 时，$FWT(A)=A$，即 $a_0^{\prime }=a_0$，底层肯定成立。

（后面我们会知道，无论什么规则，$n=1$这条永远不变，因为 $FWT$ 规则有前一半后一半 $A_0,A_1$ 也就说明 $n$ 至少为 $2$）

类似归纳法，我们只需要向上推即可。

• 我们知道 $FWT(\cdot )$ 含义，是这个下标位置的子集之和。
• 求解 $FWT(A_1)$ 时，因为区间长度为 $\frac{n}{2}$。
  • 因为递归求解，其实我们此时眼中，其下标只是 $0\sim \frac{n}{2}-1$。
  • 而在求解 $FWT(A)$ 时，其下标就变成了 $\frac{n}{2}\sim (n-1)$。
• 而 $FWT(A_0)$ 对应的下标，则没有任何变化，所以变成 $FWT(A)$ ，这部分对应的值也没有任何变化。即 $FWT(A)[0:n/2]=FWT(A_0)$

而后半部分 $FWT(A)[n/2,n]$ 则还需要再考虑：

• 对于其中的一个下标 $k$，我们知道起二进制表示为 1xxxx，他已经统计了子集形式为 1???? 的子集 $j$。
• 现在只需要确保统计进其子集形式为 0???? 的子集 $j$ 即可，而其中 ???? 显然不能大于 xxxx。
• 而这个信息，在 $FWT(A_0{)}_{(k-n/2)}$ 中恰好就有，其下标表示恰好就为 0xxxx。
• 也就是说，在位置偏移 $n/2$ 的下标，恰好就有我们要的信息。
• 直接加这个数字即可。

于是做多项式加法即可，按下标相加，即 $O(n)$ 得到了 $FWT(A)[n/2,n]=FWT(A_0)+FWT(A_1)$。

最后拼接起来有 $FWT(A)=[FWT(A_0),FWT(A_0)+FWT(A_1)]$。

总的时间复杂度是 $T(n)=2T(n/2)+O(n)$，同归并排序为 $T(n)=O(n\log n)$。

当初看这块内容时，我突然想到树状数组，然后就瞬间明白了这一块内容，有异曲同工之妙。

逆变换

将得到的 $FWT(A),FWT(B)$ 做点积即得到了 $FWT(C)$。

现在考虑如何设计逆变换 $UFWT(\cdot )$ 变回 $C$。

依旧类似的，考虑分治。

对于 $A$，求解出了 $UFWT(A_0),UFWT(A_1)$，那么如何合并求解 $UFWT(A)$​ 呢？

容易知道 $n=1$，依旧有 $UFWT(A)=A$，即 $a_0=a_0^{\prime }$。

• 是一个反演的过程，
• $UFWT(\cdot )$ 的含义，是扣除掉所有这个下标位置的真子集。
• 不难知道， $UFWT(A_0)$ 对应的下标，已经扣除完了。
• 而对于 $UFWT(A_1)$，考虑下标 $k\in [n/2,n)$ ，设 k=1xxxx，其已经扣掉了 1???? 这样的真子集了。还差 0???? 这样的真子集没扣除掉。
• 此时这些信息在 $UFWT(A_0{)}_{(k-n/2)}$ 里面有。
• 只需要扣除这个数字即可。

故，最终有 $UFWT(A)=[UFWT(A_0),UFWT(A_1)-UFWT(A_0)]$。

时间复杂度依旧为 $O(n\log n)$。

代码实现

按照上面的思路，我们可以很轻松将递归代码模拟出来。

这个递归过程是**“尾递归”，学过编译原理的同学应该知道，是可以化简为迭代**的。

迭代实现，常数更小

我们发现，

• 只有在处理 $FWT(A_1)/UFWT(A_1)$ 信息的时候，需要加上/扣掉 $FWT(A_0)/UFWT(A_0)$ 对应位置的信息。
• 而处理 $FWT(A_0)/UFWT(A_0)$，则什么都不用做。
• 所以重点应该放在前者上面。

进一步，我们解读迭代代码，更深刻理解

个人认为迭代代码实现的思路，跟狄利克雷前后缀和非常相似

• 读者若更多的空闲时间，可以学习下狄利克雷前后缀和
• 可能会对这份代码有更深的理解。

其中的 type = 1 时表示正变换，type =-1 表示逆变换。

解释如下：

• 外层 $x$ 用于表示当前 $A$ 的长度，
• $k$ 表示区间长度的一半，也就是 $A_1$ 要找的信息，扣掉这个长度，就在 $A_0$ 中找到相应的
• $i$ 表示枚举每段长度为 $x$ 的区间的左端点， $j$ 遍历这个区间的每个位置。
• 内层 (a[i + j + k] += a[i + j] * type) %= P;
  • 拆分为 (a[i + j + k] += a[i + j] * type) 表示加上/减去对应位置的元素。
  • 因为赋值号 返回的是一个左值，所以直接外面套个括号对 a[i + j + k]取模。
  • 最终展示为 (a[i + j + k] += a[i + j] * type) %= P;。

void Or(LL a[], LL type)
{ // 迭代实现，常数更小
    for (int x = 2; x <= n; x <<= 1)
    {
        int k = x >> 1;
        for (int i = 0; i < n; i += x)
            for (int j = 0; j < k; j++)
                (a[i + j + k] += a[i + j] * type) %= P;
    }
}

与运算

即求解 $c_k={\sum }_{i\cap j=k}{a}_i{b}_j$ 。

FWT变换和逆变换

其实与运算，跟或运算完全类似，因为他们是对称的。

• 因为子集和超集，是相对应的。

• $i$ 是 $k$ 的子集，那么 $k$ 就是 $i$ 的超集。

于是可以通过以下性质来构造规则：

• 集合 $k$ 的超集的超集，还是 $k$​ 的超集。
• 超集的超集之间两两乘积最后统一求和，等价于超集的和和与另一个超集的总和的乘积，即满足交换律。

并且，我们在处理 $A_0$ 和 $A_1$ 的关系时，原本是对于 $A_1$ 要运算，现在变成对于 $A_0$ 要运算。

模仿上面的过程，很容易得到：

• $FWT(A)=[FWT(A_0)+FWT(A_1),FWT(A_1)]$。

• $UFWT(A)=[UFWT(A_0)-UFWT(A_1),UFWT(A_1)]$。

此过程留给读者自行做练习，不在赘述。

复杂度依旧为 $O(n\log n)$。

代码实现

只需要将或运算实现中，a[i + j + k] 与 a[i + j] 调换，即可得到与运算的代码

void And(LL a[], LL type)
{ // 迭代实现，常数更小
    for (int x = 2; x <= n; x <<= 1)
    {
        int k = x >> 1;
        for (int i = 0; i < n; i += x)
            for (int j = 0; j < k; j++)
                
                (a[i + j] += a[i + j + k] * type) %= P;
    }
}

异或运算

即求解 $c_k={\sum }_{i\oplus j=k}{a}_i{b}_j$​ 。

也是考得最多的。

FWT变换

约定：

• $popcnt(x)$
• 表示 $x$ 二进制中 $1$ 的数量。
• $x\circ y$
• 二者按位与后，剩余的 $1$ 的数量，的奇偶性。
• 即 $x\circ y=popcnt(x\cap y)\mathrm{mod}2$

可以发现 $\circ$ 对于 $\oplus$ 有分配律。

• 即 $(x\circ y)\oplus (x\circ z)=x\circ (y\oplus z)$​
• 汉字意思是：
  • $x$ 和 $y$ 与的奇偶性，跟 $x$ 和 $z$ 与的奇偶性，的异或
  • 结果等价于 $x$ 和 $y\oplus z$ 与的奇偶性。

这个分配律就是构造 $FWT$ 规则的关键。

考虑证明该分配律 $(x\circ y)\oplus (x\circ z)=x\circ (y\oplus z)$：

• 这个式子的值域是什么？
  • 没错只有 $0$ 和 $1$，可以理解为 奇偶性相同、不同 两种情况。
• 考虑将新增一位最高位，其实 $x$ 新增的一位是什么并不重要，对 $y,z$​ 分类：
  • 如果是 $0,0$ 的话，左边两个奇偶性结果都不变，右边 $y\oplus z$ 也不变
  • 如果是 $1,1$ 的话，左边两个奇偶性同时改变，右边 $y\oplus z$ 也不变，但是左边同时改变奇偶性之后，异或结果仍然跟先前相同（奇偶性只有相同、不同）
  • 如果是 $0,1$ 或 $1,0$ 的话，左边两个奇偶性不改变，左边奇偶性最后改变，同理右边 $y\oplus z$ 改变，右边奇偶性也改变。
• Q.E.D.

于是我们可以构造 $FWT(A)$：

• $FWT(A{)}_i=a_i^{\prime }={\sum }_{i\circ j=0}{a}_j-{\sum }_{i\circ j=1}{a}_j$。
• 即跟当前位的奇偶性为 $0$ 的减去奇偶性为 $1$ 的。

证明，对于下标 $k$ 有：
$$\begin{array}{rl}FWT(A{)}_i\times FWT(B{)}_i& =a_i^{\prime }\cdot b_i^{\prime }\\ & =(\sum _{i\circ j=0}{a}_j-\sum _{i\circ j=1}{a}_j)(\sum _{i\circ k=0}{b}_k-\sum _{i\circ k=1}{b}_k)\\ & =(\sum _{i\circ j=0}{a}_j\sum _{i\circ k=0}{b}_k+\sum _{i\circ j=1}{a}_j\sum _{i\circ k=1}{b}_k)-(\sum _{i\circ j=0}{a}_j\sum _{i\circ k=1}{b}_k+\sum _{i\circ j=1}{a}_j\sum _{i\circ k=0}{b}_k)\end{array}$$
考虑此处求和符号 $\frac{+\infty }{2}$ (bushi) 的下标，

• 因为 $\circ $ 运算的结果只有 $0,1$ ，所以我们对其变形，并不容易有信息的损失，
• 如 $i\circ (j\oplus k)=0$​ 就代表了第一个括号内的所有内容。

于是有：
$$\begin{array}{rl}FWT(A{)}_i\times FWT(B{)}_i& =\sum _{i\circ (j\oplus k)=0}{a}_j{b}_k-\sum _{i\circ (j\oplus k)=1}{a}_j{b}_k\\ & =FWT(C{)}_i\end{array}$$
验证了规则成立，接下来考虑如何求解，依旧采用分治：

假设已经求解出了 $FWT(A_0),FWT(A_1)$，考虑如何合并求解 $FWT(A{)}_i={\sum }_{i\circ j=0}{a}_j-{\sum }_{i\circ j=1}{a}_j$​

依旧是归纳法，

对于区间长度为 $n=1$ 时，我们令 $FWT(A)=A$，即 $a_0^{\prime }=a_0$。

• 可能有读者就疑惑了，为什么此处不代入 ${\sum }_{i\circ j=0}{a}_j-{\sum }_{i\circ j=1}{a}_j$ 呢？
• 其实做这些变换的时候， $n=1$ 层的规则我们永远不变，就是保留数组本身。
• 因为 FWT 是线性变换，对 $n=1$ 做这些变换没有任何意义，我们最终到 $n=1$ 想要拿回来原本的信息时，还需要再乘一个线性系数，纯属浪费时间。
• 另一方面，我们前文提到的对于 $FWT(A)$ 的求解方式，前提是默认其有 $A_0,A_1$，即 $n\ge 2$ 才有那一套规则，而非 $n=1$ 时。

继续归纳法，我们向上推。

• 我们知道 $FWT(\cdot )$ 含义，加上 $i\cup j$ 奇偶性为偶的，减去奇偶性为奇的。

• 而 $FWT(A_0)$ 对应的下标 $k$，

  • 对于原本的下标，最高位新增了 $0$，因为 $0\cap 0=0$ 所以并不会改变任何奇偶性，所以这部分保留。
  • 还得计算 $[n/2,n)$ 这些下标跟它的关系，因为 $0\cap 1=0$，所以也不会改变原先的奇偶性，只需要在累加上 $FWT(A_1{)}_{(k+n/2)}$ 即可。

• 于是有 $FWT(A)[0,n/2]=FWT(A_0)+FWT(A_1)$

• 同样的考虑 $FWT(A_1)$ 对应的下标 $k$：

  • 因为 $1\cap 0=0,1\cap 1=1$。
  • 所以不会改变前半段奇偶性，会改变后半段奇偶性。

• 于是有 $FWT(A)[0,n/2]=FWT(A_0)-FWT(A_1)$

综上有：

• $FWT(A)=[FWT(A_0)+FWT(A_1),FWT(A_0)-FWT(A_1)]$

逆变换

类似模仿在“或运算”那一节中的考虑方法，其实我们可以就可以推导出 $UFWT$ 是什么样子的。

但是，笔者觉得那样太无聊了。

我们再来一点更简单、更新的方式求解逆变换。

另一个角度的 FWT

前文我已说过 FWT是一个线性变换，既然是线性变换，也就存在线性变换矩阵，也就可以表示为矩阵乘法的形式。

我们重新看或运算的公式：

• $FWT(A)=[FWT(A_0),FWT(A_0)+FWT(A_1)]$
• $UFWT(A)=[UFWT(A_0),UFWT(A_1)-UFWT(A_0)]$
  设 $f_0=FWT(A_0),f_1=FWT(A_1)$，于是正变换可以表示成：
  $$\begin{array}{r}FWT(A)=\left[\begin{array}{cc}1& 0\\ 1& 1\end{array}\right]\left[\begin{array}{c}{f}_0\\ f_1\end{array}\right]=\left[\begin{array}{c}{f}_0\\ f0+f_1\end{array}\right]=\left[\begin{array}{c}{f}_0^{\prime }\\ f_1^{\prime }\end{array}\right]\end{array}$$
  那么，考虑逆变换，我们现在知道 $f_0^{\prime },f_1^{\prime }$，如何倒推回去？
  对的，乘一个逆矩阵，就是逆变换了
  $$\begin{array}{r}UFWT(A)=\left[\begin{array}{cc}1& 0\\ -1& 1\end{array}\right]\left[\begin{array}{c}{f}_0^{\prime }\\ f_1^{\prime }\end{array}\right]=\left[\begin{array}{c}{f}_0^{\prime }\\ f{1}^{\prime }-f_0\end{array}\right]=\left[\begin{array}{c}{f}_0\\ f_1\end{array}\right]\end{array}$$
  与我们最初推出来的 $UFWT(A)=[UFWT(A_0),UFWT(A_1)-UFWT(A_0)]$ 的形式一模一样！

于是我们可以类似求解 异或运算 的逆变换：

求解如下矩阵的逆矩阵即可：
$$\left[\begin{array}{cc}1& 1\\ 1& -1\end{array}\right]$$
解得为：
$$\left[\begin{array}{cc}0.5& 0.5\\ 0.5& -0.5\end{array}\right]$$
故对于异或运算， $UFWT(A)=[\frac{UFWT(A_0)+UFWT(A_1)}{2},\frac{UFWT(A_0)-UFWT(A_1)}{2}]$

异或运算代码实现

我们接着上文，实现异或代码

正变换时 type = 1，逆变换时 type = 1 / 2，注意此处是 $2$ 的逆元。

总体实现跟上面相同，最内层循环写成这样仅为了可以通过 type 切换，当然读者也可以写 if-else 来替换增加可读性。

void Xor(LL *a, LL type)
{
    for (LL x = 2; x <= n; x <<= 1)
    {
        LL k = x >> 1;
        for (LL i = 0; i < n; i += x)
        {
            for (LL j = 0; j < k; j++)
            {
                (a[i + j] += a[i + j + k]) %= P;
                (a[i + j + k] = a[i + j] - a[i + j + k] * 2) %= P;
                (a[i + j] *= type) %= P;
                (a[i + j + k] *= type) %= P;
            }
        }
    }
}

同或运算

即求解 $c_k={\sum }_{i\odot j=k}{a}_i{b}_j$​ 。

FWT变换和逆变换

与异或运算类似的，

约定如下符号

• $x\circ y$
  • 二者按位或后，剩余的 $1$ 的数量，的奇偶性。
  • 即 $x\circ y=popcnt(\overline{x\cup y})\mathrm{mod}2$
  • 先或后非，笔者认为oi-wiki里面写错了。

同样，可以发现 $\circ$ 对于 $\odot$ 有分配律。

类似的构造 $FWT(A{)}_i={\sum }_{i\circ j=0}{a}_j-{\sum }_{i\circ j=1}{a}_j$

• 因为 $\cup$ 和 $\cap$ 是对称的
• 其实 $\oplus$ 和 $\odot$ 也是对称的
• 只需要在求解 FWT 的结果，将 $A_0$ 和 $A_1$ 对调即可，当然拼接时位置也需要对调

类似上面的过程，最终解得
$$\begin{gathered} FWT(A)=[FWT(A_1)-FWT(A_0),FWT(A_1)+FWT(A_0)] \\ UFWT(A)=[\frac{UFWT(A_1)-UFWT(A_0)}{2},\frac{UFWT(A_1)+UFWT(A_0)}{2}] \end{gathered}$$

FWT 是线性变换

前文中已经提到了 $FWT$ 是线性变换，也强调了 $n=1$ 时的特殊性。

既然是线性变换，也就满足如下性质：

• $FWT(A+B)=FWT(A)+FWT(B)$

• $FWT(c\cdot A)=c\cdot FWT(A)$

另一方面，也可以理解为 $n=1$ 时我们就算设计 $FWT$ 规则时乘以一个倍数，最终变换回来也毫无影响

例题

例题1：洛谷P4717 【模板】快速莫比乌斯/沃尔什变换 (FMT/FWT)

题目描述

给定长度为 $2^n$ 两个序列 $A,B$，设

$$C_i=\sum _{j\oplus k=i}{A}_j\times B_k$$

分别当 $\oplus$ 是 or, and, xor 时求出 $C$。

输入格式

第一行，一个整数 $n$。
第二行，$2^n$ 个数 $A_0,A_1,\dots ,A_{2^n-1}$。
第三行，$2^n$ 个数 $B_0,B_1,\dots ,B_{2^n-1}$。

输出格式

三行，每行 $2^n$ 个数，分别代表 $\oplus$ 是 or, and, xor 时 $C_0,C_1,\dots ,C_{2^n-1}$ 的值 $\mathrm{mod}998244353$。

样例输入 #1

2
2 4 6 8
1 3 5 7

样例输出 #1

2 22 46 250
88 64 112 56
100 92 68 60

提示

$1\le n\le 17$。

是 FWT 的模板题，我们代入上面的推导即可。

因为代码中其实有类似的处理流程，此时我们可以用函数指针（或者auto）来让代码优雅一点

#include <bits/stdc++.h>
using namespace std;

typedef long long LL;
constexpr int N = 2e6 + 10;
constexpr int P = 998244353;

int n, m;
LL a[N], b[N], c[N];
LL sa[N], sb[N];

LL fpow(LL a, LL b, LL md) 
{
    LL res = 1;
    while (b) 
    {
        if (b & 1)
            res = res * a % md;
        a = a * a % md;
        b >>= 1;
    }
    return res;
}

void Or(LL a[], LL type) 
{
    for (int x = 2; x <= n; x <<= 1) 
    {
        int k = x >> 1;
        for (int i = 0; i < n; i += x)
            for (int j = 0; j < k; j++)
                (a[i + j + k] += a[i + j] * type) %= P;
    }
}

void And(LL a[], LL type) 
{
    for (int x = 2; x <= n; x <<= 1) 
    {
        int k = x >> 1;
        for (int i = 0; i < n; i += x)
            for (int j = 0; j < k; j++)
                (a[i + j] += a[i + j + k] * type) %= P;
    }
}

void Xor(LL *a, LL type) 
{
    for (LL x = 2; x <= n; x <<= 1) 
    {
        LL k = x >> 1;
        for (LL i = 0; i < n; i += x) 
        {
            for (LL j = 0; j < k; j++) 
            {
                (a[i + j] += a[i + j + k]) %= P;
                (a[i + j + k] = a[i + j] - a[i + j + k] * 2) %= P;
                (a[i + j] *= type) %= P;
                (a[i + j + k] *= type) %= P;
            }
        }
    }
}

void work(auto fwt, LL inv) 
{
    for (int i = 0; i < n; i++) 
    {
        a[i] = sa[i];
        b[i] = sb[i];
    }
    fwt(a, 1);
    fwt(b, 1);
    for (int i = 0; i < n; i++) 
    {
        c[i] = a[i] * b[i] % P;
    }
    fwt(c, inv); // 逆变换
    for (int i = 0; i < n; i++) 
    {
        if (c[i] < 0) // 可能有负数
            c[i] += P;
        cout << c[i] << " \n"[i == n - 1];
    }
}

void solve() 
{
    cin >> n;
    n = 1 << n;
    for (int i = 0; i < n; i++)
        cin >> sa[i];
    for (int i = 0; i < n; i++)
        cin >> sb[i];
    work(Or, -1); // 对应的fwt函数，以及逆变换系数
    work(And, -1);
    work(Xor, fpow(2, P - 2, P));
}

int main() 
{
    cin.tie(0)->sync_with_stdio(false);
    cout.tie(0);
    solve();
}

例题2：BZOJ4589 Hard Nim

有 $n$ 堆石子，每堆石子的数量均为不超过 $M$ 的质数，进行 Nim 游戏。

求解先手必败的开局局面数量，答案取模 $10^9+7$。

范围 $N\le 10^{18},M\le 5\times 10^4$，数据组数为不超过 $80$ 组。

构造多项式 $A$，若 $i\le M$ 且 $i$ 是质数，则 $a_i=1$，否则 $a_i=0$。

于是我们需要求解的，就是 $A$ 异或卷积 $N$ 次的多项式 $C$ 即：

• $C=\underset{N\text{times}}{\underbrace{A\oplus A\oplus \cdots \oplus A}}$

• 因为每个位置都是任取一个不超过 $M$ 的质数，乘法原理选出 $N$ 个即可

根据 SG 函数，我们最终需要选择异或和为 $0$ 的。

即 $c_0$​ 即是答案。

于是，我们可以：

• 先求得 $FWT(A)$
• 然后对于 $d_i’=(a_{i}')^N $，求得多项式 $D$
• 随后做逆变换 $C=UFWT(D)$

也就是全部将这些 $A$ 都做一遍 $FWT$，然后逐个位置做点积，对于单个位置也就是快速幂了。

总时间复杂度为 $O(M\log M\log N)$​​ 。

类似 FFT ，我们上面的 $FWT$ 的区间长度都是 $2$ 的整数次幂，这边的 $M$ 可能不是。

所以我们可以找到第一个大于 $M$ 的 $len=2^i$，然后以 $len$ 进行 $FWT$，并且令 $[M+1,len)$ 均为 $0$， 效果就是等价的。

因为 BZOJ 现在已经无法提交了，选手可以试一下如下大样例，过了应该就是对了

135186 13215
3156 11536
999999999 50000
129999999 50000
1000000000 47521
1000000000000000000 50000
985000000000000000 47785

结果依次为

270253099
208387442
957671404
882780244
81076371
567652446
904847008

代码实现：

#include <bits/stdc++.h>
using namespace std;

constexpr int N = 50000 + 10, P = 1e9 + 7;
typedef long long LL;
LL n;
int m;

LL fpow(LL a, LL b, LL md) 
{
    LL res = 1;
    while (b) 
    {
        if (b & 1)
            res = res * a % md;
        a = a * a % md;
        b >>= 1;
    }
    return res;
}

void Xor(LL *a, LL type, int n) 
{
    for (LL x = 2; x <= n; x <<= 1) 
    {
        LL k = x >> 1;
        for (LL i = 0; i < n; i += x) 
        {
            for (LL j = 0; j < k; j++) 
            {
                (a[i + j] += a[i + j + k]) %= P;
                (a[i + j + k] = a[i + j] - a[i + j + k] * 2) %= P;
                (a[i + j] *= type) %= P;
                (a[i + j + k] *= type) %= P;
            }
        }
    }
}

int p[N], tot;
bool isnp[N];

void init(int n) 
{
    isnp[1] = true;
    for (int i = 2; i <= n; i++) 
    {
        if (!isnp[i])
            p[++tot] = i;
        for (int j = 1; i * p[j] <= n; j++) 
        {
            int t = i * p[j];
            isnp[t] = true;
            if (i % p[j] == 0)
                break;
        }
    }
}

LL a[N * 2];

void solve() 
{
    memset(a, 0, sizeof a);
    int len = 1;
    while (len <= m)
        len <<= 1;
    for (int i = 1; i <= tot && p[i] <= m; i++)
        a[p[i]] = 1;
    Xor(a, 1, len);
    for (int i = 0; i < len; i++)
        a[i] = fpow(a[i], n, P);
    Xor(a, fpow(2, P - 2, P), len);
    if (a[0] < 0)
        a[0] += P;
    cout << a[0] << '\n';
}

int main() 
{
    cin.tie(0)->sync_with_stdio(false);
    cout.tie(0);
    init(50000);
    while (cin >> n >> m)
        solve();
}

[例题3：2020牛客多校2 E题](E-Exclusive OR_2020牛客暑期多校训练营（第二场） (nowcoder.com))

大致题意

给出 $1\le n\le 2\times 10^5$ 个数 $a_i(0\le a_i<2^{18})$，每个数字可以重复选，问选择出 $k\in [1,n]$ 个数字时，最大异或和是多少。

输出一行，包含 $k$ 个整数，表示答案。

聪明的读者可以发现，这题跟例题2非常相似。

我们同样构造多项式 $F(x)$，若有 $a_i$ 这个数，则 $f_i=1$， 否则 $f_i=0$。

也为了避免跟题意的符号冲突，这边换成 $F$

其实也就是初始时可以取到的那些异或和。

同样求解 $F$ 异或卷积 $k$ 次的多项式 $C$ 即：

• $C=F^k=\underset{k\text{times}}{\underbrace{F\oplus F\oplus \cdots \oplus F}}$

那么最大的 $i$ 满足 $c_i>0$，即是答案（因为数值表示方案数，而我们只需要知道是否存在即可）

容易知道 $C=F^k=F^{k-1}\oplus F$

也就是做 $n-1$ 次 FWT即可，设 $V=2^{18}$，那么每次的复杂度是 $O(V\log V)$ 的，总复杂度是 $O(nV\log V)$ 的。

显然不能接受，需要优化。

设 $g_i$ 表示取 $i$ 个数的最大异或和。

• 于是显然有 $g_i\ge g_{i-2}$
• 因为我们只需要再取两个相同的数，就至少会跟原来一样。

从线性基的角度，

• 不考虑取几个数，那么要取最大异或和，取出 $18$​ 个数即可到达最大异或和。

• （你可以想一想在线性基上跑最大异或和）

• 也就是说，再取多余的数，没有任何意义了。

那最多做 $18$ 次异或卷积，在这之后都让 $g_i=g_{i-2}$ ，这样是不是对的呢？

• 答案是：错的，需要做 $19$ 次。
• 可以参考一组 hack 只求到 $18$ 的测试数据 $1,2,4,\cdots ,2^{17},0,0,\cdots$
• 其中 $g_{17}=2^{18}-1-1,g_{19}=2^{18}-1$。

于是，当 $i\le 19$ 时我们做异或卷积，当 $i>19$ 时，我们直接让 $g_i=g_{i-2}$ 即可。

时间复杂度为 $O(V{\log }^{2}V)$。

#include <bits/stdc++.h>
using namespace std;

typedef long long LL;
constexpr int N = 1 << 18;
const int P = 1e9 + 7;
int n, m;

void Xor(LL *a, LL type, int n) 
{
    for (LL x = 2; x <= n; x <<= 1) 
    {
        LL k = x >> 1;
        for (LL i = 0; i < n; i += x) 
        {
            for (LL j = 0; j < k; j++) 
            {
                (a[i + j] += a[i + j + k]) %= P;
                (a[i + j + k] = a[i + j] - a[i + j + k] * 2) %= P;
                (a[i + j] *= type) %= P;
                (a[i + j + k] *= type) %= P;
            }
        }
    }
}

LL fn[N], f[N];
int g[N];

void solve() 
{
    cin >> n;
    int x;
    int mx = 1;
    for (int i = 0; i < n; i++) 
    {
        cin >> x;
        f[x] += 1;
        mx = max(x, mx);
    }

    int len = 1;
    while (len <= mx)
        len <<= 1;

    for (int i = 0; i <= mx; i++)
        fn[i] = f[i];

    Xor(f, 1, len); // 做一次正变换，省得以后每次做

    for (int i = len - 1; i >= 0; i--) 
    {
        if (fn[i] != 0) 
        {
            g[1] = i;
            break;
        }
    }

    LL inv2 = (P + 1) / 2;
    for (int k = 2; k <= n; k++) 
    {
        if (k > 19) 
        {
            g[k] = g[k - 2];
        } 
        else 
        {
            Xor(fn, 1, len);
            for (int i = 0; i < len; i++)
                fn[i] = fn[i] * f[i] % P;
            Xor(fn, inv2, len);

            for (int i = len - 1; i >= 0; i--) 
            {
                if (fn[i] != 0) 
                {
                    g[k] = i;
                    break;
                }
            }
        }
    }

    for (int k = 1; k <= n; k++)
        cout << g[k] << " \n"[k == n];
}

int main() 
{
    cin.tie(0)->sync_with_stdio(false);
    cout.tie(0);
    solve();
}

例题4：hdu5909 Tree Cutting

题目翻译

Byteasar有一棵包含 $n$ 个顶点的树 $T$，这些顶点用 $1$ 到 $n$ 的整数标记。树中的每个顶点都有一个整数值 $v_i$。非空树 $T$ 的权值等于 $v_1\oplus v_2\oplus \dots \oplus v_n$ ，其中 $\oplus$ 表示按位异或运算。现在对于从 $[0,m)$ 的每一个整数 $k$，请计算值等于 $k$ 的非空子树的数量。$T$ 的子树是 $T$ 的子图，且它也是一棵树。

输入格式

第一行包含一个整数$T(1\le T\le 10)$，表示测试用例的数量。

在每个测试用例中，

第一行包含两个整数 $n(n\le 1000)$ 和 $m(1\le m\le 2^{10})$ ，表示树的大小和 $v$ 的上界。

第二行包含 $n$ 个整数 $v_1,v_2,\dots ,v_n(0\le v_i<m)$，表示每个顶点的值。

• 接下来的 $n-1$ 行每行包含两个整数 $a_i,b_i$ ，表示顶点 $a_i$ 和 $b_i$ 之间有一条边 $(1\le a_i,b_i\le n)$。
• 保证 $m$ 可以表示为 $2^k$ ，其中 $k$ 是非负整数。

输出格式

对于每个测试用例，输出一行包含 $m$ 个整数，第 $i$ 个数表示值等于 $i$ 的非空子树的数量。答案很大，请对 $10^9+7$ 取模。

不难看出这是一道树形 DP 的题目。

设 $f[u][i]$ 表示在 $u$ 的子树中，异或和为 $i$ 且根为 $u$ 的连通子树的个数。

考虑节点 $u$，假设已经跟前 $i-1$ 个合并信息转移了状态，接下来考虑第 $i$ 个儿子 $v$。

那么有状态转移方程：

• $f[u][i\oplus j]=f[u][i\oplus j]+f[u][i]\times f[v][j]$​
• 转移需要枚举 $i,j$，复杂度为 $O(m^2)$。
• 总时间复杂度为 $O(n{m}^2)$ 显然不行。

我们发现，转移的瓶颈在于 $f[u][i]\times f[v][j]$，而他其实就是异或卷积的形式，对于所有 $i\oplus j$ 恰好为某个值的两两乘积。

于是我们用 FWT 优化这一步，即可将转移的复杂度从 $O(m^2)$ 优化到 $O(m\log m)$。

总时间复杂度为 $O(nm\log m)$。

#include <bits/stdc++.h>
using namespace std;

typedef long long LL;
typedef pair<int, int> PII;

constexpr int N = 1000 + 10, M = 1 << 10;
const int P = 1e9 + 7;
constexpr int INF = 0x3f3f3f3f;

int n, m;

void Xor(LL *a, LL type, int n) 
{
    for (LL x = 2; x <= n; x <<= 1) 
    {
        LL k = x >> 1;
        for (LL i = 0; i < n; i += x) 
        {
            for (LL j = 0; j < k; j++) 
            {
                (a[i + j] += a[i + j + k]) %= P;
                (a[i + j + k] = a[i + j] - a[i + j + k] * 2) %= P;
                (a[i + j] *= type) %= P;
                (a[i + j + k] *= type) %= P;
            }
        }
    }
}

vector<int> g[N];
int v[N];
LL f[N][M];
LL inv2 = (P + 1) / 2; // 2的逆元

void dfs(int u, int fa) 
{
    for (auto v : g[u])
        if (v != fa)
            dfs(v, u);
    f[u][v[u]] = 1; // 只有当前点构成的子树

    for (auto v : g[u]) 
    {
        if (v != fa) 
        {
            // 保留原本的
            vector<LL> tmp(m);
            for (int i = 0; i < m; i++)
                tmp[i] = f[u][i];

            // 计算卷积部分的
            Xor(f[u], 1, m);
            Xor(f[v], 1, m);
            for (int i = 0; i < m; i++)
                f[u][i] = f[u][i] * f[v][i] % P;
            Xor(f[u], inv2, m);
            Xor(f[v], inv2, m); // 需要变回去，因为统计答案的时候还需要

            for (int i = 0; i < m; i++)
                f[u][i] = (f[u][i] + tmp[i]) % P;
        }
    }
}

void solve() 
{
    cin >> n >> m;
    for (int i = 1; i <= n; i++) 
    {
        cin >> v[i];
        g[i].clear(); // 多测清空
    }
    for (int i = 1; i <= n; i++) 
        for (int j = 0; j < m; j++)
            f[i][j] = 0;

    for (int i = 1; i < n; i++) 
    {
        int u, v;
        cin >> u >> v;
        g[u].push_back(v);
        g[v].push_back(u);
    }
    dfs(1, 0);
    vector<LL> ans(m);
    // 统计以每个点位根的答案
    for (int i = 0; i < m; i++) 
    {
        for (int j = 1; j <= n; j++)
            ans[i] = (ans[i] + f[j][i]) % P;
        if (ans[i] < 0)
            ans[i] += P;
    }
    for (int i = 0; i < m; i++)
        cout << ans[i] << " \n"[i + 1 == m];
}

int main() 
{
    cin.tie(0)->sync_with_stdio(false);
    cout.tie(0);

    int T;
    cin >> T;
    while (T--)
        solve();
}

例题5：CodeForces662C Binary Table

题目翻译

给定一个 $n$ 行 $m$ 列的表格，每个元素都是 $0/1$ 。

每次操作可以选择一行或一列，把 $0/1$ 翻转，即把 $0$ 换为 $1$ ，把 $1$ 换为 $0$ 。

请问经过若干次操作后，表格中最少有多少个 $1$​ 。

输入格式

第一行，包含两个整数 $n$ 和 $m(1\le n\le 20,1\le m\le 10^5)$。

随后的 $n$ 行，每行一个长度为 $m$ 的 $01$ 字符串。

输出格式

仅一行一个整数，表示答案。

容易发现 $n\le 20$ 很小，而 $m\le 10^5$ 很大。

不难想到枚举 $n$ 的所有状态 $2^n$。

然后再去看每一列，

• 共有 $m$ 列，每列 $n$ 个 $0/1$​​ 其实可以压缩为一个数字 $v$。
• 假设我们当前枚举的状态是
• 然后根据每一列 $1$ 的数量的多少，单独决定是否需要翻转。

这样复杂度就是 $O(2^{n}m+nm)$。 $2^{20}$ 大概是 $10^6$，这样也就是 $10^{11}$ 显然是无法接受的。

注意区分行状态、列状态。

考虑优化：

• 类似桶排序的思路，约定 $c[x]$ 表列状态 $x$ 为列的数量
• $f(x)$ 表示列状态为 $x$ 时最小 $1$ 的数量（也就是不翻转、翻转，两个结果的 $popcnt$ 取个 $\min$）

若当前行翻转的状态为 $y$，那么当前的 $1$ 的数量就是：
$$\sum _{x=0}^{2^n-1}c[x]f(x\oplus y)$$
突然发现两个下标的异或，是个固定值 $y$。

于是可以变成异或卷积的形式：
$$\sum _{i\oplus j=y}c[i]f(j)$$
而上面这个式子，其实只是 FWT 结果里面的一项，我们最终需要求解的结果是：
$$\underset{y=0}{\overset{2^n-1}{\min }}\sum _{i\oplus j=y}c[i]f(j)$$
其实只需要做一次 FWT 就可以都求出来了，FWT的数组长度为 $2^n$，所以复杂度为 $O(2^n\log (2^n))=O(2^n\cdot n)$

故总的时间复杂度为 $O(2^n\cdot n+nm)$。

#include <bits/stdc++.h>
using namespace std;

typedef long long LL;
constexpr int B = 20, N = (1 << 20), M = 1e5 + 10;
const int P = 1e9 + 7;
constexpr int INF = 0x3f3f3f3f;

void Xor(LL *a, LL type, int n) 
{
    for (LL x = 2; x <= n; x <<= 1) 
    {
        LL k = x >> 1;
        for (LL i = 0; i < n; i += x) 
        {
            for (LL j = 0; j < k; j++) 
            {
                (a[i + j] += a[i + j + k]) %= P;
                (a[i + j + k] = a[i + j] - a[i + j + k] * 2) %= P;
                (a[i + j] *= type) %= P;
                (a[i + j + k] *= type) %= P;
            }
        }
    }
}

int n, m;
string s[B];
LL c[N];
LL f[N];

void solve() 
{
    cin >> n >> m;

    int U = (1 << n) - 1;
    for (int i = 0; i <= U; i++) //预处理f[x]
        f[i] = min(__builtin_popcount(i), __builtin_popcount(i ^ U));

    for (int i = 0; i < n; i++)
        cin >> s[i];

    for (int i = 0; i < m; i++) 
    {
        int x = 0;
        for (int j = 0; j < n; j++)
            x |= (int)(s[j][i] - '0') << j;
        c[x]++;
    }

    int len = 1 << n;
    Xor(c, 1, len);
    Xor(f, 1, len);
    for (int i = 0; i < len; i++)
        c[i] = c[i] * f[i] % P;

    LL inv2 = (P + 1) / 2;
    Xor(c, inv2, len);

    LL ans = 1e18;
    for (int i = 0; i < len; i++) 
    {
        if (c[i] < 0)
            c[i] += P;
        ans = min(ans, c[i]);
    }
    cout << ans << '\n';
}

int main() 
{
    cin.tie(0)->sync_with_stdio(false);
    cout.tie(0);
    solve();
}

后言

本文的例题，更偏向于一类 FWT 的基本应用，

或者说：

• 更像推导到某一步，然后发现可以用 FWT 去优化。
• 是一类套板子的题目

虽然只是 FWT 的一角，但也是 FWT 最容易遇到的情况和题目类型。

后续我可能会出一篇从矩阵角度的 FWT 的。

• 类似 $K$ 维 FWT
• 也更贴近 FWT 本质。

参考

1. oi-wiki：快速沃尔什变换
2. 知乎：牛客竞赛：算法笔记｜快速沃尔什变换 FWT