[ZJOJ2014] 力 解题报告 (FFT)

题目链接：

https://www.luogu.org/problemnew/show/P3338

题目：

给出$n$个数$q_i$，令$F_j=sum_{i<j}frac{q_iq_j}{(i-j)^2}-sum_{i>j}frac{q_iq_j}{(i-j)^2}$，$E_i=frac{F_i}{q_i}$，求$E_i$

题解：

显然$E_j=sum_{i=1}^{j-1}frac{q_i}{(i-j)^2}-sum_{i=j+1}^{n}frac{q_i}{(i-j)^2}$

不妨设$G_i=frac{1}{i^2},G_0=0$,$T_i=q_i,T_0=0$

那么$E_j=sum_{i=1}^{j-1}T_iG_{i-j}-sum_{i=j+1}^{n}T_iG_{i-j}$

=$sum_{i=0}^{j}T_iG_{i-j}-sum_{i=j+1}^{n}T_iG_{i-j}$

显然减号前的部分是一个卷积的形式，是把$T​$和$G​$卷起来的第$i​$项

设$W_i=T_{n-i+1}(1<=i<=n),W_0=0$

减号后的$=sum_{i=j+1}^nW_{n-i+1}G_{i-j}$

$=sum_{i=1}^{n-j}W_iG_{n+1-i-j}$

=$sum_{i=0}^{n-j+1}W_iG_{n+1-i-j}$

这显然也是一个卷积的形式，是把$W$和$G$卷起来的第$n+1-j$项

$E_j=sum_{i=0}^{j}T_iG_{i-j}-sum_{i=0}^{n-j+1}W_iG_{n+1-i-j}$

 代码：

#include<algorithm>
#include<cstring>
#include<cstdio>
#include<iostream>
#include<cmath>
using namespace std;
typedef double db;

const int N=4e5+15;
const double pi=acos(-1.0);
struct complex
{
    double x,y;
    complex (double xx=0,double yy=0) {x=xx;y=yy;}
}A[N],B[N];
complex operator + (complex a,complex b) {return complex(a.x+b.x,a.y+b.y);} 
complex operator - (complex a,complex b) {return complex(a.x-b.x,a.y-b.y);} 
complex operator * (complex a,complex b) {return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);} 
int n;
int r[N<<1];
void fft(int limit,complex *a,int type)
{
    for (int i=0;i<limit;i++) if (i<r[i]) swap(a[i],a[r[i]]);
    for (int len=1;len<limit;len<<=1)
    {
        complex wn=complex(cos(pi/len),type*sin(pi/len));
        for (int k=0;k<limit;k+=(len<<1))
        {
            complex w=complex(1,0);
            for (int l=0;l<len;l++,w=w*wn)
            {
                complex Nx=a[k+l],Ny=w*a[k+l+len];
                a[k+l]=Nx+Ny;
                a[k+l+len]=Nx-Ny;
            }
        }
    }
}
void work(db *a,db *b,db *c)
{
    int limit=1,l=0;
    while (limit<n+n) limit<<=1,++l;
    for (int i=0;i<=limit;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    for (int i=0;i<=limit;i++) A[i].x=a[i],B[i].x=b[i],A[i].y=0,B[i].y=0;
    fft(limit,A,1);fft(limit,B,1);
    for (int i=0;i<=limit;i++) A[i]=A[i]*B[i];
    fft(limit,A,-1);
    for (int i=0;i<=limit;i++) c[i]=A[i].x/limit;
}
db q[N],p[N],f[N],a[N],b[N];
int main()
{    
    scanf("%d",&n);
    for (int i=1;i<=n;i++) 
    {
        scanf("%lf",&q[i]);
        p[i]=q[i];
        f[i]=(db)(1.0/i/i);
    }
    reverse(p+1,p+1+n);
    f[0]=0;q[0]=0;p[0]=0; 
    work(q,f,a);
    work(p,f,b);
    for (int i=1;i<=n;i++) printf("%lf
",a[i]-b[n-i+1]);
    return 0;
}