[BZOJ 3527] [ZJOI2014]力(FFT)
题面
求$ F_j=sum_{i<j} frac{q_i q_j}{(i-j)^2}-sum_{i>j} frac{q_i q_j}{(i-j)^2} $
题外话:这题叫”力“，大概是因为这个式子长得很像电荷间作用力的公式(F=frac{kq_1q_2}{r^2})。
分析
设(E_i=frac{F_i}{q_i})(题外话:其实就是电场强度).
[egin{aligned}E_i &=sum_{j<i} frac{q_j}{(j-i)^2}-sum_{j>i} frac{q_j}{(j-i)^2} \ E_i&=sum_{j=1}^{i-1} frac{q_j}{(j-i)^2} - sum_{j=i+1}^{n} frac{q_j}{(j-i)^2} end{aligned}
]
令(f(i)=q_i, g(i)=frac{1}{i^2})
[E_i=sum_{j=1}^{i-1} f(j)g(i-j) - sum_{j=i+1}^{n} f(j)g(j-i)
]
注意到第一个式子就是卷积的形式，但是第二个式子(j+(j−i))不是一个常数，考虑变换形式。
令(g'(i)=g(n-i)),那么(g(j-i)=g'(n-j+i))
于是式子化成:
[E_i=sum_{j=1}^{i-1} f(j)g(i-j) - sum_{j=i+1}^{n} f(j)g'(n-j+i)
]
(j+(n-j+i)=n-i),为一常数，就可以直接卷积了
但是，做两次FFT有一点麻烦，考虑将g数组开两倍，同时包含(g,g’)，因为n+i>n,把g'放在n以后的部分。即
[g(i)=egin{cases} -frac{1}{(n-i)^2},i in[1,n] \ -frac{1}{(i-n)^2},i in[n+1,2n]end{cases}
]
代码
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<complex>
#define maxn 500000
using namespace std;
typedef complex<double> com;
const double pi=acos(-1.0);
int rev[maxn+5];
com a[maxn+5],b[maxn+5],c[maxn+5];
void fft(com *x,int n,int type){
for(int i=0;i<n;i++){
if(i<rev[i]) swap(x[i],x[rev[i]]);
}
for(int len=1;len<n;len*=2){
int sz=len*2;
com wn1=com(cos(2*pi/sz),type*sin(2*pi/sz));
for(int l=0;l<n;l+=sz){
int r=l+len-1;
com wnk=com(1,0);
for(int i=l;i<=r;i++){
com tmp=x[i+len];
x[i+len]=x[i]-wnk*tmp;
x[i]=x[i]+wnk*tmp;
wnk*=wn1;
}
}
}
}
int n;
int main(){
double x;
scanf("%d",&n);
for(int i=0;i<n;i++){
scanf("%lf",&x);
a[i]=x;
}
for(int i=0;i<n;i++) b[i]=-1.0/(1.0*(n-i)*(n-i));
for(int i=n+1;i<=n*2;i++) b[i]=1.0/(1.0*(i-n)*(i-n));
int k=0,m=1;
while(m<=n*2){
m*=2;
k++;
}
for(int i=0;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
fft(a,m,1);
fft(b,m,1);
for(int i=0;i<m;i++) c[i]=a[i]*b[i];
fft(c,m,-1);
for(int i=n;i<=n*2-1;i++) printf("%.3f
",c[i].real()/m);
}