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, heta)) 。

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

[z_1(r_1, heta_1)*z_2(r_2, heta_2)=z(r_1*r_2, heta_1+ heta_2) ]

即模长相乘，幅角相加。

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

单位根

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

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

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

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

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

一、

[omega_n^k=cos(frac{2kpi}{n})+sin(frac{2kpi}{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*xdots 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) ]

设 (kle 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}}^0dots 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_jsum_{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) 了。