1.13 上午-FWT
前言
勿让将来,辜负曾经
劝退了好叭!
正文
知识点
引子
如果屏幕前的你接触过 FFT 或者 NTT 的话,可能会更好理解这一类加速多项式卷积运算的优化算法
好叭,也许你没学过,单着并不影响我们学习 FWT。这个神奇的小算法是用于优化位运算卷积的,可以将 \(O(n^2)\) 的暴力枚举做法转化为 \(O(n \log n)\) 的优秀做法
为了方便表述,我们记多项式 \(A,B,C\) 的第 \(i\) 次项的系数 分别为 \(A_i,B_i,C_i\)
什么是位运算卷积?相信各位都解过高次方程,对于两个多项式 \(A,B\),计算 \(A \times B\) 的过程被称作加法卷积。形式化地,有:
强调:上面这个东西是加法卷积,不是位运算卷积!(加法卷积是 FFT 和 NTT 用于优化的问题)
卷积形式
注意到这个求和号下方有一个神奇的小式子,叫做 \(k=i+j\),那么我们把加号改成位运算,就是位运算卷积,于是乎就出现了——
- 或运算卷积
- 与运算卷积
- 异或运算卷积
接下来我们要尝试使用 FWT 对上述卷积进行复杂度优化
写在废话之前
如果没有接触过 FFT,你可以选择略过下面三个自然段
FFT 可以用于优化加法卷积,其原理完全可以类比到 FWT 上。一个多项式存在两种表示方法——系数表示法和点值表示法。
所谓的加法卷积就是系数表示法之间的某种运算,由于多项式长度都是 \(O(n)\) 的,所以朴素暴力的做法就是 \(O(n^2)\) 的。然而,描述一个多项式只需要 \(n\) 个点去表示即可(也就是点值表示法),时间复杂度 \(O(n)\)。
所以,FFT 能够优化时间复杂度的原理是:先将多项式的系数表示法转化为点值表示法(DFT),可以借助单位根加速(FFT),然后点值表示法进行运算,最后再将点值表示法转回系数表示法(IFFT)
FWT 也可以效仿 FFT,给多项式做一些变换(FWT),使得他们之间的卷积运算不再消耗大量的时间复杂度,然后再变换回去(IFWT)。具体地,我们希望完成如下操作:
针对于不同的位运算,FWT 变换有不同的构造方式。云落能力有限,无法给出构造的详细思维过程,只能给出构造的结论以及其正确性证明
FWT_or(FWT 加速或运算卷积)
构造方式:\(FWT[P]_i = \sum_{j \cup i = i} P_j\)
说人话就是,对于多项式 \(FWT[P]\) 的第 \(i\) 项,它等于原多项式 \(P\) 中所有 \(P_j\) 的和,其中 \(j\) 满足 j|i=i
(写成代码形式兴许会好理解些?)
先不考虑如何 \(O(n \log n)\) 解决这个变换问题,先考虑该构造的正确性。简单推一下式子:
卷积的部分搞定了,那么该怎么在优秀的复杂度中实现 \(P \to FWT[P]\) 捏?
考虑分治
位运算有一些非常好的性质,叫做按位运算(废话)。也就是说,我们可以尝试,顺次取出最高位,将一个规模较大的问题划分成规模较小的问题。然后噼里啪啦的就一顿分治就 \(O(n \log n)\) 力!
这么干讲不太好懂,还是湿讲比较好
钦定 \(FWT[P]\) 有 \(2^n\) 项(长度不足可以用系数 \(0\) 补齐),那么将这 \(2^n\) 项分为 \(FWT[P]_0,...,FWT[P]_{2^{n-1}-1}\) 与 \(FWT[P]_{2^{n-1}},...,FWT[P]_{2^n-1}\) 两部分,进行一些分类讨论……
为了方便表述,我们将原序列 \(P\) 也效仿上述划分方式,并记前半部分为 \(P_0\),后半部分为 \(P_1\)
对于 \(FWT[P]\) 来说,考虑前半部分的贡献来源。由于或运算的神奇性质,显然不可能由 \(P_1\) 贡献,也就只能由 \(P_0\) 贡献。因此,\((FWT[P]_0,...,FWT[P]_{2^{n-1}-1}) = FWT[P_0]\)。同理,考虑后半部分的贡献。很显然无论这一项是什么,他都会给 \(FWT[P]\) 的后半部分造成贡献。因此,\(FWT[P]_{2^{n-1}},...,FWT[P]_{2^n-1} = FWT[P_0] + FWT[P_1]\)
结合上述分析,我们发现:在或运算卷积下,\(FWT[P]\) 总是满足:
其中 merge
表示多项式的拼接,(可以类比字符串的拼接哈!)
从代码实现上来说,递归分治的写法肯定没问题啦,就是常数有点大,所以建议写成循环的形式,外面套一层长度枚举即可
正变换 FWT 说的差不多了,接下来该说逆变换 IFWT 力!
现在相当于已知 \(FWT[P]\),求 \(IFWT[FWT[P]]\)(也就是 \(P\))。还是回到刚才多项式拆半的过程,注意到贡献给 \(FWT[P]_0,...,FWT[P]_{2^{n-1}-1}\) 的只有 \(P_0\),所以反过来也是一样的。同理,贡献给 \(FWT[P]_{2^{n-1}},...,FWT[P]_{2^n-1}\) 的既有 \(P_0\) 又有 \(P_1\),那么求 \(P_1\) 就直接把 \(P_0\) 的部分刨去即可
形式化地,变换如下:
从代码实现上来说,FWT 和 IFWT 的结构极为相似,甚至可以说只差一个正负号,所以可以传入一个参数 type
来表示进行的是正变换还是逆变换,是一种节省码量的小技巧
然后就没有然后了,进度条 \(25 \% / 100 \%\)
FWT_and(FWT 加速与运算卷积)
有了上面或运算的铺垫,与运算讲起来可能更舒服一些
构造方式:\(FWT[P]_i = \sum_{j \cap i = i} P_j\)
正确性证明:
正变换:
逆变换:
与运算和或运算换汤不换药,自己手推吧(其实是云落敲不动了,好累 qaq)!
然后就没有然后了,进度条 \(50 \% / 100 \%\)
FWT_xor(FWT 加速异或运算卷积)
这个还是要好好说一下的!
为了方便表述,定义一种运算 \(\circ\),其运算法则可以被表述为 \(x \circ y = \operatorname{popcount}(x \cap y) \bmod 2\)
说人话就是,\(x,y\) 二进制按位与后的 “\(1\)” 的个数的奇偶性
该运算满足
-
交换律:\(x \circ y = y \circ x\)
-
结合律:\((x \circ y) \circ z = x \circ (y \circ z)\)
-
分配律:\(x \circ (y \oplus z) = (x \circ y) \oplus (x \circ z)\)
请注意分配律中间的符号表示异或运算!
构造方式:\(FWT[P]_i = \sum_{i \circ j = 0} P_j - \sum_{i \circ j = 1} P_j\)
正确性证明也很抽象,诸位凑合看,要是对公式过敏的话就把这个构造方式背下来
对于异或运算的 FWT 变换,仍旧是使用分治——对多项式 \(P\) 进行拆半操作
考虑 \(FWT[P]_0,...,FWT[P]_{2^{n-1}-1}\) 的贡献来源,发现不论是 \(P_0\) 还是 \(P_1\) 都可以给它造成贡献(建议自己手模一下,建议打 “\(\circ\)” 运算的真值表)。同理,对于 \(FWT[P]_{2^{n-1}},...,FWT[P]_{2^n-1}\),发现 \(P_0\) 可能会多余贡献,且多出来那部分恰好是 \(P_1\)(因为 \(0 \cap 1 = 0,1 \cap 1 = 1\)),所以要在 \(P_0\) 中刨去 \(P_1\) 给的那一部分贡献
综上 ,正变换递归式如下:
逆变换的式子写起来也不难。在若干年以前,即便是我们并不会使用二元一次方程组,但是我们依旧可以解决这个问题。
根据我们小学二年级就学过的…… ——毕导
已知两数之和,以及两数之差,求两个数分别是——?
钦定诸位都会上面这个问题,那么逆变换的形式可以写出来了,形如:
逆变换易得:
如果你听了云落的歪门邪道码量优化小技巧,要注意逆变换的时候 type
不要传入 \(-1\),而应传入 \(\frac{1}{2}\)
然后就没有然后了,进度条 \(75 \% / 100 \%\)
子集卷积
相比于或运算卷积,多了一条限制,形如 \(i \cap j = 0\)。做一步转化,结合另一条限制 \(i \cup j = k\),显然有 \(|i| + |j| = |k|\),也就是转化为了集合大小上的关系。所以直接把原来的多项式系数按照集合大小进行分类,对于相同集合大小的多项式系数进行 FWT 变换,最后多项式 \(A,B\) 合并的时候暴力加法卷积即可
建议结合代码食用
然后就没有然后了,进度条 \(100 \% / 100 \%\)
我们终会重逢
一题一解
T1 【模板】快速莫比乌斯/沃尔什变换 (FMT/FWT)(P4717)
本来想直接放题解代码的,谁复制粘贴谁就习题棕名大礼包,想想算了,还是放自己的代码哈!
点击查看代码
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int maxn=(1<<17)+5,mod=998244353,inv=499122177;
int n,A[maxn],B[maxn];
int a[maxn],b[maxn],c[maxn];
inline void init(){
memcpy(a,A,sizeof(A));
memcpy(b,B,sizeof(B));
memset(c,0,sizeof(c));
return;
}
inline void mul(){
for(int i=0;i<n;i++){
c[i]=a[i]*b[i]%mod;
}
return;
}
inline void output(){
for(int i=0;i<n;i++){
cout<<c[i]<<' ';
}
cout<<endl;
return;
}
inline void FWT_or(int *P,int type){
for(int i=1;i<n;i<<=1){
int p=i<<1;
for(int j=0;j<n;j+=p){
for(int k=0;k<i;k++){
P[i+j+k]=(P[i+j+k]+(P[j+k]*type%mod+mod)%mod)%mod;
}
}
}
return;
}
inline void FWT_and(int *P,int type){
for(int i=1;i<n;i<<=1){
int p=i<<1;
for(int j=0;j<n;j+=p){
for(int k=0;k<i;k++){
P[j+k]=(P[j+k]+(P[i+j+k]*type%mod+mod)%mod)%mod;
}
}
}
return;
}
inline void FWT_xor(int *P,int type){
for(int i=1;i<n;i<<=1){
int p=i<<1;
for(int j=0;j<n;j+=p){
for(int k=0;k<i;k++){
int X=P[j+k],Y=P[i+j+k];
P[j+k]=(X+Y)%mod;
P[i+j+k]=(X-Y+mod)%mod;
P[j+k]=(P[j+k]*type)%mod;
P[i+j+k]=(P[i+j+k]*type)%mod;
}
}
}
return;
}
inline void OR(){
init();
FWT_or(a,1);
FWT_or(b,1);
mul();
FWT_or(c,-1);
output();
return;
}
inline void AND(){
init();
FWT_and(a,1);
FWT_and(b,1);
mul();
FWT_and(c,-1);
output();
return;
}
inline void XOR(){
init();
FWT_xor(a,1);
FWT_xor(b,1);
mul();
FWT_xor(c,inv);
output();
return;
}
signed main(){
ios::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);
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];
}
OR();
AND();
XOR();
return 0;
}
T2 【模板】子集卷积(P6097)
板,不想多说,现在这个东西不卡常数咯!
温馨提示:模数是 \(10^9 + 9\) 而不是 \(10^9 +7\)
点击查看代码
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int MAXN=22,maxn=(1<<20)+2,mod=1e9+9;
int n,a[MAXN][maxn],b[MAXN][maxn];
int len,c[MAXN][maxn];
inline int lowbit(int x){
return x&(-x);
}
inline int cal(int x){
int res=0;
while(x){
res++;
x-=lowbit(x);
}
return res;
}
inline void FWT_or(int *P,int type){
for(int i=1;i<n;i<<=1){
int p=i<<1;
for(int j=0;j<n;j+=p){
for(int k=0;k<i;k++){
P[i+j+k]=(P[i+j+k]+(P[j+k]*type%mod+mod)%mod)%mod;
}
}
}
return;
}
inline void Juan(){
for(int i=0;i<=len;i++){
for(int j=0;j<=i;j++){
for(int k=0;k<n;k++){
c[i][k]=(c[i][k]+a[j][k]*b[i-j][k]%mod)%mod;
}
}
}
return;
}
signed main(){
ios::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);
cin>>n;
len=n;
n=(1<<n);
for(int i=0;i<n;i++){
cin>>a[cal(i)][i];
}
for(int i=0;i<n;i++){
cin>>b[cal(i)][i];
}
for(int i=0;i<=len;i++){
FWT_or(a[i],1);
FWT_or(b[i],1);
}
Juan();
for(int i=0;i<=len;i++){
FWT_or(c[i],-1);
}
for(int i=0;i<n;i++){
cout<<c[cal(i)][i]<<' ';
}
return 0;
}
T3 卡牌(P8292)
不大会 FWT 的做法……只会大力容斥
如果你们还对寿司晚宴有印象的话,显然有一个比较经典的 trick,形如根号分治(云落也不知道这么说是否准确)
对于 \(\ge 43\) 的质因数,最多出现一次。所以小质数状态压缩,大质数特殊处理
计算方案,考虑容斥
显然有 \(2^n\) 减去至少不含一个小质数的,加上至少不含两个小质数……显然是对于小质数来说,正确性和复杂度都有保证
现在把大质数的贡献加进来,注意到每次询问最多需要一个大质数。也就是对于每个大质数,按需求给答案 \(\times 2^{cnt_i}\) 或者 \(\times (2^{cnt_i} - 1)\)
剩下的就交给容斥即可
点击查看代码
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int maxn=1e6+5,maxv=2e3+5,maxp=13,mod=998244353;
int n,m,cnt[maxv];
vector<int> que;
bool isprime[maxv];
int prime[maxv],rev[maxn],tot;
vector<int> fac[maxn];
int p[maxv];
int g[(1<<maxp)+5],f[(1<<maxp)+5][305];
inline void euler(){
memset(isprime,true,sizeof(isprime));
isprime[1]=false;
for(int i=2;i<maxv;i++){
if(isprime[i]){
prime[++tot]=i;
rev[i]=tot;
}
for(int j=1;j<=tot&&i*prime[j]<maxv;j++){
isprime[i*prime[j]]=false;
if(i%prime[j]==0){
break;
}
}
}
return;
}
inline int qpow(int a,int b){
int res=1;
while(b){
if(b&1){
res=res*a%mod;
}
a=a*a%mod;
b>>=1;
}
return res;
}
signed main(){
ios::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);
cin>>n;
for(int i=1;i<=n;i++){
int x;
cin>>x;
cnt[x]++;
}
euler();
for(int i=2;i<maxv;i++){
for(int j=1;j<=tot;j++){
if(i%prime[j]==0){
if(j<=13){
p[i]|=(1<<j-1);
}
fac[i].push_back(j);
}
}
}
for(int i=0;i<(1<<maxp);i++){
for(int j=2;j<maxv;j++){
if(p[j]&i){
continue;
}
int x=fac[j][fac[j].size()-1];
f[i][x]=(f[i][x]+cnt[j])%mod;
g[i]=(g[i]+cnt[j])%mod;
}
g[i]=(g[i]+cnt[1])%mod;
}
cin>>m;
while(m--){
que.clear();
int c;
cin>>c;
int state=0;
for(int i=1;i<=c;i++){
int x;
cin>>x;
que.push_back(x);
if(rev[x]<=13){
state|=(1<<rev[x]-1);
}
}
sort(que.begin(),que.end());
int res=0;
for(int i=0;i<(1<<13);i++){
if((i|state)!=state){
continue;
}
int tol=g[i],ans=1;
for(int j=0;j<c;j++){
if(que[j]<=41){
continue;
}
int x=f[i][rev[que[j]]];
ans=ans*(qpow(2,x)-1)%mod;
tol=(tol-x+mod)%mod;
}
ans=(ans*qpow(2,tol))%mod;
int sum=0;
for(int j=0;j<=13;j++){
sum=(sum+(i>>j&1))%mod;
}
if(sum&1){
res=(res-ans+mod)%mod;
}else{
res=(res+ans)%mod;
}
}
cout<<res<<endl;
}
return 0;
}
后记
T3 的 FWT 做法留坑后填(主要是云落只能理解 T3 第一篇 luogu 题解的做法,其它的太抽象了)
完结撒花!