专栏文章

快速傅里叶变换

rretep
算法·理论参与者 47已保存评论 50

文章操作

快速查看文章及其快照的属性，并进行相关操作。

当前评论
50 条
当前快照
1 份
快照标识符
@mhz5sl5v
此快照首次捕获于
2025/11/15 01:55
4 个月前
此快照最后确认于
2025/11/29 05:24
3 个月前
查看原文
快速傅里叶变换 FFT，也就是 Fast Fourier Transform，是利用离散傅里叶变换（DFT）的高效、快速计算方法的算法统称。对于一个 OIer，FFT 最大的用处是用来加速多项式乘法。其本质是以 O(nlgn)O(nlgn) 的时间复杂度将点值表示法与系数表示法相互转换。
这篇文章没有成堆的空洞乏味的公式推导，以简洁易懂的数形结合思想解释 FFT 算法。当然，没有严格的证明也是缺点之一。

多项式的表示

多项式的表示主要有两种，第一种是最常见的系数表示法。将每个项的系数写在一起组成一个向量，缺少的项就补 00，顺序最好按照次数的大小从小到大，因为这样可以与数组的下标对应起来。
例如多项式 x2+2x+4x^2+2x+4 写成数组 [4,2,1][4,2,1]。
第二种是点值表示法，任意 kk 次的多项式都可以用 k+1k+1 个点唯一确定。这是很好理解的，kk 次多项式有 k+1k+1 个系数，也就是 k+1k+1 个未知数。若有 k+1k+1 个关于他们的方程，就能够全部解出来。
例如多项式 x2+2x+4x^2+2x+4 就可以用点值 (0,4),(1,6),(2,12)(0,4),(1,6),(2,12) 表示出来。
分别考虑两种表示法的多项式乘法。对于系数表示法，多项式乘法是基于乘法分配律的 O(n2)O(n^2) 算法，很难再做优化。而对于点值表示法，假设这两个多项式是 nn 次与 mm 次的，那么结果的多项式就需要 n+m+1n+m+1 个点值来表示。如果我们已经有了前两个多项式每个的 n+m+1n+m+1 个点值（要选一样的 xx），那么最后结果的点值只需要将这些点值的值一一乘起来就行了，复杂度是 O(n)O(n) 。
例如多项式 x2+2x+4x^2+2x+4 与 2x2+x+32x^2+x+3 相乘，结果是 44 次的多项式，可以用 55 个点值表示出来。所以为 x2+2x+4x^2+2x+4 与 2x2+x+32x^2+x+3 每个式子找到 55 个点值：
  • (0,4),(1,6),(2,12),(3,19),(4,28)(0,4),(1,6),(2,12),(3,19),(4,28)
  • (0,3),(1,6),(2,13),(3,24),(4,39)(0,3),(1,6),(2,13),(3,24),(4,39)
结果就是：(0,12),(1,36),(2,156),(3,456),(4,1092)(0,12),(1,36),(2,156),(3,456),(4,1092)
虽然点值法的乘法很快，但大部分情况题目给出的和我们需要的都是系数表示法的多项式，并且系数表示到点值表示的转换是 O(n2)O(n^2) 的复杂度，和分配律是一样的。我们需要一种能快速转换点值表示法和系数表示法的算法，而这正是 FFT 的用处。

计算点值

