bzoj3527 -- FFT

代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
#define N 800010
#define PI acos(-1.0)
struct E{
    double r,i;
    E(){}
    E(double r,double i):r(r),i(i){}
    void operator /= (double a){r/=a;i/=a;}
    E operator + (E a){return E(r+a.r,i+a.i);}
    E operator - (E a){return E(r-a.r,i-a.i);}
    E operator * (E a){return E(r*a.r-i*a.i,r*a.i+i*a.r);}
}a[N],b[N],c[N];
int i,j,k,n,m,M,r[N],l;
double p[N],q[N];
inline void FFT(E* a,int n,int d){
    for(int i=0;i<n;i++)if(r[i]>i)swap(a[i],a[r[i]]);
    for(int i=1;i<n;i<<=1){
        E wn(cos(PI/i),d*sin(PI/i));
        for(int j=0;j<n;j+=i<<1){
            E w(1,0);
            for(int k=0;k<i;k++){
                E x=a[j+k],y=a[j+k+i]*w;
                a[j+k]=x+y;a[j+k+i]=x-y;
                w=w*wn;
            }
        }
    }
    if(d==-1)for(int i=0;i<n;i++)a[i]/=n;
}
int main()
{
    scanf("%d",&n);
    for(i=1;i<=n;i++)scanf("%lf",&q[i]),b[i+n].r=q[i];
    for(i=0;i<=n<<1;i++)
    if(i<n)p[i]=-1.0/(1ll*(n-i)*(n-i));else
    if(i>n)p[i]=1.0/(1ll*(n-i)*(n-i));else p[i]=0;
    n<<=1;
    for(i=0;i<=n;i++)a[i].r=p[i];
    for(M=1;M<=n<<1;M<<=1)l++;
    for(i=0;i<M;i++)r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
    FFT(a,M,1);FFT(b,M,1);
    for(i=0;i<M;i++)a[i]=a[i]*b[i];
    FFT(a,M,-1);
    for(i=n+1;i<=n+(n>>1);i++)printf("%.3lf
",a[i].r);
    return 0;
}