[BZOJ3527][ZJOI2014]力 FFT+数学

题目链接：http://www.lydsy.com/JudgeOnline/problem.php?id=3527

首先卷积的形式是$h(i)=sum_{i=0}^jf(i)g(i-j)$，如果我们可以把式子整理成这个样子再套上FFT就成功了。

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

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

令$f(i)=q_i$，$g(i)=frac{1}{i^2}$，$t(i)=q_{n-i}$，则有

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

因为$j$无法取到$i$或者$n-i$，所以令$g(0)=0$来消除最后一项的影响。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const double pi=acos(-1.0);
struct complex{
    double r,i;
    complex (double _r=0,double _i=0){
        r=_r,i=_i;
    }
    complex operator + (const complex &_){
        return complex(r+_.r,i+_.i);
    }
    complex operator - (const complex &_){
        return complex(r-_.r,i-_.i);
    }
    complex operator * (const complex &_){
        return complex(r*_.r-i*_.i,r*_.i+_.r*i);
    }
};
void reverse(complex y[],int len){
    for(int i=1,j=len>>1,k;i+1<len;i++){
        if(i<j) swap(y[i],y[j]);
        k=len>>1;
        while(j>=k){
            j-=k;
            k>>=1;
        }
        if(j<k) j+=k;
    }
}
void FFT(complex y[],int len,int on){
    reverse(y,len);
    for(int h=2;h<=len;h<<=1){
        complex wn(cos(-on*2*pi/h),sin(-on*2*pi/h));
        for(int j=0;j<len;j+=h){
            complex w(1,0);
            for(int k=j;k<j+(h>>1);k++){
                complex u=y[k],t=w*y[k+(h>>1)];
                y[k]=u+t;
                y[k+(h>>1)]=u-t;
                w=w*wn;
            }
        }
    }
    if(on==-1) for(int i=0;i<len;i++) y[i].r/=len;
}
double q[100010],ans[100010];
complex a[300010],b[300010],c[300010];
int main(){
    int n,len=1;
    scanf("%d",&n);
    while(len<n) len<<=1;len<<=1;
    for(int i=0;i<n;i++) scanf("%lf",&q[i]);
    for(int i=0;i<n;i++) a[i]=complex(q[i],0);
    for(int i=1;i<n;i++) b[i]=complex(1.0/i/i,0);
    for(int i=n;i<len;i++) a[i]=b[i]=complex();
    b[0]=complex();
    FFT(a,len,1);
    FFT(b,len,1);
    for(int i=0;i<len;i++) c[i]=a[i]*b[i];
    FFT(c,len,-1);
    for(int i=0;i<n;i++) ans[i]=c[i].r;
    for(int i=0;i<n;i++) a[i]=complex(q[n-i-1],0);
    for(int i=1;i<n;i++) b[i]=complex(1.0/i/i,0);
    for(int i=n;i<len;i++) a[i]=b[i]=complex();
    b[0]=complex();
    FFT(a,len,1);
    FFT(b,len,1);
    for(int i=0;i<len;i++) c[i]=a[i]*b[i];
    FFT(c,len,-1);
    for(int i=0;i<n;i++) ans[i]-=c[n-i-1].r;
    for(int i=0;i<n;i++) printf("%.3lf
",ans[i]);
    return 0;
}