2016-06-01 21:36:44
题目:http://www.lydsy.com/JudgeOnline/problem.php?id=3527
我就是一个大傻叉 微笑脸
1 #include<bits/stdc++.h> 2 #define inf 1000000000 3 #define ll long long 4 #define N 500005 5 using namespace std; 6 int read(){ 7 int x=0,f=1;char ch=getchar(); 8 while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();} 9 while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();} 10 return x*f; 11 } 12 const double Pi=acos(-1.0); 13 struct CD{ 14 double x,y; 15 CD(double a=0,double b=0){x=a,y=b;} 16 friend CD operator + (CD n1,CD n2){return CD(n1.x+n2.x,n1.y+n2.y);} 17 friend CD operator - (CD n1,CD n2){return CD(n1.x-n2.x,n1.y-n2.y);} 18 friend CD operator * (CD n1,CD n2){return CD(n1.x*n2.x-n1.y*n2.y,n1.x*n2.y+n1.y*n2.x);} 19 }; 20 CD a[N],b[N],c[N],d[N]; 21 int n,nn,bit; 22 double q[N]; 23 void FFT(CD *a,int n,int type){ 24 for(int i=0,j=0;i<n;i++){ 25 if(j>i)swap(a[i],a[j]); 26 int k=n; 27 while(j&(k>>=1))j&=~k; 28 j|=k; 29 } 30 for(int i=1;i<=bit;i++){ 31 CD w_n(cos(2*type*Pi/(1<<i)),sin(2*type*Pi/(1<<i))); 32 for(int j=0;j<(1<<bit);j+=(1<<i)){ 33 CD w(1,0); 34 for(int k=j;k<j+(1<<(i-1));k++){ 35 CD tmp=a[k],tt=w*a[k+(1<<(i-1))]; 36 a[k]=tmp+tt;a[k+(1<<(i-1))]=tmp-tt; 37 w=w*w_n; 38 } 39 } 40 } 41 if(type<0)for(int i=0;i<n;i++)a[i].x=a[i].x/n; 42 } 43 int main(){ 44 n=read();nn=n; 45 for(int i=0;i<n;i++)scanf("%lf",&q[i]),a[i]=CD(q[i],0); 46 for(int i=1;i<n;i++)b[i].x=1.0/(double)(i*i),b[i].y=0; 47 n=2*n-1;bit=1; 48 while((1<<bit)<n)bit++; 49 n=1<<bit; 50 b[0]=CD(); 51 for(int i=nn;i<n;i++)a[i]=CD(),b[i]=CD(); 52 53 FFT(a,n,1);FFT(b,n,1); 54 for(int i=0;i<n;i++)c[i]=a[i]*b[i]; 55 FFT(c,n,-1); 56 for(int i=0;i<nn;i++)a[i]=CD(q[nn-i-1],0); 57 for(int i=nn;i<n;i++)a[i]=CD(); 58 FFT(a,n,1); 59 for(int i=0;i<n;i++)d[i]=a[i]*b[i]; 60 FFT(d,n,-1); 61 62 for(int i=0;i<nn;i++)c[i].x-=d[nn-i-1].x; 63 for(int i=0;i<nn;i++)printf("%.5lf\n",c[i].x); 64 return 0; 65 } 66