全文为算法导论的笔记作为备忘


$Fast \ Fourier \ Transform$($ FFT$)解决的问题是

$\displaystyle A(x)=\sum_{j=0}^{n-1}a_jx^j ,\ B(x)=\sum_{j=0}^{n-1}b_jx^j$($ n$为多项式项数记$ degree(A)=k$)

任何一个大于一个多项式次数的整数都是该多项式的次数界

​$\displaystyle C(x)=A(x)×B(x)=\sum_{j=0}^{2n-2}(\sum_{k=0}^{j}a_kb_{j-k})x^j$

$ \displaystyle A(x)=a[0]+a[1]x+a[2]x^2+a[3]x^3+···+a[n-1]x^{n-1}$ $ B(x)=b[0]+b[1]x+b[2]x^2+b[3]x^3+···+a[n-1]x^{n-1}$
$\displaystyle C(x)=c[0]+c[1]x+c[2]x^2+c[3]x^3+···+c[2n-2]x^{2n-2}$

$\displaystyle c[0]=a[0]×b[0];\ c[1]=a[0]×b[1]+b[0]×a[1]; \ c[3]=a[1]×b[2]+a[2]×b[1]$

​ $\displaystyle c_j=\sum_{j=0}^{2n-2}(\sum_{k=0}^{j}a_kb_{j-k})x^j$
$c[j+k]=\sum_{j=0}^{n-1}\sum_{k=0}^{n-1}a_j b_k$

  • $ O(n^2)$
void calC(double *a,double *b,double *c,int n)
{
    for(int i=0;i<(n<<1);i++) c[i]=0;//init
    for(int i=0;i<=(n-1)<<1;i++)
      for(int j=0;j<=i&&j<n;j++)
        c[i] += a[j] * b[i-j];
    //or
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++) c[i+j] += a[i] * b[j];
}

  • 优化到$ (nlogn)$ 注:下面使用专用符号$ i$表示$ \sqrt{-1}$

多项式系数表示法通过系数向量$ \displaystyle \textbf{A}=${$ a_0,a_1,…a_{n-1}$}来确定唯一的多项式

​ $ \displaystyle A(x)=a_0+a_1x+a_2x^2+a_3x^3+···+a_{n-1}x^{n-1}$

多项式在$ x_0$处求值可系通过数表达式利用霍纳法则$ O(n)$输出$ A(x_0)$

一个次数界为$n$的多项式的点值表达式就是一个由$ n$各点值对组成的集合

​ { $ \displaystyle (x_0,y_0), (x_1,x_2),(x_3,y_3),…,(x_{n-1},y_{n-1})$ }

使得对于任意的整数$ \displaystyle k=0,1,2,···,n-1, \ x_k$各不相同,其中$ y_k=A(x_k)$

可以证明对于一组互不相同的{ $ \displaystyle x_0,x_1,x_2,x_3,…,x_{n-1}$ }(称为该种表示法的基), 且该表示的多项式唯一

​ $\displaystyle A(x)=\sum_{j=0}^{n-1} a_jx^j=a_0+a_1x+a_2x^2+···+a_{n-1}x^{n-1}$

$ \displaystyle c_j=\sum_{k=0}^{j}a_kb_{j-k}$对应的向量$ \textbf{c}$即多项式$ A×B$的系数向量称为输入向量$ \textbf{a}$和$ \textbf{b}$的卷积$ (convolution)$,表达成$ c=a\otimes b$

