求
[E_j=dfrac{F_j}{q_j}
]
原式变成
[dfrac{sum_{i=1}^{j-1}dfrac{q_i imes q_j}{(i-j)^2}-sum_{i=j+1}^ndfrac{q_i imes q_j}{(i-j)^2}}{q_j}
]
消去 (q_j)
[sum_{i=1}^{j-1}dfrac{q_i}{(i-j)^2}-sum_{i=j+1}^ndfrac{q_i}{(i-j)^2}
]
设 (g_i=dfrac{1}{i^2})
[sum_{i=1}^{j-1}q_i imes g_{j-i}-sum_{i=j+1}^n q_i imes g_{i-j}
]
令 (g_0=q_0=0)
[sum_{i=0}^{j}q_i imes g_{j-i}-sum_{i=j}^n q_i imes g_{i-j}
]
换成差值的方式表述
[sum_{i=0}^{j}q_i imes g_{j-i}-sum_{i=0}^{n-j}q_{j+i} imes g_{i}
]
令 (m=n-j)，(q'_i=q_{n-i})
[sum_{i=0}^{j}q_i imes g_{j-i}-sum_{i=0}^{m}q'_{m-i} imes g_{i}
]
现在均变成了卷积的形式，只要 FFT 后计算两个系数的差即可。
时间复杂度 (Oleft(nlog n ight))。
code:
#include<bits/stdc++.h>
using namespace std;
#define N 262145
#define Db double
#define For(i,x,y)for(i=x;i<=(y);i++)
#define Down(i,x,y)for(i=x;i>=(y);i--)
const Db PI=acos(-1.0);
struct Complex
{
Db x,y;
Complex(Db X=0.0,Db Y=0.0):x(X),y(Y){};
Complex operator+(Complex const&_)const
{
return Complex(x+_.x,y+_.y);
}
Complex operator-(Complex const&_)const
{
return Complex(x-_.x,y-_.y);
}
Complex operator*(Complex const&_)const
{
return Complex(x*_.x-y*_.y,x*_.y+y*_.x);
}
}q[N],qq[N],g[N];
int tr[N];
void FFT(Complex*f,int len,bool type)
{
int i,p,j,now;
Complex tmp,ome,buf;
For(i,0,len-1)
if(i<tr[i])swap(f[i],f[tr[i]]);
for(p=2;p<=len;p<<=1)
{
now=p>>1;
ome=Complex(cos(2.0*PI/Db(p)),sin(2.0*PI/Db(p))*(type?-1.0:1.0));
for(j=0;j<len;j+=p)
{
buf=Complex(1.0,0.0);
For(i,j,j+now-1)
{
tmp=buf*f[i+now];
f[i+now]=f[i]-tmp;
f[i]=f[i]+tmp;
buf=buf*ome;
}
}
}
if(type)
For(i,0,len-1)f[i].x=round(f[i].x/Db(len)*1000.0)/1000.0;
}
int main()
{
int n,m=1,i;
scanf("%d",&n);
For(i,1,n)
{
scanf("%lf",&q[i].x);
qq[n-i].x=q[i].x;
g[i].x=1.0/Db(i)/Db(i);
}
while(m<=n<<1)m<<=1;
For(i,0,m-1)tr[i]=tr[i>>1]>>1|(i&1?m>>1:0);
FFT(q,m,0),FFT(g,m,0),FFT(qq,m,0);
For(i,0,m-1)q[i]=q[i]*g[i],qq[i]=qq[i]*g[i];
FFT(q,m,1),FFT(qq,m,1);
For(i,1,n)printf("%.3lf
",q[i].x-qq[n-i].x);
return 0;
}
可以学到的技巧：
- 将求和符号的起始值搞成 (0)。
- 寻找和固定的下标。
- 找不到就尝试翻转下标，即将数组反过来。