FFT前言

快速傅里叶变换 (fast Fourier transform),即利用计算机计算离散傅里叶变换（DFT)的高效、快速计算方法的统称，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。

FFT（Fast Fourier Transformation） 是离散傅氏变换（DFT）的快速算法。即为快速傅氏变换。它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。

——百度百科

FFT（Fast Fourier Transformation），中文名快速傅里叶变换，用来加速多项式乘法
朴素高精度乘法时间 O ( n 2 ) O(n^2) O(n2)，但 F F T FFT FFT能 O ( n log ⁡ 2 n ) O(n\log_2 n) O(nlog2​n)的时间解决
F F T FFT FFT名字逼格高，也难懂，其他教程写得让人看不太懂，于是自己随便写一下

• 建议对复数、三角函数相关知识有所耳闻 （不会也无所谓）

• 下面难懂的点我会从网上盗

---

多项式的系数表示法和点值表示法

• F F T FFT FFT其实是一个用 O ( n log ⁡ 2 n ) O(n\log_2n) O(nlog2​n)的时间将一个用系数表示的多项式转换成它的点值表示的算法

• 多项式的系数表示和点值表示可以互相转换

系数表示法

一个n-1次n项多项式 f ( x ) f(x) f(x)可以表示为 f ( x ) = ∑ i = 0 n − 1 a i x i f(x)=\sum^{n-1}_{i=0}a_ix^i f(x)=∑i=0n−1​ai​xi
也可以用每一项的系数来表示 f ( x ) f(x) f(x)，即 f ( x ) = { a 0 , a 1 , a 2 , . . . , a n − 1 } f(x)=\{a_0,a_1,a_2,…,a_{n-1} \} f(x)={a0​,a1​,a2​,…,an−1​}
这就是系数表示法，也就是平时数学课上用的方法

点值表示法

• 把多项式放到平面直角坐标系里面，看成一个函数

• 把 n n n个不同的 x x x代入，会得出 n n n个不同的 y y y，在坐标系内就是 n n n个不同的点

• 那么这 n n n个点唯一确定该多项式，也就是有且仅有一个多项式满足 ∀ k , f ( x k ) = y k ∀k,f(x_k)=y_k ∀k,f(xk​)=yk​

• 理由很简单，把 n n n条式子联立起来成为一个有n条方程的n元方程组，每一项的系数都可以解出来

那么 f ( x ) f(x) f(x)还可以用 f ( x ) = { ( x 0 , f ( x 0 ) ) , ( x 1 , f ( x 1 ) ) , ( x 2 , f ( x 2 ) ) , . . . , ( x n − 1 , f ( x n − 1 ) ) } f(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),…,(x_{n-1},f(x_{n-1}))\} f(x)={(x0​,f(x0​)),(x1​,f(x1​)),(x2​,f(x2​)),…,(xn−1​,f(xn−1​))}来表示
这就是点值表示法

---

高精度乘法下两种多项式表示法的区别

对于两个用系数表示的多项式，我们把它们相乘
设两个多项式分别为 A ( x ) , B ( x ) A(x),B(x) A(x),B(x)
我们要枚举 A A A每一位的系数与 B B B每一位的系数相乘
那么系数表示法做多项式乘法时间复杂度 O ( n 2 ) O(n^2) O(n2)

但两个用点值表示的多项式相乘，只需要 O ( n ) O(n) O(n)的时间

什么意思呢？

