FFT入门

据说这玩意算数学?

学,我学还不行么….

基础知识

虚数初步

定义

设 \(a、b\) 为实数,\(i ^ 2 = -1\),形如 \(a + bi\) 的数叫做复数,其中\(i\)被称为虚数单位。复数域是已知最大的域。

虚数由虚部和实部组成.表示一个虚数\(x\)可以表示为\(ai+b\),这个虚数不是一个纯虚数,其中\(ai\)称为数\(x\)的虚部,\(b\)称为这个数\(x\)的实部.

Menci写的很到位QAQ我就直接引用了.

单位根

几何意义:在复平面上,以原点为圆心,\(1\) 为半径作圆,所得的圆叫做单位圆。以原点为起点,单位圆的 \(n\) 等分点为终点,作 \(n\) 个向量。设所得的幅角为正且最小的向量对应的复数为 \(\omega_n\),称为 \(n\) 次单位根。

\(n\)均为\(2\)的整数次幂

因为复数乘法的原则(模长相乘,幅角相加),这些向量的单位根分别为\(\omega_n^2,\omega_n^3,\omega_n^4,…,\omega_n^n\),特别的,\(\omega_n^n=\omega_n^0=1\)

也可以理解为,\(n\)次单位根指的是满足\(s^n=1\)的复数,这些数共\(k\)个,分别为\(\omega_n^k,k\in{0-(n-1)}\)

可以套入欧拉公式理解.

\[e^{ik}=cos(k)+isin(k)\]

至于这个欧拉公式,可以通过泰勒公式推出来

分析一下\(e^x,sin(x),cos(x)\)的公式展开

\(\begin{aligned}e^x&=1+x+\frac{1}{2!}x^2+\frac{1}{3!}x^3+……\\sin(x)&=x-\frac{1}{3!}x^3+\frac{1}{5!}x^5+……\\cos(x)&=1-\frac{1}{2!}x^2+\frac{1}{4!}x^4+……\end{aligned} \)

将\( x=i\theta \)带入\(e^x\)中

\( \begin{aligned}e^{i\theta}&=1+i\theta+\frac{(i\theta)^2}{2!}+\frac{(i\theta)^3}{3!}+\frac{(i\theta)^4}{4!}+……\\&=1+i\theta-\frac{\theta^2}{2!}-\frac{i\theta^3}{3!}+\frac{\theta^4}{4!}+\frac{i\theta^5}{5!}-……\\&=(1-\frac{\theta^2}{2!}+\frac{\theta^4}{4!}-\frac{\theta^6}{6!}+\frac{\theta^8}{8!}-……)+i(\theta-\frac{\theta^3}{3!}+\frac{\theta^5}{5!}-\frac{\theta^7}{7!}+……)\\&=cos(\theta)+isin(\theta)\end{aligned} \)

似乎偏题了不管不管

性质

性质一:\( \omega^{2k}_{2n}=\omega^k_n \)

几何意义上就是在复平面中,两个单位根表示的向量的终点相同.

性质二:\( \omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k} \)

\( \begin{aligned}\omega_{n}^{k+\frac{n}{2}} &=\omega_{n}^{k}*\omega_{n}^{\frac{n}{2}}\\&=\omega_{n}^{k}*(cos\frac{n}{2}*\frac{2\pi}{n}+i*sin\frac{n}{2}*\frac{2\pi}{n})\\&=\omega^k_n*(cos\pi+i*sin\pi)\\&=\omega_n^k*-1\\&=-\omega^k_n\end{aligned} \)

多项式乘法

给定两个函数\(A(x),B(x)\),表示关于数\(x\)的多项式

\( \begin{aligned}A(x)&=\sum_{i=0}^{n}a_ix^i\\B(x)&=\sum_{i=0}^nb_ix^i\\C(x)&=A(x)*B(x)\&=\sum_{j+k=i,0\leq j,k\leq n}a_jb_kx^i\end{aligned} \)

