FFT

FFT

Introduction

\(FFT\) 是一种利用神奇操作，在 \(nlog\) 的时间内代替 \(n^2\) 的朴素多项式乘法的算法。由 \(DFT\) 和 \(IDFT\) 两个部分组成。

原理来自于，对于一个 \(n\) 次多项式 \(f(x)\) ，它一般会被表达为：

\[f(x)=\sum_{i=0}^n a_i*x^i \]

这是表达一个多项式最常见的形式。

那么考虑一下，对于任意一个 \(n\) 次多项式 \(f(x)\) ，如果我找出上面的 \(n+1\) 个不同的点会怎么样呢？显然，我们会得到 \(n+1\) 个函数值。

而对于这 \(n+1\) 个函数值，我们可以和函数本身构成一个含 \(n+1\) 个 \(n+1\) 元一次方程的方程组，这样可以直接解方程，得到函数的每一项的系数。

换句话来说，就算只知道 \(n+1\) 个点而不清楚函数本身的长相，我们也可以通过解方程得到函数本身。

那么对于一个函数 \(F(x)=f(x)*g(x)\) 我们如何计算呢？我们完全可以先在 \(f(x),g(x)\) 上都拿出 \(2n+1\) 个点，同时对于两个函数，选取的 \(2n+1\) 个 \(x\) 要对应相同，然后再将它们乘起来，如此就得到了 \(F(x)\) 上的 \(2n+1\) 个点。

原理的话，肯定是因为 \(F(x)=f(x)*g(x)\) 啊，带入一个确切的 \(x\) ，就能通过 \(f,g\) 的函数值乘积得到 \(F\) 的函数值了。

那么得到 \(f(x),g(x)\) 上 \(2n+1\) 个值的具体过程便被称为 \(DFT\) 了。

随后将 \(2n+1\) 个值乘起来，得到了 \(F(x)\) 上的点，然后反过来解出 \(F(x)\) ，这个具体过程便被称为 \(IDFT\) 。

具体操作过程如下。

DFT

我们发现如果要求 \(f(x)\) 上的 \(2n+1\) 个点，复杂度还是 \(n^2\) 级别的，如果想要优化这个过程，应该思考如何减少运算。

\(DFT\) 的原理就是进行分治，并利用过去已有的运算结果来计算现在的。这个过程的话，是利用复数相乘时，相当于在单位圆上转了个角度，所以我们只需要特殊构造这个复数，将其作为函数自变量代入，然后让它乘着乘着乘回出发的起点即可。

那么引入一下复数和单位根的知识点。

复数

对于一个复数 \(z\) ，它通常被表达为 \(z=a+bi\) ，其中 \(a,b\) 为实数，\(i\) 为虚数，\(i^2=-1\) 。

但同时，我们可以将复数 \(z\) 表示在一个坐标轴上，横坐标为实部系数 ，纵坐标为虚部系数，则 \(z\) 在坐标轴上的坐标可以被表示为 \(z(a,b)\) 。这个坐标轴被称为复平面。

同时，\(z\) 还可以被表示为 \((\)距离原点距离，在坐标轴上的角度\()\) ，其中，到原点的距离称之为模长，即为 \(\sqrt{a^2+b^2}\) 。而角度被称之为幅角，这个角度自然是相对于 \(x\) 轴正半轴的（显然不指夹角）。

也即 \(z(r,\theta)\) 。

那么对于复数相乘，可以证明：

\[z_1(r_1,\theta_1)*z_2(r_2,\theta_2)=z(r_1*r_2,\theta_1+\theta_2) \]

即模长相乘，幅角相加。

对于两个复数，当它们的模长为 \(1\) 时，乘积里由于模长不变，模长自然就没有意义了。

单位根

我们当然要方便计算啦，所以我们用于 \(DFT\) 的复数的模长都会是 \(1\) ，也就是说我们只需要在复平面上的单位圆上取数就好了。这样一来，所有的乘法都仅仅只是角度的变换了。

那么我们的起点显然是 \((1,0)\) 这个复数，如果要经过 \(n\) 次乘法，最后让它回到起点，那么每次乘法加上的的角度就应该要是 \(\frac{2\pi}{n}\) ，换句话说，将这个圆 \(n\) 等分，上面的点就是每次乘法的结果。

对于这个用来乘上 \(n\) 次的数，我们称之为 \(n\) 次单位根，表示为 \(\omega_n\) ，所以 \((\omega_n)^n=1\) 。

