Skip to main content

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

· 19 min read
lailai
Student & Developer

原题链接

参考资料

题意简述

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

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

分别当 \oplusor,and,xor\operatorname{or},\operatorname{and},\operatorname{xor} 时求出 CC

解题思路

约定

  1. 为统一记号,设数据规模为:
m=2nm=2^n O(mlogm)=O(n2n)O(m\log m)=O(n2^n)
  1. 为避免歧义,位运算不再使用 ,&,\mid,\&,\oplus,统一记为:
  • 按位或:or\operatorname{or}(或 \cup
  • 按位与:and\operatorname{and}(或 \cap
  • 按位异或:xor\operatorname{xor}
  • 按位同或:xnor\operatorname{xnor}

引入

快速傅立叶变换(FFT,Fast Fourier Transform)可以 O(mlogm)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(mlogm)O(m\log m) 计算 或卷积与卷积

Ck=iorj=kAi×BjC_k=\sum_{i\operatorname{or} j=k}A_i\times B_j Ck=iandj=kAi×BjC_k=\sum_{i\operatorname{and} j=k}A_i\times B_j

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

Ck=ixorj=kAi×BjC_k=\sum_{i\operatorname{xor} j=k}A_i\times B_j Ck=ixnorj=kAi×BjC_k=\sum_{i\operatorname{xnor} j=k}A_i\times B_j

有些地方会将 FMT 和 FWT 统称为 FWT。

思想

FMT 和 FWT 的思想与 FFT 类似。

首先 O(mlogm)O(m\log m) 进行 正变换

AFMT/FWTAA\xrightarrow{\text{FMT/FWT}}A' BFMT/FWTBB\xrightarrow{\text{FMT/FWT}}B'

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

C=ABC'=A'\cdot B'

最后 O(mlogm)O(m\log m) 通过 逆变换 得到结果:

CIFMT/IFWTCC'\xrightarrow{\text{IFMT/IFWT}}C

这样整体时间复杂度为 O(mlogm)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}

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

TS    TST\supseteq S\iff \overline{T}\subseteq \overline{S}

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

变换 & 反演

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

fgf\to g

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

gfg\to f

在数学和算法中,我们会遇到许多变换,如傅里叶变换、数论变换、莫比乌斯变换和沃尔什变换等。

这里解释一下常见变换的命名方式:

  • 不带“快速”前缀的名称,指的是变换本身,并不涉及具体的实现方法。
  • 而带“快速”前缀的名称,表示使用一种高效实现该变换的算法。

Zeta 变换

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

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

g0(S)=TSf0(T)g_0(S)=\sum_{T\subseteq S}f_0(T) g1(S)=TSf1(T)g_1(S)=\sum_{T\supseteq S}f_1(T)

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

莫比乌斯反演

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

f0(S)=TS(1)STg0(T)f_0(S)=\sum_{T\subseteq S}(-1)^{|S|-|T|}g_0(T) f1(S)=TS(1)TSg1(T)f_1(S)=\sum_{T\supseteq S}(-1)^{|T|-|S|}g_1(T)

可以通过 容斥原理 证明:

TS(1)STg0(T)=TS(1)STVTf0(V)=VSf0(V)VTS(1)ST=VSf0(V)k=0SV(1)k(SVk)=VSf0(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)

原理

考虑如何构造变换规则。

iorj=ki\operatorname{or} j=k(并集为 kk),则 iijj 必为 kk 的子集:

iorj=kikjki\operatorname{or} j=k\Rightarrow i\subseteq k\land j\subseteq k

iandj=ki\operatorname{and} j=k(交集为 kk),则 iijj 必为 kk 的超集:

iandj=kikjki\operatorname{and} j=k\Rightarrow i\supseteq k\land j\supseteq k

我们可以构造:

FMTor[A]k=ikAi\operatorname{FMT_{or}}[A]_k=\sum_{i\subseteq k}A_{i} FMTand[A]k=ikAi\operatorname{FMT_{and}}[A]_k=\sum_{i\supseteq k}A_{i}

则有:

FMTor[A]kFMTor[B]k=(ikAi)×(jkBj)=ikjkAi×Bj=iorjkAi×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]kFMTand[B]k=(ikAi)×(jkBj)=ikjkAi×Bj=iandjkAi×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 最高位都为 00A1A_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(mlogm)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} 的递归版本:

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;
}
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} 的更高效(常数更小)的迭代版本:

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;
}
}
}
}

