快速傅里叶变换

【注】参考自算法导论。

1. 简介

快速傅里叶变换(FFT)是实现离散傅里叶变换(DFT)和离散傅里叶逆变换(IDFT)的快速算法,其时间复杂度为 O(nlogn)O(n \log n)。DFT 在实际生活中有很多应用,比如通过离散傅里叶变换,可以将系数表示的多项式转为点值表示的多项式,从而使得多项式的乘法的复杂度由 O(n2)O(n^2) 降为 O(n)O(n)

2. 实现

快速傅里叶变换涉及到离散傅里叶变换的知识,FFT 能够实现的基础在于折半定理,折半定理提供了求大整数的 DFT 的一种思路:即通过将大整数不断划分为两部分,然后分别求出两部分的 DFT,最后在合并,本质是一种分治的思想。

3. 代码

3.1 快速傅里叶变换模板

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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
#include <bits/stdc++.h>
using namespace std;

// 复数
struct Complex{
double re, im;
Complex(): re(0), im(0) {}
Complex(double _re, double _im): re(_re), im(_im) {}
Complex operator + (const Complex cp) const {
return Complex(re+cp.re, im+cp.im);
}
Complex operator - (const Complex cp) const {
return Complex(re-cp.re, im-cp.im);
}
Complex operator * (const Complex cp) const {
return Complex(re*cp.re - im*cp.im, re*cp.im + im*cp.re);
}
Complex operator / (const double div) const {
return Complex(re/div, im/div);
}
void operator = (const Complex cp) {
re = cp.re;
im = cp.im;
}
};
// 快速傅里叶变换
struct FFT{
#ifndef _FFT_
#define _FFT_
#define ll int
#define MAXN 3000005
#endif
// 待计算
ll cnt; // 位数
ll len; // 数组大小(len = 2^cnt)
ll r[MAXN]; // 位逆序置换数组
const double Pi = acos(-1.0);
FFT(): cnt(0), len(0) {}
// 构建 2 的幂次的系数数组
void build(Complex a[], ll curlen) {
// 计算 len 所占位数
cnt = 1;
len = curlen;
while((len >>= 1)) ++cnt;
// 找到最小的大于等于 curlen 的 2 的幂值
len = (ll)1 << cnt;
if(len-curlen == curlen) {
--cnt;
len = curlen;
}
// 超出 curlen 部分的系数补 0
for(ll i = curlen; i < len; ++i) a[i] = Complex(0,0);
}
// 计算位逆序置换
void bitReverse() {
for(ll i = 0; i < len; ++i) r[i] = i;
// 本质是动态规划
// 利用 i>>1 的位逆序推导出 i 的位逆序
for(ll i = 0; i < len; ++i) r[i] = (r[i>>1]>>1) | ((i&1)<<(cnt-1));
}
// 快速傅里叶变换
void DFT(Complex y[]) {
// 位逆序置换
for(ll i = 0; i < len; ++i) {
if(i < r[i]) swap(y[i], y[r[i]]);
}
// 合并 DFT
// 枚举合并长度
for(ll i = 1; i < len; i<<=1) {
// 主单位复根
Complex wn(cos(Pi/i), sin(Pi/i));
// 枚举合并组首下标
for(ll j = 0; j < len; j += (i<<1)) {
Complex w(1,0);
// 合并组 DFT
for(ll k = j; k < j+i; ++k, w=w*wn) {
Complex u = y[k];
Complex v = w*y[k+i];
y[k] = u + v;
y[k+i] = u - v;
}
}
}
}
// 快速傅里叶逆变换
void IDFT(Complex a[]) {
// 位逆序置换
for(ll i = 0; i < len; ++i) {
if(i < r[i]) swap(a[i], a[r[i]]);
}
// 合并 DFT
// 枚举合并长度
for(ll i = 1; i < len; i<<=1) {
// 主单位复根
Complex wn(cos(Pi/i), -sin(Pi/i));
// 枚举合并组首下标
for(ll j = 0; j < len; j += (i<<1)) {
Complex w(1,0);
// 合并组 DFT
for(ll k = j; k < j+i; ++k, w=w*wn) {
Complex u = a[k];
Complex v = w*a[k+i];
a[k] = u + v;
a[k+i] = u - v;
}
}
}
for(ll i = 0; i < len; ++i) {
a[i] = a[i]/len;
}
}
};

