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

fft算法

前言

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

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

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

FFT的基本运算

记 $ω_{N} \omega_N$ ωN 表示 $e^{\frac{2 π i}{N}} e^{\frac {2\pi i} {N}}$ eN2πi 。

对于给定的多项式函数

$f (x) = \sum_{i = 0}^{n} a_{i} x^{i} f(x)=\sum_{i=0}^{n}a_ix^i$ f(x)=i=0∑naixi

记 $X X$ X 为 $f f$ f 的 $D F T DFT$ DFT ，则有

$X_{i} = \sum_{j = 0}^{N - 1} a_{j} ω_{N}^{i j}, a_{i} = \frac{1}{N} \sum_{j = 0}^{N - 1} X_{j} ω_{N}^{- i j} X_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}$ Xi=j=0∑N−1ajωNij,ai=N1j=0∑N−1XjωN−ij

以及我们知道的卷积定理

$D F T (f \cdot g) = D F T (f) \cdot D F T (g) DFT(f\cdot g)=DFT(f)\cdot DFT(g)$ DFT(f⋅g)=DFT(f)⋅DFT(g) $D F T (f + g) = D F T (f) + D F T (g) DFT(f+g)=DFT(f)+DFT(g)$ DFT(f+g)=DFT(f)+DFT(g)

还有各种各样的复数运算法则，都是本文的前置技能。

其中复数运算法则中有一条性质：

$\overline{z_{1} z_{2}} = {\overset{ˉ}{z}}_{1} \cdot {\overset{ˉ}{z}}_{2} \overline {z_1z_2}=\bar z_1 \cdot \bar z_2$ z1z2=zˉ1⋅zˉ2

其中 $\overset{ˉ}{a} \bar a$ aˉ 表示 $a a$ a 的共轭复数。

这是我们优化FFT算法需要用到的重要性质。

如何优化？

考虑到 $D F T DFT$ DFT 的对称性，我们来比较一下 $X_{i} X_i$ Xi 以及 $X_{N - i} X_{N-i}$ XN−i ：

$X_{i} = \sum_{j = 0}^{N - 1} a_{j} ω_{N}^{i j} X_i=\sum_{j=0}^{N-1}a_j \omega_{N}^{ij}$ Xi=j=0∑N−1ajωNij $X_{N - i} = \sum_{j = 0}^{N - 1} a_{j} ω_{N}^{(N - i) j} = \sum_{j = 0}^{N - 1} a_{j} ω_{N}^{- i j} = \overline{X_{i}} X_{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}$ XN−i=j=0∑N−1ajωN(N−i)j=j=0∑N−1ajωN−ij=Xi

咦？？求出 $X_{i} X_i$ Xi 就知道 $X_{N - i} X_{N-i}$ XN−i 了？那用 $N N$ N 的长度做岂不是很浪费？？

为什么 $X_{n - i} = X_{i} X_{n-i}=X_i$ Xn−i=Xi 呢？实际上这是由于我们的多项式函数的系数并没有虚部造成的。我们换一个系数有虚部的函数看看。

设

$f (x) = \sum_{i = 0}^{n} (a_{i} + i b_{i}) x^{i} f(x)=\sum_{i=0}^{n}(a_i+ib_i)x^i$ f(x)=i=0∑n(ai+ibi)xi

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

$X_{i} = \sum_{j = 0}^{N - 1} (a_{j} + i b_{j}) ω_{N}^{i j} X_i=\sum_{j=0}^{N-1}(a_j+ib_j)\omega_{N}^{ij}$ Xi=j=0∑N−1(aj+ibj)ωNij $X_{N - i} = \sum_{j = 0}^{N - 1} (a_{j} + i b_{j}) ω_{N}^{(N - i) j} = \sum_{j = 0}^{N - 1} (a_{j} + i b_{j}) ω_{N}^{- i j} = \overline{X_{i}} X_{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}}$ XN−i=j=0∑N−1(aj+ibj)ωN(N−i)j=j=0∑N−1(aj+ibj)ωN−ij=Xi

这时候上面的等号就不成立了。但我们依然可以使用共轭复数的性质对比 $X_{N - i} X_{N-i}$ XN−i 和 $X_{i} X_i$ Xi ：

$\overline{X_{N - i}} = \sum_{j = 0}^{N - 1} \overline{(a_{j} + i b_{j})} ω_{N}^{i j} = \sum_{j = 0}^{N - 1} (a_{j} - i b_{j}) ω_{N}^{i j} \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}$ XN−i=j=0∑N−1(aj+ibj)ωNij=j=0∑N−1(aj−ibj)ωNij

这个结果相当美妙。我们将 $X_{i} X_i$ Xi 与 $\overline{X_{N - i}} \overline {X_{N-i}}$ XN−i 相加减，可以得到：

