FFT(快速傅里叶变换)算法

功能

一次FFT的功能

O ( n log ⁡ n ) O(n\log n) O(nlogn)的时间将一个多项式的系数表达式变成点值表达式(这两种表达方式后文会提到)。

一次IFFT的功能

O ( n log ⁡ n ) O(n\log n) O(nlogn)的时间将一个多项式的点值表达式变成系数表达式(这两种表达方式后文会提到)。

总体功能

O ( n log ⁡ n ) O(n\log n) O(nlogn)的时间求两个多项式的乘积。

一个 n n n项的整式多项式一定可以表示成: f ( x ) = ∑ i = 0 n − 1 ( a i x i ) f(x)=\sum\limits_{i=0}^{n-1}(a_ix^i) f(x)=i=0n1(aixi)
(有些 a i a_i ai可能为 0 0 0

所以,一个系数 a a a的(有序)集合可以唯一确定一个多项式。
例如: a = { 1 , 2 , 0 , 3 , 5 } a=\{1,2,0,3,5\} a={1,2,0,3,5}表示的就是多项式: f ( x ) = 1 + 2 x + 3 x 3 + 5 x 4 f(x)=1+2x+3x^3+5x^4 f(x)=1+2x+3x3+5x4
这里 a = { 1 , 2 , 0 , 3 , 5 } a=\{1,2,0,3,5\} a={1,2,0,3,5} f ( x ) f(x) f(x)系数表达式

所以,所谓求两个多项式的乘积,就是已知两个多项式的系数表达式,求它们乘积的多项式的系数表达式

前置技能

(既然大家都这样说,我也这样说吧)

多项式的阶

之前说的 f ( x ) = ∑ i = 0 n − 1 ( a i x i ) f(x)=\sum\limits_{i=0}^{n-1}(a_ix^i) f(x)=i=0n1(aixi)中的 n n n就是它的阶,注意 n n n阶多项式不一定是 n n n次的,因为 a n − 1 a_{n-1} an1可能为 0 0 0

多项式的系数表达式

前面已说。

多项式的点值表达式

一个 n n n次多项式就是一个函数,取它图像上的 n + 1 n+1 n+1个不同的点就可以确定这个多项式。

例如 f ( x ) f(x) f(x)是个 4 4 4次多项式,取它的图像上的任意 5 5 5个不同的点: ( x 1 , f ( x 1 ) ) (x_1,f(x_1)) (x1,f(x1)) ( x 2 , f ( x 2 ) ) (x_2,f(x_2)) (x2,f(x2)),…, ( x 5 , f ( x 5 ) ) (x_5,f(x_5)) (x5,f(x5)),可以确定这个多项式的系数,因为可以列出一个 5 5 5元一次方程组解出系数。

于是一个 n n n次多项式的点值表达式就是它的图像上 n + 1 n+1 n+1个不同的点。
换句话说, n + 1 n+1 n+1个不同的点可以唯一确定一个 n n n次多项式。
(以上两句话式FFT和逆FFT的核心)

复数

在实数范围内,老师会告诉你负数没有平方根,所以我们为了让负数有平方根,将数的范围从实数扩展到了复数。
(实数也属于复数,复数另外一部分是虚数)

复数的基本单位

引入一个符号 i i i i 2 = − 1 i^2=-1 i2=1
于是任意一个复数都可以表示成 a + b i a+bi a+bi,其中 a a a b b b都是实数。
例如: x 2 = − 7 x^2=-7 x2=7的根是: x 1 = 7 i x_1=\sqrt 7i x1=7 i x 2 = − 7 i x_2=-\sqrt 7i x2=7 i

复数的运算

加减乘都和整式的运算是一样的,例如两个复数 z 1 = a 1 + b 1 i z_1=a_1+b_1i z1=a1+b1i z 2 = a 2 + b 2 i z_2=a_2+b_2i z2=a2+b2i相乘: z 1 z 2 = ( a 1 + b 1 i ) ( a 2 + b 2 i ) = ( a 1 a 2 − b 1 b 2 ) + ( a 1 b 2 + a 2 b 1 ) i z_1z_2=(a_1+b_1i)(a_2+b_2i)=(a_1a_2-b_1b_2)+(a_1b_2+a_2b_1)i z1z2=(a1+b1i)(a2+b2i)=(a1a2b1b2)+(a1b2+a2b1)i
除法这里不会用到,但也不难,可以分母实数化:

复数的除法
来自百度百科:复数除法

复平面

数轴可以表示一切实数,平面可以表示一切复数:
复平面
(其实这个图没有什么用)
在复平面内的一个点 ( x , y ) (x,y) (x,y)表示一个复数 x + y i x+yi x+yi

复根

定义

定义: ω n k = c o s ( 2 π n k ) + s i n ( 2 π n k ) i \omega_n^k=cos\left(\frac{2\pi}{n}k\right)+sin\left(\frac{2\pi}{n}k\right)i ωnk=cos(n2πk)+sin(n2πk)i
ω n 1 \omega_n^1 ωn1 n n n次单位复根。

很难懂是不是?把 ω 8 0 \omega_8^0 ω80,…, ω 8 7 \omega_8^7 ω87对应的点在复平面上画出来你就知道了。
8次单位复根
(由于几何画板的公式编辑太弱,所以只标了 A B C ABC ABC
A A A B B B C C C…对应的复数分别是 ω 8 0 \omega_8^0 ω80 ω 8 1 \omega_8^1 ω81 ω 8 2 \omega_8^2 ω82,…

为什么?数学必修四的单位圆与三角函数会告诉你的。

几个性质
  • ω n i ω n j = ω n i + j \omega_n^i\omega_n^j=\omega_n^{i+j} ωniωnj=ωni+j

证明:
由于 e i θ = c o s θ + s i n θ ⋅ i e^{i\theta}=cos\theta+sin\theta\cdot i eiθ=cosθ+sinθi(我不会证,貌似要用泰勒展开)
所以 ω n i = e 2 π n i \omega_n^i=e^{\frac{2\pi}{n}i} ωni=en2πi
所以 ω n i ω n j = e 2 π n i e 2 π n j = e 2 π n ( i + j ) = ω n i + j \omega_n^i\omega_n^j=e^{\frac{2\pi}{n}i}e^{\frac{2\pi}{n}j}=e^{\frac{2\pi}{n}(i+j)}=\omega_n^{i+j} ωniωnj=en2πien2πj=en2π(i+j)=ωni+j

  • ω 2 n 2 k = ω n k \omega_{2n}^{2k}=\omega_n^k ω2n2k=ωnk
    理解:它们对应的点是一样的

证明:
ω 2 n 2 k = c o s ( 2 π 2 n 2 k ) + s i n ( 2 π 2 n 2 k ) i = c o s ( 2 π n k ) + s i n ( 2 π n k ) i = ω n k \omega_{2n}^{2k}=cos\left(\frac{2\pi}{2n}2k\right)+sin\left(\frac{2\pi}{2n}2k\right)i=cos\left(\frac{2\pi}{n}k\right)+sin\left(\frac{2\pi}{n}k\right)i=\omega_n^k ω2n2k=cos(2n2π2k)+sin(2n2π2k)i=cos(n2πk)+sin(n2πk)i=ωnk

  • ω n k + n 2 = − ω n k \omega_n^{k+\frac{n}{2}}=-\omega_n^k ωnk+2n=ωnk n n n是偶数)
    理解:它们对应的点关于原点对称。

证明: ω n k + n 2 = c o s ( 2 π n k + π ) + s i n s ( 2 π n k + π ) i = − c o s ( 2 π n k ) − s i n ( 2 π n k ) i = − ω n k \omega_{n}^{k+\frac{n}{2}}=cos\left(\frac{2\pi}{n}k+\pi\right)+sins\left(\frac{2\pi}{n}k+\pi\right)i=-cos\left(\frac{2\pi}{n}k\right)-sin\left(\frac{2\pi}{n}k\right)i=-\omega_n^k ωnk+2n=cos(n2πk+π)+sins(n2πk+π)i=cos(n2πk)sin(n2πk)i=ωnk

  • ω n k + n = ω n k \omega_n^{k+n}=\omega_n^k ωnk+n=ωnk
    理解:转了一圈回到原来的点。

证明:同上一个证明

  • ∑ i = 0 n − 1 ( ω n k ) i = 0 \sum\limits_{i=0}^{n-1}(\omega_n^k)^i=0 i=0n1(ωnk)i=0 k ≠ 0 k\neq 0 k=0)(求和引理)
    理解:每个点和它关于原点对称的点抵消了。

