数学（FFT）：BZOJ 3527 [Zjoi2014]力

 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\n",E[i]);
76     }    
77     return 0;
78 }