快速傅里叶变换(FFT)学习笔记(其一)

2023-02-27 07:58:22

再探快速傅里叶变换(FFT)学习笔记(其一)

再探快速傅里叶变换(FFT)学习笔记(其一)

写在前面

为什么写这篇博客

笔者去年暑假刚刚学习过FFT,NTT的一些基础应用。但当时对FFT和NTT的理解还不够深入。本博客参考2016年国家集训队论文中雅礼中学毛啸的《再探快速傅立叶变换》,对之前学习时的不足之处做了补充。

为了不使篇幅过长，预计将把学习笔记分为四部分:

DFT,IDFT,FFT的定义,实现与证明
NTT的实现与证明,任意模数NTT
FFT的优化技巧
多项式相关操作

注意:本文为了严谨性可能会牺牲一些可读性

一些约定

$[p(x)]=\begin{cases}1,p(x)为真 \\ 0,p(x)为假 \end{cases}$
本文中序列的下标从0开始
若$s$是一个序列,$|s|$表示$s$的长度
若大写字母如$F(x)$表示一个多项式,那么对应的小写字母如$f$表示多项式的每一项系数，即$F(x)=\sum_{i=0}^{n-1} f_ix^i$

前置知识

多项式卷积

多项式$A(x)=\sum_{i=0}^{n-1} a_ix^i$(注意最高次项的指数为$n-1$,至于为什么要这样定义见下文)

定义1.1 对于两个多项式$A(x)=\sum_{i=0}^{n-1} a_ix_i,B(x)=\sum_{i=0}^{m-1} b_ix_i$,我们定义$A(x)$和$B(x)$的卷积为:

\[A(x) \cdot B(x)=\sum_{i=0}^{n+m-2} x^i\sum_{k=0}^i a_kb_{i-k}
\]

实际上我们只要考虑系数，也就是给出两个序列$a,b$,求序列$c$使得$$c_i=\sum_{k=0}^i a_k b_{i-k} \tag{1.1}$$

直接计算卷积显然是$O(n^2)$的,FFT可以用来优化这个过程

多项式的系数表达式和点值表达式

系数表达式就是多项式中每一项的系数$a_0,a_1 \dots a_{n-1}$

我们知道,n个点可以确定一个n次多项式.那么，一个n次多项式$F(x)$也可以表示成这样的一个集合$\{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)) \dots (x_{n-1},f(x_{n-1}))\}$.比如多项式$f(x)=x^2+x+1$就可以表示成$\{(1,3),(2,7),(-1,1)\}$

这样有什么用呢?如果已知$A(x)$和$B(x)$的点值表达式,可以$O(n)$计算$C(x)=A(x)B(x)$的点值表达式

比如$A(x)=\{(x_0,A(x_0)),(x_1,A(x_1)),(x_2,A(x_2)),....,(x_n,A(x_n))\}$

$B(x)=\{(x_0,B(x_0)),(x_1,B(x_1)),(x_2,B(x_2)),....,(x_n,B(x_n))\}$

$C(x)=A(x)\times B(x)=\{(x_0,A(x_0)\times B(x_0)),(x_1,A(x_1)\times B(x_1)),(x_2,A(x_2)\times B(x_2)),....,(x_n,A(x_n)\times B(x_n))\}$

所以:我们可以在$O(n)$的时间内计算出两个点值表达式的乘积,只需要把对应点的$y$坐标相乘即可。但是,把系数表达式和点值表达式互相转化仍然需要$O(n^2)$时间.(每代入一个$x$都需要$O(n)$时间计算$A(x)$)

单位根及其性质

定义1.2:$n$次单位根为满足:$n$是最小的满足$\omega^n=1$的复数$\omega$

容易发现,$n$次单位根有$n$个,他们分布在复平面上的单位圆上。(这里我们只讨论$n$为偶数的情况)

把单位圆$n$等分,圆上的每个分割点都对应一个单位根.从点$(1,0)$开始,逆时针把这$n$个点从$0$开始编号,第$k$个点对应的复数记作$w_n^k$,上图呈现的就是$\omega_{16}^0 \text{~} \omega_{16}^{15}$

容易发现$\omega_n^k$的幅角为$\frac{2k\pi}{n}$.(一整个圆为$2\pi \ \text{rad}$,分为$n$份,再取$k$份)，因此，有:

\[\omega_n^k=\cos \frac{2k\pi}{n}+\text{i} \sin \frac{2k\pi}{n} \tag{1.2.1}
\]