证明:
∑ i = 0 n − 1 ( ω n k ) i = 1 + ω n k + ( ω n k ) 2 + . . . + ( ω n k ) n − 1 = ( ω n k ) n − 1 ω n k − 1 = ω n k n − 1 ω n k − 1 = 1 − 1 ω n k − 1 = 0 \sum\limits_{i=0}^{n-1}(\omega_n^k)^i=1+\omega_n^k+(\omega_n^k)^2+...+(\omega_n^k)^{n-1}=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}=\frac{\omega_n^{kn}-1}{\omega_n^k-1}=\frac{1-1}{\omega_n^k-1}=0 i=0n1(ωnk)i=1+ωnk+(ωnk)2+...+(ωnk)n1=ωnk1(ωnk)n1=ωnk1ωnkn1=ωnk111=0

求多项式乘积的基本步骤

已知 k 1 k_1 k1阶多项式 f ( x ) f(x) f(x) k 2 k_2 k2阶多项式 g ( x ) g(x) g(x)的系数表达式。
我们要求的是 k 1 + k 2 k_1+k_2 k1+k2阶多项式 h ( x ) = f ( x ) ⋅ g ( x ) h(x)=f(x)\cdot g(x) h(x)=f(x)g(x)的系数表达式。

  1. 统一两个多项式的阶数为 n n n,要保证 n n n 2 2 2的幂。(系数补零即可,为什么一会说)
  2. n n n个不同的值 x x x,代入 f ( x ) f(x) f(x) g ( x ) g(x) g(x)中,得到 f ( x ) f(x) f(x) g ( x ) g(x) g(x)的点值表达式: S f ( x ) = { ( x 1 , f ( x 1 ) ) , ( x 2 , f ( x 2 ) ) , . . . , ( x n , f ( x n ) ) } S_{f(x)}=\{(x_1,f(x_1)),(x_2,f(x_2)),...,(x_n,f(x_n))\} Sf(x)={(x1,f(x1)),(x2,f(x2)),...,(xn,f(xn))} S g ( x ) = { ( x 1 , g ( x 1 ) ) , ( x 2 , g ( x 2 ) ) , . . . , ( x n , g ( x n ) ) } S_{g(x)}=\{(x_1,g(x_1)),(x_2,g(x_2)),...,(x_n,g(x_n))\} Sg(x)={(x1,g(x1)),(x2,g(x2)),...,(xn,g(xn))}(两次FFT得到)
  3. 计算 h ( x ) h(x) h(x)的点值表达式: S h ( x ) = { ( x 1 , f ( x 1 ) ⋅ g ( x 1 ) ) , ( x 2 , f ( x 2 ) ⋅ g ( x 2 ) ) , . . . , ( x n , f ( x n ) ⋅ g ( x n ) ) } S_{h(x)}=\{(x_1,f(x_1)\cdot g(x_1)),(x_2,f(x_2)\cdot g(x_2)),...,(x_n,f(x_n)\cdot g(x_n))\} Sh(x)={(x1,f(x1)g(x1)),(x2,f(x2)g(x2)),...,(xn,f(xn)g(xn))}
  4. 通过 h ( x ) h(x) h(x)的点值表达式求 h ( x ) h(x) h(x)的系数表达式。
    (一次IFFT得到)

