FWT
FWT
一、用途
一般的\(FFT\)做的是形如\(F_k=\sum_{i+j=k}g_ih_j\)的卷积，但如果我们遇到的是其他的操作符呢？
比如说：
如果运算符从不同的加法变成了位运算，该怎样做？
答案是——FWT。
二、实现方式
考虑\(FFT\)，我们通过对多项式求出它在特殊位置的点值然和将点值相乘，最后复原
那么我们同样希望对于位运算卷积也能通过某种方式转化为点值相乘，也就是将\(A\)转化为\(FWT(A)\)，然后直接将\(FWT(A)\)与\(FWT(B)\)相乘得到\(FWT(C)\)，最后再还原。
那么如何得到这个\(FWT\)多项式呢，我们依次考虑这三个不同的位运算：
（一）或卷积
即\(C_k=\sum_{i|j=k}A_iB_j\)
因为\(i|j=k\)，所以说\(i\)中二进制中为\(1\)的位置与\(j\)二进制位\(1\)的位置一定是\(k\)中的位置的子集。
于是我们可以构造出:
即满足二进制中为\(1\)的位置是\(i\)的子集的所有\(j\)的集合
那么易证：
于是这个构造是可行的。
构造是有了，那么如何快速求出\(FWT(f)\)呢?
考虑将\(A\)这个区间分成，下标二进制下最高位为\(0\)的部分\(A_0\)，与最高位为\(1\)的部分\(A_1\)
那么根据定义容易得到：
其中\(merge\)是像字符串拼接一样将它拼起来的过程
还原就反过来就行了：
（二）与卷积
类比可得：
inline void FWT_and(int *f,int tp,int tot){
for(int i=1;i<tot;i<<=1)
for(int j=0;j<tot;j+=(i<<1))
for(int k=0;k<i;++k){
if(tp==1) f[j+k]=add(f[j+k],f[j+i+k]);
else f[j+k]=dec(f[j+k],f[j+i+k]);
}
}
（三）异或卷积
异或运算是最常考，也是相对较难的运算：
这个运算的公式我太不会推，这里先给出结论，然后尝试证明：
定义\(i\otimes j\)为\(i\&j\)的二进制下\(1\)的数量的奇偶性：
证明：
于是我们就有：
inline void FWT_xor(int *f,int tp,int tot){
for(int i=1;i<tot;i<<=1)
for(int j=0;j<tot;j+=(i<<1))
for(int k=0;k<i;++k){
int x=f[j+k],y=f[j+i+k];
if(tp==1) f[j+k]=add(x,y),f[j+i+k]=dec(x,y);
else f[j+k]=1ll*iv2*add(x,y)%mod,f[j+i+k]=1ll*iv2*dec(x,y)%mod;
}
}
大家可以去尝试一下这道例题，套模板就行了。
（四）子集卷积
子集卷积求的是
子集卷积相比普通的或卷积增加了一个 \(i\&j=0\) 的条件，这个条件就相当于 \(|i|+|j|=|k|\)，因此我们可以多记一维表示该集合的大小，设 \(f_{i,j}=[|j|=i]a_j,g_{i,j}=[|j|=i]b_j\)，那么原问题即相当于求：
于是就可以通过 \(\mathcal O(n)\) 次 FWT 求解了，复杂度为 \(\mathcal O(2^nn^2)\)。
#include<bits/stdc++.h>
using namespace std;
const int N=(1<<20)+10;
const int mod=1e9+9;
int n,a[21][N],b[21][N],c[21][N],siz[N],lim;
inline int add(int x,int y){return (x+y>=mod)?x+y-mod:x+y;}
inline int dec(int x,int y){return (x-y<0)?x-y+mod:x-y;}
inline void FMT(int n,int *f,int tp){
for(int i=0;i<n;++i)
for(int j=0;j<lim;++j)
if(j&(1<<i)){
if(tp==1) f[j]=add(f[j],f[j^(1<<i)]);
else f[j]=dec(f[j],f[j^(1<<i)]);
}
}
int main(){
scanf("%d",&n);lim=1<<n;
for(int i=0;i<n;++i)
for(int j=0;j<lim;++j)
if(j&(1<<i)) ++siz[j];
for(int i=0;i<lim;++i) scanf("%d",&a[siz[i]][i]);
for(int i=0;i<lim;++i) scanf("%d",&b[siz[i]][i]);
for(int i=0;i<=n;++i){
FMT(n,a[i],1);FMT(n,b[i],1);
for(int j=0;j<lim;++j)
for(int k=0;k<=i;++k)
c[i][j]=add(c[i][j],1ll*a[k][j]*b[i-k][j]%mod);
}
for(int i=0;i<=n;++i) FMT(n,c[i],-1);
for(int i=0;i<lim;++i) printf("%d ",c[siz[i]][i]);
return 0;
}
三、k 进制FWT
注意到异或看作二进制下的加法，于是异或卷积实际上是一个多元 \(FFT\) 的过程，将每个数二进制下的每一位看作不同的变元，每一轮就是对这一变元做 \(FFT\)，乘上的矩阵：
实际上就是 \(FFT\) 的范德蒙德矩阵：
\(k\) 进制 \(FWT\) ，是求解关于 \(k\) 进制不进位加法的卷积，也同样是一个多元 \(FFT\)，在第 \(i\) 轮只需要将所有数进行分组，\(k^i\) 个数为一组，取出每组下表前 \(k-1\) 位都相同的 \(k\) 个数，对它们做长度为 \(k\) 的 \(FFT\)，也就是乘上：
\(IFWT\) 也同样是乘上 \(IDFT\) 的范德蒙德矩阵，然后再除上 \(k\)。
设 \(FWT\) 的值域为 \(k^m\)，则复杂度为 \(\mathcal O(k^{m+1}m)\)。
一道不完全模板的模板题：CF1103E：
不难发现所求就是 \((\sum_ix^{a_i})^n\)。这里的乘法是十进制不进位加法卷积。
一实现你会发现有以下两个问题：
-
\(10\) 在 \(\bmod 2^{58}\) 下没有逆元。
这个过程我们实际上要对每个数除上 \(10^5\)，注意到 \(5\) 是有逆元的，因此只需要考虑 \(2^5\)，这可以先计算出 \(\bmod 2^{63}\) 下的答案，最后直接除 \(2^5\) 即可。当然为了方便代码中用 \(unsigned\ long\ long\) 直接计算了 \(\bmod 2^{64}\)。
-
\(k\) 次单位根在模意义下不存在。
对于这个问题，考虑扩域，记 \(w\) 为 \(k\) 的答案为根满足 \(w^k=1\)。将每个数即为关于 \(w\) 的 \(k-1\) 次多项式，此时数的乘法就是两个多项式 \(\bmod {w^k-1}\) 下的卷积。
但由于 \(w^k-1\) 可约，一个数可能会有多种表示方法，而我们只需要那个只有第 \(0\) 项有值的表示。这令我们想到了分圆多项式：\(\Phi_k(\omega)\)，它满足 \(\omega\) 的阶为 \(k\) 且在 \(Q\) 上不可约。于是我们只需要让所有乘法在模 \(\Phi_k(\omega)\) 下进行就行了。
但如果直接这样做，每次乘法都要跑一个多项式取模太慢了。由于 \(\Phi_k(\omega)|\omega^k-1\)，我们不妨直接在 \(\bmod \omega^k-1\) 意义下进行计算，最后再对 \(\Phi_k(\omega)\) 取模即可。复杂度就是 \(\mathcal O(k^{m+2}m)\) 的了（\(k,m\) 与上文定义相同）。
还有一个问题是分圆多项式的计算：
\[x^n-1=\prod_{d|n}\Phi_d(x) \]取 \(\ln\)：
\[\ln (x^n-1)=\sum_{d|n}\ln\Phi_d(x)\\ \]后莫比乌斯反演：
\[\ln\Phi_n(x)=\sum_{d|n}\mu(\dfrac nd)\ln(x^n-1) \]再 \(exp\) 回来：
\[\Phi_n(x)=\prod_{d|n}(x^d-1)^{\mu(\frac nd)} \]然后就可以暴力枚举 \(k\) 的约数做多项式乘法/除法即可，由于乘的多项式只有两项有值，因此这部分可以做到 \(\mathcal O(k)\)。就可以用 \(\mathcal O(2^{\omega(k)}k)\) 的复杂度求解分圆多项式了：
#include<bits/stdc++.h> using namespace std; typedef unsigned long long ull; const int N=2e5+10; const ull iv5=14757395258967641293ull; const int pw[]={1,10,100,1000,10000,100000},lim=100000; int k,n; struct node{ ull a[11]; node(){memset(a,0,sizeof(a));} inline void clr(){memset(a,0,sizeof(a));} }f[200010]; inline node operator *(const node &x,const node &y){ node ret; for(int i=0;i<k;++i) for(int j=0;j<k;++j) ret.a[(i+j)%k]+=x.a[i]*y.a[j]; return ret; } inline node operator +(const node &x,const node &y){ node ret; for(int i=0;i<k;++i) ret.a[i]=x.a[i]+y.a[i]; return ret; } inline node ksm(node x,int k){ node ret;ret.a[0]=1; for(;k;k>>=1,x=x*x) if(k&1) ret=ret*x; return ret; } inline node operator <<(const node &x,int y){ node ret; for(int i=0;i<k;++i) ret.a[(i+y)%k]+=x.a[i]; return ret; } inline node operator >>(const node &x,int y){ node ret; for(int i=0;i<k;++i) ret.a[((i-y)%k+k)%k]+=x.a[i]; return ret; } inline void FFT(node *f,int tp){ for(int i=0;i<5;++i){ for(int j=0;j<lim;j+=pw[i+1]){ for(int w=0;w<pw[i];++w){ static node tmp[11],ret[11]; for(int t=j+w,x=0;t<j+pw[i+1];t+=pw[i],++x) tmp[x]=f[t]; if(tp==1){ for(int a=0;a<k;++a){ ret[a].clr(); for(int b=0;b<k;++b) ret[a]=ret[a]+(tmp[b]<<(a*b)); } } else{ for(int a=0;a<k;++a){ ret[a].clr(); for(int b=0;b<k;++b) ret[a]=ret[a]+(tmp[b]>>(a*b)); } } for(int t=j+w,x=0;t<j+pw[i+1];t+=pw[i],++x) f[t]=ret[x]; } } } } int p[N],pd[N],mu[N],pcnt,len; node Mod; inline void init(){ mu[1]=1; for(int i=2;i<=k;++i){ if(!pd[i]) p[++pcnt]=i,mu[i]=-1; for(int j=1;j<=pcnt&&i*p[j]<=k;++j){ int num=i*p[j];pd[num]=1; if(i%p[j]) mu[num]=-mu[i]; else{ mu[num]=0; break; } } } Mod.a[0]=1; for(int i=1;i<=k;++i){ if(k%i==0){ if(mu[k/i]==1){ for(int j=k-1;j>=0;--j) Mod.a[j]=((j>=i)?Mod.a[j-i]:0)-Mod.a[j]; } else if(mu[k/i]==-1){ for(int j=0;j<k;++j) Mod.a[j]=((j>=i)?Mod.a[j-i]:0)-Mod.a[j]; } } } len=k; while(!Mod.a[len-1]) --len; } inline ull Divn(node f){ for(int i=k-1;i>=len-1;--i){ ull tmp=f.a[i]; for(int w=i,j=len-1;j>=0;--w,--j) f.a[w]-=tmp*Mod.a[j]; } ull ans=f.a[0]; for(int i=1;i<=5;++i) ans*=iv5; return (ans>>5)%(1ull<<58); } int main(){ scanf("%d",&n); for(int i=1,x;i<=n;++i) scanf("%d",&x),f[x].a[0]++; k=10; FFT(f,1); for(int i=0;i<lim;++i) f[i]=ksm(f[i],n); FFT(f,-1); init(); for(int i=0;i<n;++i) printf("%llu\n",Divn(f[i])); return 0; }
