思路
颓柿子的题目
要求求这样的一个式子
[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)
对于Ei，显然可以
[E_i=sum_{j=0}^{i-1}frac{q_j}{(i-j)^2}-sum_{j=i+1}^nfrac{q_j}{(i-j)^2}
]
前后没什么关联，可以分开考虑，首先考虑前面部分
[sum_{j=0}^{i-1}frac{q_j}{(i-j)^2}
]
设(f(i)=q_i)，(g(i)=frac{1}{i^2})
则前面一部分变为
[sum_{j=0}^{i-1}f(j)g(i-j)
]
变成了卷积的形式
后面一部分是
[sum_{j=i+1}^{n}frac{q_j}{(i-j)^2}
]
同理变成
[sum_{j=i+1}^nf(j)g(j-i)
]
似乎没什么办法，但是我们可以把(f)整个反过来（就是(f(1))与(f(n))交换，与(f(2))与(f(n-1))交换）
设交换之后的(f)为(f')，(f'_i=f_{n-i+1})
则有
[sum_{j=0}^{i-1}f'(j)g(j-i)
]
则原来的第i项和反转后的第(n-i+1)项相同
对应相减即可
代码
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#define int long long
using namespace std;
const double pi = acos(-1.0);
struct complex{//a+bi
double a,b;
complex operator + (complex &bx){
return (complex){a+bx.a,b+bx.b};
}
complex operator - (complex &bx){
return (complex){a-bx.a,b-bx.b};
}
complex operator * (complex &bx){
return (complex){a*bx.a-b*bx.b,b*bx.a+a*bx.b};
}
complex conj(void){
return (complex){a,-b};
}
};
int n;
double q[401000],ans[401000];
complex inv[401000],wnk[401000],a[401000],b[401000];
void init(int len){
for(int i=0;i<len;i++){
wnk[i]=(complex){cos(2*pi*i/len),sin(2*pi*i/len)};
inv[i]=wnk[i].conj();
}
}
void FFT(complex *a,complex *opt,int n){
int lim=0;
while((1<<lim)<n){
lim++;
}
n=(1<<lim);
for(int i=0;i<n;i++){
int t=0;
for(int j=0;j<lim;j++)
if((i>>j)&1)
t|=(1<<(lim-j-1));
if(t<i)
swap(a[i],a[t]);
}
for(int i=2;i<=n;i<<=1){
int len=i/2;
for(complex *j=a;j!=a+n;j+=i){
for(int k=0;k<len;k++){
complex t=j[k+len]*opt[n/i*k];
j[k+len]=j[k]-t;
j[k]=j[k]+t;
}
}
}
}
void Do_FFT(int n){
int lim=0;
while((1<<lim)<n){
lim++;
}
n=(1<<lim);
init(n);
FFT(a,wnk,n);
FFT(b,wnk,n);
for(int i=0;i<n;i++)
a[i]=a[i]*b[i];
FFT(a,inv,n);
for(int i=0;i<n;i++)
a[i].a/=n;
}
signed main(){
// freopen("1.in","r",stdin);
// freopen("test.out","w",stdout);
scanf("%lld",&n);
for(int i=1;i<=n;i++){
scanf("%lf",&q[i]);
}
for(int i=1;i<=n;i++){
a[i].a=q[i];
b[i].a=1.0/(i*i);
a[i].b=0;
b[i].b=0;
}
Do_FFT(2*n+4);
for(int i=1;i<=n;i++)
ans[i]=a[i].a;
memset(a,0,sizeof(a));
memset(b,0,sizeof(b));
for(int i=1;i<=n;i++){
a[i].a=q[i];
b[i].a=1.0/(i*i);
a[i].b=0;
b[i].b=0;
}
for(int i=1,j=n;i<j;i++,j--)
swap(a[i],a[j]);
Do_FFT(2*n+4);
for(int i=1;i<=n;i++)
ans[i]-=a[n-i+1].a;
for(int i=1;i<=n;i++)
printf("%.3lf
",ans[i]);
return 0;
}