【算法】快速傅里叶变换(FFT)初探

【参考】

  证明部分摘自:https://oi.men.ci/fft-notes/

【简介】

  快速傅里叶变换(FFT)是一种可以在O(nlogn)O(nlogn)时间内完成的离散傅里叶变换(DFT)算法,在OI中主要用于加速向量卷积/多项式乘法运算。

【前置技能】

【引入】

  有两个多项式A(x)A(x)B(x)B(x),求C(x)=A(x)B(x)C(x)=A(x)*B(x)

A(x)=i=0n1aixiB(x)=i=0n1bixiC(x)=A(x)B(x)=i=02n2(j=0iajbij)xi\begin{gathered} A(x)=\sum_{i=0}^{n-1}a_ix^i\\ B(x)=\sum_{i=0}^{n-1}b_ix^i\\ C(x)=A(x)*B(x)=\sum_{i=0}^{2n-2}(\sum_{j=0}^i a_j*b_{i-j})x^i\\ \end{gathered}

  用朴素的做法复杂度显然是O(n2)O(n^2)的。
  如果再推推?

C(x)=A(x)B(x)=i=02n2(j=0iajbij)xi=i=02n2j=0i(ajxj)(bijxij)=j=02n2(ajxj)i=02n2(bixi)\begin{aligned} C(x)&=A(x)*B(x)\\ &=\sum_{i=0}^{2n-2}(\sum_{j=0}^i a_j*b_{i-j})x^i\\ &=\sum_{i=0}^{2n-2}\sum_{j=0}^i(a_j*x^j)(b_{i-j}*x^{i-j})\\ &=\sum_{j=0}^{2n-2}(a_j*x^j)\sum_{i=0}^{2n-2}(b_i*x^i)\\ \end{aligned}

  式子左边和右边分别是多项式A(x),B(x)A(x),B(x)的点值表示,求一个多项式的点值表示和通过点值表示求多项式都可以通过FFT和IDFT做到O(nlogn)O(nlogn),所以我们就可以在O(nlogn)O(nlogn)的时间进行多项式乘法运算了。

【系数表示法】

  设A(x)A(x)表示一个n1n-1次多项式,所有项的系数组成的nn维向量(a0,a1,a2,...,an1)(a_0, a_1, a_2, ..., a_{n-1})唯一确定了一个多项式,这种表示法叫系数表示法。

A(x)=i=0n1aixiA(x)=\sum_{i=0}^{n-1}a_ix^i

【点值表示法】

  如果将一组互不相同的nnxx代入多项式,得到nn个值,组成的nn维向量为(x0,x1,x2,...,xn1),(y0,y1,y2,...,yn1)(x_0, x_1, x_2, ..., x_{n-1}),(y_0, y_1, y_2, ..., y_{n-1}),唯一确定了一个多项式(证明在下方),这种表示法叫点值表示法。

yi=A(xi)=j=0n1ajxijy_i=A(x_i)=\sum_{j=0}^{n-1} a_jx_i^j

  一个n1n-1次多项式在nn个不同点的取值唯一确定一个多项式的证明:
  若命题不成立,则存在两个不同的n1n-1次多项式A(x),B(x)A(x),B(x)满足i[0,n1],A(xi)=B(xi)\forall i\subseteq [0,n-1],A(x_i)=B(x_i)
  令C(x)=A(x)B(x)C(x)=A(x)-B(x),则C(x)C(x)也是一个n1n-1次多项式,且满足i[0,n1],C(xi)=0\forall i\subseteq [0,n-1],C(x_i)=0
  即C(x)C(x)nn个根,与代数基本定理(一个n1n-1次多项式在复数域上至多只有n1n-1个根)相矛盾,故原命题成立。

【卷积】

  (fg)(n)(f*g)(n)称为f,gf,g的卷积,其离散的定义为:

(fg)(n)=f(τ)g(nτ)(f*g)(n)=\sum_{-\infty }^{\infty}f(\tau)g(n-\tau)

  特征为τ+(nτ)=n\tau+(n-\tau)=n

