快速沃尔什变换（FWT）学习笔记

Part 1 简介

FFT 本质上是处理加法卷积，即 $A_iB_j$ 贡献到 $C_{i+j}$ 。而 FWT 则是处理位运算卷积，即 $A_iB_j$ 贡献到 $C_{i\star j}$ ，其中 $\star$ 是某种位运算。

FWT 的思想与 FFT 相近，也是创造一种数列到数列的线性变换 $\operatorname{FWT}(A)$ 满足它与它的逆变换都可以快速计算且 $\operatorname{FWT}(A\star B)=\operatorname{FWT}(A)\cdot\operatorname{FWT}(B)$ ，其中 $\cdot$ 表示两个数列对应位置相乘。

有了这样的线性变换就可以先求出 $\operatorname{FWT}(A)$ 和 $\operatorname{FWT}(A)$ ，算出 $\operatorname{FWT}(A\star B)=\operatorname{FWT}(A)\cdot\operatorname{FWT}(B)$ ，最后求出 $C=\operatorname{IFWT}(A\star B)$ 。

本文中假定 $|A|=|B|=n$ 且 $n$ 是 $2$ 的整次幂。

Part 2 原理

先来解决 $\operatorname{FWT}(A)$ 如何求解，不妨设 $\operatorname{FWT}(A)_i=\sum\limits_{j=0}^{n-1}c_{i,j}A_j$ ，然后探讨 $c$ 需要满足的条件。那么有：

\begin{aligned} (\text{FWT}(A)\cdot \text{FWT}(B))_i&=\sum\limits_{j=0}^{n-1} c_{i,j}A_j\sum\limits_{k=0}^{n-1}c_{i,k}B_k\\ &=\sum\limits_{j=0}^{n-1}\sum\limits_{k=0}^{n-1} c_{i,j}c_{i,k}A_jB_k\\ \text{FWT}(A\star B)_i=\text{FWT}(C)&=\sum\limits_{j=0}^{n-1} c_{i,j}C_j\\ &=\sum\limits_{j=0}^{n-1}c_{i,j}\sum\limits_{1\le k,l\le n,k\star l=j} A_kB_l\\ &=\sum\limits_{j=0}^{n-1}\sum\limits_{k=0}^{n-1} c_{i,j\star k}A_jB_k \end{aligned}

所以有：

c_{i,j}c_{i,k}=c_{i,j\star k}

注意到由于 $\star$ 是位运算，所以不妨钦定 $c$ 也可以拆位处理。令 $(a)_2=\overline{a_0a_1a_2\dots}$ 即 $a$ 的二进制表示，不妨钦定 $c$ 满足 $c_{i,j}=c_{i_0,j_0}c_{i_1,j_1}c_{i_2,j_2}\dots$ ，这样做的好处是只要知道 $c_{0/1,0/1}$ 就可以求出 $c_{i,j}$ ，而且也有：

c_{i_l,j_l}c_{i_l,k_l}=c_{i_l,j_l\star k_l}\Leftrightarrow c_{i,j}c_{i,k}=c_{i,j\star k}

但是暴力求 $\text{FWT}(A)$ 是 $O(n^2)$ 的，考虑优化，令 $a'$ 为 $a$ 去掉二进制最高位的数，按位折半即有：

\begin{aligned} \text{FWT}(A)_i&=\sum\limits_{j=0}^{\frac{n}{2}-1} c_{i_0,j_0}c_{i',j'}A_{j}+\sum\limits_{j=\frac{n}{2}}^{n-1} c_{i_0,j_0}c_{i',j'}A_{j}\\ &=c_{i_0,0}\sum\limits_{j=0}^{\frac{n}{2}-1} c_{i',j'}A_{j}+c_{i_0,1}\sum\limits_{j=\frac{n}{2}}^{n-1} c_{i',j'}A_{j} \end{aligned}

那么考虑 $i_0$ 的取值，有：