求值计算的逆称为插值 定理30.1:插值多项式的唯一性[1]:当插值多项式的次数界等于点值表达式的点值对插值才是明确的

  • 点值表达式下多项式加法运算$\displaystyle C$:{$ (x_0,y_0+{y_0}’), (x_1,y_1+{y_1}’),···,(x_k,y_{n-1}+{y_{n-1}}’) $}

  • 点值表达式下多项式乘法运算:$ degree(C)=degree(A)+degree(B)$需扩展$ A(x),B(x)$为$ 2n$项

    得出$ C$: { $ \displaystyle (x_0,y_0{y_0}’), (x_1,y_1{y_1}’),···,(x_k,y_{n-1}{y_{n-1}}’) $

目前的思路由系数表达式霍斯法则$ O(4n)$求$ 4n$个点的多项式值转换成点值表达式相乘$ O(n)$,得出的结果通过求逆$ O(n^2)$~$ O(n^3)$再转换成系数表达式

  • 此处还能进行的优化通过$ DFT$(离散傅里叶变换)在特定的单位复数根 处达到$ O(nlogn)$求值和插值

下面着重于单位复数根

$ n$次单位复数根是满足$\displaystyle {\omega}^n=1$的复数$\omega$,$n$次单位复数根恰好有$ n$个 :对于$ k=0,1,2,3···n-1$ ,$ \omega_n=e^{\frac{2\pi ik}{n}}(k=1)$称为主$ n$次单位根

  • 复数的指数定义即欧拉公式: $ \displaystyle e^{iu}=cos(u)+isin(u)$

可得这$n$个单位复数根均匀地分布在复平面的原点为圆心的单位半径的圆周上

$ n$个$ n$次单位复数根 $ \displaystyle \omega_n^0,\omega_n^1,\omega_n^2,···,\omega_n^{n-1}$,且在乘法意义下形成一个与加法群具有相同结构的群

单位根基本性质

  • $ \displaystyle \omega_n^n=\omega_n^0=1$,$ \omega_n^{-1}=\omega_n^{n-1}$, $ \omega_n^j\omega_n^k=\omega_n^{(j+k)mod \ n}$
  • 消去引理: $ \displaystyle \omega_{dn}^{dk}={(e^{\frac{2 \pi i}{dn}})}^{dk}=e^{\frac{2 \pi ik}{n}}=\omega_n^{k}$,该引理得出$ (\omega_n^k)^2=\omega_{\frac{n}{2}}^k$
  • 折半引理:对于$ (k>0的偶数)$那么这$ n$个$n$次单位复数根的平方的集合就是$ \displaystyle \frac{n}{2}$个$\frac{n}{2}$单位复数根集合。由$ \omega_n^n=1$ 易证:$\displaystyle (\omega^{\frac{n}{2}+k})^2 =(\omega_n^k)^2(k=0,1,2,3\cdot \cdot \cdot, n-1)$ 那么如果对$ n$次单位复数根进行平方将得到每个$ \frac{n}{2}$次单位根$ 2$次
  • 求和引理: $\displaystyle \sum_{j=0}^{n-1}=(\omega_n^k)^j=0$

$ DFT与FFT$(求$ A(x_0)$)

  • $ DFT$

$ \displaystyle A(x)=\sum_{j=0}^{n-1}a_jx^j$ 在 $ \omega_n^0,\omega_n^1,\omega_n^2, ···,\omega_n^{n-1}$即在$ n$个$ n$次单位复数根处的值(此处$ n$实际是$ 2n$个)组成的向量即为$ DFT_n(a)$

假设$ A(x)$的系数表示$ a=(a_0,a_1,a_2,···,a_{n-1})$ 对$ (k=0,1,2,···,n-1)$定义结果$ y_k$:

​ $ \displaystyle y_k=A(\omega_n^k)=\sum_{j=0}^{n-1}a_j\omega^{kj}$

记结果向量$ \displaystyle \textbf{y}=(y_0,y_1,···,y_{n-1})$就是系数向量 $ \displaystyle \textbf{a}=(a_0,a_1,a_2,···,a_{n-1})$的离散傅里叶变换$ (DFT)$记$ y=DFT_n(a)$

  • $FFT$

前置假设$n$是$2$的幂次 $FFT$ 采用奇偶分治

定义两个形式一致的次数界为$\frac{n}{2}$的多项式$A^{[0]}(x),A^{[1]}(x)$

​ $\displaystyle A^{[0]}(x)=a_0+a_2x+a_4x^2···+a_{n-2}x^{\frac{n}{2}-1}$

​ $\displaystyle A^{[1]}(x)=a_1+a_3x+a_5x^2···+a_{n-1}x^{\frac{n}{2}-1}$

​ 一致形式$\displaystyle A(x)=a_0+a_1x+a_2x^2+···+a_{\frac{n}{2}-1}x^{\frac{n}{2}-1}$

有$\displaystyle A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)$

