Fast Fourier Transform

快速傅里叶变换(FFT)

这几天在家没事慢慢研究了一下集训最后的一个专题FFT, 陈犇讲完还是有点不是很懂, 网上的很多资料都讲的不是很清晰, 由于这个算法太美了, 我有必要写点什么东西表示对前人智慧的赞美! (顺便练练打公式)

1.FFT这玩意是什么?有什么用?

FFT ( fast Fourier transform )的中文名是快速傅里叶变换, 在信号领域有很大的作用, 在ACM中我们常用FFT来求解多项式乘法

举个例子:

$y_1=2x^4+3x^2+5x+7$
$y_2=5x^3+7x^2+11$
FFT就是用来快速地计算$y_1\cdot y_2$所得到的多项式
($y1\cdot y2=10x^7+14x^6+15x^5+68x^4+70x^3+82x^2+55x+77$)

2.FFT预备知识

1.多项式

1.多项式定义

$A(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-2}x^{n-2}+a_{n-1}x^{n-1}$

$$A(x)=\sum_{i=0}^{n-1}a_ix^i$$

2.多项式表达方式

• 系数表达:$(a_0,a_1,a_2,\cdots,a_{n-2},a_{n-1})$
• 点值表达:选取$n$个(或以上)不同的值:$x_0,x_1,...x_{n-2},x_{n-1}$, 对多项式$A(x)$求值得到$A(x_0),A(x_1),···,A(x_{n-1})$
  那么就称$\{(x_i,A(x_i)),0\le i\lt n,i\in\mathbb Z\}$为多项式$A(x)$的点值表示

举个例子:

$f(x)=2+3x+x^2$这个多项式的系数表达式就是$(2,3,1)$, 点值表达式就是$\{(0,2),(1,6),(2,12)\}$或者$\{(1,6),(2,12),(3,20),(4,30)\}$
$n-1$次多项式只要选取任意不同的$n$个点(或以上)都是可以的

显然不管是点值表达式还是系数表达式都可以唯一确定一个$n-1$次多项式, 并且它们是可以相互推导的

3.多项式的乘法

• 系数表达式下的乘法
  $$A(x)=\sum_{i=0}^{n-1}a_ix^i,B(x)=\sum_{i=0}^{n-1}b_ix^i$$
  那么$A(x)\cdot B(x)$就可以表示为
  $$C(x)=\sum_{i=0}^{2n-2}(\sum_{j=0}^ia_jb_{i-j})x^i$$
• 点值表达式下的乘法
  $$\{x_i,C(x_i)\}=\{x_i,A(x_i)\cdot B(x_i)\}$$

举个例子:

已知$f_1(x)=x^2+1,f_2(x)=x+3$
$f_1(x)$的点值表达式是$\{(0,1),(1,2),(2,5),(3,10)\}$
$f_2(x)$的点值表达式是$\{(0,3),(1,4),(2,5),(3,6)\}$
那么$f_1(x)\cdot f_2(x)$的点值表达式就是$\{(0,3),(1,8),(2,25),(3,60)\}$

2.单位复根

0.引入

先不考虑什么叫做单位复根
先在平面直角坐标系下画一个单位圆, 即$\{(x,y)| x^2+y^2=1\}$
现在我们要在这个圆上取$n$个点$(n>0),$使得这$n$个点恰好把圆弧分成了$n$等份, 规定第$1$个点一定是$(1,0)$

举个例子:

当$n=3$时, 这些点的坐标按照逆时针顺序分别是

$$(1,0),(-\frac{1}{2},\frac{\sqrt 3}{2}),(-\frac{1}{2},-\frac{\sqrt 3}{2})$$
当 $n=4$时, 这些点的坐标按照 逆时针顺序分别是
$$(1,0),(0,1),(-1,0),(0,-1)$$
当 $n=8$时, 这些点的坐标按照 逆时针顺序分别是
$$(1,0),(\frac{\sqrt 2}{2},\frac{\sqrt 2}{2}),(0,1),(-\frac{\sqrt 2}{2},\frac{\sqrt 2}{2}),$$
$$(-1,0),(-\frac{\sqrt 2}{2},-\frac{\sqrt 2}{2}),(0,1),(\frac{\sqrt 2}{2},-\frac{\sqrt 2}{2})$$
如果我们把 平面直角坐标系看做是 复平面, 那么上面那些 点就变成了一个 复数, 它们正是 单位复根

