裸的FFT,小心i*i爆int!!!
1 #include <iostream> 2 #include <cstring> 3 #include <cstdio> 4 #include <cmath> 5 using namespace std; 6 const int maxn=300010; 7 const double PI=acos(-1.0); 8 struct complex{ 9 double r,i; 10 complex(long double r_=0.0,long double i_=0.0){ 11 r=r_;i=i_; 12 } 13 complex operator+(complex a){ 14 return complex(r+a.r,i+a.i); 15 } 16 complex operator-(complex a){ 17 return complex(r-a.r,i-a.i); 18 } 19 complex operator*(complex a){ 20 return complex(r*a.r-i*a.i,i*a.r+r*a.i); 21 } 22 }A[maxn],B[maxn],C[maxn]; 23 24 void Rader(complex *a,int len){ 25 for(int i=1,j=len>>1;i<len-1;i++){ 26 if(i<j)swap(a[i],a[j]); 27 int k=len>>1; 28 while(j>=k){ 29 j-=k; 30 k>>=1; 31 } 32 j+=k; 33 } 34 } 35 36 void FFT(complex *a,int len,int on){ 37 Rader(a,len); 38 for(int h=2;h<=len;h<<=1){ 39 complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h)); 40 for(int j=0;j<len;j+=h){ 41 complex w(1,0); 42 for(int k=j;k<j+(h>>1);k++){ 43 complex x=a[k]; 44 complex y=a[k+(h>>1)]*w; 45 a[k]=x+y; 46 a[k+(h>>1)]=x-y; 47 w=w*wn; 48 } 49 } 50 } 51 if(on==-1) 52 for(int i=0;i<len;i++) 53 a[i].r/=1.0*len; 54 } 55 56 double f[maxn],E[maxn]; 57 int main(){ 58 int n,len=1; 59 scanf("%d",&n); 60 while(len<=n*2)len<<=1; 61 for(int i=1;i<=n;i++) 62 scanf("%lf",&f[i]); 63 64 for(int i=1;i<=n;i++){ 65 A[i-1].r=f[i]; 66 B[n-i].r=f[i]; 67 C[i].r=1.0/i/i; 68 } 69 FFT(A,len,1);FFT(B,len,1);FFT(C,len,1); 70 for(int i=0;i<len;i++) 71 A[i]=A[i]*C[i],B[i]=B[i]*C[i]; 72 FFT(A,len,-1);FFT(B,len,-1); 73 for(int i=0;i<n;i++){ 74 E[i]=A[i].r-B[n-i-1].r; 75 printf("%.3f ",E[i]); 76 } 77 return 0; 78 }