好神仙啊
[F_{j}=sum_{i<j}frac{q_iq_j}{(i-j)^2}-sum_{j<i}frac{q_iq_j}{(i-j)^2}
]
求(frac{F_j}{q_j})
显然
[frac{F_j}{q_j}=sum_{i<j}frac{q_i}{(i-j)^2}-sum_{j<i}frac{q_i}{(i-j)^2}
]
先来看前面的那个柿子如何去搞
设(x=j-i)
那么
[sum_{i<j}frac{q_i}{(i-j)^2}=sum_{x=1}^{j-1}frac{q_{j-x}}{x^2}
]
我们搞出来两个多项式，(G(x)=frac{1}{x^2},H(x)=q_x)
那么就会发现
[G imes H(j)=sum_{i=1}^{j-1}G(i)H(j-i)=sum_{i=1}^{j-1}frac{q_{j-i}}{i^2}
]
哎这不就是了我们要求的东西了吗
我们发现还有后面那个东西，我们只需要把(q)反置再来一遍(FFT)就好了
代码
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#define maxn 500005
#define re register
#define LL long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
#define double long double
const double Pi=acos(-1);
inline int read()
{
char c=getchar();int x=0;while(c<'0'||c>'9') c=getchar();
while(c>='0'&&c<='9') x=(x<<3)+(x<<1)+c-48,c=getchar();
return x;
}
struct complex
{
double r,c;
complex (double a=0,double b=0) {r=a,c=b;}
}f[maxn],g[maxn],og,og1,t;
complex operator +(complex a,complex b) {return complex(a.r+b.r,a.c+b.c);}
complex operator -(complex a,complex b) {return complex(a.r-b.r,a.c-b.c);}
complex operator *(complex a,complex b) {return complex(a.r*b.r-a.c*b.c,a.r*b.c+a.c*b.r);}
int rev[maxn],len,n;
double q[maxn],ans[maxn];
inline void FFT(complex *f,int v)
{
for(re int i=0;i<=len;i++) if(i<rev[i]) std::swap(f[i],f[rev[i]]);
for(re int i=2;i<=len;i<<=1)
{
int ln=i>>1;
og1=complex(cos(Pi/ln),v*sin(Pi/ln));
for(re int l=0;l<len;l+=i)
{
og=complex(1,0);
for(re int x=l;x<l+ln;x++)
{
t=og*f[x+ln];
f[x+ln]=f[x]-t;
f[x]=f[x]+t;
og=og*og1;
}
}
}
}
int main()
{
scanf("%d",&n);
for(re int i=1;i<=n;i++) scanf("%Lf",&q[i]);
for(re int i=1;i<=n;i++) f[i].r=q[i],f[i].c=0;
for(re int i=1;i<=n;i++) g[i].r=1.0/i/i,g[i].c=0;
len=1;while(len<2*n+2) len<<=1;
for(re int i=0;i<=len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)?(len>>1):0);
FFT(f,1),FFT(g,1);
for(re int i=0;i<len;i++) f[i]=f[i]*g[i];
FFT(f,-1);
for(re int i=0;i<len;i++) ans[i]=f[i].r/len;
memset(f,0,sizeof(f)),memset(g,0,sizeof(g));
for(re int i=1;i<=n;i++) f[i].r=q[n-i+1],f[i].c=0;
for(re int i=1;i<=n;i++) g[i].r=1.0/i/i,g[i].c=0;
FFT(f,1),FFT(g,1);
for(re int i=0;i<len;i++) f[i]=f[i]*g[i];
FFT(f,-1);
for(re int i=1;i<=n;i++) printf("%.3Lf
",ans[i]-f[n-i+1].r/len);
return 0;
}