题意
首先化简式子:
(E_j=frac{F_j}{q_j})
(=frac{sumlimits_{i=1}^{j-1}frac{q_i*q_j}{(i-j)^2}-sumlimits_{i=j+1}^{n}frac{q_i*q_j}{(i-j)^2}}{q_j})
(=sumlimits_{i=1}^{j-1}frac{q_i}{(i-j)^2}-sumlimits_{i=j+1}^{n}frac{q_i}{(i-j)^2})
令(f_i=q_i,g_i=frac{1}{i^2}):
(=sumlimits_{i=1}^{j-1}f_i*g_{j-i}-sumlimits_{i=j+1}^nf_i*g_{j-i})
显然前面那个式子就是个卷积的形式,可以直接(FFT),考虑后面那个式子:
(sumlimits_{i=j+1}^nf_i*g_{j-i})
令(f'(x))表示(f(x))翻转后的多项式:
(=sumlimits_{i=1}^{j-1}f'_i*g_{j-i})
于是分别对两式做(FFT),最后将第二个翻转后相减即可。
code:
#include<bits/stdc++.h>
using namespace std;
const int maxn=1e6+10;
const double Pi=acos(-1.0);
int n,lim=1,len;
int pos[maxn];
double q[maxn];
struct cplx{double x,y;}f[maxn],g[maxn],h[maxn],a[maxn],b[maxn];
cplx operator+(cplx a,cplx b){return (cplx){a.x+b.x,a.y+b.y};}
cplx operator-(cplx a,cplx b){return (cplx){a.x-b.x,a.y-b.y};}
cplx operator*(cplx a,cplx b){return (cplx){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
inline void FFT(cplx* a,int op)
{
for(int i=0;i<lim;i++)if(i<pos[i])swap(a[i],a[pos[i]]);
for(int l=1;l<lim;l<<=1)
{
cplx wn=(cplx){cos(Pi/l),op*sin(Pi/l)};
for(int i=0;i<lim;i+=l<<1)
{
cplx w=(cplx){1,0};
for(int j=0;j<l;j++,w=w*wn)
{
cplx x=a[i+j],y=w*a[i+l+j];
a[i+j]=x+y,a[i+l+j]=x-y;
}
}
}
if(op==1)return;
for(int i=0;i<lim;i++)a[i].x/=lim;
}
int main()
{
scanf("%d",&n);
for(int i=1;i<=n;i++)scanf("%lf",&q[i]);
for(int i=1;i<=n;i++)f[i].x=h[i].x=q[i],g[i].x=1.0/(1.0*i*i);
reverse(h+1,h+n+1);
while(lim<(n<<1))lim<<=1,len++;
for(int i=0;i<lim;i++)pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
FFT(f,1);FFT(g,1);FFT(h,1);
for(int i=0;i<lim;i++)a[i]=f[i]*g[i],b[i]=h[i]*g[i];
FFT(a,-1);FFT(b,-1);
reverse(b+1,b+n+1);
for(int i=1;i<=n;i++)printf("%.3lf
",a[i].x-b[i].x);
return 0;
}