忽略极大的常数,暴力的时间复杂度是 O ( n 2 ) O(n^2) O(n2),FFT可以把它优化到 O ( n log ⁡ n ) O(n\log n) O(nlogn)

FFT(IFFT)是通过计算一半的点值表达式得到另一半,要计算的那一半重复前面的步骤。

FFT

[特别提醒]

  • 因为点值表达式中的 x x x集合取何值无影响,所以取算的越快的越好。
  • 一次FFT(快速傅里叶变换)是将一个多项式的系数表达式集合变成特定( x x x ω n 0 , . . . , ω n n − 1 \omega_n^0,...,\omega_n^{n-1} ωn0,...,ωnn1)点值表达式中的 y y y集合。
  • 一次IFFT(快速傅里叶变换逆变换)是将一个多项式的特定( x x x ω n 0 , . . . , ω n n − 1 \omega_n^0,...,\omega_n^{n-1} ωn0,...,ωnn1)点值表达式中的 y y y集合变成系数表达式集合。
  • 注意 n n n 2 2 2的幂。

递归版FFT

核心公式

我们代入求点值表达式的 n n n x x x是: ω n 0 , . . . , ω n n − 1 \omega_n^0,...,\omega_n^{n-1} ωn0,...,ωnn1

设多项式是 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+a1x+a2x2+...+an1xn1
A 0 ( x ) = a 0 + a 2 x + a 4 x 2 + . . . + a n − 2 x n 2 − 1 A_0(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{\frac{n}{2}-1} A0(x)=a0+a2x+a4x2+...+an2x2n1
A 1 ( x ) = a 1 + a 3 x + a 5 x 2 + . . . + a n − 1 x n 2 − 1 A_1(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{\frac{n}{2}-1} A1(x)=a1+a3x+a5x2+...+an1x2n1
显然 A ( x ) = A 0 ( x 2 ) + x A 1 ( x 2 ) A(x)=A_0(x^2)+xA_1(x^2) A(x)=A0(x2)+xA1(x2)(可以自己带进去算一下)

那么 A ( ω n k ) = A 0 ( ω n 2 k ) + ω n k A 1 ( ω n 2 k ) = A 0 ( ω n 2 k ) + ω n k A 1 ( ω n 2 k ) A(\omega_n^k)=A_0(\omega_n^{2k})+\omega_n^kA_1(\omega_n^{2k})=A_0(\omega_{\frac{n}{2}}^{k})+\omega_n^kA_1(\omega_{\frac{n}{2}}^{k}) A(ωnk)=A0(ωn2k)+ωnkA1(ωn2k)=A0(ω2nk)+ωnkA1(ω2nk)
A ( ω n k + n 2 ) = A 0 ( ω n 2 k + n ) + ω n k + n 2 A 1 ( ω n 2 k + n ) = A 0 ( ω n 2 k ) − ω n k A 1 ( ω n 2 k ) = A 0 ( ω n 2 k ) − ω n k A 1 ( ω n 2 k ) A(\omega_n^{k+\frac{n}{2}})=A_0(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A_1(\omega_n^{2k+n})=A_0(\omega_{n}^{2k})-\omega_n^kA_1(\omega_{n}^{2k})=A_0(\omega_{\frac{n}{2}}^{k})-\omega_n^kA_1(\omega_{\frac{n}{2}}^{k}) A(ωnk+2n)=A0(ωn2k+n)+ωnk+2nA1(ωn2k+n)=A0(ωn2k)ωnkA1(ωn2k)=A0(ω2nk)ωnkA1(ω2nk)

这两个只有正负号不同!

算法流程

我们定义FFT(A,n),表示将系数表达式的系数集合 A A A转换为点值表达式的 y y y值集合(代入的 x x x ω n 0 , ω n 1 , . . . ω n n − 1 \omega_n^0,\omega_n^1,...\omega_n^{n-1} ωn0,ωn1,...ωnn1)。
那么,FFT(A,n)中,先提出 A 0 A_0 A0 A 1 A_1 A1,再执行FFT(A0,n/2)FFT(A1,n/2)
执行完后的 A 0 A_0 A0 A 1 A_1 A1已经变成了当 “ 代入的 x x x ω n 2 0 , ω n 2 1 , . . . ω n 2 n 2 − 1 \omega_\frac{n}{2}^0,\omega_\frac{n}{2}^1,...\omega_\frac{n}{2}^{\frac{n}{2}-1} ω2n0,ω2n1,...ω2n2n1 ” 时,对应的 y y y的集合。
然后就可以循环从 0 0 0 n 2 \frac{n}{2} 2nA[i]=A0[i]+w*A1[i]A[i+n/2]=A0[i]-w*A1[i]

代码

等讲了IFFT的递归版一起看。

非递归版FFT

算法

非递归版展示
(来自不知名的大佬

这是 n = 8 n=8 n=8时的递归FFT的系数示意图。
发现最后得到的系数序列,每个数的二进制反转后恰好是最开始的序列(为什么?用心去感受一下吧)。

所以只需要把最开始的系数变成最后的样子,再一层层向上合并即可。
也就是说,把最开始的每个系数和它二进制反转后的数交换位置。

定义R[i]表示i的二进制反转后的数,Li的二进制位数(即 log ⁡ 2 n \log_2 n log2n),那么有:
R[i]=(R[i>>1]>>1)|((i&1)<<(L-1))

下面是解释:
R[i>>1]表示i除去最后一位后前面几位的反转,那么把它和i的最后一位换个位置就好了。

代码

等讲了IFFT的非递归版一起看。

IFFT

FFT的矩阵表达

( y 0 y 1 y 2 y 3 ⋮ y n − 1 ) = ( ω n 0 ω n 0 ω n 0 ω n 0 ⋯ ω n 0 ω n 0 ω n 1 ω n 2 ω n 3 ⋯ ω n n − 1 ω n 0 ω n 2 ω n 4 ω n 6 ⋯ ω n 2 ( n − 1 ) ω n 0 ω n 3 ω n 6 ω n 9 ⋯ ω n 3 ( n − 1 ) ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ ω n 0 ω n n − 1 ω n 2 ( n − 1 ) ω n 3 ( n − 1 ) ⋯ ω n ( n − 1 ) ( n − 1 ) ) × ( a 0 a 1 a 2 a 3 ⋮ a n − 1 ) \left( \begin{matrix} y_0\\ \\ y_1\\ \\ y_2\\ \\ y_3\\ \\ \vdots\\ \\ y_{n-1} \end{matrix} \right)= \left( \begin{matrix} \omega_n^0 & \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\\ \\ \omega_n^0 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n-1}\\ \\ \omega_n^0 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2(n-1)}\\ \\ \omega_n^0 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3(n-1)}\\ \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots\\ \\ \omega_n^0 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \cdots & \omega_n^{(n-1)(n-1)}\\ \end{matrix} \right)\times\\ \left( \begin{matrix} a_0\\ \\ a_1\\ \\ a_2\\ \\ a_3\\ \\ \vdots\\ \\ a_{n-1} \end{matrix} \right) y0y1y2y3yn1=ωn0ωn0ωn0ωn0ωn0ωn0ωn1ωn2ωn3ωnn1ωn0ωn2ωn4ωn6ωn2(n1)ωn0ωn3ωn6ωn9ωn3(n1)ωn0ωnn1ωn2(n1)ωn3(n1)ωn(n1)(n1)×a0a1a2a3an1

FFT可以表示为这样, { y } \{y\} {y}是点值表达式的 y y y值集合, { a } \{a\} {a}是系数表达式的系数集合。
(矩阵乘法可以自己百度一下,在这里就是 y k = ∑ i = 0 n − 1 ( a i ω n k i ) y_k=\sum\limits_{i=0}^{n-1}\left(a_i\omega_n^{ki}\right) yk=i=0n1(aiωnki)

例如: y 2 = f ( x 2 ) = f ( ω n 2 ) = a 0 ( ω n 2 ) 0 + a 1 ( ω n 2 ) 1 + a 2 ( ω n 2 ) 2 + . . . + a n − 1 ( ω n 2 ) n − 1 y_2=f(x_2)=f(\omega_n^2)=a_0(\omega_n^2)^0+a_1(\omega_n^2)^1+a_2(\omega_n^2)^2+...+a_{n-1}(\omega_n^2)^{n-1} y2=f(x2)=f(ωn2)=a0(ωn2)0+a1(ωn2)1+a2(ωn2)2+...+an1(ωn2)n1 = a 0 ω n 0 + a 1 ω n 2 + a 2 ω n 4 + . . . + a n − 1 ω n 2 ( n − 1 ) =a_0\omega_n^0+a_1\omega_n^2+a_2\omega_n^4+...+a_{n-1}\omega_n^{2(n-1)} =a0ωn0+a1ωn2+a2ωn4+...+an1ωn2(n1) = ∑ i = 0 n − 1 ( a i ω n 2 i ) =\sum\limits_{i=0}^{n-1}\left(a_i\omega_n^{2i}\right) =i=0n1(aiωn2i)

IFFT的矩阵表达

V = ( ω n 0 ω n 0 ω n 0 ω n 0 ⋯ ω n 0 ω n 0 ω n 1 ω n 2 ω n 3 ⋯ ω n n − 1 ω n 0 ω n 2 ω n 4 ω n 6 ⋯ ω n 2 ( n − 1 ) ω n 0 ω n 3 ω n 6 ω n 9 ⋯ ω n 3 ( n − 1 ) ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ ω n 0 ω n n − 1 ω n 2 ( n − 1 ) ω n 3 ( n − 1 ) ⋯ ω n ( n − 1 ) ( n − 1 ) ) V=\left( \begin{matrix} \omega_n^0 & \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\\ \\ \omega_n^0 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n-1}\\ \\ \omega_n^0 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2(n-1)}\\ \\ \omega_n^0 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3(n-1)}\\ \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots\\ \\ \omega_n^0 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \cdots & \omega_n^{(n-1)(n-1)}\\ \end{matrix} \right) V=ωn0ωn0ωn0ωn0ωn0ωn0ωn1ωn2ωn3ωnn1ωn0ωn2ωn4ωn6ωn2(n1)ωn0ωn3ωn6ωn9ωn3(n1)ωn0ωnn1ωn2(n1)ωn3(n1)ωn(n1)(n1)
( y 0 y 1 y 2 y 3 ⋮ y n − 1 ) = V × ( a 0 a 1 a 2 a 3 ⋮ a n − 1 ) \left( \begin{matrix} y_0\\ \\ y_1\\ \\ y_2\\ \\ y_3\\ \\ \vdots\\ \\ y_{n-1} \end{matrix} \right)= V\times\\ \left( \begin{matrix} a_0\\ \\ a_1\\ \\ a_2\\ \\ a_3\\ \\ \vdots\\ \\ a_{n-1} \end{matrix} \right) y0y1y2y3yn1=V×a0a1a2a3an1
那么 ( y 0 y 1 y 2 y 3 ⋮ y n − 1 ) × V − 1 = ( a 0 a 1 a 2 a 3 ⋮ a n − 1 ) \left( \begin{matrix} y_0\\ \\ y_1\\ \\ y_2\\ \\ y_3\\ \\ \vdots\\ \\ y_{n-1} \end{matrix} \right)\times V^{-1}= \left( \begin{matrix} a_0\\ \\ a_1\\ \\ a_2\\ \\ a_3\\ \\ \vdots\\ \\ a_{n-1} \end{matrix} \right) y0y1y2y3yn1×V1=a0a1a2a3an1

