[BZOJ 2194] 快速傅立叶之二 · FFT

来源：互联网发布：vr视频软件编辑：程序博客网时间：2024/04/30 13:29

求 $C(k)=\sum^{n-1}_{i=k}A_{i}B_{i-k}$

这个式子不好搞啊，然后我们就把B翻一下好了！

变成 $C(k)=\sum^{n-1}_{i=k}A_{i}B_{n-1-i+k}$

然后我们发现，i+(n-1-i+k)=n-1+k，只与k有关了！

然后我们另 $D(n-1+k)=\sum^{n-1}_{i=k}A_{i}B_{n-1-i+k}$ ，这个就很显然一个裸的卷积了，FFT模板直接上

所以最后对应的是C(0~n-1)对应D(n-1,n+n-2)。

然后是试了一下用<complex>写，代码简洁了很多吖。。。

#include <cstdio>#include <algorithm>#include <cstring>#include <cmath>#include <complex>using namespace std;#define f(i,x) for (int i=0;i<x;i++) const int N=4*1e5+10;const double pi=acos(-1);typedef complex<double> CP;int n,t,ans[N];CP a[N],b[N];void fft(CP a[],int flag){for (int i=1,j=0;i<t-1;i++){for (int s=t;(~j)&s;j^=(s>>=1));if (i<j) swap(a[i],a[j]);}for (int i=2;i<=t;i<<=1){CP wn(cos(2*pi/i),flag*sin(2*pi/i));for (int j=0;j<t;j+=i){CP w(1,0);for (int k=j;k<j+i/2;k++){CP x=a[k],y=a[k+i/2]*w;a[k]=x+y;a[k+i/2]=x-y;w*=wn;}}}if (flag==1) return ;f(i,t) a[i].real()/=t;}int main(){scanf("%d",&n);for(t=1;t<=n+n;t<<=1);f(i,n){int x,y;scanf("%d %d",&x,&y);a[i].real()=x;b[n-i-1].real()=y;}fft(a,1);fft(b,1);f(i,t) a[i]*=b[i];fft(a,-1);for (int i=n-1;i<n+n-1;i++) printf("%d\n",(int)(a[i].real()+0.5));return 0;}

0 0