【复数、复平面】

  相关知识见:http://www.ruanyifeng.com/blog/2012/09/imaginary_number.html

【单位根】

  在复平面上以原点为起点,单位圆的nn等分点为终点,作nn个向量,所得的幅角为正且最小的向量对应的复数ωn\omega _n,称为nn次单位根,则其余的n1n-1个向量对应的复数分别为ωn2,ωn3...,ωnn\omega _n^2,\omega _n^3...,\omega _n^n,显然ωn0=ωnn=1\omega _n^0=\omega _n^n=1
  因为单位根的幅角为1n\frac{1}{n},所以有ωnk=cosk2πn+isink2πn\omega _n^k=cos k\frac{2\pi}{n}+i sin k\frac{2\pi}{n}

【单位根的性质】

  ①ω2n2k=ωnk\omega _{2n}^{2k}=\omega _n^k
  ②ωnk+n2=ωnk\omega _n^{k+\frac{n}{2}}=-\omega _n^k
  结合单位圆较易理解,证明略

【快速傅里叶变换(FFT)】

  点值向量(A(ωn0),A(ωn1),A(ωn2),...,A(ωnn1))(A(\omega_n^0),A(\omega_n^1),A(\omega_n^2),...,A(\omega_n^{n-1}))被称为其系数向量(a0,a1,a2,...,an1)(a_0, a_1, a_2, ..., a_{n-1})的离散傅里叶变换。
  把A(x)A(x)按奇偶次项分组。
  令

A1(x)=(a0+a2x+a4x2+...+an2xn21)A2(x)=(a1+a3x+a5x2+...+an1xn21)A(x)=A1(x2)+xA2(x2)\begin{gathered} A_1(x)=(a_0+a_2x+a_4x^2+...+a_{n-2}x^{\frac{n}{2}-1})\\ A_2(x)=(a_1+a_3x+a_5x^2+...+a_{n-1}x^{\frac{n}{2}-1})\\ A(x)=A_1(x^2)+xA_2(x^2)\\ \end{gathered}

  若k<n2k<\frac{n}{2},求A(ωnk)A(\omega _n^k)

A(ωnk)=A1(ωn2k)+ωnkA2(ωn2k)=A1(ωn2k)+ωnkA2(ωn2k)\begin{aligned} A(\omega _n^k)&=A_1(\omega _n^{2k})+\omega _n^kA_2(\omega _n^{2k})\\ &=A_1(\omega _{\frac{n}{2}}^k)+\omega _n^kA_2(\omega _{\frac{n}{2}}^k) \end{aligned}

  求A(ωnk+n2)A(\omega _n^{k+\frac{n}{2}})

A(ωnk+n2)=A1(ωn2k+n)+ωnk+n2A2(ωn2k+n)=A1(ωn2k×ωnn)ωnkA2(ωn2k×ωnn)=A1(ωn2k)ωnkA2(ωn2k)=A1(ωn2k)ωnkA2(ωn2k)\begin{aligned} A(\omega_n^{k+\frac{n}{2}})&=A_1(\omega _n^{2k+n})+\omega _n^{k+\frac{n}{2}}A_2(\omega _n^{2k+n})\\ &=A_1(\omega _n^{2k}\times \omega_n^n)-\omega_n^kA_2(\omega_n^{2k}\times \omega_n^n)\\ &=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k})\\ &=A_1(\omega_{\frac{n}{2}}^k)-\omega_n^kA_2(\omega_{\frac{n}{2}}^k) \end{aligned}

  根据上方的推导,我们显然可以使用一个分治的算法来求得所有的点值,效率O(nlogn)O(nlogn)
  这就是Cooley-Tukey算法。

【傅里叶逆变换(IDFT)】

  将点值表示的多项式转换为系数表示,可以使用傅里叶逆变换。
  设(y0,y1,y2,...,yn1)(y_0, y_1, y_2, ..., y_{n-1})(a0,a1,a2,...,an1)(a_0, a_1, a_2, ..., a_{n-1})的离散傅里叶变换,设另一个向量(c0,c1,c2,...,cn1)(c_0, c_1, c_2, ..., c_{n-1})满足