单位根有如下的性质:

定理1.2.1 $\forall p \in \mathbb{Z}, \omega_n^k=\omega_{pn}^{pk}$

代数证明:

$\omega_{pn}^{pk}=\cos \frac{2pk \pi}{pn}+ \text{i}\sin \frac{2pk\pi}{pn}=\cos \frac{2k\pi}{n}+\text{i} \sin \frac{2k\pi}{n}=\omega_n^k$

几何证明:单位圆分成$n$份，数$k$次，和分成$pn$份，数$pk$次所代表复数是一样的

定理1.2.2 $\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^k$

代数证明:

\[\begin{aligned} \omega_{n}^{k+\frac{n}{2}} &= \cos \frac{2(k+\frac{n}{2}\pi)}{n}+ \text{i} \sin \frac{2(k+\frac{n}{2}\pi)}{n}\\ &=\cos(\frac{2k \pi}{n}+\pi)+\text{i}\sin(\frac{2k \pi}{n}+\pi) \\&= -\cos \frac{2k\pi}{n}-\text{i} \sin \frac{2k\pi}{n}=-\omega_n^k\end{aligned}
\]

倒数第二步用到了三角函数的诱导公式

几何证明:

$\frac{n}{2}$代表半圈，所以这个式子的意思是从编号为$k$的复数开始走半圈，代表的复数和之前的正好相反（实部和虚部都相反）

定理1.2.3 $\omega_n^p \times \omega_n^q=\omega_n^{p+q}$

令$\alpha=\frac{2p\pi}{n},\beta=\frac{2p\pi}{n}$

则$\omega_n^p=\cos \alpha+ \text{i} \sin \alpha,\omega_n^q=\cos \beta+ \text{i} \sin \beta$

\[\begin{aligned}
\omega_{n}^{p} \times \omega_{n}^{q}&=(\cos \alpha \cos \beta-\sin \alpha \sin \beta)+i(\sin \alpha \cos \beta+\cos \alpha \sin \beta)\\ &=\cos (\alpha+\beta)+i \sin
(\alpha+\beta)\\&=\omega_{n}^{p+q}
\end{aligned}\]

证明用到了复数乘法和三角函数的和角公式.注意到定理1.2.2实际上是定理1.2.3的推论，证明留给读者。这个公式极为重要，在FFT的公式推导中会用到。

定理1.2.3的另一个推论，$\omega_n^k=(\omega_n^1)^k$

证明显然

定理1.2.4

\[\forall v \in \mathbb{N},\frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{v k}=[v \bmod n=0]
\]

证明:

当$v \bmod n =0$时, $\omega_{n}^{vk}=\cos\frac{2vk \pi}{n}+\text{i} \sin \frac{2vk \pi}{n}$

那么$\frac{vk}{n}$一定是整数,所以$\cos\frac{2vk \pi}{n}=1,\sin \frac{2vk \pi}{n}=0,\omega_n^{vk}=1$

所以$\frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{v k}=\frac{1}{n} \times n =1$

当$v \mod n \neq 0$时,容易发现$\omega_n^v \neq 1$,

那么上式实际上是一个等比数列,根据等比数列求和公式有:

$$\begin{aligned} \frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{v k}&=\frac{1-\omega_n^{vn}}{n(1-\omega_n^v)} \ &= \frac{1-\omega_1^v}{n(1-\omega_n^v)}\end{aligned}$$

而$\omega_1^v=\cos(2v\pi)+i\sin(2v \pi)=1$

因此$\frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{v k}=0$

证毕.

这个定理会在IDFT的正确性证明中用到

DFT和IDFT

在上文中我们提到了多项式的系数表达式和点值表达式，并发现点值表达式可以简化卷积运算。

离散傅里叶变换(Discrete Fourier Transform,DFT):是把系数表达式转化为点值表达式的一种算法

逆离散傅里叶变换(Inverse Discrete Fourier Transform,IDFT):是把点值表达式转化为系数表达式的一种算法。

那么求$A$和$B$的卷积就可以先对$A$和$B$做DFT,再将$A$和$B$的点值表达式相乘得到$C$,最后再对$C$做IDFT

DFT的过程

对于多项式$A,B$,我们找到最小的大于$\max(|A|,|B|)$的2的整数幂$n$,至于为什么是2的整数幂后面会提到。给定序列$a$,对于$k \notin [0,|a|)$,我们定义$a_k=0$