直接求解\(C(x)\)的时间复杂度为\(N^2\)的,如果存在\((x,y)\)相同的两个多项式相乘可以直接变成\(O(N)\)的.

点值表示&插值

点值

点值表示的就是简单的类似于函数的对应关系.

对于每一个不同的\(x\)都有一个\(y\)与其对应,不过没有函数那么多奇怪的要求.

也可以理解为一个\(n\)维向量通过点值表示可以表达出一个新的\(n\)维向量.

\[ \begin{aligned}W_i&=A(x_i)\&=\sum_{j=0}^{n-1}a_jx_i^j\end{aligned} \]

插值

拉格朗日插值法详解

DFT&IDFT

FFT的本质是对DFT和IDFT的优化.

DFT就是系数表达转化成点值表示.

IDFT就是DFT的逆,就是讲点值表示转化成系数表达.

什么意思?

DFT:把这个多项式(可以顺带理解为连续的图像)转化成断开的点

IDFT:把这些断开的点还原成一个连续的图像(多项式

比如用FFT解决A*B问题.就是类似于求前面的系数可以用\(log n\)的复杂度解决高精度问题qwq.

然而不用FFT这样暴力求多项式复杂度是\(O(N^2)\)的.

这就需要本篇的重点

FFT

FFT

Fast-DFT

对于时间复杂度较高的DFT和IDFT而言,我们考虑从其性质出发从而优化其复杂度.

讨论性质优化

对于多项式\(A\)而言

\( A(x)=a_0+a_1x+a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6…. \)

对其系数\(a\)按奇偶分开.

\( A(x)=(a_0+a_2x^2+a_4x^4+a^6x^6+…+a_{n-1}x^{n-2})+(a_1x+a_3x^3+a_5x^5+a_7x^7+…+a_{n-1}x^{n-1}) \)

令\( B(x)=(a_0+a_2x+a_4x^2+a^6x^4+…+a_{n-2}x^{\frac{n}{2}-1}),\\C(x)=(a_1+a_3x+a_5x^3+a_7x^5+…+a_{n-1}x^{\frac{n}{2}-1}) \)

所以有\(A(x)=B(x^2)+xC(x^2)\)

我们可以尝试求一下\( A(\omega_n^k) \)的值,假设\( k\leq \frac{n}{2} \)

\( \begin{aligned}A(\omega_n^{2k})&=B(\omega_n^{2k})+\omega_n^kC(\omega_n^{2k})\\&=B(\omega_{\frac{n}{2}}^k)+\omega_n^kC(\omega_{\frac{n}{2}}^k)\end{aligned} \)

仅仅是\(\omega_n^k\)是不够的,还需要求一下\(\omega_n^{k+\frac{n}{2}}\)

\( \begin{aligned}A(\omega_n^{k+\frac{n}{2}})&=B(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}C(\omega_n^{2k+n})\\&=B(\omega_n^{2k}*\omega^{n}_{n})-\omega_{n}^{k}C(\omega_n^{2k}*\omega_n^n)\\&=B(\omega_n^{2k})-\omega_n^kC(\omega_n^{2k})\\&=B(\omega^k_{\frac{n}{2}})-\omega_n^kC(\omega_{\frac{n}{2}}^k)\end{aligned} \)

突然发现无论是\(A(\omega_n^k)\ or\ A(\omega_n^{k+\frac{n}{2}})\)最后都可以化简成\(\frac{n}{2}\).

可以通过这个优化复杂度!

刚好我们考虑的两种情况可以取遍\([0,n-1]\)对于的所有值.

确实可以让原来的取值变成一半,分析一下复杂度.

复杂度分析

因为可以直接缩减一半的缘故,所以我们可以用\(O(n)\)的时间从而求出\(A(n)\),而我们只是需要多求个\(B(\frac{n}{2})\ and\ C(\frac{n}{2})\)

所以复杂度\(T(n)\)可以用以下式子表示

\( \begin{aligned}T(n)&=2*T(\frac{n}{2})+O(n)\\&=(4*T(\frac{n}{4})+2O(\frac{n}{2}))+O(n)\\&=(8*T(\frac{n}{8})+4O(\frac{n}{4})+O(n))+O(n)\\&=2^kT(1)+O(n)+…+O(n)+O(n)+O(n)\\&=nO(1)+k*O(n)=(k+1)O(n)\\&=O(nlog_2n)\end{aligned} \)

可以理解为不断把\(T(x)\)拆成\(2*T(\frac{x}{2})+O(x)\)

其余的根据主定理(master theorem)套用即可得到复杂度\(O(nlog_2n)\)辣

Fast-IDFT

化简优化

现在我们需要把这些断开的点还原成多项式,how to do it?

考虑一下之前的那个乘法,其实就是两个大的矩阵乘起来,而我们需要分离出系数,也就需要把其中一个矩阵消掉.

之前的矩阵:

\( \begin{bmatrix}A(x_0)\\A(x_1)\\A(x_2)\\A(x_3)\\……\\A(x_{n-1})\end{bmatrix}=\begin{bmatrix}1&1&1&1&1&…&1\\1&\omega_n^1&\omega_n^2&\omega_n^3&\omega_n^4&…&\omega_n^{n-1} \\1&\omega_n^2&\omega_n^4&\omega_n^6&\omega_n^8&…&\omega_n^{2(n-1)} \\1&\omega_n^3&\omega_n^6&\omega_n^9&\omega_n^{16}&…&\omega_n^{3(n-1)} \\&&&……\\1&\omega_n^{n-1}&\omega_n^{2(n-1)}&\omega_n^{3(n-1)}&\omega_n^{4(n-1)}&…&\omega_n^{(n-1)^2} \end{bmatrix}\begin{bmatrix}a_0\\a_1\\a_2\\a_3\\……\\a_{n-1}\end{bmatrix} \)

我们需要把中间的矩阵消掉,并且这个矩形有一个非常优秀的性质我们可以利用

这个矩阵的中有一种非常特殊的性质,对该矩阵的每一项取倒数,再除以n就可以得到该矩阵的反矩阵

如何取负是一个很骚的操作.

首先我们考虑最开始的欧拉函数.已知 \( e^{i\pi}=-1 , e^{2i\pi}=-1 \) ,假设一个数\( a,a^k=1 \),可以知道这个数\( a=e^\frac{2 \pi}{k} \)(因为 \( {e^\frac{2 \pi}{k}}^k=e^{2 \pi}=1 \).我们需要\(\frac{1}{a}\)的话,可以直接把\( \pi \)取负,所以就很简单辣qwq.

其他操作与DFT相同.

代码优化

显而易见,我们可以得到一个递归版的.

显而易见,这个常数巨大.

我们还是需要优化!

考虑一下递归时的奇偶分组的情况.

考虑一组{0,1,2,3,4,5,6,7}

{0,1,2,3,4,5,6,7}

{0,2,4,6},{1,3,5,7}

{0,4},{2,6},{1,5},{3,7}

{0},{4},{2},{6},{1},{5},{3},{7}

数字 二进制 数字 二进制
0 000 0 000
1 001 4 100
2 010 2 010
3 011 6 110
4 100 1 001
5 101 5 101
6 110 3 011
7 111 7 111

刚好二进制下逆序!

所以我们运用雷德算法,二进制下翻转一下,可以对代码优化,从递归变成迭代!

可以发现改变顺序后的排列是原排列每个数的二进制位逆序置换,如果我们知道了当前位置应该填的数,可以用类似于加法进位的方法求出下一个位置应该填的数:

从最高位开始,循环执行:检查当前位,若为00,将其置为11并退出循环;否则将其置为00,并检查低位;

求得当前位置的应有的下标后,可以与对应位置的数交换,因为二进制位逆序置换的数是一一对应的,所以每个数至多只会被交换一次;

以上过程被称为雷德算法。

代码实现


By:Wahacer

2018.3.10

18:28

发表评论

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