快速傅里叶变换(FFT)与多项式乘法学习笔记 | Bill Yang's Blog
0%

快速傅里叶变换(FFT)与多项式乘法学习笔记

本文章基于Menci’s Blog进行补充,原文已经讲述得很清楚,有数学基础的同学可以直接看原文。
快速傅里叶变换(Fast Fourier Transform,简称FFT),可以在$O(n\log n)$时间内完成离散傅里叶变换(Discrete Fourier Transform,简称DFT),是在数论中常用的加速多项式乘法的计算,常用于快速计算卷积。

定义

多项式

我们将一个多项式$\sum\limits_{i=0}^{n-1}a_ix^i$记为$A(x)$。

多项式的次数

众所周知,多项式的次数为其最高项的次数。
$A(x)$的次数为$n-1$。

多项式的系数表示

对于一个多项式$A(x)$,我们将其系数写作一个向量$(a_0,a_1,a_2,\ldots,a_{n-1})$。
根据一个$n$维的系数向量可以唯一确定一个$n-1$次的多项式。
称向量$(a_0,a_1,a_2,\ldots,a_{n-1})$为多项式$A(x)$的系数表示。

多项式的点值表示

将一组互不相同的$n$个$x_i$代入多项式$A(x)$,得到$n$个值$y_i$。
称二元组向量$(\lbrace x_0,y_0\rbrace,\lbrace x_1,y_1\rbrace,\lbrace x_2,y_2\rbrace,\ldots,\lbrace x_{n-1},y_{n-1}\rbrace)$为多项式$A(x)$的点值表示。
在后文中为了叙述方便,简单地将$(y_0,y_1,y_2,\ldots,y_{n-1})$也称作点值表示。

求值与插值

定理:一个$n-1$次多项式在$n$个不同点的取值唯一确定了该多项式。

证明:假设命题不成立,存在两个不同的$n-1$次多项式$A(x),B(x)$,满足对于任何$i\in[0,n-1]$,有$A(x_i)=B(x_i)$。
令$C(x)=A(x)-B(x)$,则$C(x)$也是一个$n-1$次多项式。对于任何$i\in[0,n-1]$,有 $C(x_i)=0$。
即$C(x)$有$n$个根,这与代数基本定理(一个$n-1$次多项式在复数域上有且仅有$n-1$个根)相矛盾,故$C(x)$并不是一个$n-1$次多项式,原命题成立,证毕。

如果我们按照定义求一个多项式的点值表示,时间复杂度为$O(n^2)$。
已知多项式的点值表示,求其系数表示,可以使用插值。朴素的插值算法时间复杂度为$O(n^2)$。

多项式乘法

我们定义两个多项式$A(x)=\sum\limits_{i=0}^{n-1}a_ix^i$与$B(x)=\sum\limits_{i=0}^{n-1}b_ix^i$相乘的结果为$C(x)$(假设两个多项式次数相同,若不同可在后面补零)。

两个$n-1$次多项式相乘,得到的是一个$2n-2$次多项式,时间复杂度为$O(n^2)$。

如果使用两个多项式在$2n-1$个点处取得的点值表示,那么

时间复杂度为$O(n)$。

故我们有了一个初步的想法,能否有一种工具可以快速求出一个多项式的点值表示,然后做$O(n)$的点值乘法,快速利用点值表示求出系数表示。
当然有这样的工具,在本篇文章中使用复数转化系数表示与点值表示。

复数

设$a,b$为实数,$i^2=-1$,形如$a+bi$的数称作复数,其中$i$称作虚数单位,$a$是该复数的实部,$b$是该复数的虚部。
复数域是已知的最大的域。

复平面

在复平面中,$x$轴代表实数、$y$轴(除原点外的所有点)代表虚数。每一个复数$a+bi$对应复平面上一个从$(0,0)$指向$(a,b)$的向量。

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

复数运算遵循向量运算的规则:

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

欧拉公式

根据欧拉公式有:

该公式曾经被评选为世界上最优美的公式

单位根

