「OI」快速傅里叶变换 FFT

快速傅里叶变换 $\texttt{FFT}$ 支持在 $\Theta(n\log _ 2n)$ 的时间内计算两个 $n$ 度的多项式的乘法。由于两个整数的乘法也可以被当作多项式乘法,因此这个算法也可以用来加速大整数的乘法计算。

概述

离散傅里叶变换(Discrete Fourier Transform,缩写为 $\operatorname{DFT}$),是傅里叶变换在时域和频域上都呈离散的形式。

$\texttt{FFT}$ 是一种高效实现 $\operatorname{DFT}$ 的算法,称为快速傅立叶变换(Fast Fourier Transform,FFT)。它对傅里叶变换的理论并没有新的发现,但是对于在计算机系统或者说数字系统中应用离散傅立叶变换,可以说是进了一大步。快速数论变换 $\texttt{NTT}$ 是快速傅里叶变换 $\texttt{FFT}$ 在数论基础上的实现。

前置知识

多项式的表示

系数表示法

系数表示法就是用一个多项式的各个项系数来表达这个多项式,即使用一个系数序列来表示多项式:

$$f(x)=\sum _ {i=0}^na_ ix^i\Leftrightarrow f(x)=\{a _ 0,a _ 1,\cdots,a _ n\}$$

点值表示法

点值表示法是把这个多项式看成一个函数,从上面选取 $n+1$ 个点,从而利用这 $n+1$ 个点来唯一的表示这个函数。

通俗地说,多项式由系数表示法转为点值表示法的过程,就是 $\operatorname{DFT}$ 的过程。相对地,把一个多项式的点值表示法转化为系数表示法的过程,就是 $\operatorname{IDFT}$。而 $\texttt{FFT}$ 就是通过取某些特殊的 $x$ 的点值来加速 $\operatorname{DFT}$ 和 $\operatorname{IDFT}$ 的过程。

单位复根

$$\omega^k _ n=e^{\frac{2\pi k}{n}}=\cos(\frac{2\pi k}{n})+i\sin(\frac{2\pi k}{n})$$

简单性质

  1. $$\omega^n _ n=1$$
  2. $$\omega^k _ n=\omega^{2k} _ {2n}$$
  3. $$\omega^{k+n} _ {2n}=-\omega^k _ {2n}$$

上面三条性质带入欧拉公式即可证明。

快速傅里叶变换 FFT

$\texttt{FFT}$ 算法的基本思想是分治。分治思想体现在将多项式分为奇次项和偶次项处理。

理论基础

以一个简单函数 $f(x)$ 为例:

$$f(x)=a _ 0+a _ 1x+a _ 2x^2+a _ 3x^3+a _ 4x^4+a _ 5x^5+a _ 6x^6+a _ 7x^7$$

奇偶分类,提公因式:

$$
\begin{aligned}
f(x)&=(a _ 0+a _ 2x^2+a _ 4x^4+a _ 6x^6)+(a _ 1x+a _ 3x^3+a _ 5x^5+a _ 7x^7)\\
&=(a _ 0+a _ 2x^2+a _ 4x^4+a _ 6x^6)+x(a _ 1+a _ 3x^2+a _ 5x^4+a _ 7x^6)
\end{aligned}
$$

不妨设:

$$g(x)=a _ 0+a _ 2x+a _ 4x^2+a _ 6x^3$$
$$h(x)=a _ 1+a _ 3x+a _ 5x^2+a _ 7x^3$$

$$
f(x)=g(x^2)+x\times h(x^2)
$$

所以我们有

$$
\begin{aligned}
\operatorname{DFT}(f(\omega^k _ n))&=\operatorname{DFT}(g((\omega^k _ n)^2))+\omega^k _ n\times\operatorname{DFT}(h((\omega^k _ n)^2))\\
&=\operatorname{DFT}(g(\omega^{2k} _ n))+\omega^k _ n\times\operatorname{DFT}(h(\omega^{2k} _ n))\\
&=\operatorname{DFT}(g(\omega^k _ \frac{n}{2}))+\omega^k _ n\times\operatorname{DFT}(h(\omega^k _ \frac{n}{2}))
\end{aligned}
$$

