给定长度为 2n 的两个序列 A,B,设:
Ck=i⊕j=k∑Ai×Bj
分别当 ⊕ 是 or,and,xor 时求出 C。
- 为统一记号,设数据规模为:
m=2n
O(mlogm)=O(n2n)
- 为避免歧义,位运算不再使用 ∣,&,⊕,统一记为:
- 按位或:or(或 ∪)
- 按位与:and(或 ∩)
- 按位异或:xor
- 按位同或:xnor
快速傅立叶变换(FFT,Fast Fourier Transform)可以 O(mlogm) 计算 加法卷积(多项式乘法):
Ck=i+j=k∑Ai×Bj
若将求和符号中的 加法(+)替换为 位运算(or,and,xor,xnor),变成 位运算卷积 该如何处理呢?
直接枚举的时间复杂度为 O(m2),效率较低,因此需要使用 FMT 和 FWT。
快速莫比乌斯变换(FMT,Fast Möbius Transform)可以 O(mlogm) 计算 或卷积 和 与卷积:
Ck=iorj=k∑Ai×Bj
Ck=iandj=k∑Ai×Bj
快速沃尔什变换(FWT,Fast Walsh Transform)可以 O(mlogm) 计算 异或卷积 和 同或卷积:
Ck=ixorj=k∑Ai×Bj
Ck=ixnorj=k∑Ai×Bj
有些地方会将 FMT 和 FWT 统称为 FWT。
FMT 和 FWT 的思想与 FFT 类似。
首先 O(mlogm) 进行 正变换:
AFMT/FWTA′
BFMT/FWTB′
然后 O(m) 逐位相乘(点积):
C′=A′⋅B′
最后 O(mlogm) 通过 逆变换 得到结果:
C′IFMT/IFWTC
这样整体时间复杂度为 O(mlogm)。
每个二进制数可以视为一个 集合,二进制位为 1 的位置表示对应元素属于该集合。
(10010110)2⟺{a1,a2,a4,a7}
两个数 按位或 对应集合的 并集,而 按位与 对应 交集。
(0101)2or(1100)2=(1101)2⟺{a0,a2}∪{a2,a3}={a0,a2,a3}
(0101)2and(1100)2=(0100)2⟺{a0,a2}∩{a2,a3}={a2}
超集是补集的子集的补集,超集与子集在补集意义上互为对偶。
T⊇S⟺T⊆S
简单来说,就是将 0/1 取反,因此子集和超集是 对称 的。
变换(Transform)是指根据特定规则,将一个对象(如数、向量、函数、图形等)映射为另一个对象的过程。
f→g
变换的逆过程称为 逆变换,亦称 反演(Inversion)。
g→f
在数学和算法中,我们会遇到许多变换,如傅里叶变换、数论变换、莫比乌斯变换和沃尔什变换等。
这里解释一下常见变换的命名方式:
- 不带“快速”前缀的名称,指的是变换本身,并不涉及具体的实现方法。
- 而带“快速”前缀的名称,表示使用一种高效实现该变换的算法。
考虑一个集合函数 f:2V→Z,其中 2V 表示 V 的幂集,即 V 的所有子集构成的集合。
定义 子集和(Sum Over Subsets)和 超集和(Sum Over Supersets):
g0(S)=T⊆S∑f0(T)
g1(S)=T⊇S∑f1(T)
计算 子集和 和 超集和 的过程(f→g)称为 Zeta 变换(Zeta Transform)。
Zeta 的逆变换(g→f)称为 莫比乌斯反演(Möbius Inversion)。
f0(S)=T⊆S∑(−1)∣S∣−∣T∣g0(T)
f1(S)=T⊇S∑(−1)∣T∣−∣S∣g1(T)
可以通过 容斥原理 证明:
T⊆S∑(−1)∣S∣−∣T∣g0(T)=T⊆S∑(−1)∣S∣−∣T∣V⊆T∑f0(V)=V⊆S∑f0(V)V⊆T⊆S∑(−1)∣S∣−∣T∣=V⊆S∑f0(V)k=0∑∣S∣−∣V∣(−1)k(k∣S∣−∣V∣)=V⊆S∑f0(V)[V=S]=f0(S)
以上为子集的情形,超集的证明同理。
莫比乌斯反演也广泛出现于其他数学和算法领域,如 数论中的莫比乌斯反演。
它们的共同点是在某种偏序上逆转累加关系以恢复原函数。
考虑如何构造变换规则。
若 iorj=k(并集为 k),则 i 和 j 必为 k 的子集:
iorj=k⇒i⊆k∧j⊆k
若 iandj=k(交集为 k),则 i 和 j 必为 k 的超集:
iandj=k⇒i⊇k∧j⊇k
我们可以构造:
FMTor[A]k=i⊆k∑Ai
FMTand[A]k=i⊇k∑Ai
则有:
FMTor[A]k⋅FMTor[B]k=(i⊆k∑Ai)×j⊆k∑Bj=i⊆k∑j⊆k∑Ai×Bj=iorj⊆k∑Ai×Bj=FMTor[C]k
FMTand[A]k⋅FMTand[B]k=(i⊇k∑Ai)×j⊇k∑Bj=i⊇k∑j⊇k∑Ai×Bj=iandj⊇k∑Ai×Bj=FMTand[C]k
以上两式分别对应 子集和 和 超集和,即 Zeta 变换 和 莫比乌斯反演。
考虑如何计算 Zeta 变换。
直接枚举的时间复杂度为 O(m2),效率较低。借鉴 FFT,我们可以考虑分治。
设序列 A 的前一半为 A0,后一半为 A1。不难发现,A0 最高位都为 0,A1 的最高位都为 1。
那么,A0 的子集就是它本身的子集,而 A1 满足条件的子集是它本身和 A0 的子集:
FMTor[A]=merge(FMTor[A0],FMTor[A1]+FMTor[A0])
其中,merge 表示将两个序列首尾相连(类似字符串拼接),加法(+)表示逐项相加。
这样我们就能在 O(mlogm) 的时间复杂度得到 FMTor[A] 和 FMTand[A]。
既然 A0 的子集是它本身,A1 的子集是 FMT[A1]+FMT[A0],即可得到反演的递推式:
IFMTor[A]=merge(IFMTor[A0],IFMTor[A1]−IFMTor[A0])
同理可得:
FMTand[A]=merge(FMTand[A0]+FMTand[A1],FMTand[A1])
IFMTand[A]=merge(IFMTand[A0]−IFMTand[A1],IFMTand[A1])
正变换和逆变换可以合并为一个函数:正变换取 type=1,逆变换取 type=−1。
Tor[A]=merge(Tor[A0],Tor[A1]+Tor[A0]×type)
Tand[A]=merge(Tand[A0]+Tand[A1]×type,Tand[A1])
综上,可用分治写出 or 和 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 和 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;
}
}
}
}
一维序列 A 的 前缀和(Prefix Sum)序列 S,表示索引 不超过 它的元素之和。
Si=j=1∑iAj
假设有 n 个维度,每个维度的索引范围为 [0,p)。
扩展到 高维前缀和,就是各维度索引都 不超过 它的元素之和。
Si1,…,in=j1≤i1,…,jn≤in∑Aj1,…,jn
计算高维前缀和的常用方法有两种:容斥原理 与 逐维前缀和。
根据 容斥原理,可得:
Si1,…,in=T⊆{1,2,…,n}∑(−1)∣T∣AT
但该方法的时间复杂度为 O(2npn),效率较低。
注意到,n 维前缀和可视为 n 次逐维求和:
Si1,…,in=j1≤i1,…,jn≤in∑Aj1,…,jn=j1≤i1∑⋯jn≤in∑Aj1,…,jn
因此,每次只考虑一个维度,固定其余维度,求若干个一维前缀和。
对所有 n 个维度分别求和后,即得到 n 维前缀和。
该方法的时间复杂度为 O(npn)。
当 p=2,每维只有 0/1 两种索引时,高维前缀和 就是 子集和。
同理,也可以通过 高维后缀和(K-Dimensional Suffix Sum)得到 超集和。
所以,我们可以通过 高维前缀和 求出 Zeta 变换。
时间复杂度为 O(n2n)=O(mlogm)。
类似地,对所有 n 个维度分别做差分后,可通过 高维差分 得到莫比乌斯反演。
我们可以写出 or 和 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=0∑Ai×Bj
第一个限制 iorj=k,可以通过 或卷积 轻松处理。
(A∗orB)k=iorj=k∑Ai×Bj
第二个限制 iandj=0,即没有公共元素,等价于 ∣i∣+∣j∣=∣iorj∣。
iandj=0⟺∣i∣+∣j∣=∣iorj∣
我们可以多开一个维度,记录 集合大小。具体令:
Fi,j={Aj0∣j∣=i∣j∣=i,Gi,j={Bj0∣j∣=i∣j∣=i
则有:
Hi,j=k=0∑i(Fk∗orGi−k)j
答案为:
Ci=H∣i∣,i
时间复杂度为 O(n22n)。
考虑构造异或运算的变换。
若令 i∘j 表示 i∩j 中 1 数量的奇偶性,即:
i∘j=popcount(i∩j)mod2
其中 popcount(x) 表示 x 二进制中 1 的个数。
则有:
(i∘k)xor(j∘k)=(ixorj)∘k
我们可以构造:
FWTxor[A]k=i∘k=0∑Ai−i∘k=1∑Ai
则有:
FWTxor[A]k⋅FWTxor[B]k=(i∘k=0∑Ai−i∘k=1∑Ai)×j∘k=0∑Bj−j∘k=1∑Bj=i∘k=0∑Aij∘k=0∑Bj+i∘k=1∑Aij∘k=1∑Bj−i∘k=0∑Aij∘k=1∑Bj+i∘k=1∑Aij∘k=0∑Bj=(ixorj)∘k=0∑AiBj−(ixorj)∘k=1∑AiBj=FWTxor[C]k
而以上式子就是 沃尔什变换。
我们依然考虑分治。
设序列 A 的前一半为 A0,后一半为 A1。不难发现,A0 最高位都为 0,A1 的最高位都为 1。
当索引 i 的最高位为 0 时,无论与最高位为 0 还是 1 的 j 进行 ∘ 运算,因 0∩0=0,0∩1=0,奇偶性不变,因此这一部分为两者之和:
FWTxor[A]i=FWTxor[A0]i+FWTxor[A1]i
当索引 i 的最高位为 1 时,与最高位为 0 的 j 运算结果为 1∩0=0(奇偶性不变),与最高位为 1 的 j 运算结果为 1∩1=1(奇偶性翻转),因此这一部分为两者之差:
FWTxor[A]i=FWTxor[A0]i−FWTxor[A1]i
综上可得:
FWTxor[A]=merge(FWTxor[A0]+FWTxor[A1],FWTxor[A0]−FWTxor[A1])
逆变换易得:
IFWTxor[A]=merge(2IFWTxor[A0]+IFWTxor[A1],2IFWTxor[A0]−IFWTxor[A1])
同理,将 ∘ 定义中的 ∩ 替换为 ∪,即可得到同或:
FWTxnor[A]=merge(FWTxnor[A1]−FWTxnor[A0],FWTxnor[A1]+FWTxnor[A0])
IFWTxnor[A]=merge(2IFWTxnor[A1]−IFWTxnor[A0],2IFWTxnor[A1]+IFWTxnor[A0])
正变换和逆变换可以合并为一个函数:正变换取 type=1,t=1,逆变换取 type=−1,t=21。
Txor[A]=merge((Txor[A0]+Txor[A1])×t,(Txor[A0]−Txor[A1])×t)
Txnor[A]=merge((Txnor[A1]−Txnor[A0])×t,(Txnor[A1]+Txnor[A0])×t)
注意:模意义下的 21 指 2 的乘法逆元,即 2mod+1。
综上,可用分治写出 xor 和 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 和 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
- 数据生成:调用
std::mt19937_64
,生成 224≈1.68×107 个随机整数
- 测试方法:每个算法连续运行 100 次,记录耗时并计算平均值
类型/算法 | 倍增迭代 | 分治递归 | 高维前缀和 |
---|
或卷积(or) | 575.31 | 665.67(+15.7%) | 966.18(+67.9%) |
与卷积(and) | 559.65 | 693.12(+23.8%) | 961.38(+71.8%) |
异或卷积(xor) | 939.02 | 1167.31(+24.3%) | N/A |
同或卷积(xnor) | 937.16 | 1164.39(+24.2%) | 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;
}