首页 > 笔记 > FFT总结

FFT总结

前置知识

多项式定义

以$x$为变量的多项式是定义在一个代数域$\mathbb{F}$上,将函数$A(x)$表示为形式和

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

则我们称$a_0,a_1,a_2 \cdots a_{n-1}$为多项式的系数。所有的系数都属于$\mathbb{F}$,比如复数域$\mathbb{C}$。

如果一个多项式$A(x)$的最高次项系数为$a_k$,那我们称$A(x)$的次数为$k$,记作$degree(A)=k$

多项式运算

多项式加法

这个是非常好想的东西,就是把每一项合起来

$$C(x)=A(x)+B(x)=\sum_{i=0}^{n-1}a_ix^i+\sum_{j=0}^{n-1}b_jx^j=\sum_{k=0}^{n-1}c_kx^k$$

多项式乘法

方法是将$A(x)$中的每一项和$B(x)$中的每一项相乘,然后同类项合并。比如:

也就是

$$C(x)=\sum_{i=0}^{2n-2}c_ix^j \ and \ c_i=\sum_{j=0}^ia_jb_{i-j}$$

显然

$$degree(C)=degree(A)+degree(B)$$

多项式表示

上面的计算中多项式都是用系数来表示多项式。但是我们发现这样的方法乘法的复杂度就是$O(n^2)$的。非常的慢。

我们还可以换一种方式来考虑,用函数的方法来思考。对于一次函数,我们用两个点就可以确定它的多项式。

然后我们就可以用点值来表示一个多项式。我们用$n$个不同的点来表示的话:

$${(x_0,y_0),(x_1,y_1),\cdots ,(x_n,y_n)}$$

表示的多项式也是唯一的。

加法的时候直接把每个点加上:

$${(x_0,y_0+y'_0),(x_1,y_1+y’_1),\cdots (x_{2n-1},y_{2n-1}+y'_{2n-1})}$$

是$O(n)$的。

考虑乘法的时候,我们必须面对上面的那个问题$degree(C)=degree(A)+degree(B)$。也就是说,我们必须找$2n$个点让它构成多项式。所以我们需要对$A$,$B$的点值进行拓展。给定$A$的拓展式:

$${(x_0,y_0),(x_1,y_1),\cdots (x_{2n-1},y_{2n-1})}$$

还有$B$的拓展式:

$${(x_0,y'_0),(x_1,y’_1),\cdots (x_{2n-1},y'_{2n-1})}$$

这样就可以得到$C$:

$${(x_0,y_0y'_0),(x_1,y_1y’_1),\cdots (x_{2n-1},y_{2n-1}y'_{2n-1})}$$

也是$O(n)$的。

复数

复数的定义

复数,为实数的延伸,它使任一多项式方程式都有根。

符号表示

通常表示成$a+bi$的形式,$a,b$为实数。且$\sqrt{-1}=i$

它的实数$a$为它的实部,$b$为它的虚部。

所有复数的集合叫$\mathbb{C}$,我们可以认为实数集$\mathbb{R}$是它的一个子集。仅当$a=a+0i$

复平面

复平面,是x轴表示实数与y轴(除了原点)表示虚数,建立起来的复数的几何表示。

每一个复数$a+bi$都能表示为一个$(0,0)$指向$(a,b)$的向量。

向量的模长为$\sqrt{a^2+b^2}$就是它的长度。表示从x轴正半轴到该向量的转角的有向(以逆时针为正方向)角叫做幅角。

复数相加遵循平行四边形定则。复数相乘时,模长相乘,幅角相加。

欧拉公式

内容

它的主要内容是对于任意的x,下列等式恒成立

$$e^{ix}=\cos x+i \sin x$$

证明

对于$x \in I$,定义$f(x)=\frac{\cos x+i\sin x}{e^{ix}}$。对它求导

$$
\begin{align*}
f'(x)&= \frac{(-\sin x + i \cos x) \cdot e^{ix} - (\cos x+i \sin x) \cdot i \cdot e^{ix}}{(e^{ix})^2} \\
&= \frac{ -\sin x \cdot e^{ix} - i^2 \sin x \cdot e^{ix}}{(e^{ix})^2}\\
&= \frac{ -\sin x \cdot e^{ix} + \sin x \cdot e^{ix}}{(e^{ix})^2}\\
&= 0
\end{align*}
$$

设$[a,b] \in I,c \in (a,b)$

则由拉格朗日中值定理得:

$$f'(c)=\frac{f(b)-f(a)}{b-a}$$