也就是说,我们只要构造出 V − 1 V^{-1} V1,IFFT就解决了。


换句话说,要构造一个矩阵 V − 1 V^{-1} V1,使得 V × V − 1 = A V\times V^{-1}=A V×V1=A A A A是单位矩阵)
其中 A = ( 1 0 0 ⋯ 0 0 1 0 ⋯ 0 0 0 1 ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 0 ⋯ 1 ) A= \left( \begin{matrix} 1 & 0 & 0 & \cdots & 0\\ 0 & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & 1 \end{matrix} \right) A=1000010000100001

(计算一下会发现任何矩阵乘 A A A都是它本身)


于是傅里叶说,
V − 1 = ( 1 n ω n 0 1 n ω n 0 1 n ω n 0 1 n ω n 0 ⋯ 1 n ω n 0 1 n ω n 0 1 n ω n − 1 1 n ω n − 2 1 n ω n − 3 ⋯ 1 n ω n − ( n − 1 ) 1 n ω n 0 1 n ω n − 2 1 n ω n − 4 1 n ω n − 6 ⋯ 1 n ω n − 2 ( n − 1 ) 1 n ω n 0 1 n ω n − 3 1 n ω n − 6 1 n ω n − 9 ⋯ 1 n ω n − 3 ( n − 1 ) ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ 1 n ω n 0 1 n ω n − ( n − 1 ) 1 n ω n − 2 ( n − 1 ) 1 n ω n − 3 ( n − 1 ) ⋯ 1 n ω n − ( n − 1 ) ( n − 1 ) ) V^{-1}= \left( \begin{matrix} \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^0 & \cdots & \frac{1}{n}\omega_n^0\\ \\ \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^{-1} & \frac{1}{n}\omega_n^{-2} & \frac{1}{n}\omega_n^{-3} & \cdots & \frac{1}{n}\omega_n^{-(n-1)}\\ \\ \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^{-2} & \frac{1}{n}\omega_n^{-4} & \frac{1}{n}\omega_n^{-6} & \cdots & \frac{1}{n}\omega_n^{-2(n-1)}\\ \\ \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^{-3} & \frac{1}{n}\omega_n^{-6} & \frac{1}{n}\omega_n^{-9} & \cdots & \frac{1}{n}\omega_n^{-3(n-1)}\\ \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots\\ \\ \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^{-(n-1)} & \frac{1}{n}\omega_n^{-2(n-1)} & \frac{1}{n}\omega_n^{-3(n-1)} & \cdots & \frac{1}{n}\omega_n^{-(n-1)(n-1)}\\ \end{matrix} \right) V1=n1ωn0n1ωn0n1ωn0n1ωn0n1ωn0n1ωn0n1ωn1n1ωn2n1ωn3n1ωn(n1)n1ωn0n1ωn2n1ωn4n1ωn6n1ωn2(n1)n1ωn0n1ωn3n1ωn6n1ωn9n1ωn3(n1)n1ωn0n1ωn(n1)n1ωn2(n1)n1ωn3(n1)n1ωn(n1)(n1)

