luogu3338 [ZJOI2014]力

我发现我的构造方法好像不太一样而且比较显然？……先读入 (q) 数组（下表从零开始）。

记 (i < j) 时，(a_{i-j}=-1/i^2)；(i > j) 时，(a_{i-j}=1/i^2)；(i = j) 时，(a_{i-j}=0)。

答案 (E_i=sum_{j=0}^{n-1}a_{i-j}q_j)，可以用 FFT 优化，于是就做完了……吗？

发现 (a) 的下标可能会为负，那我们就整体平移一下，使得 (E_i=sum_{j=0}^{n-1}a_{i-j+n-1}q_j)，那么答案就是 (E) 数组的 (0+n-1 ldots n-1+n-1) 项了。（原先是 (0 ldots n-1) 项）

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
int n, lim=1, limcnt, rev[524305];
double q[100005];
const double PI=acos(-1.0);
struct Complex{
	double x, y;
	Complex(double u=0.0, double v=0.0){
		x = u; y = v;
	}
	Complex operator+(const Complex &u)const{
		return Complex(x+u.x, y+u.y);
	}
	Complex operator-(const Complex &u)const{
		return Complex(x-u.x, y-u.y);
	}
	Complex operator*(const Complex &u)const{
		return Complex(x*u.x-y*u.y, x*u.y+y*u.x);
	}
}a[524305], b[524305];
void fft(Complex a[], int opt){
	for(int i=0; i<lim; i++)
		if(i<rev[i])
			swap(a[i], a[rev[i]]);
	for(int i=2; i<=lim; i<<=1){
		int tmp=i>>1;
		Complex wn=Complex(cos(2*PI/i), opt*sin(2*PI/i));
		for(int j=0; j<lim; j+=i){
			Complex w=Complex(1.0, 0.0);
			for(int k=0; k<tmp; k++){
				Complex tmp1=a[j+k], tmp2=w*a[j+k+tmp];
				a[j+k] = tmp1 + tmp2;
				a[j+k+tmp] = tmp1 - tmp2;
				w = w * wn;
			}
		}
	}
	if(opt==-1)
		for(int i=0; i<lim; i++)
			a[i].x /= lim;
}
int main(){
	cin>>n;
	for(int i=0; i<n; i++)
		scanf("%lf", &b[i].x);
	for(int i=-n+1; i<=n-1; i++){
		if(i<0)
			a[i+n-1].x = -1.0 / i / i;
		else if(i==0)
			a[i+n-1].x = 0;
		else
			a[i+n-1].x = 1.0 / i / i;
	}
	while(lim<=3*(n-1))	lim <<= 1, limcnt++;
	for(int i=0; i<lim; i++)
		rev[i] = (rev[i>>1]>>1) | ((i&1)<<(limcnt-1));
	fft(a, 1);
	fft(b, 1);
	for(int i=0; i<lim; i++)
		a[i] = a[i] * b[i];
	fft(a, -1);
	for(int i=0; i<n; i++)
		printf("%.12f
", a[i+n-1].x);
	return 0;
}