FFT学习笔记

FFT

FFT 常用于加速多项式乘法。

点值表示法

先考虑如何表示一个多项式。

最常见的是给定长度为 \(n+1\) 的系数序列 \(a\) 来表示多项式 \(F(x)=\sum\limits_{i=0}^na_ix^i\)，做多项式乘法时直接乘法分配律，时间复杂度是 \(O(n^2)\) 的。

另一种方法叫做点值表示法，也就是用平面上 \(n+1\) 个点（横坐标两两不同）来确定一个 \(n\) 次的多项式。

假设有两个多项式 \(F(x)\) 和 \(G(x)\)，最高次分别为 \(n\) 和 \(m\)。那么计算 \(H(x)=F(x)G(x)\) 时，任取 \(x_1,x_2,\cdots,x_{n+m+1}\)，求出所有 \(F(x_i)\) 和 \(G(x_i)\)，那么 \(H(x_i)=F(x_i)+G(x_i)\)。这样子时间复杂度是 \(O(n)\) 的。

而 FFT 的作用就是：将系数表示法在 \(O(n\log n)\) 转化成点值表示法。然后在将点值表示法用 IFFT 转化回系数表示法。

快速傅里叶变换

对于一个多项式 \(F(x)=\sum\limits_{i=0}^na_ix^i\)，目的是快速求出在一些 \(x\) 上的取值。

先考虑一种特殊情况：\(F(x)=x^2\) 时，显然求出 \(F(x)\) 时也知道 \(F(-x)\) 的值，因为它是偶函数；类似的 \(F(x)=x^3\) 时，知道 \(F(x)\) 也能知道 \(F(-x)\)，因为它是奇函数。

幸运的是 \(x\) 可以任意取，所以可以取正负成对的 \(x\)，减少了一半的计算量。

对于任意一个 \(F(x)\)，可以将其分解：

\[\begin{aligned}F(x)&=a_0+a_1x^1+a_2x^2+a_3x^3+\cdots+a_nx^n\\&=(a_0+a_2x^2+a_4x^4+\cdots)+x(a_1+a_3x^2+a_5x^4+\cdots)\end{aligned} \]

令

\[F_e(x)=a_0+a_2x+a_4x^2+\cdots \]

\[F_o(x)=a_1+a_3x+a_5x^2+\cdots \]

显然 \(F_e(x^2)\) 和 \(F_o(x^2)\) 都是偶函数，那么

\[F(x)=F_e(x^2)+xF_o(x^2) \]

\[F(-x)=F_e(x^2)-xF_o(x^2) \]

而 \(F_e\) 和 \(F_o\) 都是 \(\frac{n}{2}\) 次的，而求 \(F_e(x^2)\) 和 \(F_o(x^2)\) 是原问题的子问题，可以考虑递归解决。

这里有个问题：\(x\) 是正负成对的，但是 \(x^2\) 不一定是正负成对的。上述优化只有在 \(x\) 正负成对时有效，所以考虑构造一组 \(x\)，使得 \(x^2\) 也是正负成对的。

显然实数是绝对不行的，考虑在复数域上取 \(x\)。

单位根

令 \(n=2^k\) 以下默认 \(F(x)\) 是 \(n-1\) 次的多项式。如果不是，补上系数为 \(0\) 的高次项即可。

每个复数 \(a+bi\) 都在复平面上对应一个向量 \((a,b)\)，而两个复数相乘相当于模长相乘，幅角相加。

复数 \(a+bi\) 与 \(a-bi\) 共轭，记为 \(\overline{a+bi}=a-bi\)。共轭在复平面上的表现是关于实轴对称。

在复平面上作单位圆，然后将单位圆 \(n\) 等分，显然所有 \(n\) 等分点对应的复数都是 \(x^n=1\) 的解（模长乘积为 \(1\)，幅角相加为 \(2k\pi\)）。定义其中幅角为正且最小的复数为 \(\omega_n\)，即 \(n\) 次单位根。根据定义有

\[\omega_n=\cos\frac{2\pi}{n}+i\sin\frac{2\pi}{n} \]

\[\omega_n^k=\cos\frac{2k\pi}{n}+i\sin\frac{2k\pi}{n} \]

\[\omega_n^{k+\frac{n}{2}}=-\omega_n^k \]

因为 \(n=2^k\)，那么 \((a,b)\) 与 \((-a,-b)\) 同时在圆上，满足正负成对；同时 \((\omega_n^0)^2,(\omega_n^1)^2,\cdots,(\omega_n^{\frac{n}{2}-1})^2\) 是单位圆的 \(\frac{n}{2}\) 等分点，所以也满足正负成对。

回到 FFT，我们找到了平方之后也能正负成对的 \(x=\omega_n^k\)，要求所有 \(0\le k\le\frac{n}{2}\) 的 \(F(\omega_n^k)\)。那么：

\[\begin{aligned}F(\omega_n^k)&=F_e(\omega_n^{2k})+\omega_n^kF_o(\omega_n^{2k})\\&=F_e(\omega_\frac{n}{2}^k)+\omega_n^kF_o(\omega_\frac{n}{2}^k)\end{aligned} \]

\[\begin{aligned}F(-\omega_n^k)&=F(\omega_n^{k+\frac{n}{2}})\\&=F_e(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}F_o(\omega_n^{2k+n})\\&=F_e(\omega_\frac{n}{2}^k)-\omega_n^kF_o(\omega_\frac{n}{2}^k)\end{aligned} \]

然后 \(F_e(\omega_\frac{n}{2}^k)\) 和 \(F_o(\omega_\frac{n}{2}^k)\) 就真的是子问题了，递归解决即可。时间复杂度 \(O(n\log n)\)。

快速傅里叶逆变换