设两个点值多项式分别为 f ( x ) = { ( x 0 , f ( x 0 ) ) , ( x 1 , f ( x 1 ) ) , ( x 2 , f ( x 2 ) ) , . . . , ( x n − 1 , f ( x n − 1 ) ) } f(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),…,(x_{n-1},f(x_{n-1}))\} f(x)={(x0​,f(x0​)),(x1​,f(x1​)),(x2​,f(x2​)),…,(xn−1​,f(xn−1​))} g ( x ) = { ( x 0 , g ( x 0 ) ) , ( x 1 , g ( x 1 ) ) , ( x 2 , g ( x 2 ) ) , . . . , ( x n − 1 , g ( x n − 1 ) ) } g(x)=\{(x_0,g(x_0)),(x_1,g(x_1)),(x_2,g(x_2)),…,(x_{n-1},g(x_{n-1}))\} g(x)={(x0​,g(x0​)),(x1​,g(x1​)),(x2​,g(x2​)),…,(xn−1​,g(xn−1​))}
设它们的乘积是 h ( x ) h(x) h(x)，那么 h ( x ) = { ( x 0 , f ( x 0 ) ⋅ g ( x 0 ) ) , ( x 1 , f ( x 1 ) ⋅ g ( x 1 ) ) , . . . , ( x n − 1 , f ( x n − 1 ) ⋅ g ( x n − 1 ) ) } h(x)=\{(x_0,f(x_0)·g(x_0)),(x_1,f(x_1)·g(x_1)),…,(x_{n-1},f(x_{n-1})·g(x_{n-1}))\} h(x)={(x0​,f(x0​)⋅g(x0​)),(x1​,f(x1​)⋅g(x1​)),…,(xn−1​,f(xn−1​)⋅g(xn−1​))}

所以这里的时间复杂度只有一个枚举的 O ( n ) O(n) O(n)

• 突然感觉高精度乘法能 O ( n ) O(n) O(n)暴艹一堆题？

• 但是朴素的系数表示法转点值表示法的算法还是 O ( n 2 ) O(n^2) O(n2)的，逆操作类似

• 朴素系数转点值的算法叫DFT（离散傅里叶变换），点值转系数叫IDFT（离散傅里叶逆变换）

• 难道高精度乘法只能 O ( n 2 ) O(n^2) O(n2)了吗？

---

DFT前置知识&技能

复数

毕竟高中有所以不多说

我们把形如a+bi（a,b均为实数）的数称为复数，其中a称为实部，b称为虚部，i称为虚数单位。当虚部等于零时，这个复数可以视为实数；当z的虚部不等于零时，实部等于零时，常称z为纯虚数。复数域是实数域的代数闭包，也即任何复系数多项式在复数域中总有根。 复数是由意大利米兰学者卡当在十六世纪首次引入，经过达朗贝尔、棣莫弗、欧拉、高斯等人的工作，此概念逐渐为数学家所接受。
——百度百科

初中数学老师会告诉你没有 − 1 \sqrt{-1} −1 ​，但仅限 R R R
扩展至复数集 C C C，定义 i 2 = − 1 i^2=-1 i2=−1，一个复数 z z z可以表示为 z = a + b i ( a , b ∈ R ) z=a+bi(a,b\in R) z=a+bi(a,b∈R)
其中 a a a称为实部， b b b称为虚部， i i i称为虚数单位

• 在复数集中就可以用 i i i表示负数的平方根，如 − 7 = 7 i \sqrt{-7}=\sqrt{7}i −7 ​=7 ​i

还可以把复数看成复平面直角坐标系上的一个点，比如下面

x x x轴就是实数集中的坐标轴， y y y轴就是虚数单位 i i i轴

这个点 ( 2 , 3 ) (2,3) (2,3)表示的复数就是 2 + 3 i 2+3i 2+3i，或者想象它代表的向量为 ( 2 , 3 ) (2,3) (2,3)
其实我们还可以自己想象 （其实没有这种表达方式） 它可以表示为 ( 13 , θ ) (\sqrt{13},\theta) (13 ​,θ)
一个复数 z z z的模定义为它到原点的距离，记为 ∣ z ∣ = a 2 + b 2 |z|=\sqrt{a^2+b^2} ∣z∣=a2+b2 ​
一个复数 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

复数的运算

复数不像点或向量，它和实数一样可以进行四则运算
设两个复数分别为 z 1 = a + b i , z 2 = c + d i z_1=a+bi,z_2=c+di z1​=a+bi,z2​=c+di，那么
z 1 + z 2 = ( a + c ) + ( b + d ) i z_1+z_2=(a+c)+(b+d)i z1​+z2​=(a+c)+(b+d)i z 1 z 2 = ( a c − b d ) + ( a d + b c ) i z_1z_2=(ac−bd)+(ad+bc)i z1​z2​=(ac−bd)+(ad+bc)i

