ZJOI2014力-FFT快速傅里叶变换

Luogu-ZJOI2014力

说在前面

鸣谢$hsy$神仙的word文档。
虽然您说您被zxh虐，但是您还是把我调着打/kel…

Solution

直接上式子了。。。
$$E_j=\frac{F_j}{q_j}=\sum _{i=0}^{j-1}\frac{q_i}{(i-j{)}^2}-\sum _{i=j+1}^{n-1}\frac{q_i}{(i-j{)}^2}$$
设$$f[i]=\frac{sgn(i)}{i^2}(i!=0),f[i]=0(i=0)$$
所以$$E_j=\sum _{i=0}^{n-1}q[i]f[j-i]$$
这就是一个卷积了，所以赶紧$\mathcal{F}a\mathcal{F}a\mathcal{T}a$

Warning

FFT的空间在一般情况下要开4倍，但是本体因为$j-i$可能小于$0$所以右开大了一倍

Code

#include <cstdio>
#include <algorithm>
#include <cmath>
#define N 100010

typedef long double LD;
const LD pi = std :: acos(-1);
struct cp{
	LD a, b;
	cp(LD a_ = 0.0, LD b_ = 0.0) {
		a = a_; b = b_;
	}
	inline cp operator + (const cp &nxt) const {return cp(a + nxt.a, b + nxt.b); }
	inline cp operator - (const cp &nxt) const {return cp(a - nxt.a, b - nxt.b); }
	inline cp operator * (const cp &nxt) const {return cp(a * nxt.a - b * nxt.b, a * nxt.b + b * nxt.a); }
}q[N << 3], f[N << 3];

int re[N << 3], len, l;

inline void fft(cp *a, int fg) {
	for (int i = 0; i < len; ++i) {
		if (i < re[i]) std :: swap(a[i], a[re[i]]);
	}
	for (int l = 1; l < len; l <<= 1) {
		cp bs = cp(std :: cos(pi / l), std :: sin(pi / l) * fg);
		for (int j = 0; j < len; j += (l << 1)) {
			cp now = cp(1.0, 0.0);
			for (int k = 0; k < l; ++k) {
				cp nx = a[j + k], ny = a[j + k + l] * now;
				a[j + k] = nx + ny;
				a[j + k + l] = nx - ny;
				now = now * bs;
			}
		}
	}
}

int main() {
	int n;
	scanf("%d", &n);
	for (int i = 0; i < n; ++i) {
		scanf("%Lf", &q[i].a);
	}
	for (int i = - n + 1; i < n; ++i) {
		if (!i) f[n - 1].a = 0;
		else f[i + n - 1].a = 1.0 * abs(i) / i / i / i;
	}
	for (l = 0, len = 1; len <= 3 * n - 2; len <<= 1, ++l);
	for (int i = 0; i < len; ++i)
		re[i] = (re[i >> 1] >> 1 | ((i & 1) << (l - 1)));
	fft(q, 1);
	fft(f, 1);
	for (int i = 0; i < len; ++i) q[i] = q[i] * f[i];
	fft(q, -1);
	for (int i = 0; i < n; ++i) {
		printf("%.3Lf\n", q[i + n - 1].a / len);
	}
	return 0;
}