简单点说,(下标从零开始) V i , j = ω n i j V_{i,j}=\omega_n^{ij} Vi,j=ωnij V − 1 i , j = 1 n ω n − i j V^{-1}{}_{i,j}=\frac{1}{n}\omega_n^{-ij} V1i,j=n1ωnij


接下来要证明 V × V − 1 = A = ( 1 0 0 ⋯ 0 0 1 0 ⋯ 0 0 0 1 ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 0 ⋯ 1 ) V\times V^{-1}=A=\left( \begin{matrix} 1 & 0 & 0 & \cdots & 0\\ 0 & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & 1 \end{matrix} \right) V×V1=A=1000010000100001

根据矩阵乘法的定义, A x , y = ∑ i = 0 n − 1 ( V x , i V − 1 i , y ) A_{x,y}=\sum\limits_{i=0}^{n-1}(V_{x,i}V^{-1}{}_{i,y}) Ax,y=i=0n1(Vx,iV1i,y) = ∑ i = 0 n − 1 ( ω n x i 1 n ω n − i y ) =\sum\limits_{i=0}^{n-1}\left(\omega_n^{xi}\frac{1}{n}\omega_n^{-iy}\right) =i=0n1(ωnxin1ωniy) = 1 n ∑ i = 0 n − 1 ( ω n x − y ) i =\frac{1}{n}\sum\limits_{i=0}^{n-1}(\omega_n^{x-y})^i =n1i=0n1(ωnxy)i