在讲 FFT 之前还需要思考一下其中的原理。速度不可能平白无故变快，一定是用什么方法减少了运算量。比如分治就是将本来需要计算两次的问题分解成两个只要计算一次的子问题，并且子问题是可以通过某种对称性互相转换的，这样每次递归就可以减少一半的运算量。FFT 也运用了分治思想。
既然是分治，就需要找到点值计算中的对称性。有没有哪些特别的点值我们只需要计算一个就能直接知道另一个？有，对于奇函数和偶函数，计算一个点值后，其相反数的值就可以瞬间得到。
但大部分多项式都没有奇偶性，所以我们要加工一下。分离多项式 f(x)f(x) 中的偶次项和奇次项，偶次项组成的多项式写作 even(x)even(x)，奇次项的写作 odd(x)odd(x)，提取公因式后能写成 even(x2)even(x^2) 和 x×odd(x2)x \times odd(x^2)。于是 f(x)=even(x2)+x×odd(x2)f(x)=even(x^2)+x \times odd(x^2)。
例如有多项式 f(x)=x3+2x2+3x+4f(x)=x^3+2x^2+3x+4，那么 even(x2)=2x+4even(x^2)=2x+4，odd(x2)=x+3odd(x^2)=x+3，f(x)=even(x2)+x×odd(x2)f(x)=even(x^2)+x \times odd(x^2)。
加工了以后，因为一个数的平方和其相反数的平方是一样的，所以f(−x)=even(x2)−x×odd(x2)f(-x)=even(x^2)-x \times odd(x^2)。如此，计算一个点值时，算出 even(x2)even(x^2) 和 odd(x2)odd(x^2) 的值就能得到答案了，其相反数的值也只需要将 odd(x2)odd(x^2) 取反就能得到。因为点值可以随意选择，所以我们一定会选择正负成对出现的点，因为这样我们就能用上面的公式减少一半的计算量，因为算出正的一半，负的一半也知道了。
例如计算多项式 f(x)=x3+2x2+3x+4f(x)=x^3+2x^2+3x+4 的 22 个点值，我们选择 x=±1x=\pm1。先算 x=1x=1，even(12)=6even(1^2)=6，odd(12)=4odd(1^2)=4，那么 f(1)=6+1×4=10f(1)=6+1\times4=10。再算 x=−1x=-1，这时我们不需要再算 even(x2)even(x^2) 与 odd(x2)odd(x^2) 了，只要变加为减，f(−1)=6−1×4=2f(-1)=6-1\times4=2。
不过这样也只是常数优化，没有任何用。细心的读者可能已经发现了，计算 even(x2)even(x^2) 和 odd(x2)odd(x^2) 也是求点值的问题，我们是不是可以用同样的方法优化呢？如果可以的话分治就成立了，计算点值的时间复杂度将变成 O(nlgn)O(nlgn)。可惜这个分治是不成立的，因为只有在第一次计算中我们可以随意选择点，选择那些正负成对的点，第二次计算的时候已经平方了，这时就不可能还是正负成对的了。
例如计算 x=1,2,−1,−2x=1,2,-1,-2 上的点值，这时是正负成对的，我们只需要计算 x=1,2x=1,2 就能知道 x=−1,−2x=-1,-2 的值。但继续计算 eveneven 和 oddodd 中 x=12,22x=1^2,2^2 的点值时，就不能继续分治了，因为此时的 11 和 44 不再是正负成对了。

单位根

