递归版
#include<bits/stdc++.h>
using namespace std;
#define N 4000005
int n,m;
const double Pi=acos(-1);
int tr[N];
struct CP{
double x,y;
}f[N],g[N];
CP operator + (CP a,CP b){return (CP){a.x+b.x,a.y+b.y};}
CP operator - (CP a,CP b){return (CP){a.x-b.x,a.y-b.y};}
CP operator * (CP a,CP b){return (CP){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
void FFT(CP *f,int len,bool flag){
if(len==1) return ;
FFT(f,len/2,flag);
FFT(f+len/2,len/2,flag);
CP tG=(CP){cos(2*Pi/len),sin(2*Pi/len)},buf=(CP){1,0};
if(!flag) tG.y*=-1;
for(int k=0;k<len/2;k++){
CP tt=buf*f[k+len/2];
f[k+len/2]=f[k]-tt;
f[k]=f[k]+tt;
buf=buf*tG;
}
}
int main(){
scanf("%d%d",&n,&m);
for(int i=0;i<N-3;i++) f[i].x=f[i].y=g[i].x=g[i].y=0;
for(int i=0;i<=n;i++) scanf("%lf",&f[i].x);
for(int i=0;i<=m;i++) scanf("%lf",&g[i].x);
for(m+=n,n=1;n<=m;n<<=1);
for(int i=0;i<n;i++) tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
for(int i=0;i<n;i++){
if(i<tr[i]) swap(f[i],f[tr[i]]),swap(g[i],g[tr[i]]);
}
FFT(f,n,1),FFT(g,n,1);
for(int i=0;i<n;i++) f[i]=f[i]*g[i];
for(int i=0;i<n;i++){
if(i<tr[i]) swap(f[i],f[tr[i]]);
}
FFT(f,n,0);
for(int i=0;i<=m;i++) printf("%d ",(int)(f[i].x/n+0.49));
return 0;
}
非递归
#include<bits/stdc++.h>
using namespace std;
#define N 4000005
int n,m;
const double Pi=acos(-1);
int tr[N];
struct CP{
double x,y;
}f[N],g[N];
CP operator + (CP a,CP b){return (CP){a.x+b.x,a.y+b.y};}
CP operator - (CP a,CP b){return (CP){a.x-b.x,a.y-b.y};}
CP operator * (CP a,CP b){return (CP){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
void FFT(CP *f,bool flag){
for(int len=2;len<=n;len<<=1){
CP tG=(CP){cos(2*Pi/len),sin(2*Pi/len)};
if(!flag) tG.y*=-1;
for(int qi=0;qi<n;qi+=len){
CP buf=(CP){1,0};
for(int k=qi;k<qi+len/2;k++){
CP tt=buf*f[k+len/2];
f[k+len/2]=f[k]-tt;
f[k]=f[k]+tt;
buf=buf*tG;
}
}
}
}
int main(){
scanf("%d%d",&n,&m);
for(int i=0;i<N-3;i++) f[i].x=f[i].y=g[i].x=g[i].y=0;
for(int i=0;i<=n;i++) scanf("%lf",&f[i].x);
for(int i=0;i<=m;i++) scanf("%lf",&g[i].x);
for(m+=n,n=1;n<=m;n<<=1);
for(int i=0;i<n;i++) tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
for(int i=0;i<n;i++){
if(i<tr[i]) swap(f[i],f[tr[i]]),swap(g[i],g[tr[i]]);
}
FFT(f,1),FFT(g,1);
for(int i=0;i<n;i++) f[i]=f[i]*g[i];
for(int i=0;i<n;i++){
if(i<tr[i]) swap(f[i],f[tr[i]]);
}
FFT(f,0);
for(int i=0;i<=m;i++) printf("%d ",(int)(f[i].x/n+0.49));
return 0;
}