ZJOI2014力-FFT快速傅里叶变换

版权声明:蒟蒻的博文,dalao转载标明出处就好吖 https://blog.csdn.net/jokingcoder/article/details/89364106

Luogu-ZJOI2014力

说在前面

鸣谢 h s y \mathcal{h{\color {red}sy}} 神仙的word文档。
虽然您说您被zxh虐,但是您还是把我调着打/kel…

Solution

直接上式子了。。。
E j = F j q j = i = 0 j 1 q i ( i j ) 2 i = j + 1 n 1 q i ( i j ) 2 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 ] = s g n ( i ) i 2 ( i ! = 0 ) , f [ i ] = 0 ( i = 0 ) f[i]=\frac{sgn(i)}{i^2}(i!= 0),f[i]=0(i=0)
所以 E j = i = 0 n 1 q [ i ] f [ j i ] E_j=\sum_{i=0}^{n-1}q[i]f[j-i]
这就是一个卷积了,所以赶紧 F a F a T a \mathcal{FaFaTa}

Warning

FFT的空间在一般情况下要开4倍,但是本体因为 j i j-i 可能小于 0 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;
}

猜你喜欢

转载自blog.csdn.net/jokingcoder/article/details/89364106