\text{FWT}(A)_i=\begin{cases} c_{0,0}\sum\limits_{j=0}^{\frac{n}{2}-1} c_{i',j'}A_{j}+c_{0,1}\sum\limits_{j=\frac{n}{2}}^{n-1} c_{i',j'}A_{j}&0\le i\le \frac{n}{2}-1\\ c_{1,0}\sum\limits_{j=0}^{\frac{n}{2}-1} c_{i',j'}A_{j}+c_{1,1}\sum\limits_{j=\frac{n}{2}}^{n-1} c_{i',j'}A_{j}&\frac{n}{2}\le i\le n-1 \end{cases}

设 $A0$ 为 $A$ 中下标二进制最高位为 $0$ 的部分， $A1$ 为最高位为 $1$ 的部分，那么有：

\text{FWT}(A)_i=\begin{cases} c_{0,0}\text{FWT}(A0)_i+c_{0,1}\text{FWT}(A1)_i&0\le i\le \frac{n}{2}-1\\ c_{1,0}\text{FWT}(A0)_{i-\frac{n}{2}}+c_{1,1}\text{FWT}(A1)_{i-\frac{n}{2}}&\frac{n}{2}\le i\le n-1 \end{cases}

假设 $n=2^m$ ，则可以在 $O(m2^m)$ 即 $O(n\log n)$ 的时间复杂度内求解 $\text{FWT(A)}$ 。

对于 $\text{IFWT}(A)$ ，只需要构造出 $c_{0/1,0/1}$ 的逆即可。

Part 3 具体实现

根据 $c_{i,j}c_{i,k}=c_{i,j\star k}$ 和 $c_{i_0,j_0}c_{i_1,j_1}=c_{2i_0+i_1,2j_0+j_1}$ 构造 $c_{0/1,0/1}$ 即可，称其为位矩阵。

构造过程比较人类智慧，注意矩阵必须要有逆，即每一行和每一列都有至少一个位置不为 $0$ 且不能有两行或者两列完全一样，否则就会有维度被丢失（线性代数说法）。

由于不同的位运算的 $\text{FWT}$ 本质相同，只是 $c$ 不同，所以不妨设 $\text{FWT}\left(A,\begin{bmatrix}c_{0,0}&c_{0,1}\\c_{1,0}&c_{1,1}\end{bmatrix}\right)$ 为 $A$ 在对应的 $c$ 意义下的 $\text{FWT}$ 结果，那么有 $\text{FWT}\left(\text{FWT}\left(A,c\right),c^{-1}\right)=A$ 。

3.1 $\text{OR}$ 卷积

考虑构造满足 $c_{i,j}c_{i,k}=c_{i,j|k}$ 且存在逆的位矩阵。

$c_{0,0}c_{0,0}=c_{0,0|0}=c_{0,0}\Rightarrow c_{0,0}\in\{0,1\}$ 。

同理， $c_{0/1,0/1}\in\{0,1\}$ 。

由于 $c_{0,0}c_{0,1}=c_{0,1}$ ，所以 $c_{0,0}=1,c_{0,1}=0$ 或者 $c_{0,0}=1,c_{0,1}=1$ 。

同理， $c_{1,0}=1,c_{1,1}=0$ 或者 $c_{1,0}=1,c_{1,1}=1$ 。

那么位矩阵就有两种构造方式：

\begin{bmatrix} 1&0\\1&1 \end{bmatrix} \begin{bmatrix} 1&1\\1&0 \end{bmatrix}

Tips：

观察这个位矩阵：
$\begin{bmatrix} 1&0\\1&1 \end{bmatrix}$
注意到它满足 $c_{i,j}=[i\&j=j]$ ，也就是说这种情况下 $\text{FWT}(A,c)$ 实际上相当于子集求和。

这启发我们形如 $B_i=\sum\limits_{i\star j=i}A_j$ 和 $B_i=\sum\limits_{i\star j=j}A_j$ 这样的和式（ $\star$ 是某种位运算）也可以用 $\text{FWT}$ 来快速求。

由于第一个位矩阵满足 $c_{i,j}=[i\&j=j]$ ，所以下面采用第一个位矩阵，则设 $c^{-1}=\begin{bmatrix}x&y\\z&w\end{bmatrix}$ ，则有：

