FFT学习笔记

FFT

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

点值表示法#

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

最常见的是给定长度为 $n+1$ 的系数序列 $a$ 来表示多项式 $F(x)=\sum _{i=0}^n{a}_i{x}^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 _{i=0}^n{a}_i{x}^i$，目的是快速求出在一些 $x$ 上的取值。

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

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

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

$$\begin{array}{rl}F(x)& =a_0+a_1{x}^1+a_2{x}^2+a_3{x}^3+\cdots +a_n{x}^n\\ & =(a_0+a_2{x}^2+a_4{x}^4+\cdots )+x(a_1+a_3{x}^2+a_5{x}^4+\cdots )\end{array}$$

令

$$F_e(x)=a_0+a_{2}x+a_4{x}^2+\cdots$$

$$F_o(x)=a_1+a_{3}x+a_5{x}^2+\cdots$$

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

$$F(x)=F_e(x^2)+x{F}_o(x^2)$$

$$F(-x)=F_e(x^2)-x{F}_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{array}{rl}F({\omega }_n^k)& =F_e({\omega }_n^{2k})+{\omega }_n^k{F}_o({\omega }_n^{2k})\\ & =F_e({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k{F}_o({\omega }_{\frac{n}{2}}^k)\end{array}$$

$$\begin{array}{rl}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^k{F}_o({\omega }_{\frac{n}{2}}^k)\end{array}$$

然后 $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}$，那么，

$$\left[\begin{array}{ccccc}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}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& ({\omega }_n^{n-1}{)}^1& ({\omega }_n^{n-1}{)}^2& \cdots & ({\omega }_n^{n-1}{)}^{n-1}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ ⋮\\ a_{n-1}\end{array}\right]=\left[\begin{array}{c}F(1)\\ F({\omega }_n^1)\\ F({\omega }_n^2)\\ ⋮\\ F({\omega }_n^{n-1})\end{array}\right]$$

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

$$\mathrm{C}=\frac{1}{n}\left[\begin{array}{ccccc}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}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& (\overline{{\omega }_n^{n-1}}{)}^1& (\overline{{\omega }_n^{n-1}}{)}^2& \cdots & (\overline{{\omega }_n^{n-1}}{)}^{n-1}\end{array}\right]$$

直接矩阵乘向量是 $O(n^2)$ 的。但是共轭在复平面上表现为关于实轴对称，单位根的性质不变。所以可以套用 FFT 的过程，但是 ${\omega }_n^{\prime }=\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\sim 7$。所以可以预处理出每个元素最后的位置，然后从底层做到顶层，就可以规避递归。

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