Steven_Meng's Blog

传送门

题外话，那个式子巨神ypy说是电磁公式之类的，反正我觉得和电磁力什么很像。

$F_j = \sum_{i<j}\frac{q_i q_j}{(i-j)^2 }-\sum_{i>j}\frac{q_i q_j}{(i-j)^2 }$

考虑化简：

第一步：

$E_i = \sum_{j<i}\frac{q_j}{(i-j)^2 }-\sum_{j>i}\frac{q_j}{(i-j)^2 }$

谁都会。

令$f[i]=\frac{1}{i^2},g[i]=q_i$

有$E_i=\sum _{j<i} f[j]g[i-j]-\sum _{j<i} f[j]g[j-i]$

为什么第二个式子不能化成$f[j]g[i-j]$，第一眼看上去，这样显然不对，那还让你做，第二眼看上去，其实有一个重要的性质不能忘记：多项式下标不能为负数，$i-j<0$不能作为下标。

注意到第一个式子就是卷积的形式，但是第二个式子$j+(j-i)$不是一个常数，考虑翻转。

令$g'[i]=g[n-i]$，那么我们有$g'[n-j+i]=g[j-i]$，于是式子化成$E_i=\sum _{j<i} f[j]g[i-j]-\sum _{j<i} f[j]g'[n-j+i]$

但是，做两次FFT有一点麻烦，考虑将$g$数组开两倍，同时包含$g,g'$，即对于$i<n,g[i]=\frac{1}{(n-i)^2}$，对于$i>n,g[i]=\frac{1}{(i-n)^2}$

#include <bits/stdc++.h>
#define MAXN 1000005
using namespace std;
const double PI=acos(-1.0);
struct Complex{
	double x,y;
};
inline Complex operator + (const Complex &A,const Complex &B){
	return Complex{A.x+B.x,A.y+B.y};
}
inline Complex operator - (const Complex &A,const Complex &B){
	return Complex{A.x-B.x,A.y-B.y};
}
inline Complex operator * (const Complex &A,const Complex &B){
	return Complex{A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x};
}
int r[MAXN];
Complex a[MAXN],b[MAXN];
inline void FFT(Complex *A,int n,int type){
	for (register int i=1;i<=n;++i) if (i<r[i]) swap(A[i],A[r[i]]);
	for (register int i=1;i<n;i<<=1){
		int R=(i<<1);
		Complex Wn=Complex{cos(2*PI/R),type*sin(2*PI/R)};
		for (register int j=0;j<n;j+=R){
			Complex w=Complex{1,0};
			for (register int k=0;k<i;++k,w=w*Wn){
				Complex x=A[j+k],y=A[i+j+k]*w;
				A[j+k]=x+y,A[i+j+k]=x-y;
			}
		}
	}
}
int main(){
	int n;
	scanf("%d",&n);
	for (register int i=1;i<=n;++i){
		scanf("%lf",&a[i].x);
	}
	for (register int i=0;i<n;++i){
		b[i].x=-1.0/((double)(n-i)*(n-i));
	}
	for (register int i=n+1;i<=2*n;++i){
		b[i].x=1.0/((double)(i-n)*(i-n));
	}
	int sz=2*n,m=1,L=0;
	while (m<=sz) m<<=1,L++;
	for (register int i=0;i<=m;++i){
		r[i]=(r[i>>1]>>1|((i&1)<<(L-1)));
	}
	FFT(a,m,1),FFT(b,m,1);
	for (register int i=0;i<=m;++i) a[i]=a[i]*b[i];
	FFT(a,m,-1);
	for (register int i=n+1;i<=2*n;++i) printf("%.3lf\n",a[i].x/m);
}