复数相加也满足平行四边形法则

这张是从网上盗的

即 A B + A D = A C AB+AD=AC AB+AD=AC

复数相乘还有一个值得注意的小性质
( a 1 , θ 1 ) ∗ ( a 2 , θ 2 ) = ( a 1 a 2 , θ 1 + θ 2 ) (a_1,\theta_1)*(a_2,\theta_2)=(a_1a_2,\theta_1+\theta_2) (a1​,θ1​)∗(a2​,θ2​)=(a1​a2​,θ1​+θ2​)
即模长相乘，极角相加

---

DFT（离散傅里叶变换）

• 一定注意从这里开始所有的 n n n都默认为 2 2 2的整数次幂

对于任意系数多项式转点值，当然可以随便取任意 n n n个 x x x值代入计算
但是暴力计算 x k 0 , x k 1 , . . . , x k n − 1 ( k ∈ [ 0 , n ) ) x_k^0,x_k^1,…,x_k^{n-1}(k\in[0,n)) xk0​,xk1​,…,xkn−1​(k∈[0,n))当然是 O ( n 2 ) O(n^2) O(n2)的时间
其实可以代入一组神奇的 x x x，代入以后不用做那么多的次方运算
这些 x x x当然不是乱取的，而且取这些 x x x值应该就是 傅里叶 的主意了

考虑一下，如果我们代入一些 x x x，使每个 x x x的若干次方等于 1 1 1，我们就不用做全部的次方运算了
± 1 ±1 ±1是可以的，考虑虚数的话 ± i ±i ±i也可以，但只有这四个数远远不够

• 傅里叶说：这个圆圈上面的点都可以做到

以原点为圆心，画一个半径为 1 1 1的单位圆
那么单位圆上所有的点都可以经过若干次次方得到 1 1 1
傅里叶说还要把它给 n n n等分了，比如此时 n = 8 n=8 n=8

橙色点即为 n = 8 n=8 n=8时要取的点，从 ( 1 , 0 ) (1,0) (1,0)点开始，逆时针从 0 0 0号开始标号，标到 7 7 7号
记编号为 k k k的点代表的复数值为 ω n k \omega_n^k ωnk​，那么由模长相乘，极角相加可知 ( ω n 1 ) k = ω n k (\omega_n^1)^k=\omega_n^k (ωn1​)k=ωnk​
其中 ω n 1 \omega_n^1 ωn1​称为 n n n次单位根，而且每一个 ω \omega ω都可以求出 （我三角函数不好）
ω n k = cos ⁡ k n 2 π + i sin ⁡ k n 2 π \omega_n^k=\cos{k\over n}2π+i\sin{k\over n} 2π ωnk​=cosnk​2π+isinnk​2π

或者说也可以这样解释这条式子

注意 s i n 2 θ + c o s 2 θ = 1 sin^2\theta+cos^2\theta=1 sin2θ+cos2θ=1什么的，就容易理解了

那么 ω n 0 , ω n 1 , . . . , ω n n − 1 \omega^0_n,\omega^1_n,…,\omega^{n-1}_n ωn0​,ωn1​,…,ωnn−1​即为我们要代入的 x 0 , x 1 , . . . , x n − 1 x_0,x_1,…,x_{n-1} x0​,x1​,…,xn−1​

---

单位根的一些性质

推 F F T FFT FFT的过程中需要用到 ω \omega ω的一些性质

ω n k = ω 2 n 2 k \omega^k_n=\omega^{2k}_{2n} ωnk​=ω2n2k​

• 它们表示的点（或向量）表示的复数是相同的

• 证明

• ω n k = c o s k n 2 π + i s i n k n 2 π = c o s 2 k 2 n 2 π + i s i n 2 k 2 n 2 π = ω 2 n 2 k \omega^k_n=cos{k\over n}2π+isin{k\over n} 2π=cos{2k\over 2n}2π+isin{2k\over 2n} 2π=\omega^{2k}_{2n} ωnk​=cosnk​2π+isinnk​2π=cos2n2k​2π+isin2n2k​2π=ω2n2k​