要想得到一个系数表达式的点值表达式，我们可以任意选$n$个值代入。但是傅里叶想到,单位根有着特殊的性质，可以把$n$次单位根代入多项式。形式化的说,设多项式$A(x)=\sum_{i=0}^{n-1}a_ix^i$,则$A(x)$的DFT变换结果为$a'_k=A(\omega_n^k)=\sum_{i=0}^{n-1}a_i \omega_n^{ki}$

我们后面提到的多项式与系数序列的对应关系都类似$A(x)$与$a$的对应关系，对$A(x)$做DFT就是对序列$a$做DFT;说一个东西做DFT得到了$A(x)$,就是得到了序列$a$

容易发现DFT的时间复杂度为$O(n^2)$

IDFT的过程

代入$n$次单位根的作用在IDFT中显现出来.

定理2.1 :把多项式$A(x)$的离散傅里叶变换结果$f(\omega_n^0),f(\omega_n^1) \dots f(\omega_n^{n-1})$作为另一个多项式$B(x)$的系数，取单位根的共轭复数即$\omega_n^0,\omega_n^{-1},\omega_n^{-2} \dots \omega_n^{-(n-1)}$作为$x$代入$B(x)$，得到的每个数再除以$n$，得到的是$A(x)$的各项系数

也就是说,把点值表达式中点的$y$坐标作为另一个多项式的系数后,IDFT的过程和DFT几乎一样,只是代入的是单位根的共轭复数。

我们暂时不证明定理2.1,而是证明一个更强的定理:

定理2.2:把多项式$A(x)$的DFT结果$a$,和$B(x)$的DFT结果$b$对应相乘,得到序列$t_i=a_ib_i$;再对$t_i$做IDFT,得到的序列就是$C(x)=A(x)\cdot B(x)$的各项系数$c_i$

证明:

由卷积的定义式$(1.1)$我们有:

\[c_{r}=\sum_{p, q}[(p+q) \bmod n=r] a_{p} b_{q} \tag{2.1}
\]

看到取模想到定理1.1.4

\[\forall v \in \mathbb{N},\frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{v k}=[v \bmod n=0] \tag{2.2}
\]

故:

\[\begin{aligned}
&[(p+q) \bmod n=r] \\
=&[(p+q-r) \bmod n=0] \\
=& \frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{(p+q-r) k} \\
=& \frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{-r k} \omega_n^{p k} \omega_n^{q k}
\end{aligned}\]

代入$(2.1)$

\[\begin{align*}
c_{r} &=\sum_{p, q}[(p+q) \bmod n=r] a_{p} b_{q} \\
&=\sum_{p, q} \frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{-r k} \omega_n^{p k} \omega_n^{q k} a_{p} b_{q} \\
&=\frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{-r k} \sum_{p, q} \omega_n^{p k} a_{p} \omega_n^{q k} b_{q} \\
&=\frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{-r k} \sum_{p} \omega_n^{p k} a_{p} \sum_{q} \omega_n^{q k} b_{q}\tag{2.3} \\&=\frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{-r k} a'_k b'_k \\&=\frac{1}{n} \sum_{k=0}^{n-1} \omega_n^{-r k} t_k\end{align*} \]

容易发现定理2.1是定理2.2 的推论,即$B(x)$的系数都为1时的情况.式$(2.3)$实际上指出了DFT和IDFT求卷积的过程。

和DFT类似,IDFT的时间复杂度为$O(n^2)$

FFT

上一部分我们提到了DFT的公式，但是这样直接做还是太慢了。因此我们要利用单位根的性质来加速DFT,这个算法被称为快速傅立叶变换(Fast Fourier Transform,FFT). 同样,我们也可以加速IDFT,即IFFT.

这几个算法的关系如下图所示:

其中FFT,IFFT的时间复杂度均为$O(n \log n)$

FFT的数学证明及时间复杂度分析

我们对多项式$A(x) = a_0 + a_1x + a_2x^2 +...+a_{n-1}x^{n-1}$进行FFT,仍然把$\omega_n^k$代入

把$a$按照下标奇偶性分为两部分:

$A_1(x) = a_0 + a_2x + ... + a_{n - 2}x^{\frac{n}{2} - 1}$

$A_2(x) = a_1 + a_3x + ... + a_{n - 1}x^{\frac{n}{2} - 1}$

那么我们有:

$A(x) = A_1(x^2) + xA_2(x^2)$

将$x=\omega_n^k$代入,

对于$k < \frac{n}{2}$(根据先前的定义,n为2的次幂,所以没有取整问题)

\[\begin{align*} A(\omega_n^k) &= A_1(\omega_n^{2k}) + \omega_n^kA_2(\omega_n^{2k}) \\ &= A_1(\omega_{\frac{n}{2}}^{k}) + \omega_n^kA_2(\omega_{\frac{n}{2}}^{k}) \tag{3.1}\end{align*}
\]

证明用到了定理1.2.1:$\forall p \in \mathbb{Z}, \omega_n^k=\omega_{pn}^{pk}$

对于剩下的部分,即$\omega_{n}^{k+\frac{n}{2}}(k<\frac{n}{2})$

\[\begin{align*} A(\omega_n^{k + \frac{n}{2}}) &= A_1(\omega_n^{2k + n}) + \omega_n^{k + \frac{n}{2}}A_2(\omega_n^{2k + n}) \\ &= A_1(\omega_{\frac{n}{2}}^{k} \times \omega_n^n) + \omega_n^{k + \frac{n}{2}} A_2(\omega_{\frac{n}{2}}^{k} \times \omega_n^n) \\ &= A_1(\omega_{\frac{n}{2}}^{k}) - \omega_n^kA_2(\omega_{\frac{n}{2}}^{k}) \tag{3.2} \end{align*}
\]

证明用到了定理1.2.1和定理1.2.3:$\omega_n^p \times \omega_n^q=\omega_n^{p+q}$

式$(3.1),(3.2)$也被称为“蝴蝶操作”, 通过蝴蝶操作,我们只要知道两个多项式$A_1(x),A_2(x)$在$(\omega_{\frac{n}{2}}^{0}, \omega_{\frac{n}{2}}^{1}, \omega_{\frac{n}{2}}^{2}, ... , \omega_{\frac{n}{2}}^{\frac{n}{2} - 1})$处的点值表达式,就可以求出$A(x)$在$(\omega_n^0, \omega_n^1, \omega_n^2, ..., \omega_n^{n-1})$处的点值表达式了了。每次递归下去规模缩小一半,$n=1$时终止.由于长度$n$每次除2,且要求每次的$n$均为偶数，所以我们需要$n$是2的正整数次幂

设给长度为$n$的序列做$FFT$的时间为$T(n)$,那么我们有

\[T(n)=2T(\frac{n}{2})+n
\]

根据主定理,$T(n)=\Theta(n \log n)$

FFT的递归实现

FFT递归实现很好理解.注意可以使用C++自带的复数类<complex>,也可以自己手写复数类。

<complex>重载了+,-,*,/,=,==等运算符，在使用上很方便，几乎和操作整型变量没有区别。两个成员函数.real(),.imag()分别返回复数的实部和虚部。更多它的用法请参考这个链接

#include<complex>//c++自带复数类

const double pi=acos(-1.0);

typedef complex<double> com;

com a[maxn+5],b[maxn+5],c[maxn+5];

void fft(com *x,int n,int type){//orz Fourier

	if(n==1) return;//递归边界

	com l[n>>1],r[n>>1];

	for(int i=0;i<n;i++){//按奇偶分类

		if(!(i&1)) l[i>>1]=x[i];

		else r[i>>1]=x[i];

	}

	fft(l,n>>1,type);//递归对n/2的长度做FFT

	fft(r,n>>1,type);

	com wn1=com(cos(2*pi/n),type*sin(2*pi/n))/*w_n^0*/,wnk=com(1,0);

	for(int i=0;i<(n>>1);i++){

		x[i]=l[i]+wnk*r[i];//式(3.1)

		x[i+(n>>1)]=l[i]-wnk*r[i];//式(3.2)

		wnk*=wn1;//递推算w_n^k

	}

}

int main(){

     //....

     int m=n*2;//这个数根据题目算要求来定,这里是长度为n的两个序列卷积,那就取2n

	 n=1;

	 while(n<m) n*=2;//计算出最小的2的次幂

	 fft(a,n,1);

	 fft(b,n,1);

	 for(int i=0;i<n;i++) c[i]=a[i]*b[i];

	 fft(c,n,-1);

     for(int i=0;i<m;i++) printf("%lf ",c[i].real()/n);//记得除n

    //....

}

FFT的非递归实现

