【参考】「ZJOI2014」力 - FFT by menci
【算法】FFT处理卷积
【题解】将式子代入后，化为Ej=Aj-Bj。
Aj=Σqi*[1/(i-j)^2]，i=1~j-1。
令f(i)=qi，g(i)=1/i^2，定义f(0)=g(0)=0（方便卷积）。
Aj=Σf(i)*g(j-i)，i=0~j-1，标准的卷积形式。
而对于Bj，将g反转后就是和为i+n-1的标准卷积形式了。
第一次FFT后，记得对a数组后半部分清零后再进行第二次FFT。
复杂度O(n log n)。
#include<cstdio> #include<algorithm> #include<complex> #include<cmath> using namespace std; const int maxn=300010; const double PI=acos(-1); complex<double>a1[maxn],a2[maxn]; int n; double ans[maxn],b1[maxn],b2[maxn]; namespace fft{ complex<double>o[maxn],oi[maxn]; void init(int n){ for(int k=0;k<n;k++)o[k]=complex<double>(cos(2*PI*k/n),sin(2*PI*k/n)),oi[k]=conj(o[k]); } void transform(complex<double>*a,int n,complex<double>*o){ int k=0; while((1<<k)<n)k++; for(int i=0;i<n;i++){ int t=0; for(int j=0;j<k;j++)if(i&(1<<j))t|=(1<<(k-j-1)); if(i<t)swap(a[i],a[t]); } for(int l=2;l<=n;l*=2){ int m=l/2; for(complex<double>*p=a;p!=a+n;p+=l){ for(int i=0;i<m;i++){ complex<double>t=o[n/l*i]*p[i+m]; p[i+m]=p[i]-t; p[i]+=t; } } } } void dft(complex<double>*a,int n){transform(a,n,o);} void idft(complex<double>*a,int n){transform(a,n,oi);for(int i=0;i<n;i++)a[i]/=n;} } void multply(int n){ int N=1; while(N<n+n)N*=2; for(int i=0;i<N;i++)a1[i]=a2[i]=0; for(int i=0;i<n;i++)a1[i].real(b1[i]),a2[i].real(b2[i]); fft::init(N);fft::dft(a1,N);fft::dft(a2,N); for(int i=0;i<N;i++)a1[i]*=a2[i]; fft::idft(a1,N); } int main(){ scanf("%d",&n);n++; b1[0]=b2[0]=0; for(int i=1;i<n;i++){ scanf("%lf",&b1[i]); b2[i]=1.0/i/i; } multply(n); for(int i=1;i<n;i++)ans[i]=a1[i].real(); for(int i=0;i<n/2;i++)swap(b2[i],b2[n-i-1]); multply(n); for(int i=1;i<n;i++){ ans[i]-=a1[i+n-1].real(); printf("%.3lf\n",ans[i]); } return 0; }
