快速傅里叶变换 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+42x2+x+32x^2+x+3 相乘,结果是 44 次的多项式,可以用 55 个点值表示出来。所以为 x2+2x+4x^2+2x+42x2+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+4odd(x2)=x+3odd(x^2)=x+3f(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+422 个点值,我们选择 x=±1x=\pm1。先算 x=1x=1even(12)=6even(1^2)=6odd(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)=61×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 的值。但继续计算 evenevenoddoddx=12,22x=1^2,2^2 的点值时,就不能继续分治了,因为此时的 1144 不再是正负成对了。

单位根

可以随意选择点值是幸运的,甚至可以不在实数范围内找,因为实数中没有满足平方前符号相同,平方后变成正负成对的一组数,但复数中存在

复数的一般形式为 a+bia+bi,其中 i2=1i^2=-1aabb 均为实数。当 b=0b=0 时,这个复数就是实数,当 a=0a=0 时,这个复数是纯虚数。复数支持加、减、乘、除运算,也支持实数内的运算法则,例如分配律、结合律。复数 a+bia+bi 的共轭数为 abia-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^\circkk 为任意整数,也就是 360360^\circ的倍数),因为 11 的向量幅角为 00^\circ 也就是 360360^\circ

下面举些例子。很明显,θ=0\theta=0^\circ 就是一个解,这与实际是相符的,因为 x=1x=1 确实是一个解。当 θ=180\theta=180^\circ 也就是 1-1 的向量的幅角时,只有当 nn 是偶数时才是方程的解,不然结果为 1-1,因为 nn 若是奇数那 n×θn\times\theta 就不是 360360^\circ 的倍数 了,就不能落在 11 的向量上了。这与实际也是相符的,例如 x4=1x^4=1x=1x=-1 的解,而 x3=1x^3=1 没有。

所以除 11 以外的幅角最小的那个解一定是以 360n\frac{360^\circ}{n} 为幅角, 11 为模长的向量,这个解也就是单位根 ωn1\omega_n^1。同时,单位根幅角的倍数也一定是方程的解,因为当 θ×n=360\theta \times n=360^\circ,就有 k×θ×n=k×360k\times \theta \times n=k\times360^\circ,也是 360360^\circ的倍数,也能在 nn 次幂后会落在向量 11 上。 所有以单位根幅角的倍数为幅角,模长为 11 的向量一共有 nn 个,那正是单位圆上的 nn 等分点为终点的向量。