浅谈FFT

Fast Fourier Transportation（FFT）

·多项式的表达

系数表达

对于一个次数界为n的多项式$A(x)={\sum }_{j=0}^{n-1}{a}_j{x}^j$而言，其系数表达是由一个系数组成的向量$a=(a_0,a_1,...,a_{n-1})$。

点值表达

一个次数界为n的多项式A(x)的点值表达就是一个由n个点值对所组成的集合
$$(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})$$
使得对k=0,1,…,n-1,所有$x_k$各不相同，
$$y_k=A(x_k)$$
简单的求点值运算我们可以随意代入n个不相同的数，然后得出点对，时间复杂度$\Theta (n^2)$。后面可以看到，如果我们用一点巧妙的取值，可以使时间复杂度优化到$\Theta (n{\lg }_{2}n)$。

求值计算的逆（从一个多项式的点值表达确定的系数表达形式）称为插值。

·多项式运算

$$C_j=\sum _{k=0}^j{A}_k{B}_{j-k}$$

上述是多项式的乘法，我们把C称为A和B的卷积（convolution），表示成$C=A⨂B$。

FFT的主要思路是首先把A和B转成点值表达，然后得到C的点值表达，再逆着做一遍，得到C的系数表达。

·DFT与FFT

上述做法太慢，我们要用$\Theta (n^2)$的时间把每个多项式转成点值表达，然后再用$\Theta (n^2)$的时间转回来，明显很慢，还不如暴力。

我们想，可不可以使用某些特殊的数，使得每次可以做一次运算就可以得到多个数的呢？答案是有的：单位根复数根

n次单位复数根是满足${\omega }^n=1$的复数$\omega$。n次单位复数根恰好有n个，对于k=0,1,…,n-1，这些根是$e^{\pi ik/n}$。为了解释这个表达式，我们用复数的指数形式来定义：
$$e^{iu}=cos(u)+isin(u)$$

也就是给定一个单位圆，上面均匀地分布着n个向量，如图：

·关于n次单位复根

以上图为例我们可以发现，每一个n（这里是8）次单位复根都是一个向量，他们在乘法意义下形成一个群。

引理1：(消去引理)
$$对于任意整数n⩾0,k⩾0,以及d>0,$$

$${\omega }_{dk}^{dk}={\omega }_n^k$$

证明
$${\omega }_{dk}^{dk}=(e^{2\pi i/dn}{)}^{dk}=(e^{2\pi i/n}{)}^k={\omega }_n^k$$
引理2：(折半引理)
$$如果n>0为偶数，那么n个n次单位根的平方集合就是n/2个n/2次单位根的集合$$
证明
$$({\omega }_n^{k+n/2}{)}^2={\omega }_n^{2k+n}={\omega }_n^{2k}{\omega }_n^n=({\omega }_n^k{)}^2$$
因此${\omega }_n^k$与${\omega }_n^{k+n/2}$平方相等。

引理3：(求和引理)
$$对于任意整数n⩾0和不能被n整除的非负整数k，有$$

$$\sum _{j=0}^{n-1}({\omega }_n^k{)}^j=0$$

证明
$$\sum _{j=0}^{n-1}({\omega }_n^k{)}^j=\frac{({\omega }_n^k{)}^0(1-({\omega }_n^k{)}^n)}{1-{\omega }_n^k}=\frac{({\omega }_n^n{)}^k-1}{{\omega }_n^k-1}=\frac{(1{)}^k-1}{{\omega }_n^k-1}=0$$
因为要求k不能被n整除，而且仅当k被n整除时$${\omega }_n^k=1$$成立，同时保证分母不为0。

DFT

回顾一下，我们希望计算次数界为n的多项式
$$A(x)=\sum _{j=0}^{n-1}{a}_j{x}^j$$
在$${\omega }_n^0,{\omega }_n^1,...,{\omega }_n^{n-1}$$处的值。假设A以系数形式给出，接下来定义结果$y_k$:
$$y_k=A({\omega }_n^k)=\sum _{j=0}^{n-1}{a}_j{\omega }_n^{kj}$$
向量$y=(y_0,y_1,...,y_{n-1})$就是系数向量$a=(a_0,a_1,...,a_{n-1})$的离散傅里叶变换（DFT）。我们也记作$y=DF{T}_n(a)$。

FFT

通过使用一种称为快速傅里叶变换（FFT）的方法，利用复根的特殊性质，我们就可以在$$\Theta (n\lg n)$$的时间内计算出$DF{T}_n(a)$。

注意：通篇的n我们假设是2的整数次幂。