对于这个圆上的 \(n+1\) 个点，我们从 \((1,0)\) 开始依次标号为 \(\omega_n^0,\omega_n^1\dots\omega_n^n\) 。

那么对于单位根，有如下式：

一、

\[\omega_n^k=cos(\frac{2k\pi}{n})+sin(\frac{2k\pi}{n})i \]

二、

\[\omega_n^k=\omega_{2n}^{2k} \]

三、由于 \(\omega_n^{\frac{n}{2}}=-1\\\) （这个很显然吧，单位圆的一半）：

\[\omega_n^{k+\frac{n}{2}}=-\omega_n^k \]

DFT

设 \(n+1\) 为偶数，对于函数 \(f(x)=a_0+a_1*x\dots a_n*x_n\)，我们考虑将奇偶项分开：

\[f(x)=(a_0+a_2*x^2+...+a_{n-1}*x^{n-1})+(a_1*x+...+a_n*x^n) \]

若设：

\[f_1(x)=a_0+a_2*x+...+a_{n-1}*x^{\frac{n}{2}}\\ f_2(x)=a_1+a_3*x+...+a_n*x^{\frac{n}{2}} \]

则有：

\[f(x)=f_1(x^2)+x*f_2(x^2) \]

设 \(k\le \frac{n}{2}\) ，将单位根代入函数：

\[f(\omega_n^k)=f_1(\omega_n^{2k})+\omega_n^k*f_2(\omega_n^{2k})\\ =f_1(\omega_{\frac{n}{2}}^k)+\omega_n^k*f_2(\omega_{\frac{n}{2}}^k) \]

差不多的操作：

\[f(\omega_n^{k+\frac{n}{2}})=f_1(\omega_{\frac{n}{2}}^k)-\omega_n^k*f_2(\omega_{\frac{n}{2}}^k) \]

则只要得出 \(f_1,f_2\) 在 \(\omega_{\frac{n}{2}}^0\dots \omega_{\frac{n}{2}}^{\frac{n}{2}}\) 的所有取值，\(f\) 的所有取值也就出来了。

这样的话，只需要不断对于 \(f\) 递归下去即可，边界为常数项 \(a_0\) 。

每次可以得出 \(f\) 在区间长度下的那些取值，然后回溯的时候，父区间也就算出了当前区间长度两倍的取值。

以上就是 \(DFT\) 的过程。

IDFT

设 \(S(\omega_n^k)=1+\omega_n^k+(\omega_n^k)^2+\dots+(\omega_n^k)^{n-1}\) 。

根据等差数列求和可得，当 \(k=0\) ，\(S(\omega_n^k)=n\) ，否则结果为 \(0\) 。

设 \(DFT\) 得到了 \(n\) 次多项式 \(F(x)\) 的 \(n+1\) 个点值为 \((y_0,\dots yn)\) ，若其系数分别为 \((a_0,\dots a_n)\) ，我们考虑再设一个函数 \(G(x)\) ：

\[G(x)=y_0+y_1*x+\dots y_n*x^n \]

然后代入每个单位根的倒数：\((\omega_n^0,\omega_n^{-1}\dots \omega_n^{-n})\) 。这样我们便可以再通过 \(DFT\) 得到 \(n+1\) 个值 \(c\) 。式子为：

\[c_k=\sum_{i=0}^n y_i*(\omega_n^{-k})^i\\ =\sum_{i=0}^n(\sum_{j=0}^na_j(\omega_n^i)^j)*(\omega_n^{-k})^i\\ =\sum_{i=0}^n(\sum_{j=0}^na_j(\omega_n^j)^i)*(\omega_n^{-k})^i\\ =\sum_{i=0}^n(\sum_{j=0}^na_j(\omega_n^{j-k})^i)\\ =\sum_{j=0}^na_j\sum_{i=0}^n(\omega_n^{j-k})^i\\ =\sum_{j=0}^na_j*S(\omega_n^{j-k})\\ \]

所以当且仅当 \(j=k\) 时，\(S(\omega_n^{j-k})=n\) ，其他时候皆为 \(0\) 。

所以：

\[c_k=na_k \]

所以我们只需要通过 \(DFT\) 求出 \(c\) ，将每个数除以 \(n\) 就能得到 \(F(x)\) 的系数 \(a\) 。

这就是 \(IDFT\) 了。