Skip to main content
Original3.62K words18 min

题解:P4717 【模板】快速莫比乌斯 / 沃尔什变换 (FMT / FWT)

Summary

位运算卷积模板题:给两个长度为 2 的 n 次方的序列 A、B,分别在按位或、按位与、按位异或下求 C,其中 C_k 是所有 i 运算 j 等于 k 的 A_i 乘 B_j 之和。快速莫比乌斯变换处理或卷积、与卷积,本质是子集和与超集和,即 Zeta 变换与莫比乌斯反演;快速沃尔什变换处理异或、同或卷积。三者都是正变换后逐位相乘再逆变换,时间复杂度为 O(m log m)。文中给出位运算与集合的对应、变换与反演的原理推导,以及分治递归、倍增迭代、高维前缀和三种实现的代码与性能对比,并扩展到子集卷积。

参考资料​

题意简述​

给定两个长度为 2n2^n 的序列 A,BA,B,定义:

Ck=∑i⊕j=kAi×BjC_k=\sum_{i\oplus j=k}A_i\times B_j

分别在 ⊕\oplus 是 or⁡,and⁡,xor⁡\operatorname{or},\operatorname{and},\operatorname{xor} 时求序列 CC。

解题思路​

约定​

数据规模 记为:

m=2nm=2^n O(mlog⁡m)=O(n2n)O(m\log m)=O(n2^n)

位运算符号 记为:

  • 按位或:or⁡\operatorname{or} 或 ∪\cup;
  • 按位与:and⁡\operatorname{and} 或 ∩\cap;
  • 按位异或:xor⁡\operatorname{xor};
  • 按位同或:xnor⁡\operatorname{xnor}。

引入​

快速傅里叶变换(FFT,Fast Fourier Transform)可以 O(mlog⁡m)O(m\log m) 计算 加法卷积(多项式乘法):

Ck=∑i+j=kAi×BjC_k=\sum_{i+j=k}A_i\times B_j

若将 加法(++)替换为 位运算(or⁡,and⁡,xor⁡,xnor⁡\operatorname{or},\operatorname{and},\operatorname{xor},\operatorname{xnor}),得到 位运算卷积 该如何处理呢?

显然,直接枚举的时间复杂度为 O(m2)O(m^2),因此需要引入 FMT 和 FWT。

简介​

快速莫比乌斯变换(FMT,Fast Möbius Transform)可以 O(mlog⁡m)O(m\log m) 计算 或卷积 和 与卷积:

Ck=∑ior⁡j=kAi×BjC_k=\sum_{i\operatorname{or}j=k}A_i\times B_j Ck=∑iand⁡j=kAi×BjC_k=\sum_{i\operatorname{and}j=k}A_i\times B_j

快速沃尔什变换(FWT,Fast Walsh Transform)可以 O(mlog⁡m)O(m\log m) 计算 异或卷积 和 同或卷积:

Ck=∑ixor⁡j=kAi×BjC_k=\sum_{i\operatorname{xor}j=k}A_i\times B_j Ck=∑ixnor⁡j=kAi×BjC_k=\sum_{i\operatorname{xnor}j=k}A_i\times B_j

有些资料中会将 FMT 和 FWT 统称为 FWT。

思想​

FMT 和 FWT 的思想与 FFT 类似。

首先 O(mlog⁡m)O(m\log m) 进行 正变换:

A→FMT/FWTA′A\xrightarrow{\text{FMT/FWT}}A' B→FMT/FWTB′B\xrightarrow{\text{FMT/FWT}}B'

然后 O(m)O(m) 逐位相乘(点积):

C′=A′⋅B′C'=A'\cdot B'

最后 O(mlog⁡m)O(m\log m) 进行 逆变换 得到结果:

C′→IFMT/IFWTCC'\xrightarrow{\text{IFMT/IFWT}}C

这样总时间复杂度为 O(mlog⁡m)O(m\log m)。

基础知识​

位运算 & 集合​

每个 二进制数 可以视为一个 集合,二进制位为 11 的位置表示对应元素属于该集合。

(10010110)2  ⟺  { a1,a2,a4,a7 }(10010110)_2\iff\set{a_1,a_2,a_4,a_7}

