超详细易懂FFT（快速傅里叶变换）及代码实现

昨天学了一晚上，终于搞懂了FFT。希望能写一篇清楚易懂的题解分享给大家，也进一步加深自己的理解。 FFT算是数论中比较重要的东西，听起来就很高深的亚子。但其实学会了（哪怕并不能完全理解)，会实现代码，并知道怎么灵活运用 （背板子） 就行。接下来进入正题。

FFT（Fast Fourier Transformation），中文名快速傅里叶变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。 而在信奥中，一般用来加速多项式乘法。 朴素高精度乘法的时间为 O ( n 2 ) O(n^2) O(n2)，但FFT能将时间复杂度降到 O ( n l o g 2 n ) O(nlog_2n) O(nlog2​n) 学习FFT之前，需要了解一些有关复数和多项式的知识。

F [ x ] = y = a 0 x 0 + a 1 x 1 + a 2 x 2 + . . . . . . a n x n F[x]=y=a_0x^0+a_1x^1+a_2x^2+......a_nx^n F[x]=y=a0​x0+a1​x1+a2​x2+......an​xn { a 0 , a 1 , a 2 , . . . , a n a_0,a_1,a_2,...,a_n a0​,a1​,a2​,...,an​} 是这个多项式每一项的系数，所以这是多项式的系数表示法

在函数图像中， F [ x ] F[x] F[x]这个多项式可以被n个点唯一确定，即代入n个点作为 x x x,分别解出对应的 y y y，得到n个式子。把这n条式子联立起来成为一个有n条方程的n元方程组，每一项的系数都可以解出来.（可类比二元一次方程） 也就是说，使用{ ( x 0 , f [ x 0 ] ) (x_0,f[x_0]) (x0​,f[x0​]), ( x 1 , f [ x 1 ] ) (x_1,f[x_1]) (x1​,f[x1​]),…, ( x n , f [ x n ] ) (x_n,f[x_n]) (xn​,f[xn​])}就可以完整描述出这个多项式,这就是 多项式的点值表示法

设两个多项式分别为 f ( x ) f(x) f(x), g ( x ) g(x) g(x),我们要把这两个多项式相乘 (即求卷积)。 如果用系数表示法: 我们要枚举 f f f的每一位的系数与 g g g的每一位的系数相乘，多项式乘法时间复杂度 O ( n 2 ) O(n^2) O(n2)，这也是我们所熟知的高精度乘法的原理。 如果用点值表示法: f [ x ] f[x] f[x]={ ( x 0 , f [ x 0 ] ) (x_0,f[x_0]) (x0​,f[x0​]), ( x 1 , f [ x 1 ] ) (x_1,f[x_1]) (x1​,f[x1​]),…, ( x n , f [ x n ] ) (x_n,f[x_n]) (xn​,f[xn​])} g [ x ] g[x] g[x]={ ( x 0 , g [ x 0 ] ) (x_0,g[x_0]) (x0​,g[x0​]), ( x 1 , g [ x 1 ] ) (x_1,g[x_1]) (x1​,g[x1​]),…, ( x n , g [ x n ] ) (x_n,g[x_n]) (xn​,g[xn​])} f [ x ] ∗ g [ x ] f[x]*g[x] f[x]∗g[x]={ ( x 0 , f [ x 0 ] ∗ g [ x 0 ] ) (x_0,f[x_0]*g[x_0]) (x0​,f[x0​]∗g[x0​]), ( x 1 , f [ x 1 ] ∗ g [ x 1 ] ) (x_1,f[x_1]*g[x_1]) (x1​,f[x1​]∗g[x1​]),…, ( x n , f [ x n ] ∗ g [ x n ] ) (x_n,f[x_n]*g[x_n]) (xn​,f[xn​]∗g[xn​])} 我们可以发现，如果两个多项式取相同的 x x x，得到不同的 y y y值，那么只需要 y y y值对应相乘就可以了！ 复杂度只有枚举 x x x的 O ( n ) O(n) O(n) 那么问题转换为将多项式系数表示法转化成点值表示法。 朴素系数转点值的算法叫DFT（离散傅里叶变换）,优化后为FFT（快速傅里叶变换），点值转系数的算法叫IDFT（离散傅里叶逆变换），优化后为IFFT(快速傅里叶逆变换)。之后我会分别介绍。

其实不理解卷积也没关系，但这里顺便提一下，可以跳过的 卷积与傅里叶变换有着密切的关系。利用一点性质，即两函数的傅里叶变换的乘积等于它们卷积后的傅里叶变换，能使傅里叶分析中许多问题的处理得到简化。

F ( g ( x ) ∗ f ( x ) ) = F ( g ( x ) ) F ( f ( x ) ) F(g(x)*f(x)) = F(g(x))F(f(x)) F(g(x)∗f(x))=F(g(x))F(f(x))

其中F表示的是傅里叶变换

高中数学会详细讲解，知道的可以跳过这一部分，没学过也没关系，看以下内容应该能很清楚的理解。