单位根的代数定义:$n$次单位根指$\omega^n=1$的所有复数解,其中幅角为正且最小的复数解。记作$\omega_n$。
单位根的几何定义:在复平面上,以原点为圆心,$1$为半径作圆,所得的圆叫做单位圆。以原点为起点,单位圆的$n$等分点为终点,作$n$个向量。设所得的幅角为正且最小的向量对应的复数为$\omega_n$,称为$n$次单位根。

根据欧拉公式,我们可以将单位根用具体复数表示出来:

单位根的幅角为周角的$1\over n$。

我们可以以此定义其他的复数解(向量)$\omega_n^2,\omega_n^3,\ldots,\omega_n^n$:
(代数定义)$\omega_n^k=(e^\frac{2\pi i}n)^k$。
(几何定义)单位圆上其他的向量。(可以根据代数定义与复数运算法则得到)

不难证明:$\omega_n^0=\omega_n^n=1$。

单位根的性质

消去引理:对于任意的正整数$n,d$和非负整数$k$,满足$\omega_{dn}^{dk}=\omega_n^k$

代数证明:

几何证明:
从几何意义上看,在复平面上,二者表示的向量终点相同。
更具体的,有:

折半引理:若$n$是偶数,有$\omega_n^{k+\frac n2}=-\omega_n^k$

代数证明:

几何证明:

求和引理:若对于正整数$n,k$,$n\nmid k$,有$\sum\limits_{j=0}^{n-1}(\omega_n^k)^j=0$

证明:
令$S(\omega_n^k)=\sum\limits_{j=0}^{n-1}(\omega_n^k)^j$。
因为$k\neq0$,公差不为$1$,根据等比数列求和公式可得:

因为$n\nmid k$,因此分母不为零。
当$k=0$时,不难发现$S(\omega_n^k)=n$。

快速傅里叶变换

考虑多项式$A(x)$的表示。将$n$次单位根的$0$到$n-1$次幂带入多项式的系数表示,所得点值向量$(A(\omega_n^0),A(\omega_n^1),\ldots,A(\omega_n^{n-1}))$称为其系数向量$(a_0,a_1,\ldots,a_{n-1})$的离散傅里叶变换

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

在下文中我们将多项式项数$n$视作$2$的整数幂,若不足在高位补零。

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

则有

  • 假设$k\lt \frac n2$,考虑$A(\omega_n^k)$:

    本处利用了单位根的消去引理。

  • 考虑$A(\omega_n^{k+\frac n2})$:

    本处利用了单位根的消去引理与折半引理。

这样分类我们可以遍历所有的$k\in[0,n-1]$。

也就是说,如果已知$A_1(x)$和$A_2(x)$在$\omega_\frac n2^0,\omega_\frac n2^1,\ldots,\omega_\frac n2^{\frac n2-1}$处的点值,就可以在$O(n)$的时间内求得$A(x)$在$\omega_n^0,\omega_n^1,\ldots,\omega_n^{n-1}$处的取值。而关于$A_1(x)$和$A_2(x)$的问题都是相对于原问题规模缩小了一半的子问题,分治的边界为一个常数项$a_0$。

根据主定理,该分治算法的时间复杂度为:

这就是最常用的FFT算法 —— Cooley-Tukey 算法。

傅里叶逆变换

将点值表示的多项式转化为系数表示,同样可以使用快速傅里叶变换,这个过程称为傅里叶逆变换
逆变换的推导是构造+化简的过程。

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

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

将上式展开,得:

根据单位根的求和引理,我们可知$\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i$在$j\neq k$时为$0$,但总有一次且仅有一次会出现$j=k$,当$j=k$时$\sum\limits_{i=0}^{n-1}(\omega_n^0)^i=n$。

故$c_k=na_k$。
因此$a_i=\frac1n c_i$。

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

实现

C++ 的 STL 在头文件complex中提供一个复数的模板实现complex<T>,其中T为实数类型,一般取double,在对精度要求较高的时候可以使用long double__float128(不常用)。

不难证明单位根的倒数等于其共轭复数。因此在程序实现中,为了减小误差,通常使用conj()取得idft所需的「单位根的倒数」。

