[BZOJ 3527] [ZJOI2014]力(FFT)

[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);
}