1.单位复根定义

方程$x^n=1$在复数域内一定有$n$个不同的解
那么规定: $w_n^k$是方程$x^n=1$的第$k+1$个根$(0\le k\lt n)$
也就是对应了上面我们从圆上逆时针选取的$n$个点中的第$k$个点

举个例子:

$w_2^0=1$$(x^2=1$的第$1$个根是 $1)$
$w_2^1=1$$(x^2=1$的第$2$个根是 $-1)$

$w_4^0=1$$(x^4=1$的第$1$个根是 $1)$
$w_4^1=i$$(x^4=1$的第$2$个根是 $i)$
$w_4^2=1$$(x^4=1$的第$3$个根是 $-1)$
$w_4^3=i$$(x^4=1$的第$4$个根是 $-i)$

$w_8^0=1$$(x^8=1$的第$1$个根是 $1)$
$w_8^1=i$$(x^8=1$的第$2$个根是 $\frac{\sqrt 2}{2}+\frac{\sqrt 2}{2}i)$
$w_8^2=1$$(x^8=1$的第$3$个根是 $i)$
$w_8^3=i$$(x^8=1$的第$4$个根是 $-\frac{\sqrt 2}{2}+\frac{\sqrt 2}{2}i)$
$w_8^4=1$$(x^8=1$的第$5$个根是 $-1)$
$w_8^5=i$$(x^8=1$的第$6$个根是 $-\frac{\sqrt 2}{2}-\frac{\sqrt 2}{2}i)$
$w_8^6=1$$(x^8=1$的第$7$个根是 $-i)$
$w_8^7=i$$(x^8=1$的第$8$个根是 $\frac{\sqrt 2}{2}-\frac{\sqrt 2}{2}i)$

我们不难归纳得到:

$$w_n^k=cos(\frac{2\pi k}{n})+sin(\frac{2\pi k}{n})i$$
规定: 当 $k\ge n$时, $w_n^k=w_n^{k\%n}$
例如: $w_4^6=w_4^2,w_8^{20}=w_8^4$

2.单位复根性质

下面的只讨论$n$是$2$的次方时,即: $n=2^m,(m\ge 0)$时的性质
画图我们可以很容易的证明

$$(w_n^k)^2=w_{n/2}^k$$
例如: $(w_8^3)^2=w_4^3$, $(w_8^6)^2=w_4^6=w_4^2$

3.FFT原理

用一句话来说就是:将原多项式的系数表达式转换成点值表达式$\bigl( O(nlog_2n)\bigr)$, 再将多项式的点值表达式相乘得到目标表达式的点值表达式$\bigl( O(n)\bigr)$, 再将点值表达式转换成为系数表达式$\bigl( O(nlog_2n)\bigr)$

为了能够在$O(nlog_2n)$时间内得到多项式的点值表达式, FFT采取了一些小技巧, 选择$n$个单位复根进行求值, 由于单位复根的一些性质, 我们可以在$O(nlog_2n)$求出它们的值

现在我们有一个最高次幂是$m$的多项式:

$$A(x)=a_0+a_1x+a_2x+\cdots+a_{m-1}x^{m-1}+a_mx^m$$
1.先将这个多项式的最高次数扩充到$2$的次方$n-1$

