浅谈FFT

Posted 2022-06-28 chandery

Fast Fourier Transportation（FFT）

·多项式的表达

系数表达

对于一个次数界为n的多项式\(A(x)=\sum_j=0^n-1a_jx^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_2n)?\)。

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

·多项式运算

\[ C_j=\sum_k=0^jA_kB_j-k \]

上述是多项式的乘法，我们把C称为A和B的卷积（convolution），表示成\(C=A\bigotimes 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\geqslant0,k\geqslant0,以及d>0, \]

\[ \omega^dk_dk=\omega^k_n \]

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

引理3：(求和引理)
\[ 对于任意整数n\geqslant0和不能被n整除的非负整数k，有 \]

\[ \sum_j=0^n-1(\omega^k_n)^j=0 \]

证明
\[ \sum_j=0^n-1(\omega^k_n)^j=\frac(\omega^k_n)^0(1-(\omega^k_n)^n)1-\omega^k_n=\frac(\omega^n_n)^k-1\omega^k_n-1=\frac(1)^k-1\omega^k_n-1=0 \]
因为要求k不能被n整除，而且仅当k被n整除时\[\omega^k_n=1\]成立，同时保证分母不为0。

DFT

回顾一下，我们希望计算次数界为n的多项式
\[ A(x)=\sum_j=0^n-1a_jx^j \]
在\[\omega^0_n,\omega^1_n,...,\omega^n-1_n\]处的值。假设A以系数形式给出，接下来定义结果\(y_k\):
\[ y_k=A(\omega^k_n)=\sum_j=0^n-1a_j\omega^kj_n \]
向量\(y=(y_0,y_1,...,y_n-1)\)就是系数向量\(a=(a_0,a_1,...,a_n-1)\)的离散傅里叶变换（DFT）。我们也记作\(y=DFT_n(a)\)。

FFT

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

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

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

\[ A^[0](x)=a_0+a_2x+a_4x^2+...+a_n-2x^n/2-1 \]

\[ A^[1](x)=a_1+a_3x+a_5x^2+...+a_n-1x^n/2-1 \]

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

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

所以就有代码：

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[ \beginmatrix y_0\ y_1\ y_2\ \vdots\ y_n-1 \endmatrix \right]= \left[ \beginmatrix 1 & 1 & 1 & 1 &\cdots& 1\ 1 & \omega_n & \omega^2_n & \omega^3_n &\cdots& \omega^n-1_n\ 1 & \omega^2_n & \omega^4_n & \omega^6_n &\cdots& \omega^2(n-1)_n\ 1 & \omega^3_n & \omega^6_n & \omega^9_n &\cdots& \omega^3(n-1)_n\ \vdots & \vdots& \vdots& \vdots& \ddots& \vdots\ 1 & \omega^n-1_n & \omega^2(n-1)_n & \omega^3(n-1)_n &\cdots& \omega^(n-1)(n-1)_n \endmatrix \right] \left[ \beginmatrix a_0\ a_1\ a_2\ \vdots\ a_n-1 \endmatrix \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));