BZOJ3527: [Zjoi2014]力 [FFT]

#include <bits/stdc++.h>
#define pi acos(-1.0)
const int maxn = 300000+5;
using namespace std;
struct E{
    double real,imag;
    E(double real=0,double imag=0):real(real),imag(imag){}
    friend E operator +(E A,E B){
        return E(A.real+B.real,A.imag+B.imag);
    }
    friend E operator -(E A,E B){
        return E(A.real-B.real,A.imag-B.imag);
    }
    friend E operator *(E A,E B){
        return E(A.real*B.real-A.imag*B.imag,A.imag*B.real+A.real*B.imag);
    }
};
int n,m,L,rev[maxn];
E f[maxn],rf[maxn],g[maxn],e1[maxn],e2[maxn];
void fft(E *a,int ty){
    for(int i=0;i<n;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int i=1;i<n;i<<=1){
        E wn=E(cos(pi/i),ty*sin(pi/i));
        for(int p=(i<<1),j=0;j<n;j+=p){
            E w(1,0);
            for(int k=0;k<i;++k,w=w*wn){
                E x=a[j+k],y=w*a[j+k+i];
                a[j+k]=x+y;a[j+k+i]=x-y;
            }
        }
    }
    if(ty==-1)for(int i=0;i<n;++i)a[i].real/=n;
}
int main(){
    scanf("%d",&n);--n;
    for(int i=0;i<=n;++i){double x;
        scanf("%lf",&x);
        f[i] = x;
        rf[n-i] = x;
    }
    for(int i=1;i<=n;++i)g[i]=1.0/i/i;//g[0]=0
    m=n*2;for(n=1;n<=m;n<<=1)L++;
    for(int i=0;i<n;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    fft(f,1);fft(rf,1);fft(g,1);
    for(int i=0;i<n;++i)e1[i]=f[i]*g[i];
    for(int i=0;i<n;++i)e2[i]=rf[i]*g[i];
    fft(e1,-1);fft(e2,-1);
    for(int i=0;i<=m/2;++i)
        printf("%.3f
",e1[i].real-e2[m/2-i].real);
}