$$\begin{align} A(x)&=a_0+a_1x+a_2x+\cdots+a_mx^m+0\cdot x^{m+1}+0\cdot x^{m+2}+\cdots+0\cdot x^{n-1}(n-1\ge m)\\ &=a_0+a_1x+a_2x^2+\cdots+a_{n-2}x^{n-2}+a_{n-1}x^{n-1} \end{align}$$
2.分别带入 $n$个不同的单位复根$(w_n^0,w_n^1,w_n^2,...,w_n^{n-1})$
$$\begin{cases} A(w_n^0)=a_0+a_1w_n^0+\cdots+a_{n-1}(w_n^0)^{n-1}\\ A(w_n^1)=a_0+a_1w_n^1+\cdots+a_{n-1}(w_n^1)^{n-1}\\ A(w_n^2)=a_0+a_1w_n^2+\cdots+a_{n-1}(w_n^2)^{n-1}\\ \cdots \cdots\\ A(w_n^{n-2})=a_0+a_1w_n^{n-2}\cdots+a_{n-1}(w_n^{n-2})^{n-1}\\ A(w_n^{n-1})=a_0+a_1w_n^{n-1}\cdots+a_{n-1}(w_n^{n-1})^{n-1}\\ \end{cases}$$
可以写成一个矩阵:
$$\left(\begin{matrix} A(w_n^0)\\A(w_n^1)\\A(w_n^2)\\\vdots\\A(w_n^{n-2})\\A(w_n^{n-1}) \end{matrix}\right)= \left(\begin{matrix} 1&w_n^0&(w_n^0)^2&\cdots&(w_n^0)^{n-2}&(w_n^0)^{n-1}&\\ 1&w_n^1&(w_n^1)^2&\cdots&(w_n^1)^{n-2}&(w_n^1)^{n-1}&\\ 1&w_n^2&(w_n^2)^2&\cdots&(w_n^2)^{n-2}&(w_n^2)^{n-1}&\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 1&w_n^{n-2}&(w_n^{n-2})^2&\cdots&(w_n^{n-2})^{n-2}&(w_n^{n-2})^{n-1}&\\ 1&w_n^{n-1}&(w_n^{n-1})^2&\cdots&(w_n^{n-1})^{n-2}&(w_n^{n-1})^{n-1}&\\ \end{matrix}\right) \left(\begin{matrix} a_0\\a_1\\a_2\\\vdots\\a_{n-2}\\a_{n-1} \end{matrix}\right)$$
没错我们就是要在$O(nlog_2n)$时间内得到$A(w_n^0),A(w_n^1),A(w_n^2),\cdots,A(w_n^{n-2}),A(w_n^{n-1})$这$n$个值 ! ! !
3.单独考虑其中任意一个方程
$$A(w_n^k)=a_0+a_1w_n^k+a_2(w_n^k)^2+\cdots+a_{n-1}(w_n^k)^{n-1}$$
3.将奇数项和偶数项分开
得到
$$\begin{align} A(w_n^k)&=\bigl(a_0+a_2(w_n^k)^2+\cdots+a_{n-2}(w_n^k)^{n-2}\bigr)+\bigl(a_1w_n^k+a_3(w_n^k)^3+\cdots+a_{n-1}(w_n^k)^{n-1}\bigr)\\ &=\bigl(a_0+a_2(w_n^k)^2+\cdots+a_{n-2}(w_n^k)^{n-2}\bigr)+w_n^k\bigl(a_1+a_3(w_n^k)^2+\cdots+a_{n-1}(w_n^k)^{n-2}\bigr)\\ &=\bigl(a_0+a_2w_{n/2}^k+\cdots+a_{n-2}(w_{n/2}^k)^{\frac{n}{2}-1}\bigr)+w_n^k\bigl(a_1+a_3w_{n/2}^k+\cdots+a_{n-1}(w_{n/2}^k)^{\frac{n}{2}-1}\bigr)\\ \end{align}$$
4.假设两个新的多项式:
$A_0(x)=a_0+a_2x+\cdots+a_{n-2}x^{\frac{n}{2}-1}$
$A_1(x)=a_1+a_3x+\cdots+a_{n-1}x^{\frac{n}{2}-1}$
那么$A(w_n^k)=A_0(w_{n/2}^k)+w_n^kA_1(w_{n/2}^k)$
我们只要再次用同样的方法计算$A_0(x)$以及$A_1(x)$就可以在$O(nlog_2n)$时间内求出$A(x)$的所有点值表达式了