由上文得:

$$f'(c)=0$$

$$f(a)=f(b)$$

所以此函数为常数函数

$$f(x)=f(0),\frac{\cos x+i\sin x}{e^{ix}}=\frac{\cos 0+i\sin 0}{e^{0}}=1$$

所以

$$e^{ix}=\cos x+i \sin x$$

证毕。

单位根

定义

$n$次单位根是$n$次幂为$1$的复数。它们位于复平面的单位圆上,构成正$n$边形的顶点,其中一个顶点是$1$。

设所得的幅角为正且最小的向量对应的复数为$ \omega_n$? ,称为 $n$ 次单位根。

由复数乘法的定义(模长相乘,幅角相加)可知,其与的 $n - 1$ 个向量对应的复数分别为$\omega_n^1,\omega_n^2 \cdots \omega_n^n$,其中$\omega_n^n=\omega_n^0=1$

因为这是$n$个点分布在单位圆上,且每个幅角对应的是周角的$k \over n$,由上面的欧拉公式得,

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

性质
性质1 (消去引理)

对于$\forall n,k,d \in \mathbb{Z},n,k \ge 0,d>0$,可得

$$\large \omega_{dn}^{dk}=\omega_n^k$$

证明:

$$\large \omega_{dn}^{dk} = (e^{\frac{2\pi i}{dn}})^{dk} = (e^{2\pi i/n})^{k} = \omega_n^k$$

性质2 (折半引理)

对于$\forall n,k \in \mathbb{Z},n > 0且为偶数$,可得

$$\large \omega_n ^ { k + \frac{n}{2} } = -\omega_n ^ k $$

证明

$$\large \omega_n ^ \frac{n}{2} \omega_n ^ k = (\cos \frac{n}{2} \cdot \frac{2\pi}{n} + i \sin \frac{n}{2} \cdot \frac{2\pi}{2}) \times \omega_n ^ k = (\cos \pi + i \sin \pi) \times \omega_n ^ k= -\omega_n ^ k $$

性质3 (求和引理)

对于$\forall n,k \in \mathbb{Z},n ≥ 1,k不能被n整除且k\not=0 $,可得

$$\sum_{i=0}^{n-1}(\omega_n^k)^i=0$$

证明

先证明一下几何级数求和的性质

几何级数求和:

$$\sum_{i=0}^nx^i=\frac{x^{n+1}-1}{x-1}$$

证明

$$\sum_{i=0}^nx^i=1+x+x^2+\cdots+x^n \cdots \cdots \cdots \cdots(1)$$

$$x\sum_{i=0}^nx^i=x+x^2+x^3+\cdots+x^{n+1} \cdots \cdots \cdots \cdots(2)$$

$$(2)-(1)=(x-1)\sum_{i=0}^nx^i=x^{n+1}-1$$

$$\sum_{i=0}^nx^i=\frac{x^{n+1}-1}{x-1}$$

所以对于上式

$$\sum_{i=0}^{n-1}(\omega_n^k)^i=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}=\frac{(\omega_n^n)^k-1}{\omega_n^k-1}=\frac{1^k-1}{\omega_n^k-1}=0$$

我们也可以考虑一下$k=0$的情况,此时$\omega_n^k=1$所以

$$\sum_{i=0}^{n-1}(\omega_n^k)^i=\sum_{i=0}^{n-1}1^i=n\ \ \ \ (k=0)$$

FFT思想

在上文,我们看到了假如用点值表示的多项式乘法可以在$O(n)$的时间里做完。FFT加速乘法就是我们先把多项式由系数表示转为点值表示,然后乘完再转回去。

FFT过程

DFT

我们希望计算

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

在$\omega_0,\omega_1,\cdots \omega_{n-1}$处的值($n$个$n$次单位复数根处)。假如$A$由系数形式给出$a=(a_0,a_1,\cdots a_{n-1})$,接下来对于$k=0,1,\cdots n-1$,定义结果$y_k$

$$y_k=A(\omega_n^k)=\sum_{i=0}^{n-1}a_i\omega_n^{ki}$$

则向量$y=(y_0,y_1,\cdots y_{n-1})$就是系数向量$a=(a_0,a_1,\cdots a_{n-1})$的离散傅里叶变换,记作$y=DFT_n(a)$

按照朴素算法来求离散傅里叶变换,时间复杂度仍然为$O(n^2)$

FFT

考虑将多项式按照系数下标的奇偶分为两部分

