营业日志 —— FFT-JZTXT

Q: 【模板】多项式乘法（FFT）

给定一个 \(n\) 次多项式 \(F(x)\)，和一个 \(m\) 次多项式 \(G(x)\)。

请求出 \(F(x)\) 和 \(G(x)\) 的卷积。

暴力很容易实现，但是时间复杂度为 \(O(n^2)\)

如何优化？

使用 FFT 可以有效将复杂度降至 \(O(n\log n)\)

先来一个引理：

对于 \(A(x)=a_0+a_1x+a_2x^2+...+a_nx^n\)，用函数上任意 \(n + 1\) 个互异的点均可唯一确定 \(A\).（多项式的点表示法）

\(proof.\) 对于这 \(n+1\) 个点 \((x_1, y_1), ..., (x_{n + 1}, y_{n + 1})\)，有：

\[\begin{cases} a_0 + a_1x_1 + a_2x_1^2 + ... + a_nx_1^n = y_1 \\ a_0 + a_1x_2 + a_2x_2^2 + ... + a_nx_2^n = y_2 \\ ... \\ a_0 + a_1x_n + a_2x_n^2 + ... + a_nx_n^n = y_n \end{cases} \]

方程组等价于行列式：

\[\left (\begin{array}{rrrr} 1 &x_1 & x_1^2 &... & x_1^n \\ 1 &x_2 & x_2^2 &... & x_2^n \\ & & ... \\ 1 &x_n & x_n^2 &... & x_n^n \end{array}\right) = \prod_{i, j}^{n}(x_i-x_j) \neq 0 \]

故一解！\(\square\)

注：

行列式为范德蒙行列式的转置，与原来等价
引理不仅在 \(\text{R}\) 上成立，也在 \(\text{C}\) 上成立

那么用点表示法有什么用呢？

对于 \(A(x), B(x), C(x) = A(x)B(x)\)

\(A(x), B(x)\) 取 \(x_1, x_2, ..., x_{n+m+1}\)

那么 \(C(x)\) 上取的点即为 \(\{(x_i, A(x_i)B(x_i))\}\)

时间复杂度 \(O(n)\).

取哪些点呢？

先来一手复数：\(z=a+bi\)

\[z_1+z_2=(a+bi)+(c+di)=(a+c)+(b+d)i \]

\[z_1\times z_2 = (a+bi)(c+di) = (ac-bc) + (ad+bc)i \]

取的点就是 \(n\) 次单位根 \(w_n^k\)：

单位根：\(w_n^n - 1 = 0\) 的 \(n\) 个根。

几个性质：

\(\forall i\neq j, w_n^i \neq w_n^j\)，即 \(n\) 个不同点
\(w_n^k = \cos \frac{2k\pi}{n}+i\sin \frac{2k\pi}{n}\)
\(w_n^0 = w_n^n = 1\)
\(w_{2n}^{2k} = w_n^k\)
\(w_n^{k + \frac{n}{2}} = -w_n^k\)

如何快速求 \(A(w_n^k), k = 0, 1, ..., n - 1\) 呢？

\[A(x) = a_0 + a_1x + ... + a_{n - 1}x^{n - 1} \]

令

\[\begin{cases} A_1(x) = a_0 + a_2x + a_4x^2 + ... + a_{n-2}x^{\frac{n}{2} - 1} \\ A_2(x) = a_1 + a_3x + a_5x^2 + ... + a_{n - 1}x^{\frac{n}{2} - 1} \end{cases} \]

则

\[A(x) + A_1(x^2) + xA_2(x^2) \]

将 \(x\) 代入为 \(w_n^k\)：

若 \(k\in [0, \frac{n}{2}-1], A(w_n^{k}) = A_1(w_n^{2k}) + w_n^kA_2(w_n^{2k}) = A_1(w_n^{k + \frac{n}{2}}) + w_n^kA_2(w_{\frac{n}{2}}^k)\)

若 \(k\in [\frac{n}{2}, n - 1], A(w_n^{k + \frac{n}{2}}) = A_1(w_n^{2k + n}) + w_n^{k + \frac{n}{2}}A_2(w_n^{2k + n}) = A_1(w_n^{2k}) + w_n^{k + \frac{n}{2}}A_2(w_n^{2k}) = A_1(w_n^{k + \frac{n}{2}}) - w_n^kA_2(w_{\frac{n}{2}}^k)\)

上下两式对比可得第一项相同，第二项互为相反数。

故可以分治求出。

这就是 FFT 的蝶形算法。

可发现 \(a_i\) 对应蝶形算法为其二进制下反转。

\(proof.\) 蝶形算法等价于每次奇偶反，即二进制下反转！\(\square\)

rev[i] = (r[i >> 1] >> 1) | (r[i] & 1) << (bit - 1)

#include <bits/stdc++.h>
using namespace std;

const int N = 3e6 + 10;
const double pi = acos(-1);

int n, m;
struct Complex {
	double x, y;
	Complex operator+ (const Complex& t) const {return {x + t.x, y + t.y};}
	Complex operator- (const Complex& t) const {return {x - t.x, y - t.y};}
	Complex operator* (const Complex& t) const {return {x * t.x - y * t.y, x * t.y + y * t.x};}
} a[N], b[N]; 
int rev[N], bit, tot;

void fft(Complex a[], int inv) {
	for (int i = 0; i < tot; i ++ )
		if (i < rev[i]) swap(a[i], a[rev[i]]);
	for (int mid = 1; mid < tot; mid <<= 1) {
		auto w1 = Complex({cos(pi / mid), inv * sin(pi / mid)});
		for (int i = 0; i < tot; i += mid * 2) {
			auto wk = Complex({1, 0});
			for (int j = 0; j < mid; j ++, wk = wk * w1) {
				auto x = a[i + j], y = wk * a[i + j + mid];
				a[i + j] = x + y, a[i + j + mid] = x - y;
			}
		}
	}
}

int main() {
	scanf("%d%d", &n, &m);
	for (int i = 0; i <= n; i ++ ) scanf("%lf", &a[i].x);
	for (int i = 0; i <= m; i ++ ) scanf("%lf", &b[i].x);
	while ((1 << bit) < n + m + 1) bit ++;
	tot = 1 << bit;
	for (int i = 0; i < tot; i ++) rev[i] = ((rev[i >> 1] >> 1) | ((i & 1) << (bit - 1)));
	fft(a, 1), fft(b, 1);
	for (int i = 0; i < tot; i ++ ) a[i] = a[i] * b[i];
	fft(a, -1);
	for (int i = 0; i <= n + m; i ++ ) printf("%d ", (int)(a[i].x / tot + 0.5));
	
	return 0;
}