很遗憾,FFT的递归实现常数极大(虽然非递归实现常数也很大),几乎无法通过任何OI中FFT的题目。

我们观察递归过程中根据奇偶改变系数位置的方式,例如长度为8的情况:

\[\begin{aligned}
&\left(a_{0}, a_{1}, a_{2}, a_{3}, a_{4}, a_{5}, a_{6}, a_{7}\right)\\
&\left(a_{0}, a_{2}, a_{4}, a_{6}\right)\left(a_{1}, a_{3}, a_{5}, a_{7}\right)\\
&\left(a_{0}, a_{4}\right)\left(a_{2}, a_{6}\right)\left(a_{1}, a_{5}\right)\left(a_{3}, a_{7}\right)\\
&\left(a_{0}\right)\left(a_{4}\right)\left(a_{2}\right)\left(a_{6}\right)\left(a_{1}\right)\left(a_{5}\right)\left(a_{3}\right)\left(a_{7}\right)
\end{aligned}\]

观察$0,4,2,6,1,5,3,7$的二进制表示$000,100,010,110,001,101,011,111$，

比较$0,1,2,3,4,5,6,7$的二进制表示$000,001,010,011,100,101,110,111$

我们发现,初始时位置$x$上的数会被交换到x二进制位翻转后的位置,记为$rev(x)$。注意这里的反转指的是反过来写,而不是每一位交换01。这种操作也被称为“位逆序置换”。那么只要处理出这个序列,每次FFT的都是一段连续的区间,可以非递归求解。

$rev(x)$可以根据如下代码递推求出

for(int i=0;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));

证明: 对于数$i$,我们已经求出$i>>1$的置换结果,如$k=5,i=(10111)_2,i>>1=(1011)_2$;那么我们先反转从左到右前$k-1$位,即$rev[i>>1]=rev[(01011)_2]=(11010)_2$,右移一位去掉前导0(翻转后到了末尾)就得到了前$k-1$位反转的结果; 然后再加上最后一位,最后一位的值是$(i\text{&}1)$,要左移到最前面.

int rev[maxn+5];//数x位逆序置换后变成rev[x]

int n,k;//n=2^k

void fft(com *x,int n,int type){//type=1做FFT,type=-1做IFFT

    for(int i=0;i<n;i++) if(i<rev[i]) swap(x[i],x[rev[i]]);//做置换,i<rev[i]防止一个数被换两次换回原位

    for(int len=1;len<n;len*=2){

        int sz=len*2;//当前处理的长度为len*2,len=1的情况不用计算,就是x[i]

        com wn1=com(cos(2*pi/sz),sin(2*pi/sz)*type);//计算n次单位根w_n^1,type=-1时得到共轭复数

        for(int l=0;l<n;l+=sz){//枚举一对长度为len的区间，对这两个区间做蝴蝶操作

            int r=l+len-1;//[l,l+len-1]为前半段,[l+len,r]为后半段

            com wnk=1;

            for(int i=l;i<=r;i++){

                com tmp=x[i+len];//用tmp临时存储值

                x[i+len]=x[i]-wnk*tmp;//式(3.1)

                x[i]=x[i]+wnk*tmp;//式(3.2)

                wnk*=wn1;//递推计算w_n^k

            }

        }

    }

    if(type==-1) for(int i=0;i<n;i++) x[i]/=n;//IFFT记得除n

}

int main(){

    n=1,k=0;

	while(n<m){

        n*=2;

        k++;

    }//预处理n,k

    for(int i=0;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));//预处理rev

    fft(a,n,1);

    fft(b,n,1);

    for(int i=0;i<n;i++) ans[i]=a[i]*b[i];

    fft(ans,n,-1);

}

FFT的局限性

显然FFT的过程中会出现浮点数的精度误差。比如在递推单位根$w_n^k=w_n^{k-1}\times w_n^1$时的精度误差比较大，但一般可以接受。如果追求更高精度，可以对于每个$k$重新计算$\omega_n^k=\cos \frac{2k\pi}{n}+\text{i} \sin \frac{2k\pi}{n}$,但时间开销较大。

如果有取模运算，可以使用第二节中提到的NTT算法

另外，在第三节中，我们会介绍减少FFT常数的几种方法，以及对任意长度的序列做DFT的方法。

例题

[1.BZOJ 3527] [ZJOI2014]力(FFT)

2.[Codeforces 580D]Fizzy Search(FFT)

码农公寓