两个数 按位或 对应集合的 并集,而 按位与 对应 交集。

(0101)2or⁡(1100)2=(1101)2  ⟺  { a0,a2 }∪{ a2,a3 }={ a0,a2,a3 }(0101)_2\operatorname{or}(1100)_2=(1101)_2\iff\set{a_0,a_2}\cup\set{a_2,a_3}=\set{a_0,a_2,a_3} (0101)2and⁡(1100)2=(0100)2  ⟺  { a0,a2 }∩{ a2,a3 }={ a2 }(0101)_2\operatorname{and}(1100)_2=(0100)_2\iff\set{a_0,a_2}\cap\set{a_2,a_3}=\set{a_2}

超集是补集的子集的补集,即超集和子集在补集意义上互为对偶。

T⊇S  ⟺  T‾⊆S‾T\supseteq S\iff\overline{T}\subseteq\overline{S}

简单来说,就是将 0/10/1 取反。因此子集和超集是 对称 的。

变换 & 反演​

变换(Transform)是指根据特定规则,将一个对象(例如数、向量、函数、图形等)映射为另一个对象的过程。

f→gf\to g

变换的逆过程称为 逆变换(Inverse Transform),亦称 反演(Inversion)。

g→fg\to f

变换广泛出现于数学和算法领域,例如傅里叶变换、数论变换、莫比乌斯变换和沃尔什变换。

常见变换的命名方式:

  • 不带「快速」前缀的名称,指变换本身,不涉及具体的实现算法。
  • 带有「快速」前缀的名称,表示一种高效实现该变换的算法。

Zeta 变换​

考虑一个集合函数 f:2V→Zf:2^V\to\mathbb{Z},其中 2V2^V 表示 VV 的幂集,即 VV 的所有子集构成的集合。

定义 子集和(Sum Over Subsets)和 超集和(Sum Over Supersets):

g0(S)=∑T⊆Sf0(T)g_0(S)=\sum_{T\subseteq S}f_0(T) g1(S)=∑T⊇Sf1(T)g_1(S)=\sum_{T\supseteq S}f_1(T)

计算 子集和 和 超集和 的过程(f→gf\to g)称为 Zeta 变换(Zeta Transform)。

莫比乌斯反演​

Zeta 的逆变换(g→fg\to f)称为 莫比乌斯反演(Möbius Inversion)。

f0(S)=∑T⊆S(−1)∣S∣−∣T∣g0(T)f_0(S)=\sum_{T\subseteq S}(-1)^{|S|-|T|}g_0(T) f1(S)=∑T⊇S(−1)∣T∣−∣S∣g1(T)f_1(S)=\sum_{T\supseteq S}(-1)^{|T|-|S|}g_1(T)

可以通过 容斥原理 证明:

∑T⊆S(−1)∣S∣−∣T∣g0(T)=∑T⊆S(−1)∣S∣−∣T∣∑V⊆Tf0(V)=∑V⊆Sf0(V)∑V⊆T⊆S(−1)∣S∣−∣T∣=∑V⊆Sf0(V)∑k=0∣S∣−∣V∣(−1)k(∣S∣−∣V∣k)=∑V⊆Sf0(V)[V=S]=f0(S)\begin{aligned} \sum_{T\subseteq S}(-1)^{|S|-|T|}g_0(T) & =\sum_{T\subseteq S}(-1)^{|S|-|T|}\sum_{V\subseteq T}f_0(V) \\ & =\sum_{V\subseteq S}f_0(V)\sum_{V\subseteq T\subseteq S}(-1)^{|S|-|T|} \\ & =\sum_{V\subseteq S}f_0(V)\sum_{k=0}^{|S|-|V|}(-1)^k\binom{|S|-|V|}{k} \\ & =\sum_{V\subseteq S}f_0(V)[V=S] \\ & =f_0(S) \end{aligned}

以上为子集的情形,超集的证明同理。

莫比乌斯反演也广泛出现于其他数学和算法领域,例如 数论中的莫比乌斯反演。

它们的共同点是在某种偏序上逆转累加关系以恢复原函数。

快速莫比乌斯变换(FMT)​

原理​

考虑如何构造变换规则。