举个例子:

多项式$f(x)=a_0+a_1x+a_2x^2+a_3x^3$, 我们要计算$A(w_4^0),A(w_4^1),A(w_4^2),A(w_4^3)$这4个点的值
1.
$A(w_4^0)=A_0(w_2^0)+w_4^0A_1(w_2^0)$
$A(w_4^1)=A_0(w_2^1)+w_4^1A_1(w_2^1)$
$A(w_4^2)=A_0(w_2^0)+w_4^2A_1(w_2^0)$
$A(w_4^3)=A_0(w_2^1)+w_4^3A_1(w_2^1)$
2.
$A(w_2^0)=A_0(w_1^0)+w_2^0A_1(w_1^0)$
$A(w_2^1)=A_0(w_1^0)+w_2^1A_1(w_1^0)$
3.
因为$w_1^0=1$所以可以直接得到它的值, 也就是当前的系数

上述过程就称为 离散傅里叶变换(Discrete Fourier Transform)

我们将多个多项式的点值表达式相乘, 获得了目标表达式的点值表达式
那么: 怎么获得目标表达式的系数表达式呢?

下面我们将要介绍 逆离散傅里叶变换(Inverse Discrete Fourier Transform)

$$\left(\begin{matrix} A(x_0)\\A(x_1)\\A(x_2)\\\vdots\\A(x_{n-2})\\A(x_{n-1}) \end{matrix}\right)= \left(\begin{matrix} 1&x_0&(x_0)^2&\cdots&(x_0)^{n-2}&(x_0)^{n-1}&\\ 1&x_1&(x_1)^2&\cdots&(x_1)^{n-2}&(x_1)^{n-1}&\\ 1&x_2&(x_2)^2&\cdots&(x_2)^{n-2}&(x_2)^{n-1}&\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 1&x_{n-2}&(x_{n-2})^2&\cdots&(x_{n-2})^{n-2}&(x_{n-2})^{n-1}&\\ 1&x_{n-1}&(x_{n-1})^2&\cdots&(x_{n-1})^{n-2}&(x_{n-1})^{n-1}&\\ \end{matrix}\right) \left(\begin{matrix} a_0\\a_1\\a_2\\\vdots\\a_{n-2}\\a_{n-1} \end{matrix}\right)$$

上面的矩阵就是$DFT$的过程, 下面的矩阵就是$IDFT$的过程

$$\left(\begin{matrix} a_0\\a_1\\a_2\\\vdots\\a_{n-2}\\a_{n-1} \end{matrix}\right)=\frac{1}{n} \left(\begin{matrix} 1&x_0^{-1}&(x_0)^{-2}&\cdots&(x_0)^{-n+2}&(x_0)^{-n+1}&\\ 1&x_1^{-1}&(x_1)^{-2}&\cdots&(x_1)^{-n+2}&(x_1)^{-n+1}&\\ 1&x_2^{-1}&(x_2)^{-2}&\cdots&(x_2)^{-n+2}&(x_2)^{-n+1}&\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 1&x_{n-2}^{-1}&(x_{n-2})^{-2}&\cdots&(x_{n-2})^{-n+2}&(x_{n-2})^{-n+1}&\\ 1&x_{n-1}^{-1}&(x_{n-1})^{-2}&\cdots&(x_{n-1})^{-n+2}&(x_{n-1})^{-n+1}&\\ \end{matrix}\right) \left(\begin{matrix} A(x_0)\\A(x_1)\\A(x_2)\\\vdots\\A(x_{n-2})\\A(x_{n-1}) \end{matrix}\right)$$
其实差异并不大, 我们只需要将原来DFT中$w_n^k$换成$w_n^{-k}$,结果再乘上$\frac{1}{n}$就从$DFT$变成了$IDFT$