ω n k + n 2 = − ω n k \omega^{k+{n \over 2}}_n=-\omega_n^k ωnk+2n​​=−ωnk​

• 它们表示的点关于原点对称，所表示的复数实部相反，所表示的向量等大反向

• 证明

• ω n n 2 = c o s n 2 n 2 π + i s i n n 2 n 2 π = c o s π + i s i n π = − 1 \omega^{n\over 2}_n=cos{{n\over 2}\over n}2\pi+isin{{n\over 2}\over n}2\pi=cos\pi+isin\pi=-1 ωn2n​​=cosn2n​​2π+isinn2n​​2π=cosπ+isinπ=−1

• （这个东西和 e i x = c o s x + i s i n x e^{ix}=cosx+isinx eix=cosx+isinx与 e i π + 1 = 0 e^{i\pi}+1=0 eiπ+1=0有点关系，我不会就不讲了）

ω n 0 = ω n n \omega^0_n=\omega^n_n ωn0​=ωnn​

• 都等于 1 1 1，或 1 + 0 i 1+0i 1+0i

---

FFT（快速傅里叶变换）

虽然 D F T DFT DFT搞出来一堆很牛逼的 ω \omega ω作为代入多项式的 x x x值
但只是代入这类特殊 x x x值法的变换叫做 D F T DFT DFT而已，还是要代入单位根暴力计算

• DFT还是暴力 O ( n 2 ) O(n^2) O(n2)…

但 D F T DFT DFT可以分治来做，于是 FFT（快速傅里叶变换） 就出来了
首先设一个多项式 A ( x ) A(x) A(x)
A ( x ) = ∑ i = 0 n − 1 a i x i = a 0 + a 1 x + a 2 x 2 + . . . + a n − 1 x n − 1 A(x)=\sum^{n-1}_{i=0}a_ix^i=a_0+a_1x+a_2x^2+…+a_{n-1}x^{n-1} A(x)=i=0∑n−1​ai​xi=a0​+a1​x+a2​x2+…+an−1​xn−1

按 A ( x ) A(x) A(x)下标的奇偶性把 A ( x ) A(x) A(x)分成两半，右边再提一个 x x x

A ( x ) = ( a 0 + a 2 x 2 + . . . + a n − 2 x n − 2 ) + ( a 1 x + a 3 x 3 + . . . + a n − 1 x n − 1 ) A(x)=(a_0+a_2x^2+…+a_{n-2}x^{n-2})+(a_1x+a_3x^3+…+a_{n-1}x^{n-1}) A(x)=(a0​+a2​x2+…+an−2​xn−2)+(a1​x+a3​x3+…+an−1​xn−1)

= ( 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_0+a_2x^2+…+a_{n-2}x^{n-2})+x(a_1+a_3x^2+…+a_{n-1}x^{n-2}) =(a0​+a2​x2+…+an−2​xn−2)+x(a1​+a3​x2+…+an−1​xn−2)

两边好像非常相似，于是再设两个多项式 A 1 ( x ) , A 2 ( x ) A_1(x),A_2(x) A1​(x),A2​(x)，令

A 1 ( x ) = a 0 + a 2 x + a 4 x 2 + . . . + a n − 2 x n 2 − 1 A_1(x)=a_0+a_2x+a_4x^2+…+a_{n-2}x^{{n\over 2}-1} A1​(x)=a0​+a2​x+a4​x2+…+an−2​x2n​−1 A 2 ( x ) = a 1 + a 3 x + a 5 x 2 + . . . + a n − 1 x n 2 − 1 A_2(x)=a_1+a_3x+a_5x^2+…+a_{n-1}x^{{n \over 2}-1} A2​(x)=a1​+a3​x+a5​x2+…+an−1​x2n​−1

比较明显得出
A ( x ) = A 1 ( x 2 ) + x A 2 ( x 2 ) A(x)=A_1(x^2)+xA_2(x^2) A(x)=A1​(x2)+xA2​(x2)

再设 k < n 2 k<{n\over 2} k<2n​，把 ω n k \omega^k_n ωnk​作为 x x x代入 A ( x ) A(x) A(x)（接下来几步变换要多想想单位根的性质）