$$A(x)=(a_0+a_2x^2+a_4x^4+\cdots +a_{n-2}x^{n-2})+(a_1x+a_3x^3+a_5x^5+\cdots +a_{n-1}x^{n-1})$$

$$
A_1(x)=a_0+a_2x+a_4x^2+\cdots +a_{n-2}x^{\frac{n}{2}-1}\\
A_2(x)=a_1+a_3x+a_5x^2+\cdots +a_{n-1}x^{\frac{n}{2}-1}
$$

显然

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

假如$k<\frac{n}{2}$,我们要求$A(\omega_n^k)$

$$
\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)
\end{align*}
$$

这一步用到了消去引理

假如$k>\frac{n}{2}$,我们把$k$分一下$A(\omega_n^{k+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_n^{2k} \times \omega_n^n)-\omega_n^kA_2(\omega_n^{2k} \times \omega_n^n)\\
&= 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)
\end{align*}
$$

这一步我们用到了折半引理,消去引理,还有$\omega_n^n=1$

通过这样的转化,我们发现了当$k \in [0,\frac{n}{2}-1]$的时候,$k,k+\frac{n}{2}$取遍了$[0,n-1]$。

也就是说,如果已知 $A_1(x)$ 和 $A_2(x)$ 在 $\omega_\frac{n}{2}^0,\omega_\frac{n}{2}^1,\cdots \omega_\frac{n}{2}^{\frac{n}{2}-1}$ 处的点值,就可以在$O(n)$的时间内求得$ A(x) $在$\omega_n^0,\omega_n^1,\cdots \omega_n^{n-1}$处的取值。而关于$ A_1(x) $和$ A_2(x) $的问题都是相对于原问题规模缩小了一半的子问题。我们可以递归得跑下去。

然后复杂度就是

$$T(n)=2T(\frac{n}{2})+O(n)=O(n \log n)$$

标准的主定理。。。

逆FFT

我们考虑一个$(y_0,y_1,y_2,\cdots,y_{n-1})$,为$(a_0,a_1,a_2,\cdots,a_{n-1})$的离散傅里叶变换。再考虑另一个向量$(c_0,c_1,c_2,\cdots,c_{n-1})$满足

$$c_k=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i$$

即多项式$B(x)=y_0+y_1x+y_2x+\cdots+y_{n-1}x$在$\omega_n^0,\omega_n^{-1},\omega_n^{-2},\cdots,\omega_n^{-{n-1}}$处的点值表示。

展开一下:
$$
\begin{align*}
c_k &= \sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i \\
&= \sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^i)^j)(\omega_n^{-k})^i\\
&= \sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^j)^i)(\omega_n^{-k})^i\\
&= \sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^j)^i(\omega_n^{-k})^i)\\
&= \sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^{j-k})^i\\
&= \sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i\\
\end{align*}
$$
根据求和引理,我们设

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

上式就可以化为

$$c_k=\sum_{j=0}^{n-1}a_jS(\omega_n^{j-k})$$

显然,当$k=j$时$S(\omega_n^{j-k})=n$,否则$S(\omega_n^{j-k})=0$。即

$$c_i=na_i,a_i=\frac{1}{n}c_i$$

所以,使用单位根的倒数代替单位根,做一次类似快速傅里叶变换的过程,再将结果每个数除以$n$,即为傅里叶逆变换的结果。

FFT实现

递归实现

例题

蝴蝶操作

~~其实我不知道有啥用~~

就是我们看一下合并的时候

发现w*a2[i]计算了两次,于是

实际优化其实也不算多。。

反正算导上写了。。。。

迭代实现

我们观察一下递归的序列

发现递归序列相当于原序列的二进制翻转。

我们就可以模拟一个二进制加法,但是是反向的二进制加法。

考虑正向二进制加法的过程,相当于从最低位开始,找到第一个 0,然后把这个 0 改成 1,之前的 1 全部变成 0。那么反向二进制加法就是从最高位开始,找到第一个 0,然后把这个 0 改成 1,前面的 1 全部改成 0。

然后我们操作的时候,就直接把原序列变成递归序列,然后更新就可以了。

例题

效率变化非常明显:

例题

大整数乘法

cv3123

其实多项式乘法就是模拟了一个竖式乘法的过程。改动是需要处理一下进位。

喂丸待


4 COMMENTS

  1. attack2017-12-09 上午8:59

    厉害厉害,请问什么时候能补完剩下的呢?

如果你觉的这篇文章不错,分享给朋友吧!

打开微信“扫一扫”,打开网页后点击屏幕右上角分享按钮

×