高维前缀和

前缀和

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

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

高维前缀和

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

扩展到 高维前缀和,就是各维度索引都 不超过 它的元素之和。

Si1,,in=j1i1,,jninAj1,,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)TATS_{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=j1i1,,jninAj1,,jn=j1i1jninAj1,,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(mlogm)O(n2^n)=O(m\log m)

逆变换

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

代码实现

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

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;
}
}
}
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=iorj=kiandj=0Ai×BjC_k=\sum_{\substack{i\operatorname{or} j=k\\i\operatorname{and} j=0}}A_i\times B_j

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

(AorB)k=iorj=kAi×Bj(A*_{\operatorname{or}}B)_k=\sum_{i\operatorname{or} j=k}A_i\times B_j

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

iandj=0    i+j=iorji\operatorname{and} j=0\iff |i|+|j|=|i\operatorname{or} j|

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

Fi,j={Ajj=i0ji,Gi,j={Bjj=i0jiF_{i,j}=\begin{cases} A_j & |j|=i \\ 0 & |j|\ne i \end{cases}, G_{i,j}=\begin{cases} B_j & |j|=i \\ 0 & |j|\ne i \end{cases}

则有:

Hi,j=k=0i(FkorGik)jH_{i,j}=\sum_{k=0}^i(F_{k}*_{\operatorname{or}}G_{i-k})_j

答案为:

Ci=Hi,iC_i=H_{|i|,i}

时间复杂度为 O(n22n)O(n^22^n)

快速沃尔什变换(FWT)

原理

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

若令 iji\circ j 表示 iji\cap j11 数量的奇偶性,即:

ij=popcount(ij)mod2i\circ j=\operatorname{popcount}(i\cap j)\bmod 2

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

则有:

(ik)xor(jk)=(ixorj)k(i\circ k)\operatorname{xor} (j\circ k)=(i\operatorname{xor} j)\circ k

我们可以构造:

FWTxor[A]k=ik=0Aiik=1Ai\operatorname{FWT_{xor}}[A]_k=\sum_{i\circ k=0}A_i-\sum_{i\circ k=1}A_i

则有:

FWTxor[A]kFWTxor[B]k=(ik=0Aiik=1Ai)×(jk=0Bjjk=1Bj)=(ik=0Aijk=0Bj+ik=1Aijk=1Bj)(ik=0Aijk=1Bj+ik=1Aijk=0Bj)=(ixorj)k=0AiBj(ixorj)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 最高位都为 00A1A_1 的最高位都为 11

当索引 ii 的最高位为 00 时,无论与最高位为 00 还是 11jj 进行 \circ 运算,因 00=0,01=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 时,与最高位为 00jj 运算结果为 10=01\cap 0=0(奇偶性不变),与最高位为 11jj 运算结果为 11=11\cap 1=1(奇偶性翻转),因此这一部分为两者之差:

FWTxor[A]i=FWTxor[A0]iFWTxor[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 的乘法逆元,即 mod+12\frac{mod+1}{2}

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

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);
for(int i=0;i<mid;i++)
{
ll x=a0[i],y=a1[i],t=type==1?1:mod+1>>1;
a0[i]=(x+y)*t%mod;
a1[i]=(x-y+mod)*t%mod;
}
}
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);
for(int i=0;i<mid;i++)
{
ll x=a0[i],y=a1[i],t=type==1?1:mod+1>>1;
a0[i]=(y-x+mod)*t%mod;
a1[i]=(y+x)*t%mod;
}
}

倍增迭代

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

void fwt_xor(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++)
{
ll x=a[i+j],y=a[i+j+k],t=type==1?1:mod+1>>1;
a[i+j]=(x+y)*t%mod;
a[i+j+k]=(x-y+mod)*t%mod;
}
}
}
}
void fwt_xnor(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++)
{
ll x=a[i+j],y=a[i+j+k],t=type==1?1:mod+1>>1;
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,生成 2241.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\%)N/A\text{N/A}
同或卷积(xnor\operatorname{xnor}937.16937.161164.39(+24.2%)1164.39(+24.2\%)N/A\text{N/A}

参考代码

#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)
{
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],t=type==1?1:mod+1>>1;
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];
template<typename F>
void calc(F 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;
}