BZOJ3527: [Zjoi2014]力

传送门

首先显然 $E[j]=sum_{i=1}^{j-1}frac{q[i]}{(i-j)^2}-sum_{i=j+1}^{n}frac{q[i]}{(i-j)^2}$

考虑怎么 $FFT$，设 $g[i]=frac{sgn(i)}{i^2}$

则 $E[j]=sum_{i=1}^{n}q[i]g[j-i]$，显然的 $FFT$ 式子

但是发现 $g$ 有负数下标，不好处理

所以设 $h[i]$ 使得 $g[i]=h[i+n-1]$，那么 $E[j]=sum_{i=1}^{n}q[i]h[j-i+n-1]$

为了使 $q,h$ 下标之和与左边下标相同，设 $S[n-1+j]=E[j]$，那么 $S[n-1+j]=sum_{i=1}^{n}q[i]h[j-i+n-1]$

然后就可以 $FFT$ 算 $S$，最后输出 $E[i]$ 即输出 $S[n-1+i]$

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;
typedef long long ll;
typedef long double ldb;
inline int read()
{
    int x=0,f=1; char ch=getchar();
    while(ch<'0'||ch>'9') { if(ch=='-') f=-1; ch=getchar(); }
    while(ch>='0'&&ch<='9') { x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
    return x*f;
}
const int N=6e5+7;
const ldb pi=acos(-1.0);
struct CP {
    ldb x,y;
    CP (ldb xx=0,ldb yy=0) { x=xx,y=yy; }
    inline CP operator + (const CP &tmp) const {
        return CP(x+tmp.x,y+tmp.y);
    }
    inline CP operator - (const CP &tmp) const {
        return CP(x-tmp.x,y-tmp.y);
    }
    inline CP operator * (const CP &tmp) const {
        return CP(x*tmp.x-y*tmp.y,x*tmp.y+y*tmp.x);
    }
}Q[N],H[N];
int n,p[N];
void FFT(CP *A,int len,int type)
{
    for(int i=0;i<len;i++) if(i<p[i]) swap(A[i],A[p[i]]);
    for(int mid=1;mid<len;mid<<=1)
    {
        CP wn(cos(pi/mid),type*sin(pi/mid));
        for(int R=mid<<1,j=0;j<len;j+=R)
        {
            CP w(1,0);
            for(int k=0;k<mid;k++,w=w*wn)
            {
                CP x=A[j+k],y=w*A[j+mid+k];
                A[j+k]=x+y;
                A[j+mid+k]=x-y;
            }
        }
    }
}
int main()
{
    n=read();
    for(int i=0;i<n;i++) scanf("%Lf",&Q[i].x);
    for(int i=-n+1;i<n;i++) H[n-1+i].x=i ? (i<0 ? -1.0/i/i : 1.0/i/i) : 0;
    int len=1,tot=0;
    while(len<=3*(n-1)) len<<=1,tot++;
    for(int i=0;i<len;i++) p[i]=(p[i>>1]>>1) | ((i&1)<<(tot-1));
    FFT(Q,len,1); FFT(H,len,1);
    for(int i=0;i<=len;i++) Q[i]=Q[i]*H[i];
    FFT(Q,len,-1);
    for(int i=0;i<n;i++) printf("%.12Lf
",Q[n-1+i].x/len);
    return 0;
}