【BZOJ-3527】力 FFT

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
using namespace std;
struct Complex{
	double r,i;
	Complex(double R=0.0,double I=0.0) {r=R; i=I;}
	Complex operator + (const Complex & A) const {return Complex(r+A.r,i+A.i);}
	Complex operator - (const Complex & A) const {return Complex(r-A.r,i-A.i);}
	Complex operator * (const Complex & A) const {return Complex(r*A.r-i*A.i,r*A.i+i*A.r);}
};
#define MAXN 600010
#define Pai acos(-1.0)
Complex A[MAXN],B[MAXN],C[MAXN],D[MAXN];
int N,len;
double a[MAXN]; 
inline void Prework()
{
	len=1;
	while (len<((N-1)<<1)) len<<=1;
	for (int i=0; i<=N-1; i++) A[i]=Complex(a[i],0);
	for (int i=N; i<len; i++) A[i]=Complex(0,0);
	for (int i=0; i<=N-1; i++) B[i]=Complex(a[N-1-i],0);
	for (int i=N; i<len; i++) B[i]=Complex(0,0);
	for (int i=0; i<=N-1; i++) if (i) D[i]=C[i]=Complex(1.0/i/i,0);
	for (int i=N; i<len; i++) D[i]=C[i]=Complex(0,0);
}
inline void Rader(Complex *x)
{
    for (int i=1,j=len>>1,k; i<len-1; i++)
        {
            if (i<j) swap(x[i],x[j]);
            k=len>>1;
            while (j>=k) j-=k,k>>=1;
            if (j<k) j+=k;
        }
}
inline void DFT(Complex *x,int opt)
{
    Rader(x);
    for (int h=2; h<=len; h<<=1)
        {
            Complex Wn( cos(opt*2*Pai/h),sin(opt*2*Pai/h) );
            for (int i=0; i<len; i+=h)
                {
                    Complex W(1,0);
                    for (int j=i; j<i+h/2; j++)
                        {
                            Complex u=x[j],t=W*x[j+h/2];
                            x[j]=u+t; x[j+h/2]=u-t;
                            W=W*Wn;
                        }
                }
        }
    if (opt==-1)
        for (int i=0; i<len; i++) x[i].r/=len;
}
inline void FFT(Complex *x,Complex *y)
{
	DFT(x,1); DFT(y,1);
	for (int i=0; i<len; i++) x[i]=x[i]*y[i];
	DFT(x,-1);
}
int main()
{
	scanf("%d",&N);
	for (int i=0; i<N; i++) scanf("%lf",&a[i]);
	Prework();
//	for (int i=0; i<len; i++) printf("%.6lf\n",A[i].r); puts("=================");
//	for (int i=0; i<len; i++) printf("%.6lf\n",B[i].r); puts("=================");
//	for (int i=0; i<len; i++) printf("%.6lf\n",C[i].r); puts("=================");
	FFT(A,C); FFT(B,D);
	for (int i=0; i<N; i++) printf("%.3lf\n",A[i].r-B[N-1-i].r);
	return 0;
}