【jzoj3617】【ZJOI2014】【力】【fft】

题目大意

解题思路

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

我们设f[i]=1/i/i，$g_i=\sum_{j<i}q_j/(i-j)^2=\sum_{j<i}q_jf_{i-j}$

我们发现这个数一个卷积，我们把q和f当作两个从左到右是从低到高位的数做乘法（在最前面各加一个零位，消去当前位的影响），由于乘法的定义，乘后的第i位即为g[i]。对于后面的部分可以类似处理。

code

#include<cmath>
#include<cstdio>
#include<algorithm>
#define LF double
#define LL long long
#define min(a,b) ((a<b)?a:b)
#define max(a,b) ((a>b)?a:b)
#define fo(i,j,k) for(int i=j;i<=k;i++)
#define fd(i,j,k) for(int i=j;i>=k;i--)
using namespace std;
int const mxn=1e5;
int n,N,M;LF pi=acos(-1);
struct rec{
    LF x,y;
    rec(LF X=0,LF Y=0){x=X;y=Y;}
};
rec operator+(rec x,rec y){return rec(x.x+y.x,x.y+y.y);}
rec operator-(rec x,rec y){return rec(x.x-y.x,x.y-y.y);}
rec operator*(rec x,rec y){return rec(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);}
rec a[mxn*4+10],b[mxn*4+10],c[mxn*4+10],d[mxn*4+10],t[mxn*4+10];
void DFT(rec *a,int tag){
    fo(i,0,N-1){
        int pos=0;
        for(int ii=i,bit=1;bit<=M;pos=(pos<<1)+(ii&1),ii=ii>>1,bit++);
        t[pos]=a[i];
    }
    for(int size=2;size<=N;size=size<<1){
        int half=size>>1;
        fo(i,0,half-1){
            rec w=rec(cos(tag*pi*i/half),sin(tag*pi*i/half));
            for(int j=i;j<N;j+=size){
                rec u=t[j],v=w*t[j+half];
                t[j]=u+v;
                t[j+half]=u-v;
            }
        }
    }
    fo(i,0,N-1)a[i]=t[i];
}
void FFT(rec *a,rec *b){
    DFT(a,1);
    fo(i,0,N-1)a[i]=a[i]*b[i];
    DFT(a,-1);
    fo(i,0,N-1)a[i].x=a[i].x/N;
}
int up(LF x){
    return int(x)+(int(x)!=x);
}
int main(){
    freopen("d.in","r",stdin);
    freopen("d.out","w",stdout);
    scanf("%d",&n);
    M=up(log(n*2+2)/log(2)),N=pow(2,M);
    fo(i,1,n)scanf("%lf",&a[i].x),b[n-i+1].x=a[i].x;
    fo(i,1,n)c[i].x=1.0/i/i;
    DFT(c,1);
    FFT(a,c);
    FFT(b,c);
    fo(i,1,n)printf("%lf\n",a[i].x-b[n-i+1].x);
    return 0;
}