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

前言

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

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

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

FFT的基本运算

记 ${\omega }_N$ 表示 $e^{\frac{2\pi i}{N}}$ 。
对于给定的多项式函数

$$f(x)=\sum _{i=0}^n{a}_i{x}^i$$

记 $X$ 为 $f$ 的 $DFT$ ，则有

$$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}$$

以及我们知道的卷积定理

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

还有各种各样的复数运算法则，都是本文的前置技能。
其中复数运算法则中有一条性质：
$$\overline{z_1{z}_2}={\bar{z}}_1\cdot {\bar{z}}_2$$
其中 $\bar{a}$ 表示 $a$ 的共轭复数。
这是我们优化FFT算法需要用到的重要性质。

如何优化？

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

$$X_i=\sum _{j=0}^{N-1}{a}_j{\omega }_N^{ij}$$$$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}$$

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

$$f(x)=\sum _{i=0}^n(a_i+i{b}_i)x^i$$

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

$$X_i=\sum _{j=0}^{N-1}(a_j+i{b}_j){\omega }_N^{ij}$$$$X_{N-i}=\sum _{j=0}^{N-1}(a_j+i{b}_j){\omega }_N^{(N-i)j}=\sum _{j=0}^{N-1}(a_j+i{b}_j){\omega }_N^{-ij}\sout{=\overline{X_i}}$$

这时候上面的等号就不成立了。但我们依然可以使用共轭复数的性质对比 $X_{N-i}$ 和 $X_i$ ：
$$\overline{X_{N-i}}=\sum _{j=0}^{N-1}\overline{(a_j+i{b}_j)}{\omega }_N^{ij}=\sum _{j=0}^{N-1}(a_j-i{b}_j){\omega }_N^{ij}$$

这个结果相当美妙。我们将 $X_i$ 与 $\overline{X_{N-i}}$ 相加减，可以得到：
$$X_i+\overline{X_{N-i}}=2\sum _{j=0}^{N-1}{a}_j{\omega }_N^{ij}$$$$X_i-\overline{X_{N-i}}=2i\sum _{j=0}^{N-1}{b}_j{\omega }_N^{ij}$$

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

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

$$f(x)=\sum _{i=0}^n{a}_i{x}^i$$

我们可以将其改写为
$$f(x)=\sum _{i=0}^{\lfloor \frac{n}{2}\rfloor }(a_{2i}+i{a}_{2i+1})x^i$$

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

在多项式乘法上的应用

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

普通的多项式乘法

对于两个函数

$$f(x)=\sum _{i=0}^n{a}_i{x}^i,g(x)=\sum _{i=0}^m{b}_i{x}^i$$

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

我们已经将 $f,g$ 改写成了以下形式：
$$f(x)=\sum _{i=0}^{\lfloor \frac{n}{2}\rfloor }(a_{2i}+i{a}_{2i+1})x^i,g(x)=\sum _{i=0}^{\lfloor \frac{m}{2}\rfloor }(b_{2i}+i{b}_{2i+1})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$$

记 $X,Y,Z$ 为 $f,g,h$ 的 $DFT$ 。
$$X_i=\sum _{j=0}^{N-1}(a_{2j}+i{a}_{2j+1}){\omega }_N^{ij},Y_i=\sum _{j=0}^{N-1}(b_{2j}+i{b}_{2j+1}){\omega }_N^{ij}$$$$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}$$$$X_i{Y}_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}$$

对比括号里的四项我们发现只有一项不同，我们将两式相减：
$$\begin{array}{rl}{Z}_i-X_i{Y}_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}^j{a}_{2k+1}{b}_{2(j-k)-1}{\omega }_N^{ij}+\sum _{j=0}^{N-1}\sum _{k=0}^j{a}_{2k+1}{b}_{2(j-k)+1}{\omega }_N^{ij}\\ & =\sum _{j=0}^{N-1}\sum _{k=0}^j{a}_{2k+1}{b}_{2(j-k)+1}{\omega }_N^{i(j+1)}+\sum _{j=0}^{N-1}\sum _{k=0}^j{a}_{2k+1}{b}_{2(j-k)+1}{\omega }_N^{ij}\\ & =(1+{\omega }_N^i)\sum _{j=0}^{N-1}\sum _{k=0}^j{a}_{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{array}$$

这里利用到了上面的结论
$$X_i-\overline{X_{N-i}}=2i\sum _{j=0}^{N-1}{a}_{2j+1}{\omega }_N^{ij}$$

因此

$$Z_i=\frac{4X_i{Y}_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次 $DFT$ 的做法，直接套上上面的做法就可以优化到4次或3.5次。此外还有一种更好写的4次做法：
将 $a$ 拆成 $A_1,A_2$ ， $b$ 拆成 $B_1,B_2$ ，并构造
$$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$$
分别做 $DFT$ 之后算出 $A1B1,A1B2,A2B1,A2B2$ 的$DFT$ 分别存进两个数组的虚部和实部最后合并。

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

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