x = y x=y x=y时, ω n x − y = 1 \omega_n^{x-y}=1 ωnxy=1 A i , j = 1 A_{i,j}=1 Ai,j=1
x ≠ y x\neq y x=y时,根据求和引理, A i , j = 0 A_{i,j}=0 Ai,j=0

得证。


综上所述:

  • FFT
    ( y 0 y 1 y 2 y 3 ⋮ y n − 1 ) = ( a 0 a 1 a 2 a 3 ⋮ a n − 1 ) × ( ω n 0 ω n 0 ω n 0 ω n 0 ⋯ ω n 0 ω n 0 ω n 1 ω n 2 ω n 3 ⋯ ω n n − 1 ω n 0 ω n 2 ω n 4 ω n 6 ⋯ ω n 2 ( n − 1 ) ω n 0 ω n 3 ω n 6 ω n 9 ⋯ ω n 3 ( n − 1 ) ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ ω n 0 ω n n − 1 ω n 2 ( n − 1 ) ω n 3 ( n − 1 ) ⋯ ω n ( n − 1 ) ( n − 1 ) ) \left( \begin{matrix} y_0\\ \\ y_1\\ \\ y_2\\ \\ y_3\\ \\ \vdots\\ \\ y_{n-1} \end{matrix} \right)= \left( \begin{matrix} a_0\\ \\ a_1\\ \\ a_2\\ \\ a_3\\ \\ \vdots\\ \\ a_{n-1} \end{matrix} \right)\times\\ \left( \begin{matrix} \omega_n^0 & \omega_n^0 & \omega_n^0 & \omega_n^0 & \cdots & \omega_n^0\\ \\ \omega_n^0 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n-1}\\ \\ \omega_n^0 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2(n-1)}\\ \\ \omega_n^0 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3(n-1)}\\ \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots\\ \\ \omega_n^0 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \cdots & \omega_n^{(n-1)(n-1)}\\ \end{matrix} \right) y0y1y2y3yn1=a0a1a2a3an1×ωn0ωn0ωn0ωn0ωn0ωn0ωn1ωn2ωn3ωnn1ωn0ωn2ωn4ωn6ωn2(n1)ωn0ωn3ωn6ωn9ωn3(n1)ωn0ωnn1ωn2(n1)ωn3(n1)ωn(n1)(n1)
  • IFFT
    ( a 0 a 1 a 2 a 3 ⋮ a n − 1 ) = ( y 0 y 1 y 2 y 3 ⋮ y n − 1 ) × ( 1 n ω n 0 1 n ω n 0 1 n ω n 0 1 n ω n 0 ⋯ 1 n ω n 0 1 n ω n 0 1 n ω n − 1 1 n ω n − 2 1 n ω n − 3 ⋯ 1 n ω n − ( n − 1 ) 1 n ω n 0 1 n ω n − 2 1 n ω n − 4 1 n ω n − 6 ⋯ 1 n ω n − 2 ( n − 1 ) 1 n ω n 0 1 n ω n − 3 1 n ω n − 6 1 n ω n − 9 ⋯ 1 n ω n − 3 ( n − 1 ) ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ 1 n ω n 0 1 n ω n − ( n − 1 ) 1 n ω n − 2 ( n − 1 ) 1 n ω n − 3 ( n − 1 ) ⋯ 1 n ω n − ( n − 1 ) ( n − 1 ) ) \left( \begin{matrix} a_0\\ \\ a_1\\ \\ a_2\\ \\ a_3\\ \\ \vdots\\ \\ a_{n-1} \end{matrix} \right)= \left( \begin{matrix} y_0\\ \\ y_1\\ \\ y_2\\ \\ y_3\\ \\ \vdots\\ \\ y_{n-1} \end{matrix} \right)\times\left( \begin{matrix} \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^0 & \cdots & \frac{1}{n}\omega_n^0\\ \\ \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^{-1} & \frac{1}{n}\omega_n^{-2} & \frac{1}{n}\omega_n^{-3} & \cdots & \frac{1}{n}\omega_n^{-(n-1)}\\ \\ \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^{-2} & \frac{1}{n}\omega_n^{-4} & \frac{1}{n}\omega_n^{-6} & \cdots & \frac{1}{n}\omega_n^{-2(n-1)}\\ \\ \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^{-3} & \frac{1}{n}\omega_n^{-6} & \frac{1}{n}\omega_n^{-9} & \cdots & \frac{1}{n}\omega_n^{-3(n-1)}\\ \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots\\ \\ \frac{1}{n}\omega_n^0 & \frac{1}{n}\omega_n^{-(n-1)} & \frac{1}{n}\omega_n^{-2(n-1)} & \frac{1}{n}\omega_n^{-3(n-1)} & \cdots & \frac{1}{n}\omega_n^{-(n-1)(n-1)}\\ \end{matrix} \right) a0a1a2a3an1=y0y1y2y3yn1×n1ωn0n1ωn0n1ωn0n1ωn0n1ωn0n1ωn0n1ωn1n1ωn2n1ωn3n1ωn(n1)n1ωn0n1ωn2n1ωn4n1ωn6n1ωn2(n1)n1ωn0n1ωn3n1ωn6n1ωn9n1ωn3(n1)n1ωn0n1ωn(n1)n1ωn2(n1)n1ωn3(n1)n1ωn(n1)(n1)

