学习笔记 - FFT(快速傅里叶变换) | Lucky_Glass's Blog
0%

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

其实是复习笔记了…… Tab. 第 233 次复习


Part0. FFT 用来干什么

众所周知,用朴素算法实现 多项式乘法 $f(x)\cdot g(x)$:分别枚举 $f(x)$ 和 $g(x)$ 中的一个项,系数相乘、指数相加,得到答案中的一个项(的一部分)。

不妨设 $f(x),g(x)$ 的次数都是 $n$,那么上述的朴素算法就是 $O(n^2)$ 的。这样的时间复杂度是多项式乘法的一个瓶颈……

而 FFT 就是用来解决多项式乘法的问题的,可以把它的时间复杂度优化到 $O(n \log n)$。

实际上,FFT 并不是直接计算多项式乘法,而是把原来的多项式 $f(x),g(x)$ 在 $O(n\log n)$ 的复杂度内转换为它的 点值表示(后面会讲),而点值表示的多项式相乘的时间复杂度是 $O(n)$ 的。最后再用 $O(n\log n)$ 的时间复杂度把所得多项式的点值表示转化为一般形式。


Part1. 点值表示

如果完全没有了解过 FFT 的 reader 们看完 Part0 都有一点懵。那么现在就来科普一下什么是点值表示。

Part1/1. 什么是点值表示

对于一个 $n$ 次的多项式 $f(x)$,我们将 $n+1$ 个不同的 $x_i$ 代入这个 $f(x)$,可以得到 $n+1$ 个多项式的值 $y_i$。根据这 $n+1$ 个数对 $(x_i,y_i)$ 可以唯一确定一个 $n$ 次的多项式。

也就是说一个多项式与一个点值表示是一一对应的

那么 FFT 完成的操作就是:

  • 把已知的一个多项式转化成对应的点值表示;
  • 把已知的点值表示转换成对应的多项式。

复杂度都是 $O(n\log n)$。

Part1/2. 点值表示的乘法

假设我们用 $\{x_0,x_1,x_2,…,x_n\}$ 来表示 $f(x)$ 为 $\{(x_0,f(x_0)),(x_1,f(x_1)),…,(x_n,f(x_n))\}$、表示 $g(x)$ 为 $\{(x_0,g(x_0)),(x_1,g(x_1)),…,(x_n,g(x_n))\}$。

那么 $f(x)\cdot g(x)$ 的点值表达式就是 $\{(x_0,f(x_0)\cdot g(x_0)),(x_1,f(x_1)\cdot g(x_1)),…,(x_n,f(x_n)\cdot g(x_n))\}$ ,其实就是对应项相乘,可见它的时间复杂度是 $O(n)$ 的。


Part2. 算法思想

Part2/1. 小科普:复数

是复数(complex)不是负数(negative)

首先了解到虚数 $i=\sqrt{-1}$。

然后一个复数 $\mathbf a$ 可以表示为实部 $x$ 和虚部 $y\times i$ 之和 —— $\mathbf a=x+y\times i$。对于一个实数,我们可以把它放在“一维数轴”上;那么对于一个复数,我们需要把它放在“二维数轴”,也就是直角坐标系上。

png1

对于复数的乘法,比如 $\mathbf a=x+yi,\mathbf b=x’+y’i$ ,计算 $\mathbf {ab}$:

  1. 可以展开直接相乘:$\mathbf{ab}=(x+yi)(x’+y’i)=(xx’-yy’)+(xy’+x’y)i$ ;
  2. 可以从几何角度理解:
    对于一个“向量”表示的复数,有如下定义(实际上不满足向量的运算法则):
    png2
    那么两个复数相乘等于“幅角相加、模长相乘”

Part2/2. 单位复根

Hint. 复数的模长 $|\mathbf a|$(这个不是绝对值),表示把它表示成向量后的长度,即 $|\mathbf a|=\sqrt{x^2+y^2}$,

定义 n 次单位复根 $\mathbf \omega_n^i$ 为满足 $\mathbf z^n=1$ 的复数 $\mathbf z$。比如 $n=3$ 时,单位复根表现在直角坐标系上为:

png3

可以发现 $n$ 次单位复根满足以下性质($\mathbf \omega_n^i$ 表示第 $i$ 个单位复根):

  • $|\mathbf \omega_n^i|=1$(都在单位圆上);
  • $\mathbf \omega_n^i$ 的幅角为 $\frac {2\pi i}{n}$;Tab. 这里的 $i$ 是下标不是虚数,角度为弧度制($360^\circ$ 弧度角为 $2\pi$)
  • 总共有 $n$ 个单位复根,记为 $\omega_n^0,\omega_n^1,…,\omega_n^{n-1}$。

另外,还有一些比较明显的性质:

  • $\omega_n^i\cdot\omega_n^j=\omega_n^{i+j}$;
  • $\omega_n^{i+n}=\omega_n^i$;
  • $\omega_n^i=\omega_{kn}^{ki}$;
  • $\omega_n^i=-\omega_n^{i+n/2}$ (这里的“负”表示大小相等、方向相反)

以上这些结论都可以通过在单位圆上画出单位根来证明。

单位复根有什么用呢?因为 $n$ 次单位复根恰好有 $n$ 个,就可以把这 $n$ 个单位复根分别代入 $n-1$ 次多项式 $f(x)$ —— 得到点值表达 $\{(\omega_n^0,f(\omega_n^0)),(\omega_n^1,f(\omega_n^1)),…,(\omega_n^{n-1},f(\omega_n^{n-1}))\}$

Part2/3. 傅里叶正变换(一般形式 to 点值表达)