$X_{i} + \overline{X_{N - i}} = 2 \sum_{j = 0}^{N - 1} a_{j} ω_{N}^{i j} X_i+\overline {X_{N-i}}=2\sum_{j=0}^{N-1}a_j\omega_{N}^{ij}$ Xi+XN−i=2j=0∑N−1ajωNij $X_{i} - \overline{X_{N - i}} = 2 i \sum_{j = 0}^{N - 1} b_{j} ω_{N}^{i j} X_i-\overline {X_{N-i}}=2i\sum_{j=0}^{N-1}b_j\omega_{N}^{ij}$ Xi−XN−i=2ij=0∑N−1bjωNij

也就是说，我们只要对 $f f$ f 做了 $D F T DFT$ DFT ，就可以在 $O (n) O(n)$ O(n) 时间内得到其实部和虚部分别做 $D F T DFT$ DFT 的结果。

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

$f (x) = \sum_{i = 0}^{n} a_{i} x^{i} f(x)=\sum_{i=0}^{n}a_ix^i$ f(x)=i=0∑naixi

我们可以将其改写为

$f (x) = \sum_{i = 0}^{⌊ \frac{n}{2} ⌋} (a_{2 i} + i a_{2 i + 1}) x^{i} f(x)=\sum_{i=0}^{ \lfloor \frac{n} {2}\rfloor}(a_{2i}+ia_{2i+1})x^i$ f(x)=i=0∑⌊2n⌋(a2i+ia2i+1)xi

以提高我们的时空效率。

当然也可以不按这种方法拆开，但这种拆分方式有利于之后的多项式乘法优化。

在多项式乘法上的应用

$F F T FFT$ FFT 最大的应用是解决多项式乘法的问题，我们将这个优化回归应用到多项式乘法上。

普通的多项式乘法

对于两个函数

$f (x) = \sum_{i = 0}^{n} a_{i} x^{i}, g (x) = \sum_{i = 0}^{m} b_{i} x^{i} f(x)=\sum_{i=0}^{n}a_ix^i,g(x)=\sum_{i=0}^{m}b_ix^i$ f(x)=i=0∑naixi,g(x)=i=0∑mbixi

我们当然可以用上面的方式把 $f f$ f 和 $g g$ g 的奇偶项的 $D F T DFT$ DFT 分别乘起来再分别做 $I D F T IDFT$ IDFT 得到 $f \cdot g f\cdot g$ f⋅g ，但这样我们的优化就失去了实用价值。我们考虑怎样使用更少的 $I D F T IDFT$ IDFT 次数得到原函数。

我们已经将 $f, g f,g$ f,g 改写成了以下形式：

$f (x) = \sum_{i = 0}^{⌊ \frac{n}{2} ⌋} (a_{2 i} + i a_{2 i + 1}) x^{i}, g (x) = \sum_{i = 0}^{⌊ \frac{m}{2} ⌋} (b_{2 i} + i b_{2 i + 1}) x^{i} f(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$ f(x)=i=0∑⌊2n⌋(a2i+ia2i+1)xi,g(x)=i=0∑⌊2m⌋(b2i+ib2i+1)xi

我们的目标函数为

$h (x) = \sum_{i = 0}^{⌊ \frac{n + m}{2} ⌋} \sum_{j = 0}^{i} (a_{2 j} b_{2 (i - j)} + a_{2 j + 1} b_{2 (i - j) - 1} + i (a_{2 j} b_{2 (i - j) + 1} + a_{2 j + 1} b_{2 (i - j)})) x^{i} h(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$ h(x)=i=0∑⌊2n+m⌋j=0∑i(a2jb2(i−j)+a2j+1b2(i−j)−1+i(a2jb2(i−j)+1+a2j+1b2(i−j)))xi

记 $X, Y, Z X,Y,Z$ X,Y,Z 为 $f, g, h f,g,h$ f,g,h 的 $D F T DFT$ DFT 。

$X_{i} = \sum_{j = 0}^{N - 1} (a_{2 j} + i a_{2 j + 1}) ω_{N}^{i j}, Y_{i} = \sum_{j = 0}^{N - 1} (b_{2 j} + i b_{2 j + 1}) ω_{N}^{i j} X_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}$ Xi=j=0∑N−1(a2j+ia2j+1)ωNij,Yi=j=0∑N−1(b2j+ib2j+1)ωNij $Z_{i} = \sum_{j = 0}^{N - 1} \sum_{k = 0}^{j} (a_{2 k} b_{2 (j - k)} + a_{2 k + 1} b_{2 (j - k) - 1} + i (a_{2 k} b_{2 (j - k) + 1} + a_{2 k + 1} b_{2 (j - k)})) ω_{N}^{i j} Z_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}$ Zi=j=0∑N−1k=0∑j(a2kb2(j−k)+a2k+1b2(j−k)−1+i(a2kb2(j−k)+1+a2k+1b2(j−k)))ωNij $X_{i} Y_{i} = \sum_{j = 0}^{N - 1} \sum_{k = 0}^{j} (a_{2 k} b_{2 (j - k)} - a_{2 k + 1} b_{2 (j - k) + 1} + i (a_{2 k} b_{2 (j - k) + 1} + a_{2 k + 1} b_{2 (j - k)})) ω_{N}^{i j} X_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}$ XiYi=j=0∑N−1k=0∑j(a2kb2(j−k)−a2k+1b2(j−k)+1+i(a2kb2(j−k)+1+a2k+1b2(j−k)))ωNij