\begin{cases} x+0z=1\\ y+0w=0\\ x+z=0\\ y+w=1 \end{cases}

解得：

\begin{cases} x=1\\ y=0\\ z=-1\\ w=1 \end{cases}

所以 $c^{-1}=\begin{bmatrix}1&0\\-1&1\end{bmatrix}$ 。

3.2 $\text{AND}$ 卷积

$c_{0,0}c_{0,0}=c_{0,0\&0}=c_{0,0}\Rightarrow c_{0,0}\in\{0,1\}$ 。

同理， $c_{0/1,0/1}\in\{0,1\}$ 。

由于 $c_{0,0}c_{0,1}=c_{0,0}$ ，所以 $c_{0,0}=0,c_{0,1}=1$ 或 $c_{0,0}=1,c_{0,1}=1$ 。

同理， $c_{1,0}=0,c_{1,1}=1$ 或 $c_{1,0}=1,c_{1,1}=1$ 。

那么位矩阵就有两种构造方式：

\begin{bmatrix} 0&1\\1&1 \end{bmatrix} \begin{bmatrix} 1&1\\0&1 \end{bmatrix}

由于第一个位矩阵满足 $c_{i,j}=[i|j=2^k-1]$ ，所以采用第一个位矩阵，同理，待定系数法求逆得 $c^{-1}=\begin{bmatrix}-1&1\\1&0\end{bmatrix}$ 。

3.3 $\text{XOR}$ 卷积

由于对于任意的 $x,y$ ，均有 $c_{0,0}c_{x,y}=c_{x,y}$ ，所以 $c_{0,0}=1$ 。

根据 $c_{1,1}c_{1,1}=c_{1,0}$ 且矩阵不存在为 $0$ 的行，所以 $c_{1,0}$ 与 $c_{1,1}$ 均非 $0$ 。

根据 $c_{1,0}c_{1,0}=c_{1,0}$ 且 $c_{1,0}\not=0$ 可得 $c_{1,0}=1$ 。

根据， $c_{0,1}c_{0,1}=c_{1,0}$ ，可得 $c_{0,1}=-1$ 或 $c_{0,1}=1$ 。

同理， $c_{1,1}c_{1,1}=c_{1,0}$ ， $c_{1,1}=-1$ 或 $c_{1,1}=1$ 。

那么位矩阵就有两种构造方式：

\begin{bmatrix} 1&-1\\1&1 \end{bmatrix} \begin{bmatrix} 1&1\\1&-1 \end{bmatrix}

同样的，由于第二个位矩阵满足 $c_{i,j}=(-1)^{|i\&j|}$ （ $|a|$ 为 $a$ 二进制表示中 $1$ 的个数），所以采用第二个位矩阵，求逆得 $c^{-1}=\begin{bmatrix}\frac{1}{2}&\frac{1}{2}\\\frac{1}{2}&-\frac{1}{2}\end{bmatrix}$ 。

3.4 代码实现

直接套

\text{FWT}(A)_i=\begin{cases} c_{0,0}\text{FWT}(A0)_i+c_{0,1}\text{FWT}(A1)_i&0\le i\le \frac{n}{2}-1\\ c_{1,0}\text{FWT}(A0)_{i-\frac{n}{2}}+c_{1,1}\text{FWT}(A1)_{i-\frac{n}{2}}&\frac{n}{2}\le i\le n-1 \end{cases}

即可，P4717 【模板】快速莫比乌斯/沃尔什变换 (FMT/FWT) 代码如下：

#include <iostream>
#include <cstdio>

using namespace std;

const int p=998244353,inv2=499122177;

inline int getlen(int n)
{
	int res=1;
	while(res<n) res<<=1;
	return res;
}

const int ORC[2][2]={{1,0},{1,1}},IORC[2][2]={{1,0},{p-1,1}};
const int ANDC[2][2]={{0,1},{1,1}},IANDC[2][2]={{p-1,1},{1,0}};
const int XORC[2][2]={{1,1},{1,p-1}},IXORC[2][2]={{inv2,inv2},{inv2,p-inv2}};