可以随意选择点值是幸运的，甚至可以不在实数范围内找，因为实数中没有满足平方前符号相同，平方后变成正负成对的一组数，但复数中存在。
复数的一般形式为 a+bia+bi，其中 i2=−1i^2=-1，aa 和 bb 均为实数。当 b=0b=0 时，这个复数就是实数，当 a=0a=0 时，这个复数是纯虚数。复数支持加、减、乘、除运算，也支持实数内的运算法则，例如分配律、结合律。复数 a+bia+bi 的共轭数为 a−bia-bi，也就是含 ii 的项取反。
复平面对于复数就如数轴对于实数，不同的是这是个二维平面，平面上的任意坐标 (a,b)(a,b) 代表了复数 a+bia+bi。复平面实际上就是数轴多了个与其垂直的虚轴，是实数领域的拓展。
有了复平面以后，我们可以将任意复数表示成从原点指向其坐标的一个向量。向量的幅角为该向量与正实轴之间划过的角度。向量的模长是坐标到原点的距离。
这样，复数的加法和乘法就能用向量的运算表示出其几何意义了。加法就是向量中的平行四边形法则，头尾相接然后连线。乘法则是幅角相加，模长相乘。
如坐标系中一样，复平面上也有单位圆——以原点为圆心，11 为半径。
在圆上作 nn 个点将单位圆等分成 nn 份（以单位圆与实轴正半轴的交点为一个等分点）。然后以原点为起点，圆上的这 nn 个等分点为终点，作出 nn 个向量。其中幅角最小的那个的向量被称为 nn 次单位根，记做 ωn1\omega_n^1，其余的向量分别为 ωn2\omega_n^2，ωn3\omega_n^3，…\dots，ωnn\omega_n^n。
单位根的代数意义是 xn=1x^n=1 的解，一定会有 nn 个解，这 nn 个解正是 ωn1\omega_n^1，ωn2\omega_n^2，ωn3\omega_n^3，…\dots，ωnn\omega_n^n，其中实数解只有 11，有时有 −1-1，其余都是复数解。
为什么 xn=1x^n=1 的解会这么规律的排布在单位圆上呢？回忆我们刚才提过的复数乘法的几何意义，又因为单位圆上的向量模长都为 11，所以模长相乘后不会变，只有幅角相加了。因为幂运算来自于乘法所以也可以这样解释，一个单位圆上的向量的 nn 次方就是将其幅角乘 nn。于是求解 xn=1x^n=1 就转换成了找到角度 θ\theta，使 θ×n=k×360∘\theta \times n=k\times360^\circ（kk 为任意整数，也就是 360∘360^\circ的倍数），因为 11 的向量幅角为 0∘0^\circ 也就是 360∘360^\circ。
下面举些例子。很明显，θ=0∘\theta=0^\circ 就是一个解，这与实际是相符的，因为 x=1x=1 确实是一个解。当 θ=180∘\theta=180^\circ 也就是 −1-1 的向量的幅角时，只有当 nn 是偶数时才是方程的解，不然结果为 −1-1，因为 nn 若是奇数那 n×θn\times\theta 就不是 360∘360^\circ 的倍数 了，就不能落在 11 的向量上了。这与实际也是相符的，例如 x4=1x^4=1 有 x=−1x=-1 的解，而 x3=1x^3=1 没有。
所以除 11 以外的幅角最小的那个解一定是以 360∘n\frac{360^\circ}{n} 为幅角， 11 为模长的向量，这个解也就是单位根 ωn1\omega_n^1。同时，单位根幅角的倍数也一定是方程的解，因为当 θ×n=360∘\theta \times n=360^\circ，就有 k×θ×n=k×360∘k\times \theta \times n=k\times360^\circ，也是 360∘360^\circ的倍数，也能在 nn 次幂后会落在向量 11 上。 所有以单位根幅角的倍数为幅角，模长为 11 的向量一共有 nn 个，那正是单位圆上的 nn 等分点为终点的向量。
有了幅角和模长，怎么具体的把这个复数写出来呢？三角函数是个很有用的公式，高中知识告诉我们 sinθsin\theta 表示当幅角为 θ\theta， 模长为 11 时的 yy 坐标，cosθcos\theta 表示 xx 坐标。复平面中 xx 和 yy 分别表示 aa 和 bb 的值，所以我们就可以用三角函数将 nn 次单位根表示出来，ωn1=cos360∘n+i×sin360∘n\omega_n^1=cos\frac{360^\circ}{n}+i\times sin\frac{360^\circ}{n}，同理 ωnk=cos360∘×kn+i×sin360∘×kn\omega_n^k=cos\frac{360^\circ\times k}{n}+i\times sin\frac{360^\circ\times k}{n}，因为之前说过 kk 次方就是将幅角乘 kk。我们还需要把公式换成弧度制，因为 C++ 的内置函数是弧度制的。弧度制很简单，只需要记住 π=180∘\pi=180^\circ 即可。于是得到 ωnk=cos2kπn+i×sin2kπn\omega_n^k=cos\frac{2k\pi}{n}+i\times sin\frac{2k\pi}{n}，这个式子叫做欧拉公式，好像也可以用泰勒展开得到，让人感叹数学的魔幻。
上图中，OMOM 为 coscos 值，ONON 为 sinsin 值。
因为 FFT 中需要正负成对出现的值，所以我们还得看看单位根的正负关系。复数 a+bia+bi 的相反数是 −a−bi-a-bi，其几何意义是 xyxy 坐标同时取反，也就是以原点为中心，旋转 180∘180^\circ。用向量表示就是模长相同，方向相反。这很重要，因为我们发现单位圆 nn 等分，当 nn 为偶数时，就能得到 n2\frac{n}{2} 对这样互为相反数的点，而它们都可以用单位根表示出来。
如上图，ω81\omega_8^1，ω82\omega_8^2，ω83\omega_8^3 对应的相反数就是 ω85\omega_8^5，ω86\omega_8^6，ω87\omega_8^7。
所以我们可以得到当 nn 为偶数时，−ωnk=ωnk+n2-\omega^k_n=\omega^{k+\frac{n}{2}}_{n} （划重点）。