算法流程

既然构造出来了 V − 1 V^{-1} V1,我们就知道要代入计算的 x x x集合就是:
{ ω n 0 , ω n − 1 , . . . , ω n − ( n − 1 ) } \{\omega_n^0,\omega_n^{-1},...,\omega_n^{-(n-1)}\} {ωn0,ωn1,...,ωn(n1)}

为了减少误差,最后统一乘 1 n \frac{1}{n} n1

把 “ 求多项式乘积的基本步骤 ” 中的第3步得到的 h ( x ) h(x) h(x)的点值表达式的 y y y集合当做一个多项式的系数表达式集合,用一次FFT求出 x x x集合是 { ω n 0 , ω n − 1 , . . . , ω n − ( n − 1 ) } \{\omega_n^0,\omega_n^{-1},...,\omega_n^{-(n-1)}\} {ωn0,ωn1,...,ωn(n1)}时,这个多项式的点值表达式的 y y y集合,就是 h ( x ) h(x) h(x)系数集合的 n n n

还记得FFT算法流程的这一步吗:
A[i]=A0[i]+w*A1[i]A[i+n/2]=A0[i]-w*A1[i]

将它变为这样就是IFFT了:
A[i]=A0[i]-w*A1[i]A[i+n/2]=A0[i]+w*A1[i]

现在知道为什么代码要一起展示了,因为只需在FFT函数中加一个参数就可以实现逆变换了。

代码

看下面。

代码

大整数乘法一般的时间复杂度是 O ( n 2 ) O(n^2) O(n2),把两个大整数看成两个多项式的系数,求它们的乘积,最后处理进位,就可以做到 O ( n log ⁡ n ) O(n\log n) O(nlogn)了。

递归版

据说递归版要爆空间(我这道题和洛谷上的那道都用递归的过了),所以大家都写非递归版。

#include<cmath>
#include<cstdio>
#include<cstring>

#define LOG 17
#define MAXN (1<<LOG)
int res[MAXN+5];
char x[MAXN+5],y[MAXN+5];

struct Complex{
    double x,y;//x+yi
    Complex(){x=y=0;}
    Complex(double a,double b):x(a),y(b){}
    Complex operator + (Complex b)const{return Complex(x+b.x,y+b.y);}
    Complex operator - (Complex b)const{return Complex(x-b.x,y-b.y);}
    Complex operator * (Complex b)const{return Complex(x*b.x-y*b.y,y*b.x+x*b.y);}
}A0[MAXN+5],B0[MAXN+5];
//复数类,好像C++有一个叫complex的库

#define PI acos(-1)
void FFT(Complex *A,int len,int sign){
    if(len==1)
        return;
    int mid=len>>1;
    Complex A1[mid+5],A2[mid+5];//按奇偶分离系数到A1和A2
    for(int i=0;i<mid;i++)
        A1[i]=A[i<<1],A2[i]=A[i<<1|1];
    FFT(A1,mid,sign),FFT(A2,mid,sign);
    for(int i=0;i<mid;i++){
        Complex w(cos(2.0*PI/len*i),sign*sin(2.0*PI/len*i));
        //w是len次单位复根的i次方
        //sign=-1是逆变换
        A[i]=A1[i]+w*A2[i],A[i+mid]=A1[i]-w*A2[i];
    }
}