递归实现

直接按照上面得到的结论来实现即可,比较直观。
Menci的代码:

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
const double PI = acos(-1);

bool inversed = false;

inline std::complex<double> omega(const int n, const int k) {
if (!inversed) return std::complex<double>(cos(2 * PI / n * k), sin(2 * PI / n * k));
else return std::conj(std::complex<double>(cos(2 * PI / n * k), sin(2 * PI / n * k)));
}

inline void transform(std::complex<double> *a, const int n) {
if (n == 1) return;

static std::complex<double> buf[MAXN];
const int m = n / 2;
// 按照系数奇偶划分为两半
for (int i = 0; i < m; i++) {
buf[i] = a[i * 2];
buf[i + m] = a[i * 2 + 1];
}
std::copy(buf, buf + n, a);

// 分治
std::complex<double> *a1 = a, *a2 = a + m;
fft(a1, m);
fft(a2, m);

// 合并两个子问题
for (int i = 0; i < m; i++) {
std::complex<double> x = omega(n, i);
buf[i] = a1[i] + x * a2[i];
buf[i + m] = a1[i] - x * a2[i];
}

std::copy(buf, buf + n, a);
}

递归实现

递归实现的FFT效率不高,实际中一般采用迭代实现。

二进制位翻转

考虑FFT到达分治递归边界时,每个数所在的位置及其二进制位:

发现规律,分治到边界后的每个数下标的等于原下标的二进制位翻转
有了这个规律,在进行二进制位反转后,就可以采用迭代的方法替代分治了(后面将讲的蝴蝶操作)。

1
2
3
4
5
6
int k=log2(n);
for(int i=0; i<n; i++) { //二进制位翻转
int t=0;
for(int j=0; j<k; j++)if(i&(1<<j))t|=(1<<(k-j-1));
if(i<t)swap(a[i],a[t]); //避免多次翻转
}

蝴蝶操作

考虑合并两个子问题的过程,假设$A_1(\omega_\frac n2^k)$和$A_2(\omega_\frac n2^k)$分别储存在$a(k)$和$a(\frac n2+k)$中,$A(\omega_n^k)$和$A(\omega_n^{k+\frac n2})$将要被存放在$b(k)$和$b(\frac n2+k)$中,合并的单位操作可表示为:

考虑加入一个临时变量$t$,使得这个过程可以在原地完成,不需要另一个数组$b$,也就是说,将$A(\omega_n^k)$和$A(\omega_n^{k+\frac n2})$存放在$a(k)$和$a(\frac n2+k)$中,覆盖原来的值。

因为每个位置只被调用一次,且在被覆盖的值刚好被调用,因此不会出错。

这一过程被称为蝴蝶操作
为什么叫蝴蝶操作呢?
因为这张图看起来像蝴蝶。。。(我觉得不像啊)

有了二进制位翻转与蝴蝶操作,我们就可以枚举区间长度$l$两两合并即可。

单位根

我们可以事先调用实部与虚部预处理单位根,在蝴蝶操作需要用到的$\omega_l^k$可以通过等式$\omega_l^k=\omega_n^{\frac nlk}$转化为预处理过的单位根(消去引理)。

二进制位翻转+蝴蝶操作 代码

其中cpcomplex<double>omega[i]是预处理过的单位根$\omega_n^i$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
void transform(cp* a,cp* omega) {
int k=log2(n);
for(int i=0; i<n; i++) { //二进制位翻转
int t=0;
for(int j=0; j<k; j++)if(i&(1<<j))t|=(1<<(k-j-1));
if(i<t)swap(a[i],a[t]); //避免翻转两次
}
for(int len=2; len<=n; len*=2) {
int mid=len>>1;
for(cp* p=a; p!=a+n; p+=len)
for(int i=0; i<mid; i++) {
cp t=omega[n/len*i]*p[mid+i];
p[mid+i]=p[i]-t;
p[i]+=t;
}
}
}

