【bzoj2194】快速傅立叶之二 FFT

时间:2023-03-08 16:55:43

题意:给定序列a,b,求序列c,\(c(k)=\sum_{i=k}^{n-1}a(i)b(i-k)\)

Solution:

\[c(k)=\sum_{i=k}^{n-1}a(i)b(i-k)\\
c(k)=\sum_{i=0}^{n-k-1}a(i+k)b(i)\\
设ar(i)=a(n-i-1)\\
c(k)=\sum_{i=0}^{n-k-1}ar(n-i-k-1)b(i)\\
可以看出这是个卷积的形式,直接裸套fft
\]

Code:

#include<bits/stdc++.h>
#define ll long long
#define Pi acos(-1.0)
using namespace std;
const int N=262144;
int n,len=1,tim,rtt[N];
struct cp{double x,y;}aa[N],bb[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.y*b.x+a.x*b.y};}
void FFT(cp *a,int flag){
for(int i=0;i<len;i++)
if(i<rtt[i]) swap(a[i],a[rtt[i]]);
for(int l=2;l<=len;l<<=1){
cp wn=(cp){cos(flag*2*Pi/l),sin(flag*2*Pi/l)};
for(int st=0;st<len;st+=l){
cp w=(cp){1,0};
for(int u=st;u<st+(l>>1);u++,w=w*wn){
cp x=a[u],y=w*a[u+(l>>1)];
a[u]=x+y,a[u+(l>>1)]=x-y;
}
}
}
}
int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=-f;ch=getchar();}
while(isdigit(ch)){x=x*10+ch-48;ch=getchar();}
return x*f;
}
int main(){
n=read();
for(int i=0;i<n;i++){
scanf("%lf",&aa[n-i-1].x);
scanf("%lf",&bb[i].x);
}
while(len<=(n<<1)) len<<=1,++tim;
for(int i=0;i<len;i++)
rtt[i]=(rtt[i>>1]>>1)|((i&1)<<(tim-1));
FFT(aa,1);FFT(bb,1);
for(int i=0;i<len;i++) aa[i]=aa[i]*bb[i];
FFT(aa,-1);for(int i=0;i<len;i++) aa[i].x/=len;
for(int i=0;i<n;i++)
printf("%d\n",(int){aa[n-i-1].x+0.5});
return 0;
}