int main(){
    while(~scanf("%s%s",x+1,y+1)){
        int len1=strlen(x+1),len2=strlen(y+1),N=1;
        while(N<len1+len2)
            N<<=1;
        //统一长度为2的幂
        for(int i=0;i<len1;i++)
            A0[i]=Complex(x[len1-i]-'0',0);
        for(int i=0;i<len2;i++)
            B0[i]=Complex(y[len2-i]-'0',0);
        //倒着存可以避免很多恶心的问题
        FFT(A0,N,1);
        FFT(B0,N,1);
        for(int i=0;i<N;i++)
            A0[i]=A0[i]*B0[i];
        //得到乘积多项式的点值表达式
        FFT(A0,N,-1);
        for(int i=0;i<N;i++)
            res[i+1]=int(A0[i].x/N+0.5);
        //注意最后统一除以n,为避免精度误差,+0.5再取整
        for(int i=1;i<=N;i++)
            res[i+1]+=res[i]/10,res[i]%=10;
        while(N>1&&res[N]==0)
            N--;
        //处理进位和前导零
        for(int i=N;i>=1;i--)
            putchar(res[i]+'0');
        putchar('\n');
        memset(A0,0,sizeof A0);
        memset(B0,0,sizeof B0);
        memset(res,0,sizeof res);
    }
}

非递归版

建议写这种。

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;

#define LOG 17
#define MAXN (1<<LOG)
int res[MAXN+5];
char x[MAXN+5],y[MAXN+5];

struct Complex{
    double x,y;//x+yi
    Complex(){x=y=0;}
    Complex(double a,double b):x(a),y(b){}
    Complex operator + (Complex b)const{return Complex(x+b.x,y+b.y);}
    Complex operator - (Complex b)const{return Complex(x-b.x,y-b.y);}
    Complex operator * (Complex b)const{return Complex(x*b.x-y*b.y,y*b.x+x*b.y);}
}A0[MAXN+5],B0[MAXN+5];

#define PI acos(-1)
int R[MAXN+5];
void FFT(Complex *A,int len,int sign){
    for(int i=0;i<len;i++)
        if(i<R[i])//不加的话会swap两次,最后什么都没有发生
            swap(A[i],A[R[i]]);
    //调整成最后的顺序
    for(int mid=1;mid<len;mid<<=1)//枚举当前要处理的区间长度的一半
        for(int i=0;i<len;i+=mid<<1)//每2*mid的长度为一段,分开处理
            for(int j=0;j<mid;j++){//处理当前2*mid长度的区间
                Complex w(cos(PI/mid*j),sign*sin(PI/mid*j));//2.0被mid消掉了
                Complex tmp1=A[i+j],tmp2=A[i+j+mid];
                A[i+j]=tmp1+w*tmp2;
                A[i+j+mid]=tmp1-w*tmp2;
            }
}

int main(){
    while(~scanf("%s%s",x+1,y+1)){
        int len1=strlen(x+1),len2=strlen(y+1),N=1,Log=0;
        while(N<len1+len2)
            N<<=1,Log++;
        for(int i=0;i<len1;i++)
            A0[i]=Complex(x[len1-i]-'0',0);
        for(int i=0;i<len2;i++)
            B0[i]=Complex(y[len2-i]-'0',0);
        for(int i=0;i<N;i++)
            R[i]=(R[i>>1]>>1)|((i&1)<<(Log-1));
        //dp得到R数组
        FFT(A0,N,1),FFT(B0,N,1);
        for(int i=0;i<N;i++)
            A0[i]=A0[i]*B0[i];
        FFT(A0,N,-1);
        for(int i=0;i<N;i++)
            res[i+1]=int(A0[i].x/N+0.5);
        for(int i=1;i<=N;i++)
            res[i+1]+=res[i]/10,res[i]%=10;
        while(N>1&&res[N]==0)
            N--;
        for(int i=N;i>=1;i--)
            putchar(res[i]+'0');
        putchar('\n');
        memset(A0,0,sizeof A0);
        memset(B0,0,sizeof B0);
        memset(res,0,sizeof res);
    }
}

后记

感谢这些大佬:

CXH大佬
十分简明易懂的FFT(快速傅里叶变换)
小学生都能看懂的FFT!!!
快速傅里叶变换(FFT)详解

  • 7
    点赞
  • 17
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
FFT(Fast Fourier Transform)是一种基于蝶形算法快速傅里叶变换方法,用于将时域信号转换为频域表示。它是一种高效计算离散傅里叶变换(DFT)的算法。 蝶形算法FFT的关键部分,它通过分治策略将DFT问题分解为规模更小的子问题。具体而言,蝶形算法将DFT的计算任务分为奇偶两组,每组包含一半的信号点。然后,它通过递归调用自身来分别计算这两组的DFT。通过这种方式,FFT可以显著减少计算复杂性。 蝶形算法的核心是蝶形运算。它是一种简单的计算,涉及两个复数的乘法和加法。蝶形运算需要两个输入值(称为蝶形输入),分别乘以一个旋转系数(通常通过旋转因子计算得到),然后相加得到两个输出值(称为蝶形输出)。这个运算重复进行,直到计算出所有输出值。 FFT的过程可以用一个蝶形图来表示,其中每个节点代表一次蝶形运算。蝶形图按照二叉树的形式组织,上层节点代表更高频率的DFT,下层节点代表更低频率的DFT。通过在蝶形图上按照特定顺序进行计算,FFT可以以线性时间复杂度完成整个DFT过程。 总结起来,FFT快速傅里叶变换的核心是蝶形算法,它通过分治策略将DFT问题分解为规模更小的子问题,然后通过蝶形运算来计算这些子问题。FFT以其高效的计算能力,广泛应用于信号处理、图像处理、通信等领域。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值