快速傅里叶变换

直接地说，快速傅里叶变换就是将单位根与其幂带入多项式求点值，并且利用单位根的性质让正负成对的性质一直保持下去，让分治一直成立，从而达到 O(nlgn)O(nlgn) 的时间复杂度。
我们已经知道了当 nn 为偶数时，−ωnk=ωnk+n2-\omega^k_n=\omega_{n}^{k+\frac{n}{2}}，所以带入 f(−x)=even(x2)−x×odd(x2)f(-x)=even(x^2)-x \times odd(x^2) 得到 f(ωnk+n2)=even(ωn2k)−ωnk×odd(ωn2k)f(\omega_{n}^{k+\frac{n}{2}})=even(\omega_n^{2k})-\omega_n^{k}\times odd(\omega_n^{2k})。翻译一下就是求出前 n2\frac{n}{2} 个根的 eveneven 和 oddodd 值后，后一半就直接知道了。这个分治是否能一直成立？会不会像实数一样递归一次后就不再是正负成对了？答案是不会，假如我们算的是正的一半（算负的一半也行）：
ωn1\omega_n^{1}，ωn2\omega_n^{2}，…\dots，ωnn2\omega_n^{\frac{n}{2}}
它们的幅角分别为：
2πn\frac{2\pi}{n}，2πn×2\frac{2\pi}{n}\times2，…\dots，2πn×n2\frac{2\pi}{n}\times \frac{n}{2}（弧度制，2π=360∘2\pi=360^\circ）
带入 eveneven 和 oddodd 求点值要先平方，平方相当于幅角翻倍：
2πn×2\frac{2\pi}{n}\times 2，2πn×4\frac{2\pi}{n}\times4，…\dots，2πn×2×n2\frac{2\pi}{n}\times 2\times\frac{n}{2}
我们惊奇地发现这 n2\frac{n}{2} 个向量平方以后重新将圆等分为了 n2\frac{n}{2} 份，当然前提是 nn 为偶数。如何理解？因为 nn 为偶数，所以这前 n2\frac{n}{2} 个向量本身平分了 xx 轴上方的半个单位圆。那么平方以后，幅角翻倍，就变成平分整个单位圆了。所以如果 n2\frac{n}{2} 还是偶数，就能继续分治下去，因为既然是平分单位圆，那就还是正负成对的。
因此，我们可以将多项式变成 22 的幂个项，多的项就当系数为 00，因为这样就能分治到底了，实现起来非常方便并且时间复杂度也是不变的。
例如多项式 x3+3x2+2x+4x^3+3x^2+2x+4 应补成 0x4+x3+3x2+2x+40x^4+x^3+3x^2+2x+4，因为 4=224=2^2。

快速傅里叶逆变换