ck=i=0n1yi(ωnk)i=i=0n1j=0n1(aj(ωni)j)(ωnk)i=i=0n1j=0n1(aj(ωnj)i(ωnk)i)=i=0n1j=0n1aj(ωnjk)i=j=0n1aj(i=0n1(ωnjk)i)\begin{aligned} c_k&=\sum_{i=0}^{n-1}y_i(\omega _n^{-k})^i\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}(a_j(\omega _n^i)^j)(\omega_n^{-k})^i\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}(a_j(\omega _n^j)^i(\omega_n^{-k})^i)\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^{j-k})^i\\ &=\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i)\\ \end{aligned}

  设

S(ωnk)=1+ωnk+(ωnk)2+...+(ωnk)n1S(\omega_n^k)=1+\omega_n^k+(\omega_n^k)^2+...+(\omega_n^k)^{n-1}

  k=0k=0时,显然S(ωnk)=nS(\omega_n^k)=n
  k0k\neq 0时,则等式两边乘上ωnk\omega_n^k

ωnkS(ωnk)=ωnk+(ωnk)2+(ωnk)3+...+(ωnk)n\omega_n^kS(\omega_n^k)=\omega_n^k+(\omega_n^k)^2+(\omega_n^k)^3+...+(\omega_n^k)^{n}

  两式相减,得

ωnkS(ωnk)S(ωnk)=(ωnk)n1S(ωnk)=(ωnk)n1ωnk1\begin{gathered} \omega_n^kS(\omega_n^k)-S(\omega_n^k)=(\omega_n^k)^n-1\\ S(\omega_n^k)=\frac{(\omega_n^k)^n-1}{\omega_n^k-1} \end{gathered}

  分子为0,分母不为0,故

S(ωnk)=0S(\omega_n^k)=0

  考虑上式

ck=j=0n1aj(i=0n1(ωnjk)i)=j=0n1ajS(ωnjk)\begin{aligned} c_k&=\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i)\\ &=\sum_{j=0}^{n-1}a_jS(\omega_n^{j-k})\\ \end{aligned}

  当j=kj=k时,S(ωnjk)=nS(\omega_n^{j-k})=n,否则为0,故

ci=naiai=1nci\begin{gathered} c_i=na_i\\ a_i=\frac{1}{n}c_i \end{gathered}

  所以使用单位根的倒数作为xx,做一次FFT,对得到的每个数除以nn,即可得到点值表示对应的系数表示。

【实现】

【迭代FFT】

  递归版FFT速度十分慢,所以一般使用迭代版的FFT。
  观察规律可知FFT分治到终止时,某个编号的初始下标是终止下标二进制位翻转后的值,然后就可以模拟分治进行迭代FFT了。

【蝴蝶操作】

  运用一个临时变量省去一个中转数组。
  若此时a(k)a(k)A1(ωn2k)A_1(\omega_{\frac{n}{2}}^k)a(k+n2)a(k+\frac{n}{2})A2(ωn2k)A_2(\omega_{\frac{n}{2}}^k),用a(k)+ωnk×a(n2+k)a(k)+\omega_n^k\times a(\frac{n}{2}+k)更新a(k)a(k),用a(k)ωnk×a(n2+k)a(k)-\omega_n^k\times a(\frac{n}{2}+k)更新a(n2+k)a(\frac{n}{2}+k)时,用一个临时变量tt进行以下操作。

t=ωnk×a(n2+k)a(n2+k)=a(k)ta(k)=a(k)+t\begin{gathered} t=\omega_n^k\times a(\frac{n}{2}+k)\\ a(\frac{n}{2}+k)=a(k)-t\\ a(k)=a(k)+t \end{gathered}

【注意】

  最后求系数取整时为减少误差一般+0.5以四舍五入

【流程】

  将数从初始位置换到终止位置。
  模拟分治进行迭代FFT。