若 ior⁡j=ki\operatorname{or}j=k(并集为 kk),则 ii 和 jj 必为 kk 的子集:

ior⁡j=k  ⟹  i⊆k∧j⊆ki\operatorname{or}j=k\implies i\subseteq k\land j\subseteq k

若 iand⁡j=ki\operatorname{and}j=k(交集为 kk),则 ii 和 jj 必为 kk 的超集:

iand⁡j=k  ⟹  i⊇k∧j⊇ki\operatorname{and}j=k\implies i\supseteq k\land j\supseteq k

我们可以构造:

FMTor⁡[A]k=∑i⊆kAi\operatorname{FMT_{or}}[A]_k=\sum_{i\subseteq k}A_i FMTand⁡[A]k=∑i⊇kAi\operatorname{FMT_{and}}[A]_k=\sum_{i\supseteq k}A_i

则有:

FMTor⁡[A]k⋅FMTor⁡[B]k=(∑i⊆kAi)×(∑j⊆kBj)=∑i⊆k∑j⊆kAi×Bj=∑ior⁡j⊆kAi×Bj=FMTor⁡[C]k\begin{aligned} \operatorname{FMT_{or}}[A]_k\cdot \operatorname{FMT_{or}}[B]_k & =\left(\sum_{i\subseteq k} A_i\right)\times\left(\sum_{j\subseteq k} B_j\right) \\ & =\sum_{i\subseteq k}\sum_{j\subseteq k}A_i\times B_j \\ & =\sum_{i\operatorname{or}j\subseteq k}A_i\times B_j \\ & =\operatorname{FMT_{or}}[C]_k \end{aligned} FMTand⁡[A]k⋅FMTand⁡[B]k=(∑i⊇kAi)×(∑j⊇kBj)=∑i⊇k∑j⊇kAi×Bj=∑iand⁡j⊇kAi×Bj=FMTand⁡[C]k\begin{aligned} \operatorname{FMT_{and}}[A]_k\cdot \operatorname{FMT_{and}}[B]_k & =\left(\sum_{i\supseteq k} A_i\right)\times\left(\sum_{j\supseteq k} B_j\right) \\ & =\sum_{i\supseteq k}\sum_{j\supseteq k}A_i\times B_j \\ & =\sum_{i\operatorname{and}j\supseteq k}A_i\times B_j \\ & =\operatorname{FMT_{and}}[C]_k \end{aligned}

以上两式分别对应 子集和 和 超集和,即 Zeta 变换 和 莫比乌斯反演。

分治递归​

正变换​

考虑如何计算 Zeta 变换。

显然,直接枚举的时间复杂度为 O(m2)O(m^2)。类似 FFT,可以考虑分治。

设序列 AA 的前一半为 A0A_0,后一半为 A1A_1。不难发现,A0A_0 最高位都为 00,A1A_1 的最高位都为 11。

那么,A0A_0 的子集就是它本身的子集,而 A1A_1 满足条件的子集是它本身和 A0A_0 的子集:

FMTor⁡[A]=merge⁡(FMTor⁡[A0],FMTor⁡[A1]+FMTor⁡[A0])\operatorname{FMT_{or}}[A]=\operatorname{merge}(\operatorname{FMT_{or}}[A_0],\operatorname{FMT_{or}}[A_1]+\operatorname{FMT_{or}}[A_0])

其中,merge⁡\operatorname{merge} 表示将两个序列首尾相连(类似字符串拼接),加法(++)表示逐项相加。

这样我们就能在 O(mlog⁡m)O(m\log m) 的时间复杂度得到 FMTor⁡[A]\operatorname{FMT_{or}}[A] 和 FMTand⁡[A]\operatorname{FMT_{and}}[A]。

逆变换​

既然 A0A_0 的子集是它本身,A1A_1 的子集是 FMT⁡[A1]+FMT⁡[A0]\operatorname{FMT}[A_1]+\operatorname{FMT}[A_0],即可得到反演的递推式:

IFMTor⁡[A]=merge⁡(IFMTor⁡[A0],IFMTor⁡[A1]−IFMTor⁡[A0])\operatorname{IFMT_{or}}[A]=\operatorname{merge}(\operatorname{IFMT_{or}}[A_0],\operatorname{IFMT_{or}}[A_1]-\operatorname{IFMT_{or}}[A_0])