3.2 大整数乘法模板

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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
#include <bits/stdc++.h>
using namespace std;

// 复数
struct Complex{
double re, im;
Complex(): re(0), im(0) {}
Complex(double _re, double _im): re(_re), im(_im) {}
Complex operator + (const Complex cp) const {
return Complex(re+cp.re, im+cp.im);
}
Complex operator - (const Complex cp) const {
return Complex(re-cp.re, im-cp.im);
}
Complex operator * (const Complex cp) const {
return Complex(re*cp.re - im*cp.im, re*cp.im + im*cp.re);
}
Complex operator / (const double div) const {
return Complex(re/div, im/div);
}
void operator = (const Complex cp) {
re = cp.re;
im = cp.im;
}
};
// 大整数乘法(快速傅里叶变换实现)
struct BigIntegersMultiply{
#ifndef _BIGINTEGERSMULTIPLY_
#define _BIGINTEGERSMULTIPLY_
#define ll int
#define MAXN 3000005
#endif
const double Pi = acos(-1.0);
// 待计算
ll cnt; // 位数
ll len; // 最终乘积数组大小(len = 2^cnt)
ll r[MAXN]; // 位逆序置换数组
ll ans[MAXN]; // 最终答案(逆序存储)
// 待输入
ll len1, len2; // 初始系数数组大小
Complex x1[MAXN], x2[MAXN]; // 两个大整数乘数(逆序储存)

BigIntegersMultiply(): cnt(0), len(0) {}
// 构建 2 的幂次的系数数组
void build() {
// 计算两个乘数最小的公共 2 的幂次 len 及所占位数 cnt
cnt = 0, len = 1;
while(len < (len1<<1) || len < (len2<<1)) ++cnt, len <<= 1;
// 超出 curlen 部分的系数补 0
for(ll i = len1; i < len; ++i) x1[i] = Complex(0,0);
for(ll i = len2; i < len; ++i) x2[i] = Complex(0,0);
}
// 计算位逆序置换
void bitReverse() {
for(ll i = 0; i < len; ++i) r[i] = i;
// 本质是动态规划
// 利用 i>>1 的位逆序推导出 i 的位逆序
for(ll i = 0; i < len; ++i) r[i] = (r[i>>1]>>1) | ((i&1)<<(cnt-1));
}
// 快速傅里叶变换
void DFT(Complex x[]) {
// 位逆序置换
for(ll i = 0; i < len; ++i) {
if(i < r[i]) swap(x[i], x[r[i]]);
}
// 合并 DFT
// 枚举合并长度
for(ll i = 1; i < len; i<<=1) {
// 主单位复根
Complex wn(cos(Pi/i), sin(Pi/i));
// 枚举合并组首下标
for(ll j = 0; j < len; j += (i<<1)) {
Complex w(1,0);
// 合并组 DFT
for(ll k = j; k < j+i; ++k, w=w*wn) {
Complex u = x[k];
Complex v = w*x[k+i];
x[k] = u + v;
x[k+i] = u - v;
}
}
}
}
// 快速傅里叶逆变换
void IDFT(Complex x[]) {
// 位逆序置换
for(ll i = 0; i < len; ++i) {
if(i < r[i]) swap(x[i], x[r[i]]);
}
// 合并 DFT
// 枚举合并长度
for(ll i = 1; i < len; i<<=1) {
// 主单位复根
Complex wn(cos(Pi/i), -sin(Pi/i));
// 枚举合并组首下标
for(ll j = 0; j < len; j += (i<<1)) {
Complex w(1,0);
// 合并组 DFT
for(ll k = j; k < j+i; ++k, w=w*wn) {
Complex u = x[k];
Complex v = w*x[k+i];
x[k] = u + v;
x[k+i] = u - v;
}
}
}
for(ll i = 0; i < len; ++i) {
x[i] = x[i]/len;
}
}
// 大整数乘法
void bigMultiply() {
build();
bitReverse();
DFT(x1), DFT(x2);
for(ll i = 0; i < len; ++i) x1[i] = x1[i]*x2[i];
IDFT(x1);
for(ll i = 0; i < len; ++i) ans[i] = ll(x1[i].re + 0.5);
for(ll i = 0; i < len; ++i) ans[i+1] += ans[i]/10, ans[i] %= 10;
len = len1 + len2 - 1;
while(len && !ans[len]) --len;
++len;
}
};