FFT利用分治策略，采用A(x)中偶数下标的系数与奇数下标的系数，分别定义两个新的次数界为n/2的多项式

$$A^{[0]}(x)=a_0+a_{2}x+a_4{x}^2+...+a_{n-2}{x}^{n/2-1}$$

$$A^{[1]}(x)=a_1+a_{3}x+a_5{x}^2+...+a_{n-1}{x}^{n/2-1}$$

可以很容易发现
$$A(x)=A^{[0]}(x^2)+x{A}^{[1]}(x^2)$$
所以原问题转化为求两个次数界为n/2的多项式$$A^{[0]}(x)$$和$$A^{[1]}(x)$$在点$$({\omega }_n^0{)}^2,({\omega }_n^1{)}^2,...,({\omega }_n^{n-1}{)}^2$$的取值。

所以我们可以发现在求出$$A^{[0]}(x^2)$$和$$A^{[1]}(x^2)$$以后，可以算出两个复根的结果：
$$y_k=y_k^{[0]}+{\omega }_n^k{y}_k^{[1]}=A^{[0]}({\omega }_n^{2k})+{\omega }_n^k{A}^{[1]}({\omega }_n^{2k})=A({\omega }_n^k)$$
还有
$$y_{k+(n/2)}=y_k^{[0]}-{\omega }_n^k{y}_k^{[1]}=y_k^{[0]}+{\omega }_n^{k+(n/2)}{y}_k^{[1]}=A^{[0]}({\omega }_n^{2k})+{\omega }^{k+(n/2)}{A}^{[1]}({\omega }_n^{2k})$$
$$=A^{[0]}({\omega }_n^{2k+n})+{\omega }^{k+(n/2)}{A}^{[1]}({\omega }_n^{2k+n})=A({\omega }_n^{k+(n/2)})$$

所以就有代码：

void FFT(comp *a,int n,int inv){
	if(n==1) return;
	int mid=n/2;
	for (int i=0;i<mid;++i) c[i]=a[i*2],c[i+mid]=a[i*2+1];
	for	(int i=0;i<n;++i) a[i]=c[i];
	FFT(a,mid,inv);
	FFT(a+mid,mid,inv);
	comp wn={cos(2.0*pi/n),inv*sin(2.0*pi/n)},w={1,0};
	for (int i=0;i<mid;++i,w=w*wn){
		c[i]=a[i]+w*a[i+mid];
		c[i+mid]=a[i]-w*a[i+mid];
	} 
	for	(int i=0;i<n;++i) a[i]=c[i];
}

·在单位复数根的插值

现在我们展示如何在单位复数根处插值来完成多项式乘法方案,使得我们把一个多项式从点值表达转换回系数表达。
我们可以把DFT写成矩阵乘积
$$\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ ⋮\\ y_{n-1}\end{array}\right]=\left[\begin{array}{cccccc}1& 1& 1& 1& \cdots & 1\\ 1& {\omega }_n& {\omega }_n^2& {\omega }_n^3& \cdots & {\omega }_n^{n-1}\\ 1& {\omega }_n^2& {\omega }_n^4& {\omega }_n^6& \cdots & {\omega }_n^{2(n-1)}\\ 1& {\omega }_n^3& {\omega }_n^6& {\omega }_n^9& \cdots & {\omega }_n^{3(n-1)}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& {\omega }_n^{3(n-1)}& \cdots & {\omega }_n^{(n-1)(n-1)}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ ⋮\\ a_{n-1}\end{array}\right]$$
尴尬的是跑得贼慢：

随便卡卡就爆了…
分治难免递归，递归常数大。
于是，考虑改进。

·蝴蝶变换

盗图一张
可以发现，每个下标的二进制形式反过来就是它们最后在序列中的位置。
于是就有了迭代打法。

void FFT(Moon *a,int inv){
    int i,j,len;
    for (i=0;i<n;++i)
        if(i<R[i])swap(a[i],a[R[i]]);
    for (len=2;len<=n;len<<=1){
        int half=len/2;
        Moon w,wn={cos(Pi/half),inv*sin(Pi/half)};
        for (j=0;j<n-i;j+=len,w={1,0}){
        	for (i=0;i<half;++i,w=w*wn){
                Moon q=w*a[j+half+i],Q=a[j+i];
                a[j+half+i]=Q-q;
                a[j+i]=Q+q;
            }
        }
    }
}
int main(){
    for (i=0;i<n;++i) R[i]=(R[i>>1]>>1)|((i&1)<<(p-1));
}