A ( ω n k ) = A 1 ( ( ω n k ) 2 ) + ω n k A 2 ( ( ω n k ) 2 ) A(\omega^k_n)=A_1((\omega^k_n)^2)+\omega^k_nA_2((\omega^k_n)^2) A(ωnk​)=A1​((ωnk​)2)+ωnk​A2​((ωnk​)2) = A 1 ( ω n 2 k ) + ω n k A 2 ( ω n 2 k ) = A 1 ( ω n 2 k ) + ω n k A 2 ( ω n 2 k ) =A_1(\omega^{2k}_n)+\omega^k_nA_2(\omega^{2k}_n)=A_1(\omega^k_{n\over2})+\omega^k_nA_2(\omega^k_{n\over 2}) =A1​(ωn2k​)+ωnk​A2​(ωn2k​)=A1​(ω2n​k​)+ωnk​A2​(ω2n​k​)

那么对于 A ( ω n k + n 2 ) A(\omega^{k+{n\over2}}_n) A(ωnk+2n​​)的话，代入 ω n k + n 2 \omega^{k+{n \over 2}}_n ωnk+2n​​
A ( ω n k + n 2 ) = A 1 ( ω n 2 k + n ) + ω n k + n 2 A 2 ( ω n 2 k + n ) A(\omega^{k+{n\over 2}}_n)=A_1(\omega^{2k+n}_n)+\omega^{k+{n\over 2}}_nA_2(\omega^{2k+n}_n) A(ωnk+2n​​)=A1​(ωn2k+n​)+ωnk+2n​​A2​(ωn2k+n​) = A 1 ( ω n 2 k ω n n ) − ω n k A 2 ( ω n 2 k ω n n ) =A_1(\omega^{2k}_n\omega^n_n)-\omega^k_nA_2(\omega^{2k}_n\omega^n_n) =A1​(ωn2k​ωnn​)−ωnk​A2​(ωn2k​ωnn​) = A 1 ( ω n 2 k ) − ω n k A 2 ( ω n 2 k ) = A 1 ( ω n 2 k ) − ω n k A 2 ( ω n 2 k ) =A_1(\omega^{2k}_n)-\omega^k_nA_2(\omega^{2k}_n)=A_1(\omega^k_{n\over2})-\omega^k_nA_2(\omega^k_{n\over2}) =A1​(ωn2k​)−ωnk​A2​(ωn2k​)=A1​(ω2n​k​)−ωnk​A2​(ω2n​k​)

• 发现了什么？

A ( ω n k ) A(\omega^k_n) A(ωnk​)和 A ( ω n k + n 2 ) A(\omega^{k+{n\over2}}_n) A(ωnk+2n​​)两个多项式后面一坨东西只有符号不同
就是说，如果已知 A 1 ( ω n 2 k ) A_1(\omega^k_{n\over 2}) A1​(ω2n​k​)和 A 2 ( ω n 2 k ) A_2(\omega^k_{n\over 2}) A2​(ω2n​k​)的值，我们就可以同时知道 A ( ω n k ) A(\omega^k_n) A(ωnk​)和 A ( ω n k + n 2 ) A(\omega^{k+{n\over2}}_n) A(ωnk+2n​​)的值
现在我们就可以递归分治来搞 F F T FFT FFT了

每一次回溯时只扫当前前面一半的序列，即可得出后面一半序列的答案
n = = 1 n==1 n==1时只有一个常数项，直接 r e t u r n return return
时间复杂度 O ( n log ⁡ 2 n ) O(n\log_2n) O(nlog2​n)

---

IFFT（快速傅里叶逆变换）

想一下，我们不仅要会 F F T FFT FFT，还要会IFFT（快速傅里叶逆变换）
我们把两个多项式相乘 （也叫求卷积），做完两遍 F F T FFT FFT也知道了积的多项式的点值表示
可我们平时用系数表示的多项式，点值表示没有意义

• 怎么把点值表示的多项式快速转回系数表示法？

