ZJOI2014力-FFT快速傅里叶变换

最新推荐文章于 2021-08-06 11:35:04 发布

原创最新推荐文章于 2021-08-06 11:35:04 发布 · 340 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#FFT #ZJOI2014

本文详细解析了Luogu-ZJOI2014力题目的解决方案，通过引入卷积概念和快速傅立叶变换(FFT)算法，提供了一个高效求解策略。文章首先定义了特定函数f[i]，然后将原问题转化为卷积形式，最后使用FFT加速计算过程。

Luogu-ZJOI2014力

说在前面

鸣谢 $\mathcal{h{\color {red}sy}}$ 神仙的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{FaFaTa}$

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;
}