那么求$A(x)$在$\displaystyle \omega_n^0,\omega_n^1,\omega_n^2, ···,\omega_n^{n-1}$ 处的值就转换成:

  1. 求次数界为$\frac{n}{2}$的多项式在$ A^{[0]}(x)和A^{[1]}(x)$ 在

    ​ $\displaystyle (\omega_n^0)^2,(\omega_n^1)^2,(\omega_n^2)^2, ···,(\omega_n^{n-1})^2$ 的取值

  2. 折半引理 可知上式由$\frac{n}{2}$个不同的数,现在已经成功把规模为$ n$的$ DFT_n$ 划分成两个规模为$ \frac{n}{2}$个元素的$ DFT_{\frac{n}{2}}$

$FFT$递归求解$n$个元素构成的 $\textbf{a}=(a_0,a_1,a_2,···,a_{n-1})$的$ DFT$

  • 小结:给出系数向量$\textbf{a}$ $DTF_n(a)$过程为先求出$n$个$ n$次单位复数根的平方。如果只有一个元素$ a_0$ 那么单位根为$ \omega_n^0=1$,$ y_0=a_0w_n^0=a0$ 将$ \frac{n}{2}$个元素带入一致形式,转换成程序语言就是不断对$ \textbf{a}$向量划分奇偶,每次划分求值。
  1. 对于$\displaystyle k=0,1,2,3···,\frac{n}{2}-1$我们定义$ y$为分治后的向量$ y_{k}^{[0]}=A^{[0]}(\omega_{\frac{n}{2}}^k)$ , $ y_k^{[1]}=A^{[1]}(\omega_{\frac{n}{2}}^k)$

​ 那么$\displaystyle y_k=y_k^{[0]}+w_n^ky_k^{[1]} =A^{[0]}(\omega_{\frac{n}{2}}^k)+w_n^kA^{[1]}(\omega_{\frac{n}{2}}^k)=A^{[0]}(\omega_n^k)+w_n^kA^{[1]}(\omega_n^k)=A(\omega_n^k)$

  1. 对于$\displaystyle k=\frac{n}{2},\frac{n}{2}+1,\frac{n}{2}+2,\frac{n}{2}+3···,n-1$那么

    $ \displaystyle y_{k+\frac{n}{2}}=y_k^{[0]}-w_n^ky_k^{[1]} =A^{[0]}(\omega_n^{2k})+w_n^{k+\frac{n}{2}}A^{[1]}(\omega_n^{2k})$

    ​ $ \displaystyle =A^{[0]}(\omega_n^{2k+n})+w_n^{2k+n}A^{[1]}\omega_n^{2k+n}=A(\omega_n^{2k+n})=A(\omega_n^{k+\frac{n}{2}})$

    ( $ \displaystyle \omega_{\frac{n}{2}}^k=\omega_{n}^{2k}$ , $ \omega_n^{2k+n}=-\omega_n^k$ )

//FFT求值C++
vector<Complex> recursive-FFT(vector<Complex>const &a) {
    int n = a.size();
    if(n == 1) return a;
    Complex w_n = complex(cos(2<em>pi/n), sin(2</em>pi/n));//构造主n次单位根 
    Complex w(1,0);
    vector<Complex> a0,a1,y;
    a0.reserve(n/2); a1.reserve(n/2); y.reserve(n)
    for(int i=0;i<n/2;i++) {
        a0.push_back(a[i<em>2]);//make_vector(a0,a2,a4,...,a n-2);//划分奇偶
        a1.push_back(a[i</em>2+1]);//make_vector(a1,a3,a5,...,a n-1);
    }
    vector<Complex> y0 = recursive-FFT(a0);//分治递归成下标为偶数的向量
    vector<Complex> y1 = recursive-FFT(a1);//奇数
    for(int k=0; k<n/2; k++) {//0~n/2-1
        y[k] = y0[k] + w * y1[k];
        y[k+n/2] = y0[k] - w * y1[k];
        w = w * w_n;
    }
    return y;
}

发表评论

电子邮件地址不会被公开。 必填项已用*标注

This site uses Akismet to reduce spam. Learn how your comment data is processed.