数集拓展到实数范围内，仍有些运算无法进行。比如判别式小于0的一元二次方程仍无解，因此将数集再次扩充，达到复数范围。 复数 z z z被定义为二元有序实数对 ( a , b ) (a,b) (a,b),记为 z = a + b i z=a+bi z=a+bi,这里 a a a和 b b b是实数，规定 i i i是虚数单位。 ( i 2 = − 1 i^2=-1 i2=−1 即 i = − 1 i=\sqrt{-1} i=−1 ​) 对于复数 z = a + b i z=a+bi z=a+bi。实数 a a a称为复数z的实部(real part),记作 r e z = a rez=a rez=a.实数 b b b称为复数z的虚部(imaginary part)记作 Imz=b. 当虚部等于零时，这个复数可以视为实数。 即当 b = 0 b=0 b=0时， z = a z=a z=a,这时复数成为实数；当且仅当 a = b = 0 a=b=0 a=b=0时，它是实数0； 当z的虚部不等于零时，实部等于零时，常称z为纯虚数。 即当 a = 0 a=0 a=0且 b ≠ 0 b≠0 b​=0时， z = b i z=bi z=bi，我们就将其称为纯虚数。 将复数的实部与虚部的平方和的正的平方根的值称为该复数的模，记作 ∣ z ∣ ∣z∣ ∣z∣。 即对于复数z=a+bi，它的模为 ∣ z ∣ = ( a 2 + b 2 ) ∣z∣=\sqrt{(a^2+b^2)} ∣z∣=(a2+b2) ​

直接两张图搞定√ (应该可以一目了然)

加法法则: ( a + b i ) + ( c + d i ) = ( a + c ) + ( b + d ) i ; (a+bi)+(c+di)=(a+c)+(b+d)i; (a+bi)+(c+di)=(a+c)+(b+d)i; 减法法则: ( a + b i ) − ( c + d i ) = ( a − c ) + ( b − d ) i ; (a+bi)-(c+di)=(a-c)+(b-d)i; (a+bi)−(c+di)=(a−c)+(b−d)i; 注：复数加减满足平行四边形法则：

乘法法则: ( a + b i ) ⋅ ( c + d i ) = ( a c − b d ) + ( b c + a d ) i (a+bi)·(c+di)=(ac-bd)+(bc+ad)i (a+bi)⋅(c+di)=(ac−bd)+(bc+ad)i 复数相乘一个重要法则：模长相乘，幅角相加。（这个定理很重要） 模长：这个向量的模长，即这个点到原点的距离。（不懂的可再看下向量的几何意义）。 幅角: 从原点出发、指向x轴正半轴的射线绕原点逆时针旋转至过这个点所经过的角。 在极坐标（可看成平面直角坐标系）下，复数可用模长r与幅角θ表示为 ( r , θ ) (r,θ) (r,θ)。对于复数 a + b i a+bi a+bi, r = ( a ² + b ² ) r=\sqrt{(a²+b²)} r=(a²+b²) ​， θ = a r c t a n ( b / a ) θ=arctan(b/a) θ=arctan(b/a)。此时，复数相乘表现为模长相乘，幅角相加。 除法法则： ( a + b i ) ÷ ( c + d i ) = [ ( a c + b d ) / ( c ² + d ² ) ] + [ ( b c − a d ) / ( c ² + d ² ) ] i (a+bi)÷(c+di)=[(ac+bd)/(c²+d²)]+[(bc-ad)/(c²+d²)]i (a+bi)÷(c+di)=[(ac+bd)/(c²+d²)]+[(bc−ad)/(c²+d²)]i

一个复数 z = a + b i z=a+bi z=a+bi的共轭复数为 a − b i a−bi a−bi（实部不变，虚部取反），记为 z ‾ = a − b i \overline{z}=a-bi z=a−bi 当复数模为1时（即|z|=1），与共轭复数互为倒数 证明： z ∗ z ‾ = a 2 − b 2 ∗ i 2 = a 2 + b 2 = ∣ z ∣ 2 = 1 z*\overline{z}=a^2-b^2*i^2=a^2+b^2=|z|^2=1 z∗z=a2−b2∗i2=a2+b2=∣z∣2=1

由于多项式乘法用点值表示比用系数表示快的多，所以我们先要将系数表示法转化成点值表示法相乘，再将结果的点值表示法转化为系数表示法的过程。 第一个过程叫做FFT（快速傅里叶变换)，第二个过程叫IFFT（快速傅里叶逆变换） 在讲这两个过程之前，首先了解一个概念：

复数 ω \omega ω满足 ω n = 1 \omega^n=1 ωn=1，称 ω \omega ω是n次单位根

