FFT算法的优化——提高FFT算法一倍的时空效率

前言

由于本人学习FFT时第一份代码并不是从普通的FFT学起的,而是直接从一个经过优化的版本开始学。

这是一个利用到共轭复数性质的优化,出自myy的论文,但在网上缺乏关于这个优化的文章,有写的都是大神,没有很细致的讲解。因此我决定自己总结一下。

本文是在默认读者了解并掌握了FFT算法的情况下写的。没有学过FFT的同学可以先去学一下。

FFT的基本运算

ωN\omega_N 表示 e2πiNe^{\frac {2\pi i} {N}}
对于给定的多项式函数

f(x)=i=0naixif(x)=\sum_{i=0}^{n}a_ix^i

XXffDFTDFT ,则有

Xi=j=0N1ajωNij,ai=1Nj=0N1XjωNijX_i=\sum_{j=0}^{N-1}a_j \omega_{N}^{ij},a_i=\frac{1}{N}\sum_{j=0}^{N-1}X_j \omega_{N}^{-ij}

以及我们知道的卷积定理

DFT(fg)=DFT(f)DFT(g)DFT(f\cdot g)=DFT(f)\cdot DFT(g)DFT(f+g)=DFT(f)+DFT(g)DFT(f+g)=DFT(f)+DFT(g)

还有各种各样的复数运算法则,都是本文的前置技能。
其中复数运算法则中有一条性质:
z1z2=zˉ1zˉ2 \overline {z_1z_2}=\bar z_1 \cdot \bar z_2
其中 aˉ\bar a 表示 aa 的共轭复数。
这是我们优化FFT算法需要用到的重要性质。

如何优化?

考虑到 DFTDFT 的对称性,我们来比较一下 XiX_i 以及 XNiX_{N-i}

Xi=j=0N1ajωNijX_i=\sum_{j=0}^{N-1}a_j \omega_{N}^{ij}XNi=j=0N1ajωN(Ni)j=j=0N1ajωNij=XiX_{N-i}=\sum_{j=0}^{N-1}a_j \omega_{N}^{(N-i)j}=\sum_{j=0}^{N-1}a_j \omega_{N}^{-ij}=\overline {X_i}

咦??求出 XiX_i 就知道 XNiX_{N-i} 了?那用 NN 的长度做岂不是很浪费??
为什么 Xni=XiX_{n-i}=X_i 呢?实际上这是由于我们的多项式函数的系数并没有虚部造成的。我们换一个系数有虚部的函数看看。

f(x)=i=0n(ai+ibi)xif(x)=\sum_{i=0}^{n}(a_i+ib_i)x^i

这里可能变量名重了= =,考虑到惯用的表达方式,我们约定此文以下部分中上下标的 ii 表示变量,其余位置表示虚数单位。

Xi=j=0N1(aj+ibj)ωNijX_i=\sum_{j=0}^{N-1}(a_j+ib_j)\omega_{N}^{ij}XNi=j=0N1(aj+ibj)ωN(Ni)j=j=0N1(aj+ibj)ωNij=XiX_{N-i}=\sum_{j=0}^{N-1}(a_j+ib_j) \omega_{N}^{(N-i)j}=\sum_{j=0}^{N-1}(a_j+ib_j) \omega_{N}^{-ij} \sout{= \overline {X_i}}

这时候上面的等号就不成立了。但我们依然可以使用共轭复数的性质对比 XNiX_{N-i}XiX_i
XNi=j=0N1(aj+ibj)ωNij=j=0N1(ajibj)ωNij\overline {X_{N-i}}=\sum_{j=0}^{N-1}\overline{(a_j+ib_j)} \omega_{N}^{ij} =\sum_{j=0}^{N-1}{(a_j-ib_j)} \omega_{N}^{ij}

这个结果相当美妙。我们将 XiX_iXNi\overline {X_{N-i}} 相加减,可以得到:
Xi+XNi=2j=0N1ajωNijX_i+\overline {X_{N-i}}=2\sum_{j=0}^{N-1}a_j\omega_{N}^{ij}XiXNi=2ij=0N1bjωNijX_i-\overline {X_{N-i}}=2i\sum_{j=0}^{N-1}b_j\omega_{N}^{ij}

也就是说,我们只要对 ff 做了 DFTDFT ,就可以在 O(n)O(n) 时间内得到其实部和虚部分别做 DFTDFT 的结果。

因此对于任意一个系数为实数的多项式

f(x)=i=0naixif(x)=\sum_{i=0}^{n}a_ix^i

我们可以将其改写为
f(x)=i=0n2(a2i+ia2i+1)xif(x)=\sum_{i=0}^{ \lfloor \frac{n} {2}\rfloor}(a_{2i}+ia_{2i+1})x^i

以提高我们的时空效率。
当然也可以不按这种方法拆开,但这种拆分方式有利于之后的多项式乘法优化。

在多项式乘法上的应用

FFTFFT 最大的应用是解决多项式乘法的问题,我们将这个优化回归应用到多项式乘法上。

普通的多项式乘法

对于两个函数

f(x)=i=0naixi,g(x)=i=0mbixif(x)=\sum_{i=0}^{n}a_ix^i,g(x)=\sum_{i=0}^{m}b_ix^i

我们当然可以用上面的方式把 ffgg 的奇偶项的 DFTDFT 分别乘起来再分别做 IDFTIDFT 得到 fgf\cdot g ,但这样我们的优化就失去了实用价值。我们考虑怎样使用更少的 IDFTIDFT 次数得到原函数。

