[ZJOI2014]力【FFT】

[ZJOI2014]力
已知$$E(j)=sum_{i < j} frac{q_i}{(i-j)^2}-sum_{i > j} frac{q_i}{(i-j)^2}$$， 求(E)
令(a(i)=q_i ,b(i)=i^{-2})
对于减号左边，$$L(j)=sum_{i < j}a(i)b(j-i)$$
对于减号右边，$$R(j)=-sum_{i > j}a(i)b(i-j)$$
令$$c(i)=b(i-n)=(i-n)^{-2}, i > n$$， 则$$L(j)=sum_{i < j}a(i)b(j-i+n-n)=sum_{i < j} a(i)c(j-i+n)=(ac)(n+j)$$
令$$c(i)=-b(n-i)=-(n-i)^{-2}, i < n$$，即$$c(n-i)=-b(i)=-i^{-2}$$,则$$R(j)=sum_{i > j}a(i)b(i-j)=sum_{i > j}a(i)c(n-i+j)=(ac)(n+j)$$
特殊的，令(c(n)=0)
则$$E(j)=(a*c)(n+j)$$
求出(a)和(c)的卷积就好了

const double pi=acos(-1.0);
const int Maxn=2e6+10;
int n, l, r[Maxn], lim=1;

struct complex{double a, b; complex(double a=0, double b=0):a(a), b(b){}} f[Maxn], g[Maxn];
complex operator +(complex &A, complex &B){return complex(A.a+B.a, A.b+B.b);}
complex operator -(complex &A, complex &B){return complex(A.a-B.a, A.b-B.b);}
complex operator *(complex &A, complex &B){return complex(A.a*B.a-A.b*B.b, A.b*B.a+A.a*B.b);}

void FFT(complex *A, int tp){
	for(int i=0; i < lim; i++) if(i < r[i]) swap(A[i], A[r[i]]);
	complex wn, w, x, y;
	for(int mid=1; mid < lim; mid<<=1){ wn=complex(cos(pi/mid), tp*sin(pi/mid));
		for(int j=0, R=mid<<1; j < lim; j+=R){ w=complex(1, 0);
			for(int k=0; k < mid; k++, w=w*wn) x=A[j+k], y=w*A[j+k+mid], A[j+k]=x+y, A[j+k+mid]=x-y;
		}
	}
}
void solve(){
	n=read();
	for(int i=1; i <= n; i++) scanf("%lf", &f[i].a); g[n].a=0;
	for(int i=0; i < n; i++) g[i].a=-1.0/(n-i)/(n-i); 
	for(int i=n+1; i <= n*2; i++) g[i].a=1.0/(i-n)/(i-n);
	while(lim <= n*3) lim <<= 1, l++; 
	for(int i=0; i < lim; i++) r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
	FFT(f, 1); FFT(g, 1); for(int i=0; i < lim; i++) f[i]=f[i]*g[i];
	FFT(f, -1); for(int i=n+1; i <= n+n; i++) printf("%lf
", f[i].a/lim);
}