隐藏
快速傅里叶变换（FFT）与多项式乘法学习笔记 | Bill Yang's Blog

路终会有尽头，但视野总能看到更远的地方。

0%

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

发表于 更新于 分类于 Valine：

本文章基于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}$转化为预处理过的单位根（消去引理）。

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

其中cp为complex<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$。

参考资料

相关文章
姥爷们赏瓶冰阔落吧~