题意
给出 (n) 个数 (q_i),
定义
[F_j = sum_{i=1}^{j-1} frac{q_i imes q_j}{(i-j)^2} - sum_{i=j+1}^{n} frac{q_i imes q_j}{(i-j)^2}
]
[E_i = frac{F_i}{q_i}
]
求 (E_i) ((i in [1,n])).
思路
首先, 提公因式, 化简可以得到,
[ E_j = sum_{i=1}^{j-1} frac{q_i}{(j-i)^2} - sum_{i=j+1}^{n} frac{q_i}{(j-i)^2} $$ (把 $q_j$ 提出来后约掉)
注意到式子可以分为两个部分, 这两个部分中的因子可以分成 $i$ 和 $j-i$ 两个部分, 即这两个部分的编号和为 $j$, 那么我们就找到加法卷积了.
设 $F[i] = frac{1}{i^2}$, $G[i] = q_i$.
则
$$ E_j = sum_{i=1}^{j-1} F[j-i] imes G[i] - sum_{i=j+1}^{n} F[j-i] imes G[i] ]
这样, 我们就可以分两步愉快地 FFT 了.
需要注意的是, 后半部分的 (j-i) 为负数, 所以需要将整个数组往前移动 (n) 位, 最后统计答案的时候也要往前移动 (n) 位 (没听懂的话可以看一下代码).
代码
#include<bits/stdc++.h>
#define _USE_MATH_DEFINES
#define ll long long
#define db double
using namespace std;
const int _=1e5+7;
const int __=3e5+7;
const db Pi=M_PI;
struct cn{
db a,b;
cn operator + (const cn &x) const {
return (cn){a+x.a,b+x.b};
}
cn operator - (const cn &x) const {
return (cn){a-x.a,b-x.b};
}
cn operator * (const cn &x) const {
return (cn){a*x.a-b*x.b,a*x.b+b*x.a};
}
};
int N,n,num[__];
db q[_],E[_];
cn f[__],g[__];
void FFT(cn *f,int id){
for(int i=0;i<n;i++)
if(i<num[i]) swap(f[i],f[num[i]]);
for(int len=2;len<=n;len<<=1){
int gap=len>>1;
cn w1=(cn){cos(2*Pi/len),sin(2*Pi/len)*id};
for(int k=0;k<n;k+=len){
cn w=(cn){1,0};
for(int i=k;i<k+gap;i++,w=w*w1){
cn tmp=w*f[i+gap];
f[i+gap]=f[i]-tmp;
f[i]=f[i]+tmp;
}
}
}
}
int main(){
//freopen("force.in","r",stdin);
cin>>N;
for(int i=1;i<=N;i++) scanf("%lf",&q[i]);
n=1; while(n<=2*N) n<<=1;
for(ll i=1;i<=N;i++){ // 这里需要开 long long
f[i].a=(db)1ll/(i*i);
g[i].a=q[i];
}
for(int i=0;i<n;i++)
num[i]=(num[i>>1]>>1)|((i&1) ?n>>1 :0);
FFT(f,1);
FFT(g,1);
for(int i=0;i<n;i++) f[i]=f[i]*g[i];
FFT(f,-1);
for(int i=1;i<=N;i++) E[i]=f[i].a/n;
memset(f,0,sizeof(f));
memset(g,0,sizeof(g));
for(ll i=0;i<=N;i++){
f[i].a=(db)1ll/((i-N)*(i-N)); // 往前移动 N 位
g[i].a=q[i];
}
f[N].a=0;
FFT(f,1);
FFT(g,1);
for(int i=0;i<n;i++) f[i]=f[i]*g[i];
FFT(f,-1);
for(int i=1;i<=N;i++) printf("%.3lf
",E[i]-f[i+N].a/n); // 统计答案时也往前移动 N 位
return 0;
}