同理

$$
\operatorname{DFT}(f(\omega^{k+\frac{n}{2}} _ n))=\operatorname{DFT}(g(\omega^k _ \frac{n}{2}))-\omega^k _ n\times\operatorname{DFT}(h(\omega^k _ \frac{n}{2}))
$$

根据上面的推导,我们得出了两个式子:

$$
\begin{aligned}
\operatorname{DFT}(f(\omega^k _ n))=\operatorname{DFT}(g(\omega^k _ \frac{n}{2}))+\omega^k _ n\times\operatorname{DFT}(h(\omega^k _ \frac{n}{2}))\\
\operatorname{DFT}(f(\omega^{k+\frac{n}{2}} _ n))=\operatorname{DFT}(g(\omega^k _ \frac{n}{2}))-\omega^k _ n\times\operatorname{DFT}(h(\omega^k _ \frac{n}{2}))
\end{aligned}
$$

不难发现,如果我们求出了 $\operatorname{DFT}(g(\omega^k _ \frac{n}{2}))$ 和 $\operatorname{DFT}(h(\omega^k _ \frac{n}{2}))$,就能同时求出 $\operatorname{DFT}(f(\omega^k _ n))$ 和 $\operatorname{DFT}(f(\omega^{k+\frac{n}{2}} _ n))$。于是递归计算即可,这样的做法是 $\Theta(n\log _ 2n)$ 的,其中 $n$ 只需满足 $n=2^m,m\in\mathbb{N _ +}$ 且 $n$ 大于 $f$ 的次数。

蝴蝶变换,进一步的优化

其实就是递归太慢了,有个非递归的做法。

这个算法还可以从“分治”的角度继续优化。我们每一次都会把整个多项式的奇数次项和偶数次项系数分开,一只分到只剩下一个系数。但是,这个递归的过程需要更多的内存。因此,我们可以先“模仿递归”把这些系数在原数组中“拆分”,然后再“倍增”地去合并这些算出来的值。然而我们又要如何去拆分这些数呢?

找个规律即可发现:其实就是原来的那个序列,每个数用二进制表示,然后把二进制翻转对称一下,就是最终那个位置的下标。我们称这个变换为蝴蝶变换。根据它的定义,我们可以在 $\Theta(n)$ 的时间内求出每个数蝴蝶变换的结果。

for(reg int i=0;i<limit;++i)
    r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));

快速傅里叶逆变换

$\operatorname{IDFT}$(傅里叶反变换)就是把上一步获得的目标多项式的点值形式转换成系数形式。

其实不考虑推导过程的话也很简单。我们可以定义一个函数,在里面加一个参数 $1$ 或者是 $-1$,然后把它乘到 $\pi$ 的身上。传入 $1$ 就是 $\operatorname{DFT}$,传入 $-1$ 就是 $\operatorname{IDFT}$,$\operatorname{IDFT}$ 最后每一个数还要除掉 $n$。

简单代码实现

inline void DFT(reg Complex a[],reg int f){
    for(reg int i=0;i<limit;++i)
        if(i<r[i])
            swap(a[i],a[r[i]]);
    for(reg int i=1;i<limit;i<<=1){
        Complex w(cos(pi/i),f*sin(pi/i));
        for(reg int j=0;j<limit;j+=(i<<1)){
            Complex e(1,0);
            for(reg int k=0;k<i;++k,e=e*w){
                static Complex x,y;
                x=a[j+k],y=a[j+k+i]*e;
                a[j+k]=x+y,a[j+k+i]=x-y;
            }
        }
    }
    if(f==-1)
        for(reg int i=0;i<limit;++i)
            a[i].x/=1.0*limit;
    return;
}

简单例题

进阶

参考资料