对比括号里的四项我们发现只有一项不同，我们将两式相减：

$\begin{matrix} Z_{i} - X_{i} Y_{i} & = \sum_{j = 0}^{N - 1} \sum_{k = 0}^{j} (a_{2 k + 1} b_{2 (j - k) - 1} + a_{2 k + 1} b_{2 (j - k) + 1}) ω_{N}^{i j} \\ = \sum_{j = 0}^{N - 1} \sum_{k = 0}^{j} a_{2 k + 1} b_{2 (j - k) - 1} ω_{N}^{i j} + \sum_{j = 0}^{N - 1} \sum_{k = 0}^{j} a_{2 k + 1} b_{2 (j - k) + 1} ω_{N}^{i j} \\ = \sum_{j = 0}^{N - 1} \sum_{k = 0}^{j} a_{2 k + 1} b_{2 (j - k) + 1} ω_{N}^{i (j + 1)} + \sum_{j = 0}^{N - 1} \sum_{k = 0}^{j} a_{2 k + 1} b_{2 (j - k) + 1} ω_{N}^{i j} \\ = (1 + ω_{N}^{i}) \sum_{j = 0}^{N - 1} \sum_{k = 0}^{j} a_{2 k + 1} b_{2 (j - k) + 1} ω_{N}^{i j} \\ = - \frac{(1 + ω_{N}^{i}) (X_{i} - \overline{X_{N - i}}) (Y_{i} - \overline{Y_{N - i}})}{4} \end{matrix} \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}$ Zi−XiYi=j=0∑N−1k=0∑j(a2k+1b2(j−k)−1+a2k+1b2(j−k)+1)ωNij=j=0∑N−1k=0∑ja2k+1b2(j−k)−1ωNij+j=0∑N−1k=0∑ja2k+1b2(j−k)+1ωNij=j=0∑N−1k=0∑ja2k+1b2(j−k)+1ωNi(j+1)+j=0∑N−1k=0∑ja2k+1b2(j−k)+1ωNij=(1+ωNi)j=0∑N−1k=0∑ja2k+1b2(j−k)+1ωNij=−4(1+ωNi)(Xi−XN−i)(Yi−YN−i)

这里利用到了上面的结论

$X_{i} - \overline{X_{N - i}} = 2 i \sum_{j = 0}^{N - 1} a_{2 j + 1} ω_{N}^{i j} X_i-\overline {X_{N-i}}=2i\sum_{j=0}^{N-1}a_{2j+1}\omega_{N}^{ij}$ Xi−XN−i=2ij=0∑N−1a2j+1ωNij

因此

$Z_{i} = \frac{4 X_{i} Y_{i} - (1 + ω_{N}^{i}) (X_{i} - \overline{X_{N - i}}) (Y_{i} - \overline{Y_{N - i}})}{4} Z_i=\frac{4X_iY_i-(1+\omega_N^{\ i})(X_i-\overline{X_{N-i}})(Y_i-\overline{Y_{N-i}})} {4}$ Zi=44XiYi−(1+ωNi)(Xi−XN−i)(Yi−YN−i)

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

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

我们知道任意模多项式乘法有8次和7次 $D F T DFT$ DFT 的做法，直接套上上面的做法就可以优化到4次或3.5次。此外还有一种更好写的4次做法：

将 $a a$ a 拆成 $A_{1}, A_{2} A_1,A_2$ A1,A2 ， $b b$ b 拆成 $B_{1}, B_{2} B_1,B_2$ B1,B2 ，并构造

$f (x) = \sum_{i = 0}^{n} (A_{1, i} + i A_{2, i}) x^{i}, g (x) = \sum_{i = 0}^{m} (B_{1, i} + i B_{2, i}) x^{i} f(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$ f(x)=i=0∑n(A1,i+iA2,i)xi,g(x)=i=0∑m(B1,i+iB2,i)xi

分别做 $D F T DFT$ DFT 之后算出 $A 1 B 1, A 1 B 2, A 2 B 1, A 2 B 2 A1B1,A1B2,A2B1,A2B2$ A1B1,A1B2,A2B1,A2B2 的 $D F T DFT$ DFT 分别存进两个数组的虚部和实部最后合并。

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

（第一次写blog，写得不好，可能还有写挂的地方，求轻喷）

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

fft算法

前言

FFT的基本运算

如何优化？

在多项式乘法上的应用

普通的多项式乘法

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

相关阅读

栏目导航

推荐阅读

热门阅读