按位与​

同理可得:

FMTand⁡[A]=merge⁡(FMTand⁡[A0]+FMTand⁡[A1],FMTand⁡[A1])\operatorname{FMT_{and}}[A]=\operatorname{merge}(\operatorname{FMT_{and}}[A_0]+\operatorname{FMT_{and}}[A_1],\operatorname{FMT_{and}}[A_1]) IFMTand⁡[A]=merge⁡(IFMTand⁡[A0]−IFMTand⁡[A1],IFMTand⁡[A1])\operatorname{IFMT_{and}}[A]=\operatorname{merge}(\operatorname{IFMT_{and}}[A_0]-\operatorname{IFMT_{and}}[A_1],\operatorname{IFMT_{and}}[A_1])

代码实现​

正变换和逆变换可以合并为一个函数:正变换取 type=1\text{type}=1,逆变换取 type=−1\text{type}=-1。

Tor⁡[A]=merge⁡(Tor⁡[A0],Tor⁡[A1]+Tor⁡[A0]×type)\operatorname{T_{or}}[A]=\operatorname{merge}(\operatorname{T_{or}}[A_0],\operatorname{T_{or}}[A_1]+\operatorname{T_{or}}[A_0]\times\text{type}) Tand⁡[A]=merge⁡(Tand⁡[A0]+Tand⁡[A1]×type,Tand⁡[A1])\operatorname{T_{and}}[A]=\operatorname{merge}(\operatorname{T_{and}}[A_0]+\operatorname{T_{and}}[A_1]\times\text{type},\operatorname{T_{and}}[A_1])

综上,可用分治写出 or⁡\operatorname{or} 和 and⁡\operatorname{and} 的递归版本:

191 Bcpp
void fmt_or(ll *a,int n,int type)
{
if(n==1)return;
int mid=n>>1;
ll *a0=a,*a1=a+mid;
fmt_or(a0,mid,type);
fmt_or(a1,mid,type);
for(int i=0;i<mid;i++)a1[i]=(a1[i]+a0[i]*type+mod)%mod;
}
194 Bcpp
void fmt_and(ll *a,int n,int type)
{
if(n==1)return;
int mid=n>>1;
ll *a0=a,*a1=a+mid;
fmt_and(a0,mid,type);
fmt_and(a1,mid,type);
for(int i=0;i<mid;i++)a0[i]=(a0[i]+a1[i]*type+mod)%mod;
}

倍增迭代​

类似 FFT,可用倍增写出 or⁡\operatorname{or} 和 and⁡\operatorname{and} 的更高效(常数更小)的迭代版本:

181 Bcpp
void fmt_or(ll *a,int n,int type)
{
for(int k=1;k<n;k<<=1)
{
for(int i=0;i<n;i+=k<<1)
{
for(int j=0;j<k;j++)
{
a[i+j+k]=(a[i+j+k]+a[i+j]*type+mod)%mod;
}
}
}
}
180 Bcpp
void fmt_and(ll *a,int n,int type)
{
for(int k=1;k<n;k<<=1)
{
for(int i=0;i<n;i+=k<<1)
{
for(int j=0;j<k;j++)
{
a[i+j]=(a[i+j]+a[i+j+k]*type+mod)%mod;
}
}
}
}

高维前缀和​

前缀和​

一维序列 AA 的 前缀和(Prefix Sum)序列 SS,表示索引 不超过 它的元素之和。

Si=∑j=1iAjS_i=\sum_{j=1}^i A_j

高维前缀和​

假设有 nn 个维度,每个维度的索引范围为 [0,p)[0,p)。

扩展到 高维前缀和(K-Dimensional Prefix Sum),就是各维度索引都 不超过 它的元素之和。

Si1,…,in=∑j1≤i1,…,jn≤inAj1,…,jnS_{i_1,\dots,i_n}=\sum_{j_1\le i_1,\dots,j_n\le i_n}A_{j_1,\dots,j_n}

计算高维前缀和的常用方法有两种:容斥原理 和 逐维前缀和。