模板

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
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
#include<algorithm>
#include<iostream>
#include<iomanip>
#include<cstring>
#include<complex>
#include<cstdlib>
#include<climits>
#include<complex>
#include<vector>
#include<cstdio>
#include<cmath>
#include<queue>
using namespace std;

#define cp complex<double>

inline const int Get_Int() {
int num=0,bj=1;
char x=getchar();
while(x<'0'||x>'9') {
if(x=='-')bj=-1;
x=getchar();
}
while(x>='0'&&x<='9') {
num=num*10+x-'0';
x=getchar();
}
return num*bj;
}

const int maxn=131072+5;
const double pi=acos(-1);

struct FastFourierTransform {
int n,rev[maxn];
cp omega[maxn],iomega[maxn];

void init(int n) {
this->n=n;
for(int i=0; i<n; i++) {
omega[i]=cp(cos(2*pi/n*i),sin(2*pi/n*i));
iomega[i]=conj(omega[i]); //共轭复数
}
int k=log2(n);
for(int i=0; i<n; i++) { //二进制位翻转
int t=0;
for(int j=0; j<k; j++)if(i&(1<<j))t|=(1<<(k-j-1));
rev[i]=t;
}
}

void transform(cp* a,cp* omega) {
for(int i=0; i<n; i++)if(i<rev[i])swap(a[i],a[rev[i]]); //避免翻转两次
for(int len=2; len<=n; len*=2) {
int mid=len>>1;
for(cp* p=a; p!=a+n; p+=len)
for(int i=0; i<mid; i++) {
cp t=omega[n/len*i]*p[mid+i];
p[mid+i]=p[i]-t;
p[i]+=t;
}
}
}

void dft(cp* a) {
transform(a,omega);
}

void idft(cp* a) {
transform(a,iomega);
for(int i=0; i<n; i++)a[i]/=n;
}
} fft;

void Multiply(const int* a1,const int n1,const int* a2,const int n2,int* ans) {
int n=1;
while(n<n1+n2)n<<=1; //化成整数
static cp c1[maxn],c2[maxn];
for(int i=0; i<n1; i++)c1[i].real(a1[i]);
for(int i=0; i<n2; i++)c2[i].real(a2[i]);
fft.init(n);
fft.dft(c1);
fft.dft(c2);
for(int i=0; i<n; i++)c1[i]*=c2[i];
fft.idft(c1);
for(int i=0; i<n1+n2-1; i++)ans[i]=round(c1[i].real());
}

一些卡常技巧

2018/2/26 根据经验update

  1. 少使用Multiply()封装,根据情况FFT。
  2. 当每次init后FFT次数很少,可以不预处理$\omega$,改为每次重新计算。(时间效率提高$10\%\sim 20\%$)
  3. 手写complex<double>类。(时间效率提高约$20\%$)
  4. 可以构造共轭复数使得卷积时只需要一次dft,但相对的精度误差会变大。(时间效率提供约$20\%\sim 30\%$)
    修改上述模板代码:
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    void idft(cp* a) {
    transform(a,iomega);
    for(int i=0; i<n; i++)a[i]=cp(a[i].real()/(n<<2),a[i].imag()/n);
    }

    void Multiply(const int* a1,const int n1,const int* a2,const int n2,int* ans) {int n=1;
    while(n<n1+n2)n<<=1; //化成整数
    static cp c[maxn];
    for(int i=0; i<max(n1,n2); i++)c[i].real(a1[i]+a2[i],a1[i]-a2[i]);
    fft.init(n);
    fft.dft(c);
    for(int i=0; i<n; i++)c[i]*=c[i];
    fft.idft(c);
    for(int i=0; i<n1+n2-1; i++)ans[i]=round(c[i].real());
    }

注意事项

  • 因为FFT使用到了多次单位根折半引理,因此我们需要在高位补上足够的$0$使得项数成为$2$的整数幂。
  • 使用FFT时可能面临精度问题,根据题目模数的分解选择使用FFT还是NTT
  • 若使用FFT/NTT计算高精度乘法,注意去掉高位的$0$。

参考资料

姥爷们赏瓶冰阔落吧~