快速傅里叶变换(FFT)相关内容汇总
分类:
IT文章
•
2025-02-04 20:23:07
快速傅里叶变换，是求两个多项式卷积的算法，其时间复杂度为$O(nlog n)$，优于普通卷积求法，且根据有关证明，快速傅里叶变换是基于变换求卷积的理论最快算法。
关于FFT的介绍，最详细易懂的是《算法导论》上的内容。
其大致介绍与代码在这里：http://www.cnblogs.com/rvalue/p/7351400.html.
1.FFT&NTT模板
![]()
1 #include<cmath>
2 #include<cstdio>
3 #include<algorithm>
4 #define rep(i,l,r) for (int i=l; i<=r; i++)
5 using namespace std;
6
7 const int N=300010;
8 const double pi=acos(-1.);
9 int n,m,L,x,rev[N];
10
11 struct C{
12 double x,y;
13 C(double _x=0,double _y=0):x(_x),y(_y){}
14 C operator +(C &b){ return C(x+b.x,y+b.y); }
15 C operator -(C &b){ return C(x-b.x,y-b.y); }
16 C operator *(C &b){ return C(x*b.x-y*b.y,x*b.y+b.x*y); }
17 }a[N],b[N];
18
19 void DFT(C a[],int n,bool f){
20 for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
21 for (int i=1; i<n; i<<=1){
22 C wn=C(cos(pi/i),(f?1:-1)*sin(pi/i));
23 for (int p=i<<1,j=0; j<n; j+=p){
24 C w(1,0);
25 for (int k=0; k<i; k++,w=w*wn){
26 C x=a[j+k],y=w*a[i+j+k]; a[j+k]=x+y; a[i+j+k]=x-y;
27 }
28 }
29 }
30 if (!f) for (int i=0; i<n; i++) a[i].x/=n;
31 }
32
33 int main(){
34 freopen("Dft.in","r",stdin);
35 freopen("Dft.out","w",stdout);
36 scanf("%d%d",&n,&m);
37 rep(i,0,n) scanf("%d",&x),a[i]=x;
38 rep(i,0,m) scanf("%d",&x),b[i]=x;
39 rep(i,0,max(n,m)) a[i]=C(a[i].x+b[i].x,a[i].x-b[i].x);
40 m=n+m; for (n=1; n<=m; n<<=1) L++;
41 rep(i,0,n-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
42 DFT(a,n,1); rep(i,0,n-1) a[i]=a[i]*a[i]; DFT(a,n,0);
43 rep(i,0,m) printf("%d ",(int)(a[i].x/4+0.5));
44 return 0;
45 }
FFT
![]()
1 void NTT(int a[],int n,bool f){
2 for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
3 for (int i=1; i<n; i<<=1){
4 int wn=ksm(3,f ? (mod-1)/(i<<1) : (mod-1)-(mod-1)/(i<<1));
5 for (int p=i<<1,j=0; j<n; j+=p){
6 int w=1;
7 for (int k=0; k<i; k++,w=1ll*w*wn%mod){
8 int x=a[j+k],y=1ll*w*a[i+j+k]%mod;
9 a[j+k]=(x+y)%mod; a[i+j+k]=(x-y+mod)%mod;
10 }
11 }
12 }
13 if (f) return;
14 int inv=ksm(n,mod-2);
15 for (int i=0; i<n; i++) a[i]=1ll*a[i]*inv%mod;
16 }
NTT
![]()
1 #include<cmath>
2 #include<cstdio>
3 #include<algorithm>
4 #include<cstring>
5 #define mem(a) memset(a,0,sizeof(a))
6 #define rep(i,l,r) for (int i=(l); i<=(r); i++)
7 using namespace std;
8
9 const int N=500010,mod=998244353,inv2=(mod+1)/2;
10 int n,k,rev[N],inv[N],X[N],Y[N],A[N],B[N],C[N],D[N],E[N],F[N],G[N];
11
12 void Print(int a[],int n=::n){ for (int i=0; i<n; i++) printf("%d ",a[i]); puts(""); }
13
14 int ksm(int a,int b){
15 int res=1;
16 for (; b; a=1ll*a*a%mod,b>>=1)
17 if (b & 1) res=1ll*res*a%mod;
18 return res;
19 }
20
21 void NTT(int a[],int n,bool f){
22 for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
23 for (int i=1; i<n; i<<=1){
24 int wn=ksm(3,f ? (mod-1)/(i<<1) : (mod-1)-(mod-1)/(i<<1));
25 for (int p=i<<1,j=0; j<n; j+=p){
26 int w=1;
27 for (int k=0; k<i; k++,w=1ll*w*wn%mod){
28 int x=a[j+k],y=1ll*w*a[i+j+k]%mod;
29 a[j+k]=(x+y)%mod; a[i+j+k]=(x-y+mod)%mod;
30 }
31 }
32 }
33 if (f) return;
34 int inv=ksm(n,mod-2);
35 for (int i=0; i<n; i++) a[i]=1ll*a[i]*inv%mod;
36 }
37
38 void mul(int a[],int b[],int l){
39 int n=1,L=0;
40 for (; n<(l<<1); n<<=1) L++;
41 for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
42 NTT(a,n,1); NTT(b,n,1);
43 for (int i=0; i<n; i++) a[i]=1ll*a[i]*b[i]%mod;
44 NTT(a,n,0); NTT(b,n,0);
45 }
46
47 void Inv(int a[],int b[],int l){
48 if (l==1){ b[0]=ksm(a[0],mod-2); return; }
49 Inv(a,b,l>>1); int n=1,L=0;
50 for (; n<(l<<1); n<<=1) L++;
51 for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
52 for (int i=0; i<l; i++) A[i]=a[i];
53 NTT(A,n,1); NTT(b,n,1);
54 for (int i=0; i<n; i++) b[i]=1ll*b[i]*(2-1ll*A[i]*b[i]%mod+mod)%mod;
55 NTT(b,n,0);
56 for (int i=l; i<n; i++) b[i]=0;
57 for (int i=0; i<n; i++) A[i]=0;
58 }
59
60 void Sqrt(int a[],int b[],int l){
61 if (l==1){ b[0]=sqrt(a[0]); return; }
62 Sqrt(a,b,l>>1); Inv(b,B,l); int n=1,L=0;
63 for (; n<(l<<1); n<<=1) L++;
64 for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
65 for (int i=0; i<l; i++) C[i]=a[i];
66 NTT(C,n,1); NTT(B,n,1); NTT(b,n,1);
67 for (int i=0; i<n; i++) b[i]=1ll*inv2*(b[i]+1ll*C[i]*B[i]%mod)%mod;
68 NTT(b,n,0);
69 for (int i=l; i<n; i++) b[i]=0;
70 for (int i=0; i<n; i++) C[i]=B[i]=0;
71 }
72
73 void Deri(int a[],int b[],int l){
74 for (int i=1; i<l; i++) b[i-1]=1ll*i*a[i]%mod;
75 }
76
77 void Inte(int a[],int b[],int l){
78 for (int i=1; i<l; i++) b[i]=1ll*a[i-1]*inv[i]%mod; b[0]=0;
79 }
80
81 void Ln(int a[],int b[],int l){
82 Deri(a,D,l); Inv(a,E,l); mul(D,E,l); Inte(D,b,l);
83 for (int i=0; i<(l<<1); i++) D[i]=E[i]=0;
84 }
85
86 void Exp(int a[],int b[],int l){
87 if (l==1){ b[0]=1; return; }
88 Exp(a,b,l>>1); Ln(b,F,l); int n=1,L=0;
89 for (; n<(l<<1); n<<=1) L++;
90 for (int i=0; i<n; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
91 for (int i=0; i<l; i++) F[i]=(-F[i]+a[i]+mod)%mod; F[0]=(F[0]+1)%mod;
92 NTT(F,n,1); NTT(b,n,1);
93 for (int i=0; i<n; i++) b[i]=1ll*b[i]*F[i]%mod;
94 NTT(b,n,0);
95 for (int i=l; i<n; i++) b[i]=0;
96 for (int i=0; i<n; i++) F[i]=0;
97 }
98
99 void Ksm(int a[],int b[],int k,int l){
100 Ln(a,G,l);
101 for (int i=0; i<l; i++) G[i]=1ll*k*G[i]%mod;
102 Exp(G,b,l);
103 }
104
105 int main(){
106 freopen("polynomial.in","r",stdin);
107 freopen("polynomial.out","w",stdout);
108 scanf("%d%d",&n,&k);
109 for (int i=0; i<n; i++) scanf("%d",&X[i]);
110 int l=1; for (; l<=n; l<<=1); inv[0]=inv[1]=1;
111 rep(i,2,l) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
112 Sqrt(X,Y,l); mem(X); //Print(Y);
113 Inv(Y,X,l); mem(Y); //Print(X);
114 Inte(X,Y,l); mem(X); //Print(Y);
115 Exp(Y,X,l); mem(Y); //Print(X);
116 Inv(X,Y,l); Y[0]=(Y[0]+1)%mod; mem(X); //Print(Y);
117 Ln(Y,X,l); X[0]=(X[0]+1)%mod; mem(Y); //Print(X);
118 Ksm(X,Y,k,l); mem(X); //Print(Y);
119 Deri(Y,X,n); mem(Y); Print(X);
120 return 0;
121 }