我们已经讲完了如何使用单位根使分治成立，从而在 O(nlgn)O(nlgn) 内将系数表示法转换成点值表示法。但反过来怎么办呢？如何在 O(nlgn)O(nlgn) 内将点值表示法转换成系数表示法？
点值的计算可以写成矩阵乘法的形式，矩阵乘法很简单，不会可以点这里。
其中最左边的矩阵是带入的单位根，每个项括号外有裹了个次方是因为多项式每项变量的次数是不一样的。a0a_0 到 an−1a_{n-1} 是多项式的系数。y0y_0 到 yn−1y_{n-1} 是算出来的点值。
假如我们已经通过点值表示法的多项式乘法得到了结果多项式的点值，也就是说 y0y_0 到 yn−1y_{n-1} 是已知的。需要求结果多项式的系数，也就是说 a0a_0 到 an−1a_{n-1} 是未知的。这里要用一些线性代数的知识。设 AA，BB 和 CC 都是矩阵，且 A∗B=CA * B=C（矩阵乘法记做 ∗*），那么就有 A−1∗C=BA^{-1}*C=B，其中矩阵 A−1A^{-1} 为矩阵 AA 的逆矩阵。注意这里的 −1-1 不是次方的意思，只是一种标记方法，表示逆矩阵。
所以我们只需要求出最左边的矩阵的逆矩阵就能通过乘以点值矩阵得到系数矩阵了。而这个工作数学家们已经帮我们完成了，涉及亿点线性代数，我就不展开了，想看完整证明的同学可以点这里。
求其逆矩阵最后的结果是：
上加一线表示共轭，也就是复数含 ii 的项取反，前面的 1n\frac{1}{n} 表示矩阵内的每个项都除以 nn。总体来说，其逆矩阵就是每项取共轭，再除以 nn。
点值矩阵乘以这个逆矩阵得到的矩阵每项再除以 nn 就是系数矩阵，但这个工作也是 O(n2)O(n^2) 的，如何加速？细心的读者可能已经发现了，这个过程和傅里叶快速变换是一样的，完全可以用相同的分治来加速，因为共轭实质上就是按 xx 轴翻折，单位圆上单位根的性质不变。代码实现中两者只有一行语句的差别，也就是单位根含 ii 项的正负号。

代码实现

洛谷 P3803 FFT
CPP
#include <bits/stdc++.h>
#define N 5000000
#define err 1e-6
using namespace std;

int n, m;
const double pi = acos(-1.0);
complex<double> a[N], b[N]; 

void fft(complex<double> *a,int len,int inv){
	if(len==1)return;
	complex<double> e[len/2],o[len/2];
	for(int i=0;i<len;i+=2)
		e[i/2]=a[i],o[i/2]=a[i+1];
	fft(e,len/2,inv); fft(o,len/2,inv);
	complex<double> wn(cos(2*pi/len),sin(2*pi/len)*inv),wk(1,0);
	for(int i=0;i<len/2;i++,wk*=wn)
		a[i]=e[i]+wk*o[i],
		a[i+len/2]=e[i]-wk*o[i];
}

int main(){
	cin>>n>>m;
	for(int i=0,buf;i<=n;i++)cin>>buf,a[i].real(buf);
	for(int i=0,buf;i<=m;i++)cin>>buf,b[i].real(buf);
	int 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]*=b[i];
	fft(a,len,-1);
	for (int i=0;i<=n+m;++i)
		printf("%.0f ",a[i].real()/len);
	return 0;
}

参考资料

Fast Fourier Transform (2020), Youtube [online]. Available from: https://www.youtube.com/watch?v=h7apO7q16V0&feature=emb_logo [Accessed on June 1 2022]
一小时学会快速傅里叶变换 (2022), zhihu [online]. Available from: https://zhuanlan.zhihu.com/p/31584464 [Accessed on May 31 2022]
快速傅里叶变换 (2022), zhihu [online]. Available from: https://zhuanlan.zhihu.com/p/347091298 [Accessed on May 31 2022]
图片均来源于网络，侵权必改！

评论

共 50 条评论，欢迎与作者交流。

正在加载评论...