我们已经将 f,gf,g 改写成了以下形式:
f(x)=i=0n2(a2i+ia2i+1)xi,g(x)=i=0m2(b2i+ib2i+1)xif(x)=\sum_{i=0}^{ \lfloor \frac{n} {2}\rfloor}(a_{2i}+ia_{2i+1})x^i,g(x)=\sum_{i=0}^{ \lfloor \frac{m} {2}\rfloor}(b_{2i}+ib_{2i+1})x^i

我们的目标函数为
h(x)=i=0n+m2j=0i(a2jb2(ij)+a2j+1b2(ij)1+i(a2jb2(ij)+1+a2j+1b2(ij)))xih(x)=\sum_{i=0}^{ \lfloor \frac{n+m} {2}\rfloor}\sum_{j=0}^i(a_{2j}b_{2(i-j)}+a_{2j+1}b_{2(i-j)-1}+i(a_{2j}b_{2(i-j)+1}+a_{2j+1}b_{2(i-j)}))x^i

X,Y,ZX,Y,Zf,g,hf,g,hDFTDFT
Xi=j=0N1(a2j+ia2j+1)ωNij,Yi=j=0N1(b2j+ib2j+1)ωNijX_i=\sum_{j=0}^{N-1}(a_{2j}+ia_{2j+1})\omega_N^{ij},Y_i =\sum_{j=0}^{N-1}(b_{2j}+ib_{2j+1})\omega_N^{ij}Zi=j=0N1k=0j(a2kb2(jk)+a2k+1b2(jk)1+i(a2kb2(jk)+1+a2k+1b2(jk)))ωNijZ_i=\sum_{j=0}^{N-1}\sum_{k=0}^j(a_{2k}b_{2(j-k)}+a_{2k+1}b_{2(j-k)-1}+i(a_{2k}b_{2(j-k)+1}+a_{2k+1}b_{2(j-k)}))\omega_N^{ij}XiYi=j=0N1k=0j(a2kb2(jk)a2k+1b2(jk)+1+i(a2kb2(jk)+1+a2k+1b2(jk)))ωNijX_iY_i=\sum_{j=0}^{N-1}\sum_{k=0}^j(a_{2k}b_{2(j-k)}-a_{2k+1}b_{2(j-k)+1}+i(a_{2k}b_{2(j-k)+1}+a_{2k+1}b_{2(j-k)}))\omega_N^{ij}

对比括号里的四项我们发现只有一项不同,我们将两式相减:
ZiXiYi=j=0N1k=0j(a2k+1b2(jk)1+a2k+1b2(jk)+1)ωNij=j=0N1k=0ja2k+1b2(jk)1ωNij+j=0N1k=0ja2k+1b2(jk)+1ωNij=j=0N1k=0ja2k+1b2(jk)+1ωNi(j+1)+j=0N1k=0ja2k+1b2(jk)+1ωNij=(1+ωN i)j=0N1k=0ja2k+1b2(jk)+1ωNij=(1+ωN i)(XiXNi)(YiYNi)4\begin{aligned} Z_i-X_iY_i&=\sum_{j=0}^{N-1}\sum_{k=0}^j(a_{2k+1}b_{2(j-k)-1}+a_{2k+1}b_{2(j-k)+1})\omega_N^{ij} \\&=\sum_{j=0}^{N-1}\sum_{k=0}^ja_{2k+1}b_{2(j-k)-1}\omega_N^{ij}+\sum_{j=0}^{N-1}\sum_{k=0}^ja_{2k+1}b_{2(j-k)+1}\omega_N^{ij}\\ &=\sum_{j=0}^{N-1}\sum_{k=0}^ja_{2k+1}b_{2(j-k)+1}\omega_N^{i(j+1)}+\sum_{j=0}^{N-1}\sum_{k=0}^ja_{2k+1}b_{2(j-k)+1}\omega_N^{ij} \\&=(1+\omega_N^{\ i})\sum_{j=0}^{N-1}\sum_{k=0}^ja_{2k+1}b_{2(j-k)+1}\omega_N^{ij} \\ &=-\frac{(1+\omega_N^{\ i})(X_i-\overline{X_{N-i}})(Y_i-\overline{Y_{N-i}})} {4}\end{aligned}

这里利用到了上面的结论
XiXNi=2ij=0N1a2j+1ωNijX_i-\overline {X_{N-i}}=2i\sum_{j=0}^{N-1}a_{2j+1}\omega_{N}^{ij}

因此

Zi=4XiYi(1+ωN i)(XiXNi)(YiYNi)4Z_i=\frac{4X_iY_i-(1+\omega_N^{\ i})(X_i-\overline{X_{N-i}})(Y_i-\overline{Y_{N-i}})} {4}

这样我们做多项式乘法就可以把时空效率都提高一倍了,从原先的3次FFT优化到1.5次FFT,并且该做法的精度上优于普通的FFT。

对任意模多项式乘法的优化

我们知道任意模多项式乘法有8次和7次 DFTDFT 的做法,直接套上上面的做法就可以优化到4次或3.5次。此外还有一种更好写的4次做法:
aa 拆成 A1,A2A_1,A_2bb 拆成 B1,B2B_1,B_2 ,并构造
f(x)=i=0n(A1,i+iA2,i)xi,g(x)=i=0m(B1,i+iB2,i)xif(x)=\sum_{i=0}^{n} (A_{1,i}+iA_{2,i})x^i,g(x)=\sum_{i=0}^{m} (B_{1,i}+iB_{2,i})x^i
分别做 DFTDFT 之后算出 A1B1,A1B2,A2B1,A2B2A1B1,A1B2,A2B1,A2B2DFTDFT 分别存进两个数组的虚部和实部最后合并。

作为黑科技还是很值得一用的。

(第一次写blog,写得不好,可能还有写挂的地方,求轻喷)

展开阅读全文

没有更多推荐了,返回首页