[BZOJ3527][ZJOI2014]力（快速傅利叶变换FFT）

Description
给出n个数qi，给出Fj的定义如下：

令$E_i=F_i/q_i$，求Ei

Input
第一行一个整数n。
接下来n行每行输入一个数，第i行表示qi。
n≤100000，0 < qi < 1000000000

Output
n行，第i行输出Ei。与标准答案误差不超过1e-2即可。

Sample Input
5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880

Sample Output
-16838672.693
3439.793
7509018.566
4595686.886
10903040.872

一开始以为要求F吓哭了..

把式子化简可得：

$E_i=\sum_{j<i}{\frac{q_j}{(i-j)^2}}-\sum_{j>i}{\frac{q_j}{(i-j)^2}}$

这两个式子都好像卷积啊，第一个式子中$\sum_{j<i}$其实等效于$\sum_{j=0}^{i-1}$对于第二个式子把他倒过来就可以满足卷积啦！
那么之后就是这样:

$E_i=\sum_{j=0}^{i-1}{\frac{q_j}{(i-j)^2}}-\sum_{j=0}^{n-i-1}{\frac{q_{n-j}}{(n-i-j)^2}}$

把他化成卷积的形式

$E_i=\sum_{j=0}^{i-1}{{q_j}\frac{1}{(i-j)^2}}-\sum_{j=0}^{n-i-1}{{q_{n-j}}\frac{1}{(n-i-j)^2}}$

那么我们令

$f(i)=q_i，g(i)=\frac{1}{i^2}，f'(i)=q_{n-i}$

那么就可以把他们化成卷积的基本形式
首先

$E_i=\sum_{j=0}^{i-1}{f(j)g(i-j)}-\sum_{j=0}^{n-i-1}{f(n-j)g(n-i-j)}$

然后

$E_i=\sum_{j=0}^{i-1}{f(j)g(i-j)}-\sum_{j=0}^{n-i-1}{f'(j)g(n-i-j)}$

两次fft求卷积即可
（貌似有一次fft的做法？我去研究一下…..）

code：

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<iostream>
using namespace std;
const int maxn=100010;
const double pi=acos(-1.0);
struct complex
{
    double r,i;
    complex(){}
    complex(double _r,double _i){r=_r,i=_i;}
    friend complex operator + (const complex &x,const complex &y){return complex(x.r+y.r,x.i+y.i);}
    friend complex operator - (const complex &x,const complex &y){return complex(x.r-y.r,x.i-y.i);}
    friend complex operator * (const complex &x,const complex &y){return complex(x.r*y.r-x.i*y.i,x.r*y.i+x.i*y.r);}
}a[maxn*4],b[maxn*4];
double h1[maxn*4],h2[maxn*4],q[maxn*4];
int R[maxn*4];
void fft(complex *y,int len,int on)
{
    for(int i=0;i<len;i++)if(i<R[i])swap(y[i],y[R[i]]);
    for(int i=1;i<len;i<<=1) 
    {
        complex wn(cos(pi/i),sin(on*pi/i)); 
        for(int j=0;j<len;j+=(i<<1)) 
        {
            complex w(1,0);
            for(int k=0;k<i;k++,w=w*wn) 
            {
                complex u=y[j+k]; 
                complex v=w*y[j+k+i];
                y[j+k]=u+v;
                y[j+k+i]=u-v;
            }
        }
    }
    if(on==-1)for(int i=0;i<len;i++)y[i].r/=len; 
}
void calc(double *h,int n)
{
    int m=n*2; int L=0;
    for(n=1;n<=m;n*=2) L++;
    memset(R,0,sizeof(R));
    for(int i=0;i<n;i++) R[i]=(R[i>>1]>>1)|(i&1)<<(L-1);
    fft(a,n,1);
    fft(b,n,1);
    for(int i=0;i<=n;i++) a[i]=a[i]*b[i];
    fft(a,n,-1);
    for(int i=0;i<=n;i++) h[i]=(a[i].r+0.5);
}
int main()
{
    int n;scanf("%d",&n);
    for(int i=1;i<=n;i++) scanf("%lf",&q[i]);
    for(int i=1;i<=n;i++)
    {
        a[i].r=q[i];
        b[i].r=1/((double)i*(double)i);
    }
    calc(h1,n);
    memset(a,0,sizeof(a));
    memset(b,0,sizeof(b));
    a[0].r=q[n];
    for(int i=1;i<=n;i++)
    {
        a[i].r=q[n-i];
        b[i].r=1/((double)i*(double)i);
    }
    calc(h2,n);
    for(int i=1;i<=n;i++)
    {
        printf("%lf\n",h1[i]-h2[n-i]);
    }
    return 0;
}