单位圆:圆心为原点、1为半径的圆 把单位圆n等分，取这n个点（或点表示的向量）所表示的复数（即分别以这n个点的横坐标为实部、纵坐标为虚部，所构成的虚数），即为n次单位根。 下图包含了当n=8时，所有的8次单位根，分别记为 ω 8 1 , ω 8 2 . . . . . , ω 8 8 \omega_8^1,\omega_8^2.....,\omega_8^8 ω81​,ω82​.....,ω88​ (图中圆的半径是1，w表示 ω \omega ω,且下标8已省略） 图是我自己画的，可能有点丑QWQ 由此我们知道如何找单位根啦 从点(1,0)开始（即 ω n 1 \omega_n^1 ωn1​）,逆时针将这n个点从0开始编号，第k个点对应的虚数记作 ω n k \omega_n^k ωnk​ 由复数相乘法则：模长相乘幅角相加​ 可得： ( ω n 1 ) k = ω n k (\omega_n^1)^k=\omega_n^k (ωn1​)k=ωnk​ 根据每个复数的幅角，可以计算出所对应的点/向量。 ω n k \omega_n^k ωnk​ 对应的点/向量是 ( c o s ⁡ k n 2 π , s i n ⁡ k n 2 π ) (cos⁡\frac kn2π,sin⁡\frac kn2π) (cos⁡nk​2π,sin⁡nk​2π)，即为复数 c o s ⁡ k n 2 π + i ∗ s i n ⁡ k n 2 π cos⁡\frac kn2π+i *sin⁡\frac kn2π cos⁡nk​2π+i∗sin⁡nk​2π

建议记住，因为对之后的分析很重要！！

至于怎么证明，就是复数相乘时模长相乘幅角相加的原则。或者你直接观察图也可以很显然的得出结论。​

对于任意多项式系数表示转点值表示，例如 F [ x ] = y = a 0 x 0 + a 1 x 1 + a 2 x 2 + . . . . . . + a n x n F[x]=y=a_0x^0+a_1x^1+a_2x^2+......+a_nx^n F[x]=y=a0​x0+a1​x1+a2​x2+......+an​xn ，可以随便取任意n个 x x x值代入计算,但这样时间复杂度是 O ( n 2 ) O(n^2) O(n2) 所以伟大数学家傅里叶取了一些特殊的点代入，从而进行优化。 他规定了点值表示中的 n n n个 x x x为 n n n个模长为1的复数。这 n n n个复数不是随机的，而是单位根。 把上述的n个复数（单位根） ω n 0 , ω n 1 . . . . . , ω n n − 1 \omega_n^0,\omega_n^1.....,\omega_n^{n-1} ωn0​,ωn1​.....,ωnn−1​代入多项式，能得到一种特殊的点值表示，这种点值表示就叫DFT（离散傅里叶变换）。

虽然DFT能把多项式转换成点值但它仍然是暴力代入 n n n个数，复杂度仍然是O(n2)，所以它只是快速傅里叶变换的朴素版。 所以我们要考虑利用单位根的性质，加速我们的运算，得到FFT（快速傅里叶变换） 对于多项式 A ( x ) = a 0 + a 1 x + a 2 x 2 + . . . + a n − 1 x n − 1 A(x)=a_0+a_1x+a_2x^2+...+a_{n−1}x^{n−1} A(x)=a0​+a1​x+a2​x2+...+an−1​xn−1 将A(x)的每一项按照下标的奇偶分成两部分： A ( x ) = a 0 + a 2 x 2 + . . . + a n − 2 x n − 2 + x ∗ ( a 1 + a 3 x 2 + . . . + a n − 1 x n − 2 ) A(x)=a_0+a_2x^2+...+a_{n−2}x^{n−2}+x*(a_1+a_3x^2+...+a_{n−1}x^{n−2}) A(x)=a0​+a2​x2+...+an−2​xn−2+x∗(a1​+a3​x2+...+an−1​xn−2) 设两个多项式 A 0 ( x ) A_0(x) A0​(x)和 A 1 ( x ) A_1(x) A1​(x)，令： A 0 ( x ) = a 0 x 0 + a 2 x 1 + . . . + a n − 2 x n / 2 − 1 A_0(x)=a_0x^0+a_2x^1+...+a_{n-2}x^{n/2-1} A0​(x)=a0​x0+a2​x1+...+an−2​xn/2−1 A 1 ( x ) = a 1 x 0 + a 3 x 1 + . . . + a n − 1 x n / 2 − 1 A_1(x)=a_1x^0+a_3x^1+...+a_{n-1}x^{n/2-1} A1​(x)=a1​x0+a3​x1+...+an−1​xn/2−1 显然， A ( x ) = A 0 ( x 2 ) + x ∗ A 1 ( x 2 ) A(x)=A_0(x^2)+x*A_1(x^2) A(x)=A0​(x2)+x∗A1​(x2) 假设 k < n kif(ch=='-')f=-1;ch=getchar();} while(ch>='0'&&ch int len=1>1)|((i&1) //i小于rev[i]时才交换，防止同一个元素交换两次，回到它原来的位置。 if(i cp w(1,0); for(int k=j;k n=read(); scanf("%s%s",s1,s2); //读入的数的每一位看成多项式的一项，保存在复数的实部 for(int i=0;i