• I D F T IDFT IDFT暴力 O ( n 2 ) O(n^2) O(n2)做？其实也可以用 F F T FFT FFT用 O ( n log ⁡ 2 n ) O(n\log_2n) O(nlog2​n)的时间搞

你有没有想过为什么傅里叶是把 ω n k \omega^k_n ωnk​作为 x x x代入而不是别的什么数？
原因是有的但是有我也看不懂
由于我是沙雕所以只用记住一个结论

• 一个多项式在分治的过程中乘上单位根的共轭复数，分治完的每一项除以 n n n即为原多项式的每一项系数

意思就是说 F F T FFT FFT和 I F F T IFFT IFFT可以一起搞

---

朴素版FFT板子

c + + c++ c++有自带的复数模板 c o m p l e x complex complex库
a . r e a l ( ) a.real() a.real()即表示复数 a a a的实部

#include<complex>
#define cp complex<double>

void fft(cp *a,int n,int inv)//inv是取共轭复数的符号
{
    if (n==1)return;
    int mid=n/2;
    static cp b[MAXN];
    fo(i,0,mid-1)b[i]=a[i*2],b[i+mid]=a[i*2+1];
    fo(i,0,n-1)a[i]=b[i];
    fft(a,mid,inv),fft(a+mid,mid,inv);//分治
    fo(i,0,mid-1)
    {
        cp x(cos(2*pi*i/n),inv*sin(2*pi*i/n));//inv取决是否取共轭复数
        b[i]=a[i]+x*a[i+mid],b[i+mid]=a[i]-x*a[i+mid];
    }
    fo(i,0,n-1)a[i]=b[i];
}

两个多项式 a , b a,b a,b相乘再转系数多项式 c c c，通常只用打这么一小段

    cp a[MAXN],b[MAXN];
    int c[MAXN];
    fft(a,n,1),fft(b,n,1);//1系数转点值
    fo(i,0,n-1)a[i]*=b[i];
    fft(a,n,-1);//-1点值转系数
    fo(i,0,n-1)c[i]=(int)(a[i].real()/n+0.5);//注意精度

很明显，FFT只能处理 n n n为 2 2 2的整数次幂的多项式
所以在 F F T FFT FFT前一定要把 n n n调到 2 2 2的次幂

这个板子看着好像很优美，但是

递归常数太大，要考虑优化…

---

FFTの优化——迭代版FFT

这个图也是盗的

这个很容易发现点什么吧？

• 每个位置分治后的最终位置为其二进制翻转后得到的位置

这样的话我们可以先把原序列变换好，把每个数放在最终的位置上，再一步一步向上合并
一句话就可以 O ( n ) O(n) O(n)预处理第 i i i位最终的位置 r e v [ i ] rev[i] rev[i]

fo(i,0,n-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));

至于蝴蝶变换它死了其实是我不会

---

真·FFT板子

void fft(cp *a,int n,int inv)
{
    int bit=0;
    while ((1<<bit)<n)bit++;
    fo(i,0,n-1)
    {
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
        if (i<rev[i])swap(a[i],a[rev[i]]);//不加这条if会交换两次（就是没交换）
    }
    for (int mid=1;mid<n;mid*=2)//mid是准备合并序列的长度的二分之一
    {
        cp temp(cos(pi/mid),inv*sin(pi/mid));//单位根，pi的系数2已经约掉了
        for (int i=0;i<n;i+=mid*2)//mid*2是准备合并序列的长度，i是合并到了哪一位
        {
            cp omega(1,0);
            for (int j=0;j<mid;j++,omega*=temp)//只扫左半部分，得到右半部分的答案
            {
                cp x=a[i+j],y=omega*a[i+j+mid];
                a[i+j]=x+y,a[i+j+mid]=x-y;//这个就是蝴蝶变换什么的
            }
        }
    }
}

这个板子好像不是那么好背
至少这个板子已经很优美了

---

FFT后记

本人版权意识薄弱……

• 本博客部分知识学习于

• https://www.cnblogs.com/RabbitHu/p/FFT.html

• https://www.cnblogs.com/zwfymqz/p/8244902.html?mType=Group#_label3

• https://blog.csdn.net/ggn_2015/article/details/68922404

N T T NTT NTT我来了