inline void FWT(int n,int a[],const int c[2][2])
{
	for(int len=2;len<=n;len<<=1)
	{
		int mid=len>>1;
		for(int l=0;l<=n-len;l+=len)
		{
			for(int k=0;k<mid;k++)
			{
				int x=a[l+k],y=a[l+mid+k];
				a[l+k]=(1ll*c[0][0]*x%p+1ll*c[0][1]*y%p)%p;
				a[l+mid+k]=(1ll*c[1][0]*x%p+1ll*c[1][1]*y%p)%p;
			}
		}
	}
}

int n;
int a[1<<17],b[1<<17],c[1<<17];

int main()
{
	scanf("%d",&n);
	n=1<<n;
	for(int i=0;i<n;i++) scanf("%d",&a[i]);
	for(int i=0;i<n;i++) scanf("%d",&b[i]);
	FWT(n,a,ORC),FWT(n,b,ORC);
	for(int i=0;i<n;i++) c[i]=1ll*a[i]*b[i]%p;
	FWT(n,c,IORC);
	for(int i=0;i<n;i++) printf("%d ",c[i]);
	printf("\n");
	FWT(n,a,IORC),FWT(n,b,IORC);
	
	FWT(n,a,ANDC),FWT(n,b,ANDC);
	for(int i=0;i<n;i++) c[i]=1ll*a[i]*b[i]%p;
	FWT(n,c,IANDC);
	for(int i=0;i<n;i++) printf("%d ",c[i]);
	printf("\n");
	FWT(n,a,IANDC),FWT(n,b,IANDC);
	
	FWT(n,a,XORC),FWT(n,b,XORC);
	for(int i=0;i<n;i++) c[i]=1ll*a[i]*b[i]%p;
	FWT(n,c,IXORC);
	for(int i=0;i<n;i++) printf("%d ",c[i]);
	printf("\n");
	return 0;
}

Part 4 更多拓展

有些时候考的往往不是裸的 $\text{FWT}$ 。

下文中若 $A$ 与 $B$ 为数列， $\star$ 为某种位运算，那么 $A\star B$ 表示 $A$ 与 $B$ 在 $\star$ 运算下的卷积结果，即 $(A\star B)_i=\sum\limits_{j\star k=i}A_jB_k$ 。

$\text{FWT}$ 应用时往往要利用它是线性变换来优化，即 $\text{FWT}(A)+\text{FWT}(B)=\text{FWT(A+B)}$ 且 $\text{FWT}(aA)=a\text{FWT}(A)$ 。

若 $A$ 只有少数项非 $0$ 则可能有分类讨论优化时间复杂度的做法。

一些例题：

CF1119H Triple，题解

4.1 离线子集卷积

C_{k}=\sum\limits_{i|j=k,i\&j=0}A_iB_j=\sum\limits_{i\subseteq k}A_iB_{k-i}

发现 $i\&j=0$ 很烦，但是不难发现它等价于 $|i|+|j|=|k|$ （ $|a|$ 表示 $a$ 二进制表示中的 $1$ 的个数），所以可以令 $SA_{i,j}=[|j|=i]A_j,SB_{i,j}=[|j|=i]B_j$ ，那么有：

R_i=\sum\limits_{j=0}^iSA_j|SB_{i-j}

由于 $\text{FWT}$ 是线性变换，所以有：

R_i=\text{IFWT}\left(\sum\limits_{j=0}^i\text{FWT}(SA_j)\cdot \text{FWT}(SB_{i-j})\right)

答案即为 $R_{|i|,i}$ ，时间复杂度 $O(m^22^m)$ ，参考代码：（P6097 【模板】子集卷积）