Tab. 以下思路来自 Tiw_Air_OAO 的 cnblogs 博客,但不尽相同:https://www.cnblogs.com/Tiw-Air-OAO/p/10162034.html

FFT 的实现方法是 分治:记多项式 $f(x)=a_0+a_1x+a_2x^2+…+a_{n-1}x^{n-1}$。

接下来把 $f(x)$ 的偶数次项、奇数次项分开:

那么可以得到 $f(x)=x\cdot f_o(x^2)+f_e(x^2)$ (直接展开这个式子就可以证明了)。

令 $p=n/2$,然后有以下结论:

所以,如果已经知道 $f_o(\omega_p^i),f_e(\omega_p^i)$,就可以 $O(1)$ 求出 $f(\omega_n^i)$ 和 $f(\omega_n^{i+p})$。实现方法就是递归:每次从求解 $f(\omega_n^i)$ 的单位复根到求解 $f_o(\omega_{n/2}^i)$ 和 $f_e(\omega_{n/2}^i)$ ;最后当多项式只剩两项(一次项和常数项)时直接代入计算。

由于 $n$ 每次会除以二,而要计算 $n$ 次,所以时间复杂度为 $O(n\log n)$。

Part2/4. 傅里叶逆变换(点值表达 to 一般形式)

其实傅里叶正变换是作的这样一个矩阵乘法:

(矩阵乘法是第一个矩阵的第 i 行依次乘第二个矩阵的第 j 列,得到结果的第 i 行、第 j 列)

记上面的第一个矩阵(系数矩阵)为 $V$。我们再定义一个矩阵 $D$:

于是计算 $D\times V$ 的第 $(i,j)$ 项:

① 当 $i=j$ 时:

② 当 $i\not=j$ 时:

Hint. 用等比数列……

因为 $\omega_n^{(j-i)n}=\omega_n^0=1$,所以 $(D\times V)_{i,j}=0$,于是发现:

所以我们把点值表达再乘上这个矩阵 $D$(就像正变换一样),最后再给每个值除以 $n$。


Part3. 离散傅里叶变换实现

虽然之前讲的算法思想是递归求解,然而实测的效率很低……于是改成 迭代(递推)

先来研究一下递归的时候系数是被怎么分的,比如 n=8 时(下面的数字代表的是多项式 $f(x)$ 的系数,如 “2” 代表的是 $f(x)$ 的第 2 项的系数 $a_2$):

png4

现在还看不出规律,换成二进制再看看:

png5

原来的系数编号的顺序是 {000,001,010,011,100,101,110,111},最后得到的编号顺序为 {000,100,010,110,001,101,011,111}……惊奇地发现每个位置上的数字的编号都反过来了(比如原来的第 4 个系数 011 变成了 110)。所以只需要做一个二进制翻转即可。


Part3. 源代码

版子题:UOJ #34

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
/*Lucky_Glass*/
#include<cmath>
#include<vector>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;

struct CComplex{
double Vr,Vi;
CComplex(double Ur=0,double Ui=0):Vr(Ur),Vi(Ui){}
};
CComplex operator +(CComplex Va,CComplex Vb){return CComplex(Va.Vr+Vb.Vr,Va.Vi+Vb.Vi);}
CComplex operator -(CComplex Va,CComplex Vb){return CComplex(Va.Vr-Vb.Vr,Va.Vi-Vb.Vi);}
CComplex operator *(CComplex Va,CComplex Vb){return CComplex(Va.Vr*Vb.Vr-Va.Vi*Vb.Vi,Va.Vr*Vb.Vi+Va.Vi*Vb.Vr);}
CComplex operator /(CComplex Va,double Vb){return CComplex(Va.Vr/Vb,Va.Vi/Vb);}

typedef CComplex cplx;
const double CPI=acos(-1);

namespace FFT{
void Change(cplx *Ci,int len,int tmp){
int lglen=log2(len);
for(int i=0;i<len;i++){
int Revi=0;
for(int j=0;j<lglen;j++) Revi|=(((i>>j)&1)<<(lglen-1-j));
if(Revi<i) swap(Ci[Revi],Ci[i]);
}
for(int s=2;s<=len;s<<=1){
int t=s/2;
cplx Vper{cos(tmp*CPI/t),sin(tmp*CPI/t)};
for(int i=0;i<len;i+=s){
cplx Vr{1,0};
for(int j=0;j<t;j++,Vr=Vr*Vper){
cplx Ve=Ci[i+j],Vo=Ci[i+j+t];
Ci[i+j]=Ve+Vo*Vr;
Ci[i+j+t]=Ve-Vo*Vr;
}
}
}
if(tmp==-1)
for(int i=0;i<=len;i++)
Ci[i]=Ci[i]/(double)len;
}
}

const int N=1e6;
int Na,Nb,Nc,len;
double Ii;
cplx Ca[N+3],Cb[N+3],Cc[N+3];

int main(){
scanf("%d%d",&Na,&Nb);
for(int i=0;i<=Na;i++) scanf("%lf",&Ca[i].Vr);
for(int i=0;i<=Nb;i++) scanf("%lf",&Cb[i].Vr);
Nc=Na+Nb;
len=1<<(int)ceil(log2(Nc));
if(len==Nc) len<<=1;
FFT::Change(Ca,len,1);FFT::Change(Cb,len,1);
for(int i=0;i<=len;i++) Cc[i]=Ca[i]*Cb[i];
FFT::Change(Cc,len,-1);
for(int i=0;i<=Nc;i++){
printf("%d",(int)round(Cc[i].Vr));
if(i==Nc) printf("\n");
else printf(" ");
}
return 0;
}

The End

Thanks for reading~

Email: lucky_glass@foxmail.com,欢迎提问~