容斥原理​

根据 容斥原理,可得:

Si1,…,in=∑T⊆{ 1,2,…,n }(−1)∣T∣ATS_{i_1,\dots,i_n}=\sum_{T\subseteq\set{1,2,\dots,n}}(-1)^{|T|}A_T

但该方法的时间复杂度为 O(2npn)O(2^np^n)。

逐维前缀和​

注意到,nn 维前缀和可视为 nn 次逐维求和:

Si1,…,in=∑j1≤i1,…,jn≤inAj1,…,jn=∑j1≤i1⋯∑jn≤inAj1,…,jnS_{i_1,\dots,i_n}=\sum_{j_1\le i_1,\dots,j_n\le i_n}A_{j_1,\dots,j_n}=\sum_{j_1\le i_1}\dots\sum_{j_n\le i_n} A_{j_1,\dots,j_n}

因此,每次只考虑一个维度,固定其余维度,求若干个一维前缀和。

对所有 nn 个维度分别求和后,即得到 nn 维前缀和。

该方法的时间复杂度为 O(npn)O(np^n)。

正变换​

当 p=2p=2,每维只有 0/10/1 两种索引时,高维前缀和 就是 子集和。

同理,可以通过 高维后缀和(K-Dimensional Suffix Sum)得到 超集和。

所以,我们也可以通过 高维前/后缀和 求出 Zeta 变换。

时间复杂度为 O(n2n)=O(mlog⁡m)O(n2^n)=O(m\log m)。

逆变换​

类似地,对所有 nn 个维度分别做差分后,可通过 高维差分 得到莫比乌斯反演。

代码实现​

我们可以写出 or⁡\operatorname{or} 和 and⁡\operatorname{and} 的高维前缀和版本:

141 Bcpp
void fmt_or(ll *a,int n,int type)
{
for(int i=1;i<n;i<<=1)
{
for(int j=0;j<n;j++)
{
if(i&j)a[j]=(a[j]+a[i^j]*type+mod)%mod;
}
}
}
145 Bcpp
void fmt_and(ll *a,int n,int type)
{
for(int i=1;i<n;i<<=1)
{
for(int j=0;j<n;j++)
{
if(!(i&j))a[j]=(a[j]+a[i^j]*type+mod)%mod;
}
}
}

子集卷积​

详见 洛谷 P6097 【模板】子集卷积。

FMT 还可以用于计算 子集卷积:

Ck=∑ior⁡j=kiand⁡j=0Ai×BjC_k=\sum_{\substack{i\operatorname{or}j=k\\i\operatorname{and}j=0}}A_i\times B_j

第一个限制 ior⁡j=ki\operatorname{or}j=k,可以通过 或卷积 轻松处理。

(A∗or⁡B)k=∑ior⁡j=kAi×Bj(A*_{\operatorname{or}}B)_k=\sum_{i\operatorname{or}j=k}A_i\times B_j

第二个限制 iand⁡j=0i\operatorname{and}j=0,即没有公共元素,等价于 ∣i∣+∣j∣=∣ior⁡j∣|i|+|j|=|i\operatorname{or}j|。

iand⁡j=0  ⟺  ∣i∣+∣j∣=∣ior⁡j∣i\operatorname{and}j=0\iff |i|+|j|=|i\operatorname{or}j|

我们可以多开一个维度,记录 集合大小。具体令:

Fi,j={Aj∣j∣=i0∣j∣≠iF_{i,j}= \begin{cases} A_j & |j|=i \\ 0 & |j|\ne i \end{cases} Gi,j={Bj∣j∣=i0∣j∣≠iG_{i,j}= \begin{cases} B_j & |j|=i \\ 0 & |j|\ne i \end{cases}

则有:

Hi,j=∑k=0i(Fk∗or⁡Gi−k)jH_{i,j}=\sum_{k=0}^i(F_k*_{\operatorname{or}}G_{i-k})_j

答案为:

Ci=H∣i∣,iC_i=H_{|i|,i}

时间复杂度为 O(n22n)=O(mlog⁡2m)O(n^22^n)=O(m\log^2 m)。

快速沃尔什变换(FWT)​

原理​

考虑构造异或运算的变换。

若令 i∘ji\circ j 表示 i∩ji\cap j 中 11 数量的奇偶性,即:

i∘j=popcount⁡(i∩j) mod 2i\circ j=\operatorname{popcount}(i\cap j)\bmod 2

其中 popcount⁡(x)\operatorname{popcount}(x) 表示 xx 二进制中 11 的个数。

则有:

(i∘k)xor⁡(j∘k)=(ixor⁡j)∘k(i\circ k)\operatorname{xor}(j\circ k)=(i\operatorname{xor}j)\circ k

我们可以构造:

FWTxor⁡[A]k=∑i∘k=0Ai−∑i∘k=1Ai\operatorname{FWT_{xor}}[A]_k=\sum_{i\circ k=0}A_i-\sum_{i\circ k=1}A_i

则有:

FWTxor⁡[A]k⋅FWTxor⁡[B]k=(∑i∘k=0Ai−∑i∘k=1Ai)×(∑j∘k=0Bj−∑j∘k=1Bj)=(∑i∘k=0Ai∑j∘k=0Bj+∑i∘k=1Ai∑j∘k=1Bj)−(∑i∘k=0Ai∑j∘k=1Bj+∑i∘k=1Ai∑j∘k=0Bj)=∑(ixor⁡j)∘k=0AiBj−∑(ixor⁡j)∘k=1AiBj=FWTxor⁡[C]k\begin{aligned} \operatorname{FWT_{xor}}[A]_k\cdot \operatorname{FWT_{xor}}[B]_k & =\left(\sum_{i\circ k=0}A_i-\sum_{i\circ k=1}A_i\right)\times\left(\sum_{j\circ k=0}B_j-\sum_{j\circ k=1}B_j\right) \\ & =\left(\sum_{i\circ k=0}A_i\sum_{j\circ k=0}B_j+\sum_{i\circ k=1}A_i\sum_{j\circ k=1}B_j\right) \\ & -\left(\sum_{i\circ k=0}A_i\sum_{j\circ k=1}B_j+\sum_{i\circ k=1}A_i\sum_{j\circ k=0}B_j\right) \\ & =\sum_{(i\operatorname{xor}j)\circ k=0}A_iB_j-\sum_{(i\operatorname{xor}j)\circ k=1}A_iB_j \\ & =\operatorname{FWT_{xor}}[C]_k \end{aligned}

而以上式子就是 沃尔什变换。

分治递归​

正变换​

我们依然考虑分治。

设序列 AA 的前一半为 A0A_0,后一半为 A1A_1。不难发现,A0A_0 最高位都为 00,A1A_1 的最高位都为 11。

当索引 ii 的最高位为 00 时,无论与最高位为 00 还是 11 的 jj 进行 ∘\circ 运算,由于 0∩0=0,0∩1=00\cap 0=0,0\cap 1=0,奇偶性不变,因此这一部分为两者之和:

FWTxor[A]i=FWTxor[A0]i+FWTxor[A1]i\mathrm{FWT_{xor}}[A]_i=\mathrm{FWT_{xor}}[A_0]_i+\mathrm{FWT_{xor}}[A_1]_i

当索引 ii 的最高位为 11 时,与最高位为 00 的 jj 运算结果为 1∩0=01\cap 0=0(奇偶性不变),与最高位为 11 的 jj 运算结果为 1∩1=11\cap 1=1(奇偶性翻转),因此这一部分为两者之差:

FWTxor[A]i=FWTxor[A0]i−FWTxor[A1]i\mathrm{FWT_{xor}}[A]_i=\mathrm{FWT_{xor}}[A_0]_i-\mathrm{FWT_{xor}}[A_1]_i

综上,可得:

FWTxor⁡[A]=merge⁡(FWTxor⁡[A0]+FWTxor⁡[A1],FWTxor⁡[A0]−FWTxor⁡[A1])\operatorname{FWT_{xor}}[A]=\operatorname{merge}(\operatorname{FWT_{xor}}[A_0]+\operatorname{FWT_{xor}}[A_1],\operatorname{FWT_{xor}}[A_0]-\operatorname{FWT_{xor}}[A_1])

逆变换​

逆变换易得:

IFWTxor⁡[A]=merge⁡(IFWTxor⁡[A0]+IFWTxor⁡[A1]2,IFWTxor⁡[A0]−IFWTxor⁡[A1]2)\operatorname{IFWT_{xor}}[A]=\operatorname{merge}\left(\frac{\operatorname{IFWT_{xor}}[A_0]+\operatorname{IFWT_{xor}}[A_1]}{2},\frac{\operatorname{IFWT_{xor}}[A_0]-\operatorname{IFWT_{xor}}[A_1]}{2}\right)

按位同或​

同理,将 ∘\circ 定义中的 ∩\cap 替换为 ∪\cup,即可得到同或:

FWTxnor⁡[A]=merge⁡(FWTxnor⁡[A1]−FWTxnor⁡[A0],FWTxnor⁡[A1]+FWTxnor⁡[A0])\operatorname{FWT_{xnor}}[A]=\operatorname{merge}(\operatorname{FWT_{xnor}}[A_1]-\operatorname{FWT_{xnor}}[A_0],\operatorname{FWT_{xnor}}[A_1]+\operatorname{FWT_{xnor}}[A_0]) IFWTxnor⁡[A]=merge⁡(IFWTxnor⁡[A1]−IFWTxnor⁡[A0]2,IFWTxnor⁡[A1]+IFWTxnor⁡[A0]2)\operatorname{IFWT_{xnor}}[A]=\operatorname{merge}\left(\frac{\operatorname{IFWT_{xnor}}[A_1]-\operatorname{IFWT_{xnor}}[A_0]}{2},\frac{\operatorname{IFWT_{xnor}}[A_1]+\operatorname{IFWT_{xnor}}[A_0]}{2}\right)

代码实现​

正变换和逆变换可以合并为一个函数:正变换取 type=1,t=1\text{type}=1,t=1,逆变换取 type=−1,t=12\text{type}=-1,t=\frac{1}{2}。

Txor⁡[A]=merge⁡((Txor⁡[A0]+Txor⁡[A1])×t,(Txor⁡[A0]−Txor⁡[A1])×t)\operatorname{T_{xor}}[A]=\operatorname{merge}((\operatorname{T_{xor}}[A_0]+\operatorname{T_{xor}}[A_1])\times t,(\operatorname{T_{xor}}[A_0]-\operatorname{T_{xor}}[A_1])\times t) Txnor⁡[A]=merge⁡((Txnor⁡[A1]−Txnor⁡[A0])×t,(Txnor⁡[A1]+Txnor⁡[A0])×t)\operatorname{T_{xnor}}[A]=\operatorname{merge}((\operatorname{T_{xnor}}[A_1]-\operatorname{T_{xnor}}[A_0])\times t,(\operatorname{T_{xnor}}[A_1]+\operatorname{T_{xnor}}[A_0])\times t)

注意:在模意义下,12\frac{1}{2} 表示 22 的乘法逆元,即 2−1=p+12(modp)2^{-1}=\frac{p+1}{2}\pmod p。

综上,可用分治写出 xor⁡\operatorname{xor} 和 xnor⁡\operatorname{xnor} 的递归版本:

261 Bcpp
void fwt_xor(ll *a,int n,int type)
{
if(n==1)return;
int mid=n>>1;
ll *a0=a,*a1=a+mid;
fwt_xor(a0,mid,type);
fwt_xor(a1,mid,type);
ll t=type==1?1:mod+1>>1;
for(int i=0;i<mid;i++)
{
ll x=a0[i],y=a1[i];
a0[i]=(x+y)*t%mod;
a1[i]=(x-y+mod)*t%mod;
}
}
264 Bcpp
void fwt_xnor(ll *a,int n,int type)
{
if(n==1)return;
int mid=n>>1;
ll *a0=a,*a1=a+mid;
fwt_xnor(a0,mid,type);
fwt_xnor(a1,mid,type);
ll t=type==1?1:mod+1>>1;
for(int i=0;i<mid;i++)
{
ll x=a0[i],y=a1[i];
a0[i]=(y-x+mod)*t%mod;
a1[i]=(y+x)*t%mod;
}
}

倍增迭代​

同理,可用倍增写出 xor⁡\operatorname{xor} 和 xnor⁡\operatorname{xnor} 更高效(常数更小)的迭代版本:

245 Bcpp
void fwt_xor(ll *a,int n,int type)
{
ll t=type==1?1:mod+1>>1;
for(int k=1;k<n;k<<=1)
{
for(int i=0;i<n;i+=k<<1)
{
for(int j=0;j<k;j++)
{
ll x=a[i+j],y=a[i+j+k];
a[i+j]=(x+y)*t%mod;
a[i+j+k]=(x-y+mod)*t%mod;
}
}
}
}
246 Bcpp
void fwt_xnor(ll *a,int n,int type)
{
ll t=type==1?1:mod+1>>1;
for(int k=1;k<n;k<<=1)
{
for(int i=0;i<n;i+=k<<1)
{
for(int j=0;j<k;j++)
{
ll x=a[i+j],y=a[i+j+k];
a[i+j]=(y-x+mod)*t%mod;
a[i+j+k]=(y+x)*t%mod;
}
}
}
}

性能分析​

  • 测试机器:CPU Apple M3 Max,RAM 36GB;
  • 测试环境:macOS Sequoia 15.6,clang 17.0.0;
  • 编译选项:-std=c++17 -O2;
  • 计时方法:调用 std::chrono::high_resolution_clock,单位为 ms\mathrm{ms};
  • 数据生成:调用 std::mt19937_64,生成 224≈1.68×1072^{24}\approx 1.68\times 10^7 个随机整数;
  • 测试方法:每个算法连续运行 100100 次,记录耗时并计算平均值。
类型/算法倍增迭代分治递归高维前缀和
或卷积(or⁡\operatorname{or})575.31575.31665.67(+15.7%)665.67(+15.7\%)966.18(+67.9%)966.18(+67.9\%)
与卷积(and⁡\operatorname{and})559.65559.65693.12(+23.8%)693.12(+23.8\%)961.38(+71.8%)961.38(+71.8\%)
异或卷积(xor⁡\operatorname{xor})939.02939.021167.31(+24.3%)1167.31(+24.3\%)—
同或卷积(xnor⁡\operatorname{xnor})937.16937.161164.39(+24.2%)1164.39(+24.2\%)—

参考代码​

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

using ll=long long;
const int mod=998244353;
const int N=20;
void fmt_or(ll *a,int n,int type)
{
for(int k=1;k<n;k<<=1)
{
for(int i=0;i<n;i+=k<<1)
{
for(int j=0;j<k;j++)
{
a[i+j+k]=(a[i+j+k]+a[i+j]*type+mod)%mod;
}
}
}
}
void fmt_and(ll *a,int n,int type)
{
for(int k=1;k<n;k<<=1)
{
for(int i=0;i<n;i+=k<<1)
{
for(int j=0;j<k;j++)
{
a[i+j]=(a[i+j]+a[i+j+k]*type+mod)%mod;
}
}
}
}
void fwt_xor(ll *a,int n,int type)
{
ll t=type==1?1:mod+1>>1;
for(int k=1;k<n;k<<=1)
{
for(int i=0;i<n;i+=k<<1)
{
for(int j=0;j<k;j++)
{
ll x=a[i+j],y=a[i+j+k];
a[i+j]=(x+y)*t%mod;
a[i+j+k]=(x-y+mod)*t%mod;
}
}
}
}
ll A[1<<N],B[1<<N],a[1<<N],b[1<<N];
void calc(auto f,int n)
{
for(int i=0;i<n;i++){a[i]=A[i];b[i]=B[i];}
f(a,n,1);
f(b,n,1);
for(int i=0;i<n;i++)a[i]=a[i]*b[i]%mod;
f(a,n,-1);
for(int i=0;i<n;i++)cout<<a[i]<<' ';
cout<<'\n';
}
int main()
{
ios::sync_with_stdio(false);
cin.tie(nullptr);
int n;
cin>>n;
n=1<<n;
for(int i=0;i<n;i++)cin>>A[i];
for(int i=0;i<n;i++)cin>>B[i];
calc(fmt_or,n);
calc(fmt_and,n);
calc(fwt_xor,n);
return 0;
}