FFT/NTT/FWT重學or新學？
阿新 • • 發佈：2022-01-02
不知不覺，就開始學習多項式了，好快啊
重新拾起我之前學的百嘛不是的\(FFT\)，新學\(NTT/FWT\)
開始吧學習部落格
FFT:快速傅立葉變換
總的來說就是快速求取兩個多項式的乘積
把板子粘在這裡Luogu3803 多項式乘法
code
#include<bits/stdc++.h> using namespace std; #define fo(i,x,y) for(int i=(x);i<=(y);i++) #define fu(i,x,y) for(int i=(x);i>=(y);i--) const double pi=acos(-1.0); const int N=1<<21; struct coex{ double r,i; coex(){} coex(double x,double y){r=x;i=y;} coex operator + (coex a){return coex(r+a.r,i+a.i);} coex operator - (coex a){return coex(r-a.r,i-a.i);} coex operator * (coex a){return coex(r*a.r-i*a.i,r*a.i+i*a.r);} }a[N],b[N],w[N]; int af[N],la,lb,lim,len; void fft(coex *a,int lim){ fo(i,0,lim-1)if(af[i]>i)swap(a[i],a[af[i]]); for(int t=lim>>1,d=1;d<lim;d<<=1,t>>=1) for(int i=0;i<lim;i+=(d<<1)) fo(j,0,d-1){ coex tmp=w[t*j]*a[i+j+d]; a[i+j+d]=a[i+j]-tmp; a[i+j]=a[i+j]+tmp; } } signed main(){ scanf("%d%d",&la,&lb); fo(i,0,la)scanf("%lf",&a[i].r); fo(i,0,lb)scanf("%lf",&b[i].r); for(lim=1,len=0;lim<=la+lb;lim<<=1,len++); fo(i,0,lim-1){ af[i]=(af[i>>1]>>1)|((i&1)<<(len-1)); w[i]=coex(cos(2.0*i*pi/lim),sin(2.0*i*pi/lim)); } fft(a,lim);fft(b,lim); fo(i,0,lim-1)a[i]=a[i]*b[i],w[i].i=-w[i].i; fft(a,lim); fo(i,0,la+lb)printf("%d ",(int)(a[i].r/lim+0.5)); }
需要注意這麼幾點:
1、\(lim\)要大於多項式次數，嚴格大於
2、搞清楚\(d\)和\(t\)的關係，以及各種2倍
3、背板子，背就完事了。
4、最後答案別忘了除以\(lim\)，預處理單位根別忘了\(2*i*pi/lim\)，別忘了除以\(lim\)
NTT:快速數論變換
code
void ntt(int *a,int lim){ fo(i,0,lim-1)if(af[i]>i)swap(a[i],a[af[i]]); for(int t=lim>>1,d=1;d<lim;d<<=1,t>>=1) for(int i=0;i<lim;i+=(d<<1)) fo(j,0,d-1){ int tmp=g[t*j]*a[i+j+d]%mod; a[i+j+d]=(a[i+j]-tmp+mod)%mod; a[i+j]=(a[i+j]+tmp)%mod; } } int bas[M],ans[M]; int a[M],b[M]; void mul(int *x,int *y,int *z){ fo(i,0,lim-1)a[i]=y[i]; fo(i,0,lim-1)b[i]=z[i]; g[0]=1;g[1]=ksm(3,(mod-1)/lim,mod); fo(i,2,lim-1)g[i]=g[i-1]*g[1]%mod; ntt(a,lim);ntt(b,lim); g[0]=1;g[1]=ksm(g[1],mod-2,mod); fo(i,2,lim-1)g[i]=g[i-1]*g[1]%mod; fo(i,0,lim-1)a[i]=a[i]*b[i]%mod; ntt(a,lim);int inv=ksm(lim,mod-2,mod); fo(i,0,lim-1)a[i]=a[i]*inv%mod; fo(i,0,m-2)a[i]=(a[i]+a[i+m-1])%mod; fo(i,0,m-2)x[i]=a[i]; }
注意：
1、注意模數，以及自己找原根
2、注意最後出來的時候除以\(lim\)
3、別忘了變換原根
多項式求逆
就是求某一個多項式的逆元\(mod\ x^n\)意義下
主要吧，這個板子還挺難背的
遞推上來的證明過程
所以結論就是\(B=2B'-AB'^2\)(\(B'\)表示\(mod\ x^{\lceil\frac{n}{2}\rceil}\)意義下的逆元)
code
#include<bits/stdc++.h> using namespace std; #define int long long #define fo(i,x,y) for(int i=(x);i<=(y);i++) #define fu(i,x,y) for(int i=(x);i>=(y);i--) int read(){ int s=0,t=1;char ch=getchar(); while(!isdigit(ch)){if(ch=='-')t=-1;ch=getchar();} while(isdigit(ch)){s=s*10+ch-'0';ch=getchar();} return s*t; } const int mod=998244353; const int N=1<<18; int ksm(int x,int y){ int ret=1; while(y){ if(y&1)ret=ret*x%mod; x=x*x%mod;y>>=1; }return ret; } int af[N],g[N],lim,len; void ntt(int *a,int lim){ fo(i,0,lim-1)if(af[i]>i)swap(a[i],a[af[i]]); for(int t=lim>>1,d=1;d<lim;d<<=1,t>>=1) for(int i=0;i<lim;i+=(d<<1)) fo(j,0,d-1){ int tmp=g[t*j]*a[i+j+d]%mod; a[i+j+d]=(a[i+j]-tmp+mod)%mod; a[i+j]=(a[i+j]+tmp)%mod; } } int a[N],b[N],c[N]; void calc(int deg){//deg表示當前mod的是x^deg if(deg==1){return b[0]=ksm(a[0],mod-2),void();} calc((deg+1)>>1); for(lim=1,len=0;lim<deg*2;lim<<=1,len++); fo(i,0,lim-1)af[i]=(af[i>>1]>>1)|((i&1)<<(len-1)); g[0]=1;g[1]=ksm(3,(mod-1)/lim); fo(i,2,lim-1)g[i]=g[i-1]*g[1]%mod; fo(i,0,deg-1)c[i]=a[i];//注意這裡只賦值前deg-1位，因為這次乘法只能用到這些位 //這裡a有deg-1位，而b只有deg+1>>1位，因為b是要平方的，並且這裡是求a的deg-1位的逆元 ntt(c,lim);ntt(b,lim); g[0]=1;g[1]=ksm(g[1],mod-2); fo(i,2,lim-1)g[i]=g[i-1]*g[1]%mod; fo(i,0,lim-1)b[i]=(2-c[i]*b[i]%mod+mod)*b[i]%mod; ntt(b,lim);int iv=ksm(lim,mod-2); fo(i,0,lim-1)b[i]=b[i]*iv%mod; fo(i,deg,lim-1)b[i]=0;//把後邊的清空，防止影響下一次運算 fo(i,0,lim-1)c[i]=0;//清空 } int n; signed main(){ n=read(); fo(i,0,n-1)a[i]=read(); calc(n); fo(i,0,n-1)printf("%lld ",b[i]); }
FWT:快速沃爾什變換
發現我們所做的卷積都是\(c=\sum_{i=1}^{n}\sum_{j=1}^{i}a_j*b_{i-j}\)
也就是\(FFT\)和\(NTT\)能處理的加法卷積，如果是與？或？異或呢？
就要用到\(FWT\)，也就是快速位運算卷積
|:或卷積
這個要是直接卷的話，好像有點困難，但是高維字首和就是或卷積，所以我們可以用高維字首和解決或卷積的問題
然而高維字首和還是沒有\(FWT\)來的直接
\(FWT\)也和\(FFT\)一樣有變換和逆變換，這個直接看最上面我給出的那個部落格吧
這次變換不能叫點值表示式了，應該叫它子集和表示式
相當於是列舉每一個二進位制位來求子集和，逆變換拆回去的時候是一樣的
code
void ort(int *a,int lim,int tp){//tp=1表示順變換,tp=2表示逆變換
for(int d=1;d<lim;d<<=1)
for(int i=0;i<lim;i+=(d<<1))
fo(j,0,d-1){
a[i+j+d]=(a[i+j+d]+a[i+j]*tp)%mod;
}
}
&:與卷積
這應該叫做高維字尾和吧，我們可以用高維字尾和解決與卷積問題
變換就是字尾和，逆變換一樣
code
void andt(int *a,int lim,int tp){//tp=1表示順變換,tp=2表示逆變換
for(int d=1;d<lim;d<<=1)
for(int i=0;i<lim;i+=(d<<1))
fo(j,0,d-1){
a[i+j]=(a[i+j]+a[i+j+d]*tp)%mod;
}
}
^:異或卷積
這個應該叫高維字尾差分
好像有點難以概括，額，不概括了，直接看碼吧
code
void xort(int *a,int lim,int tp){//tp=1表示順變換,tp=0.5表示逆變換
for(int d=1;d<lim;d<<=1)
for(int i=0;i<lim;i+=(d<<1))
fo(j,0,d-1){
int tmp=a[i+j+d];
a[i+j+d]=(a[i+j]-tmp+mod)*tp%mod;
a[i+j]=(a[i+j]+tmp)*tp%mod;
}
}
高維字首和的求法，列舉每一維，做一遍字首和
在位運算上也就是列舉每一個二進位制位，將這一位為\(0\)的加到這一位為\(1\)的上面
這樣做高維字首和雖然樸素，但是有正確性的保證，保證了不重不漏
code
#include<bits/stdc++.h>
using namespace std;
#define int long long
#define fo(i,x,y) for(int i=(x);i<=(y);i++)
#define fu(i,x,y) for(int i=(x);i>=(y);i--)
int read(){
int s=0,t=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')t=-1;ch=getchar();}
while(isdigit(ch)){s=s*10+ch-'0';ch=getchar();}
return s*t;
}
const int mod=998244353;
const int N=1<<18;
int ksm(int x,int y){
int ret=1;
while(y){
if(y&1)ret=ret*x%mod;
x=x*x%mod;y>>=1;
}return ret;
}
int n,lim,a[N],b[N],c[N],aa[N],bb[N];
void init(int lim){fo(i,0,lim-1)a[i]=aa[i],b[i]=bb[i];}
void ort(int *a,int lim,int tp){
for(int d=1;d<lim;d<<=1)
for(int i=0;i<lim;i+=(d<<1))
fo(j,0,d-1){
a[i+j+d]=(a[i+j+d]+a[i+j]*tp)%mod;
}
}
void andt(int *a,int lim,int tp){
for(int d=1;d<lim;d<<=1)
for(int i=0;i<lim;i+=(d<<1))
fo(j,0,d-1){
a[i+j]=(a[i+j]+a[i+j+d]*tp)%mod;
}
}
void xort(int *a,int lim,int tp){
for(int d=1;d<lim;d<<=1)
for(int i=0;i<lim;i+=(d<<1))
fo(j,0,d-1){
int tmp=a[i+j+d];
a[i+j+d]=(a[i+j]-tmp+mod)*tp%mod;
a[i+j]=(a[i+j]+tmp)*tp%mod;
}
}
void mix(int lim){fo(i,0,lim-1)c[i]=a[i]*b[i]%mod;}
void pt(int *a,int lim){fo(i,0,lim-1)printf("%lld ",a[i]);printf("\n");}
signed main(){
n=read();lim=1<<n;
fo(i,0,lim-1)aa[i]=read();
fo(i,0,lim-1)bb[i]=read();
init(lim);ort(a,lim,1);ort(b,lim,1);mix(lim);ort(c,lim,mod-1);pt(c,lim);
init(lim);andt(a,lim,1);andt(b,lim,1);mix(lim);andt(c,lim,mod-1);pt(c,lim);
init(lim);xort(a,lim,1);xort(b,lim,1);mix(lim);xort(c,lim,ksm(2,mod-2));pt(c,lim);
}
技巧
1、生成函式時，有關於對稱的問題，如果對稱的話，那麼兩個點的位置相加是對稱中心的兩倍，意思就是多項式卷積之後對稱軸是同一條的都在一起
2、對於生成函式的利用，我們只需要係數，而並不是真正的把\(x\)帶進去找值，我只需要係數就行了
比如：\(\sum_{i=0}^{n}\sum_{j=0}^{i}j!*2^{i-j}\)，這個東西階乘和冪可以分別看作是兩個多項式的係數
這樣我們可以通過卷積快速求出後面的和，也就是卷積後多項式的第\(i\)項係數，符合多項式乘法的規則，然後帶入求和
QQ:2953174821