P3338 [ZJOI2014]力 题解

题面

看这个题目的式子，好像就是最近学的库仑定律+电场强度？？做法推式子FFT。

考虑首先通过 \(E_i=\frac{F_i}{q_i}\) 消掉原式中的一些量。得到：

\[E_i=\sum_{i=0}^j\frac{q_i}{(i-j)^2}-\sum_{i=j}^n\frac{q_i}{(i-j)^2} \]

考虑怎么把它化成卷积形式。令 \(f_i=q_i,g_i=\frac{1}{i^2}\)，可得

\[\begin{aligned} E_i&=\sum_{i=0}^{j}g[j-i]f[i]-\sum_{i=j}^{n}g[j-i]f[i]\\ &=\sum_{i=0}^{j}g[j-i]f[i]-\sum_{i=0}^{n-j}f[j+i]g[i] \end{aligned} \]

典型的 \(\min\) 卷积形式，一般做法是考虑翻转。令 \(f'[i]=f[n-i]\)，可得

\[\begin{aligned} E_i&=\sum_{i=0}^{j}g[j-i]f[i]-\sum_{i=0}^{n-j}f'[n-j-i]g[i]\\ &=\sum_{i=0}^{j}g[j-i]f[i]-\sum_{i=0}^{t}f'[t-i]g[i] \end{aligned} \]

这样就化成两个卷积的形式了，最后就把两个多项式卷一下然后用每一项系数求答案即可。\(ans_i=l_i-r_{n-i}\)。

点击查看代码

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int N=4e5+13;
const double Pi=acos(-1); 
struct complex{
	double x,y;
	complex(double xx=0,double yy=0){x=xx,y=yy;}
	complex operator +(const complex &a)const{return complex(x+a.x,y+a.y);}
	complex operator -(const complex &a)const{return complex(x-a.x,y-a.y);}
	complex operator *(const complex &a)const{return complex(x*a.x-y*a.y,x*a.y+y*a.x);}
}a[N],b[N];
int n,r[N];
double A[N],B[N],C[N],L[N],R[N];
inline void fft(complex *f,int limit,int type){
	for(int i=0;i<limit;++i)
		if(i<r[i]) swap(f[i],f[r[i]]);
	for(int mid=1;mid<limit;mid<<=1){
		complex Wn(cos(Pi/mid),type*sin(Pi/mid));
		for(int j=0;j<limit;j+=(mid<<1)){
			complex w(1,0);
			for(int k=0;k<mid;++k,w=w*Wn){
				complex x=f[j+k],y=w*f[j+k+mid];
				f[j+k]=x+y;f[j+k+mid]=x-y;
			}
		}
	}
}
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;++i) scanf("%lf",&A[i]);
	for(int i=1;i<=n;++i) B[i]=1.0/i/i;
	for(int i=0;i<=n;++i) C[i]=A[n-i];
	int limit=1,l=0;
	while(limit<=(n<<1)) limit<<=1,++l;
	for(int i=0;i<limit;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	for(int i=0;i<limit;++i) a[i].x=A[i],b[i].x=B[i];
	fft(a,limit,1),fft(b,limit,1);
	for(int i=0;i<limit;++i) a[i]=a[i]*b[i];
	fft(a,limit,-1);
	for(int i=0;i<=n;++i) L[i]=a[i].x/limit;
	for(int i=0;i<limit;++i) a[i].x=C[i],b[i].x=B[i],a[i].y=b[i].y=0;
	fft(a,limit,1),fft(b,limit,1);
	for(int i=0;i<limit;++i) a[i]=a[i]*b[i];
	fft(a,limit,-1);
	for(int i=0;i<=n;++i) R[i]=a[i].x/limit;
	for(int i=1;i<=n;++i) printf("%.3lf\n",L[i]-R[n-i]);
	return 0;
}