已经完成系数表示法到点值表示法的转换了，现在要转回去。

设系数为 \(a_0,a_1,\cdots,a_{n-1}\)，那么，

\[\begin{bmatrix}1&1&1&\cdots&1\\1&(\omega_n^1)^1&(\omega_n^1)^2&\cdots&(\omega_n^1)^{n-1}\\1&(\omega_n^2)^1&(\omega_n^2)^2&\cdots&(\omega_n^2)^{n-1}\\\vdots&\vdots&\vdots&\ddots&\vdots\\1&(\omega_n^{n-1})^1&(\omega_n^{n-1})^2&\cdots&(\omega_n^{n-1})^{n-1}\end{bmatrix}\begin{bmatrix}a_0\\a_1\\a_2\\\vdots\\a_{n-1}\end{bmatrix}=\begin{bmatrix}F(1)\\F(\omega_n^1)\\F(\omega_n^2)\\\vdots\\F(\omega_n^{n-1})\end{bmatrix} \]

而从 \(F(\omega_n^k)\) 转回 \(a_k\)，只需要将左边的矩阵求逆即可。求逆的结果为：

\[\operatorname{C}=\frac{1}{n}\begin{bmatrix}1&1&1&\cdots&1\\1&(\overline{\omega_n^1})^1&(\overline{\omega_n^1})^2&\cdots&(\overline{\omega_n^1})^{n-1}\\1&(\overline{\omega_n^2})^1&(\overline{\omega_n^2})^2&\cdots&(\overline{\omega_n^2})^{n-1}\\\vdots&\vdots&\vdots&\ddots&\vdots\\1&(\overline{\omega_n^{n-1}})^1&(\overline{\omega_n^{n-1}})^2&\cdots&(\overline{\omega_n^{n-1}})^{n-1}\\\end{bmatrix} \]

直接矩阵乘向量是 \(O(n^2)\) 的。但是共轭在复平面上表现为关于实轴对称，单位根的性质不变。所以可以套用 FFT 的过程，但是 \(\omega_n'=\cos \frac{2\pi}{n}-i\sin\frac{2\pi}{n}\)。

Code

#include <bits/stdc++.h>
const int M = 5e6 + 10;
const double pi = acos(-1);

struct Complex{
	double real, imag;
	Complex(double _r = 0, double _i = 0){
		real = _r, imag = _i;
	}
};
inline Complex operator + (Complex A, Complex B){
	return Complex(A.real + B.real, A.imag + B.imag);
}
inline Complex operator - (Complex A, Complex B){
	return Complex(A.real - B.real, A.imag - B.imag);
}
inline Complex operator * (Complex A, Complex B){
	return Complex(A.real * B.real - A.imag * B.imag, A.real * B.imag + B.real * A.imag);
}

int n, m, len;
Complex a[M], b[M];

inline void FFT(Complex* a, int len, int v) {
	if (len == 1) return;
	Complex e[len >> 1], o[len >> 1];
	for (int i = 0; i < len; i += 2)
		e[i >> 1] = a[i], o[i >> 1] = a[i + 1];
	FFT(o, len >> 1, v), FFT(e, len >> 1, v);
	Complex Wn(cos(2 * pi / len), v * sin(2 * pi / len)), w(1, 0);
	for (int i = 0; i < (len >> 1); i++, w = w * Wn)
		a[i] = e[i] + w * o[i], a[i + (len >> 1)] = e[i] - w * o[i];
}

int main() {
	scanf("%d%d", &n, &m);
	for (int i = 0; i <= n; i++) scanf("%lf", &a[i].real);
	for (int i = 0; i <= m; i++) scanf("%lf", &b[i].real);
	
	len = 1;
	while (len < n + m + 1) len <<= 1;
	FFT(a, len, 1), FFT(b, len, 1);
	for (int i = 0; i <= len; i++) a[i] = a[i] * b[i];
	FFT(a, len, -1);
	for (int i = 0; i <= n + m; i++)
		printf("%d ", (int)(a[i].real / len + 0.5));
	return 0;
}

蝴蝶变换

递归的开销太大了，考虑迭代实现 FFT。

假设一个八阶 FFT，在递归的过程中，考虑每个元素的位置：

第一层：\([a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7]\)；

第二层：\([a_0,a_2,a_4,a_6][a_1,a_3,a_5,a_7]\)；

第三层：\([a_0,a_4][a_2,a_6][a_1,a_5][a_3,a_7]\)；

第四层：\([a_0][a_4][a_2][a_6][a_1][a_5][a_3][a_7]\)。

数字	二进制
\(0\)	\(000\)
\(4\)	\(100\)
\(2\)	\(010\)
\(6\)	\(110\)
\(1\)	\(001\)
\(5\)	\(101\)
\(3\)	\(011\)
\(7\)	\(111\)

容易发现，如果将二进制表示倒过来，就是 \(0\sim7\)。所以可以预处理出每个元素最后的位置，然后从底层做到顶层，就可以规避递归。

Code

inline void FFT(Complex* a, int v) {
	for (int i = 0; i < n; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (n >> 1) : 0);
	for (int i = 0; i < n; i++)
		if (rev[i] > i) std::swap(a[i], a[rev[i]]);
	for (int k = 2, l = 1; k <= n; k <<= 1, l <<= 1) {
		Complex Wn = Complex(cos(2 * pii / k), v * sin(2 * pii / k));
		for (int i = 0; i < n; i += k) {
			Complex w = Complex(1, 0);
			for (int j = 0; j < l; j++, w = w * Wn) {
				Complex tmp = w * a[i + j + l];
				a[i + j + l] = a[i + j] - tmp;
				a[i + j] = a[i + j] + tmp;
			}
		}
	}
}