快速傅里叶变换

简介

快速傅里叶变换(FFT)

什么是FFT?FFT(fast Fourier transform)是
用来计算离散傅里叶变换(discrete Fourier transform,FFT)及其逆变换(IDFT)的快速算法。如果没有数学分析基础,这里很难用严密的语言把它解释清楚;但如果你只需要一个直观感受的话,不妨记住这句话:DFT把时域信号变换为频域信号。

注意,时域和频域只是两种信号分析方法,并不是两种不同类别的信号。在谈论语言信号的时候,“音量小”指的是时域,而“声音尖”指的是频域。在录音的时候,音频通常按照时域形式保存,即各个时刻采样到的振幅;如果需要各个频率的数据,则需要DFT。

DFT有一个很有意思的性质:时域卷积,频域卷积;频域卷积,时域卷积。如果你不知道什么是卷积(convolution),也没关系,请再记住一句话:多项式乘法实际上是多项式系数向量的卷积。

多项式乘法

给定两个单变量多项式$A(x)$和$B(x)$,次数均不超过$n$,如果快速计算两者的乘积呢?最简单的方法就是把系数两两相乘,再相加。这样计算的时间复杂度为$O(n^2)$,当$n$很大时速度将会很慢。注意,高精度整数的乘法是多项式乘法的特殊情况(想一想,为什么),所以多项式乘法的快速乘法也可以用来做高精度乘法。

上述算法之所以慢,是因为表示方法不好。虽然“系数序列”是最自然的表示方法,但却不适合用来计算乘法。在多项式快速乘法中,需要用点值法来表示一个多项式。点值表示是一个“点-值”对的序列$\lbrace(x_0,y_0),(x_1,y_1)…(x_{n-1},y_{n-1})\rbrace$,它代表一个满足$y_k=A(x_k)$的多项式$A$。可以证明:恰好有一个次数小于$n$的多项式满足这个条件。

点值表示法非常适合做乘法:只要两个多项式的点集$\lbrace x_i \rbrace$相同,则只需要把对应的值乘起来就可以了,只需$O(n)$时间。但问题在于输入输出的时候仍需要采用传统的系数表示法,因此需要快速地在系数表示和点值表示之间转换。还记得刚才那句话吗?时域卷积,频域卷积。也就是说,多项式的系数表示法对应于时域,而点值法对应于频域,因此只需要用FFT计算出一个DFT就可以完成转换。

单位根

这样算出来的点值表示法,对应的求值点是哪些呢?答案是$2n$次单位根。所谓“$n$次单位根”是满足$x^{n}=1$ 的复数。$x^{n}=1$ 的话,$x$岂不是就等于1?很遗憾,那只是实数域中的情况。在复数域中,1恰好有$n$个单位根,$e^{2k\pi i/n}$,其中$i$是虚数单位$(i^2=-1)$, 而$e^{iu}=cos\ u+i\ sin\ u$。单位根非常特殊,因此FFT才有办法在更短的时间内求出多项式在这些点的取值。

利用FFT进行快速多项式乘法的具体步骤

  1. (补0)在两个多项式的最前面补0,得到两个$2n$次多项式,设系数向量分别为$v_1$和$v_2$。
  2. (求值)用FFT计算$f_1=DFT(v_1)$和$f_2=DFT(v_2)$。这里得到的$f_1$和$f_2$分别是两个输入多项式在$2n$次单位根出的各个取值(即点值表示)。
  3. (乘法)把两个向量$f_1$和$f_2$对应的每一维对应相乘,得带向量$f$。它对应输入多项式乘积的点值表示。
  4. (插值)用FFT计算$v=IDFT(f)$,其中$v$就是乘积的系数向量。

例题

题目来源 SHUOJ 1966

Description

给出$A_0,A_1,…A_{n-1}$和$B_0,B_1,…B_{n-1}$,求$C_0,C_1,…C_{n-1}$
其中$$C_{i+j}=\sum A_i B_j$$

Input

多组数据,第一行为一个数字$n$。
之后有两行数列,分别表示$A_0,A_1,…A_{n-1}$和$B_0,B_1,…B_{n-1}$。
$1≤ n ≤ 10^5$
$1≤ A_i,B_i ≤ 100$

Output

对于每组数据,输出一行共$n$个数字。

Sample Input

1
2
3
3
1 2 3
1 2 3

Sample Output

1
1 4 10

代码

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;

const double PI = acos(-1.0);
//复数结构体
struct Complex
{
double x, y; //实部和虚部 x+yi
Complex(double tx = 0.0, double ty = 0.0)
{
x = tx;
y = ty;
}
Complex operator-(const Complex &b) const
{
return Complex(x - b.x, y - b.y);
}
Complex operator+(const Complex &b) const
{
return Complex(x + b.x, y + b.y);
}
Complex operator*(const Complex &b) const
{
return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
}
};
/*
* 进行FFT和IFFT前的反转变换。
* 位置i和 (i二进制反转后位置)互换
* len必须取2的幂
*/
void change(Complex y[], int len)
{
int i, j, k;
for (i = 1, j = len / 2; i < len - 1; i++)
{
if (i < j)
swap(y[i], y[j]);
//交换互为小标反转的元素,i<j保证交换一次
//i做正常的+1,j左反转类型的+1,始终保持i和j是反转的
k = len / 2;
while (j >= k)
{
j -= k;
k /= 2;
}
if (j < k)
j += k;
}
}
/*
* 做FFT
* len必须为2^k形式,
* on==1时是DFT,on==-1时是IDFT
*/
void fft(Complex y[], int len, int on)
{
change(y, len);
for (int h = 2; h <= len; h <<= 1)
{
Complex wn(cos(-on * 2 * PI / h), sin(-on * 2 * PI / h));
for (int j = 0; j < len; j += h)
{
Complex w(1, 0);
for (int k = j; k < j + h / 2; k++)
{
Complex u = y[k];
Complex t = w * y[k + h / 2];
y[k] = u + t;
y[k + h / 2] = u - t;
w = w * wn;
}
}
}
if (on == -1)
for (int i = 0; i < len; i++)
y[i].x /= len;
}

const int MAXN = 400010;
Complex a[MAXN], b[MAXN], c[MAXN];
int x[MAXN], y[MAXN];
int sum[MAXN];

int main()
{
int n;
while (cin >> n)
{
int t = 1;
for (int i = 0; i < n; i++)
cin >> x[i];
for (int i = 0; i < n; i++)
cin >> y[i];
for (int i = 0; i < n; i++)
{
a[i] = Complex(x[i], 0);
b[i] = Complex(y[i], 0);
}
while (t < n * 2)
t <<= 1;
for (int i = n; i < t; i++)
{
a[i] = Complex(0, 0);
b[i] = Complex(0, 0);
}
fft(a, t, 1);
fft(b, t, 1);
for (int i = 0; i < t; i++)
c[i] = a[i] * b[i];
fft(c, t, -1);
for (int i = 0; i < t; i++)
sum[i] = (int)(c[i].x + 0.5);
for (int i = 0; i < n; i++)
{
if (i)
cout << " " << sum[i];
else
cout << sum[i];
}
cout << endl;
}
return 0;
}
捐助作者