int main()
{
	scanf("%d",&n);
	n=1<<n;
	for(int i=0;i<n;i++) scanf("%d",&a[i]);
	for(int i=0;i<n;i++) scanf("%d",&b[i]);
	for(int i=0;i<(1<<20);i++) for(int j=0;j<20;j++) popc[i]+=i>>j&1;
	for(int i=0;i<=20;i++) for(int j=0;j<n;j++) sa[i][j]=(popc[j]==i)*a[j],sb[i][j]=(popc[j]==i)*b[j];
	for(int i=0;i<=20;i++) FWT(n,sa[i],ORC),FWT(n,sb[i],ORC);
	for(int i=0;i<=20;i++) for(int j=0;j<=i;j++) for(int k=0;k<(1<<20);k++) r[i][k]=(r[i][k]+1ll*sa[j][k]*sb[i-j][k]%p)%p;
	for(int i=0;i<=20;i++) FWT(n,r[i],IORC);
	for(int i=0;i<n;i++) printf("%d ",r[popc[i]][i]);
	printf("\n");
	return 0;
}

4.2 半在线子集卷积

C_k=B_k\sum\limits_{i|j=k,i\&j=0,i\not=k}C_iA_j=B_k\sum\limits_{i\subset k}C_iA_{k-i}

和离线子集卷积类似，令 $SA_{i,j}=[|j|=i]A_j,SC_{i,j}=[|j|=i]C_j$ ，那么有：

\begin{aligned} SC_{i}&=B\cdot\sum\limits_{k=0}^{i-1}SC_k|SA_{i-k}\\ &=B\cdot\text{IFWT}\left(\sum\limits_{j=0}^i\text{FWT}(SC_j)\cdot \text{FWT}(SA_{i-j})\right) \end{aligned}

那么从小到大枚举 $i$ 计算即可。

4.3 每一位运算法则不同

给定一个长 $\log n$ 的字符串，字符集为 |&^，表示每一位要进行的位运算。

依旧是考虑：

\text{FWT}(A)_i=\begin{cases} c_{0,0}\text{FWT}(A0)_i+c_{0,1}\text{FWT}(A1)_i&0\le i\le \frac{n}{2}-1\\ c_{1,0}\text{FWT}(A0)_{i-\frac{n}{2}}+c_{1,1}\text{FWT}(A1)_{i-\frac{n}{2}}&\frac{n}{2}\le i\le n-1 \end{cases}

只不过此时 $c$ 取字符串第 $\log n$ 位对应运算的那个矩阵。

inline void FWT(int n,int a[],int w)
{
	for(int len=2,pos=1;len<=n;len<<=1,pos++)
	{
		int mid=len>>1;
		int c[2][2];
		if(w==1)
		{
			if(str[pos]=='|') memcpy(c,ORC,sizeof(ORC));
			if(str[pos]=='&') memcpy(c,ANDC,sizeof(ANDC));
			if(str[pos]=='^') memcpy(c,XORC,sizeof(XORC));
		}
		else
		{
			if(str[pos]=='|') memcpy(c,IORC,sizeof(IORC));
			if(str[pos]=='&') memcpy(c,IANDC,sizeof(IANDC));
			if(str[pos]=='^') memcpy(c,IXORC,sizeof(IXORC));
		}
		for(int l=0;l<=n-len;l+=len)
		{
			for(int k=0;k<mid;k++)
			{
				int x=a[l+k],y=a[l+mid+k];
				a[l+k]=(1ll*c[0][0]*x%p+1ll*c[0][1]*y%p)%p;
				a[l+mid+k]=(1ll*c[1][0]*x%p+1ll*c[1][1]*y%p)%p;
			}
		}
	}
}

快速沃尔什变换（FWT）学习笔记

Part 1 简介

Part 2 原理

Part 3 具体实现

3.1 OR\text{OR}OR 卷积

3.2 AND\text{AND}AND 卷积

3.3 XOR\text{XOR}XOR 卷积

3.4 代码实现

Part 4 更多拓展

4.1 离线子集卷积

4.2 半在线子集卷积

4.3 每一位运算法则不同

感谢您的支持，我会继续努力的!

3.1 $\text{OR}$ 卷积

3.2 $\text{AND}$ 卷积

3.3 $\text{XOR}$ 卷积