快速傅里叶变换（FFT）
一、简介
好像已经在 ABC 场上由于不会 FFT 被 G 题虐了至少两次了，全场就我不会 FFT。
所以我又来学 FFT 了。
$$\color{red}{ \Rightarrow 建议先学：FWT \Leftarrow }$$
了解大致思想即可。
给定一个 \(n\) 次多项式 \(F(x)\)，和一个 \(m\) 次多项式 \(G(x)\)。
请求出 \(F(x)\) 和 \(G(x)\) 的乘积。
不难看出这是一个高精度乘法板子，复杂度 \(O(n^2)\)。
FFT 基于分治，可以优化到 \(O(n \log n)\)。
二、前置知识
- 圆的弧度制
- 向量及其运算
- 复数及其运算
这些都比较基础。
\(n\) 次单位根
钦定 \(n\) 是二的正整数次幂。
定义 \(n\) 次单位根 \(\omega_n^i\) 为满足 \(z^n=1\) 的复数 \(z\)。
把单位元按角度平均分成 \(n\) 份，则 \(\omega_n^i\) 就是 \(i\) 个这样的份拼在一起，特别地有 \(\omega_n^0=1\)。
注意这些根都是向量，由此可以推导出单位根的性质：
- \(\omega_n^i = \omega_n^{i+n}\)
- \(\omega_n^i \omega_n^j = \omega_n^{i+j}\)
- \(\omega_{dn}^{di} = \omega_{n}^i\)
- \(\omega_n^i = -\omega_n^{i+\frac{n}{2}}\)
三、多项式的表示法
常见的表示方法都是系数表示法，给出多项式每一个次数对应的系数就可以唯一表示一个多项式。
用系数表示法直接做多项式乘法是没有任何前途的，因为它必须将两个多项式的每一位分别相乘，复杂度 \(O(n^2)\)。
考虑换一种表示法能不能快速完成乘法运算？
可以用点值表示法。
类似于表示一次函数 \(y=kx+b\)，还有另一种方法：两点确定一条直线，只要知道这两点 \((x_1,y_1), (x_2,y_2)\) 也可以表示出这个一次函数。
类似地，对于一个 \(n\) 次多项式 \(f\)，只需要 \(n+1\) 个点 \((x_0,f(x_0)), \dots, (x_i, f(x_i)), \dots, (x_n, f(x_n))\)。
对于多项式乘法 \(f \times g\)，就相当于每一个点的 \(y\) 值分别相乘：\((x_i, f(x_i) g(x_i) )\)。
这样做的复杂度是 \(O(n)\) 的！所以瓶颈转为：如何快速将一个多项式的 系数表示法 和 点值表示法 互相转化？
这个“转化”其实就是“变换”，也就是快速傅里叶变换（FFT）。
四、快速傅里叶正变换
大致框架
由于两个 \(n\) 次多项式乘法之后多项式有 \(2n+1\) 项，要再把它逆变换回系数表示法，就要求正变换选取至少 \(2n+1\) 个点值。
当然可以直接选取 \(2n+1\) 个 \(x\) 并计算它的 \(y\) 值，但这样的复杂度是 \(O(n^2)\) 的。
为什么这么计算低效？因为没有很好地运用多项式本身的性质。
考虑一个函数 \(y=x^2\)，需要选取 \(2\) 个点，最简单的做法就是暴力选两个点计算，但是还有更优的选择。
由于 \(y=x^2\) 是偶函数，只需要选一个 \(x \not= 0\) 计算，那么 \(-x\) 的函数值也就知道了。
对于 \(y=x^3\) 同理，它是奇函数，所以前面加个负号就行了。
扩展到一般多项式，上面的特例启发我们把多项式按照次数奇偶分组。
计算多项式 \(f(x) = 3x^5 + 2x^4 + x^3 + 7x^2 + 5x + 1\)，可以分组为 \(f(x) = \big( 2x^4 + 7x^2 + 1 \big) + x \big( 3x^4 + x^2 + 5 \big)\)，对于奇函数部分提出来一个 \(x\)。
设偶函数部分为 \(f_0(x)\)，奇函数部分为 \(f_1(x)\)，而且我们发现提了一个 \(x\) 出来之后所有次数都是偶数，所以这又是一个关于 \(x^2\) 的函数。
因此 \(f(x) = f_0(x^2) + x f_1(x^2)\)，那么 \(f(-x) = f_0(x^2) - x f_1(x^2)\)。
这件事很牛，因为只需要计算 \(f(x)\) 就可以知道 \(f(-x)\) 了，而且我们发现次数减半，\(f_0, f_1\) 是子问题，又可以递归求解了。
计算一下时间复杂度：每一层需要 \(O(n)\) 计算多项式的值，然后递归成两部分，每部分数量减半，最多递归 \(\log n\) 层，所以复杂度和归并排序一样是 \(O(n \log n)\)。
修正错误
但仔细观察会发现还存在一些问题。
在求 \(f\) 的过程中已经钦定了这样一件事：我们可以通过求 \(f(x)\)，推导出 \(f(-x)\)。
可由于上述算法中递归到下一层求的是 \(f'(x^2)\)，然而 \(x^2 = (-x)^2\)，所以这并不成立，也就是求出来的是相同 \(x\) 的值。
所以我们需要找一些数，使得它们平方之后仍然互为相反数。这就很复数了。
考虑 \(3\) 次多项式的正变换：我们需要选取 \(3+1=4\) 个点并使得它们组成两对相反数，假设选取的是 \(x_1,-x_1,x_2,-x_2\)。
把“递归”的过程反过来想成“合并”。
那么合并一层就会得到 \(x_1^2, x_2^2\)；合并两层得到 \(x_1^4\)。
由于 \(x\) 是可以任选的，不妨钦定 \(x_1=1\)，那么 \(-x_1=-1\)。
第二层有 \(x_1^2=1,x_2^2=-1\)，第三层 \(x_1^4=1\)。
由于 \(x_2^2=-1\)，可以推出 \(x_2=i,-x_2=-i\)，所以选取的四个数是 \(1,-1,i,-i\)。
整体上看，实际就是在解方程 \(x^4=1\)，求出它的四个解。
单位圆的妙用
由于递归每次要分成两半，不妨先把 \(n\) 补成 \(2\) 的幂次。
那么选取的这些 \(x_i,-x_i\) 就需要满足 \(x_i^n=1\)，即它们是 \(x^n=1\) 的 \(n\) 个解，是 \(1\) 的 \(n\) 个 \(n\) 次方根。
把这些点转移到复平面的单位圆上。则任意相邻两点和圆心连线的夹角是 \(\theta = \frac{2 \pi}{n}\)。
可以用欧拉公式表示出这些点：\(e^{i \theta} = \cos \theta + i \sin\theta\)。
定义 \(\omega_n = e^{\frac{2 \pi i}{n}}\)，则这些点就是 \(\omega_n^0, \omega_n^1, \dots\)。
回到之前单位根的性质：由于 \(\omega_n^i = -\omega_n^{i+\frac{n}{2}}\)，所以每个根都有对应的相反数。
而且在平方之后，就变成了 \(\frac{n}{2}\) 次单位根，依然是正负配对的。
解决了复数这一问题，回来修改表达式：
由于 \(x_i = \omega^i\)，所以有：
五、快速傅里叶逆变换
考虑把正变换刻画为矩阵、向量的乘积（\(p\) 表示多项式的系数）：
由于 \(x_i = \omega^i\)，可以把整个矩阵改写为：
这个矩阵被称为离散傅里叶变换（DFT）矩阵。
已知 DFT 矩阵和 \(f\) 的值，要求系数 \(p\) 相当于给 DFT 矩阵求逆再和 \(f\) 相乘。
此处省略一大堆代数推导
事实上矩阵的逆就长这样：
所以只需要把原来的 \(\omega\) 替换成 \(\omega^{-1}\) 即可！最后再乘上 \(\frac{1}{n}\)。
六、优化
蝴蝶变换
盗一张图。
这张图模拟了递归 FFT 中不断奇偶分组的过程。
通过并不强大的注意力，我们可以发现如下事实：
0 1 2 3 4 5 6 7
0 4 2 6 1 5 3 7
变换只不过是把上面数的二进制位反过来，就成了下面这些数字。
因此可以 \(O(n \log n)\) 先求出下面这行数字，记为 id 数组，然后迭代直接把下面的两两合并。
由于这样蝴蝶变换两次之后又反转回去了，所以不需要做额外的复原操作。
至于 \(O(n)\) 把原本的 \(f\) 数组按顺序排好，且不浪费多余的 \(O(n)\) 空间，可以这么做：
for (ri int i = 0; i < lim; i++)
if (i < id[i]) swap(f[i], f[id[i]]);
三步并两步
把 \(b\) 的元素放在 \(a\) 的虚部上，然后对 \(a\) 平方：
所以平方后的虚部就是答案的两倍。
这样只需要对 \(a\) 正变换一次，反过来 \(a^2\) 逆变换一次即可。
手写复数运算 + inline
一定要加 inline！！
七、代码
递归求解
#include <bits/stdc++.h>
using namespace std;
const int N = 1e6 + 5;
const double PI = acos(-1);
struct Complex {
double a, b;
Complex(double x = 0, double y = 0) { a = x, b = y; }
} ; // a + b * i
Complex operator + (Complex x, Complex y) { return (Complex){x.a + y.a, x.b + y.b}; }
Complex operator - (Complex x, Complex y) { return (Complex){x.a - y.a, x.b - y.b}; }
Complex operator * (Complex x, Complex y) { return (Complex){x.a * y.a - x.b * y.b, x.a * y.b + x.b * y.a}; }
int n, m, lim;
Complex a[N], b[N], c[N];
void FFT(int lim, Complex *f, int opt) {
if (lim == 1) return;
int mid = lim >> 1;
Complex f0[mid + 2], f1[mid + 2];
for (int i = 0; i < mid; i++) f0[i] = f[i << 1], f1[i] = f[i << 1 | 1];
FFT(mid, f0, opt), FFT(mid, f1, opt);
Complex omega(cos(2 * PI / lim), opt * sin(2 * PI / lim) ), now(1, 0);
for (int i = 0; i < mid; i++, now = now * omega) {
f[i] = f0[i] + (f1[i] * now);
f[i + mid] = f0[i] - (f1[i] * now);
}
}
int main() {
scanf("%d%d", &n, &m);
for (int i = 0; i <= n; i++) scanf("%lf", &a[i].a);
for (int i = 0; i <= m; i++) scanf("%lf", &b[i].a);
lim = 1; while (lim <= n + m) lim <<= 1;
FFT(lim, a, 1), FFT(lim, b, 1);
for (int i = 0; i < lim; i++) c[i] = a[i] * b[i];
// for (int i = 0; i < lim; i++) cout << c[i].a << ' ' << c[i].b << endl;
FFT(lim, c, -1);
for (int i = 0; i <= n + m; i++) printf("%.0lf ", fabs(c[i].a) / lim);
return 0;
}
蝴蝶变换
#include <bits/stdc++.h>
#define ri register
using namespace std;
const int N = 4e6 + 5;
const double PI = acos(-1);
struct Complex {
double a, b;
Complex(double x = 0, double y = 0) { a = x, b = y; }
} ; // a + b * i
inline Complex operator + (Complex x, Complex y) { return (Complex){x.a + y.a, x.b + y.b}; }
inline Complex operator - (Complex x, Complex y) { return (Complex){x.a - y.a, x.b - y.b}; }
inline Complex operator * (Complex x, Complex y) { return (Complex){x.a * y.a - x.b * y.b, x.a * y.b + x.b * y.a}; }
int n, m, lim;
Complex a[N], b[N], c[N];
int id[N], p[N];
inline int rev(int x, int lim) {
int y = 0; while (lim) y = (y << 1) | (x & 1), x >>= 1, lim >>= 1;
return y;
}
inline void init(int n) { for (ri int i = 0; i < n; i++) id[i] = rev(i, n - 1); }
void FFT(int lim, Complex *f, int opt) {
for (ri int i = 0; i < lim; i++)
if (i < id[i]) swap(f[i], f[id[i]]);
for (ri int len = 2; len <= lim; len <<= 1) {
ri int mid = len >> 1;
for (int i = 0; i < lim; i += len) {
ri Complex omega(cos(2 * PI / len), opt * sin(2 * PI / len)), now(1, 0);
for (ri int j = i; j < i + mid; j++, now = now * omega) {
Complex f0 = f[j], f1 = f[j + mid];
f[j] = f0 + (f1 * now);
f[j + mid] = f0 - (f1 * now);
}
}
}
}
int main() {
scanf("%d%d", &n, &m);
for (ri int i = 0; i <= n; i++) scanf("%lf", &a[i].a);
for (ri int i = 0; i <= m; i++) scanf("%lf", &b[i].a);
lim = 1; while (lim <= n + m) lim <<= 1;
init(lim);
for (ri int i = 0; i < lim; i++) a[i].b = b[i].a;
FFT(lim, a, 1);
for (ri int i = 0; i < lim; i++) c[i] = a[i] * a[i];
// for (ri int i = 0; i < lim; i++) cout << c[i].a << ' ' << c[i].b << endl;
FFT(lim, c, -1);
for (ri int i = 0; i <= n + m; i++) printf("%.0lf ", fabs(c[i].b) / lim / 2);
return 0;
}
浙公网安备 33010602011771号