[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}\)(题外话:其实就是电场强度).

\[\begin{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)=\begin{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\n",c[i].real()/m);
}