复变、控制论、ACM,学了三次忘了三次……

在这里只是整理一下提纲,详细推导与代码推荐参考该知乎文章

傅里叶变换

傅里叶级数的三角形式

\(f(t)\)是以\(T\)为周期的实值函数,且在\([-\frac{T}{2},\frac{T}{2}]\)上满足迪利克雷条件,

则在\(f(t)\)的连续点处,有

\[ f(t)=\frac{a_0}{2}+\sum_{n=1}^{+\infty}{(a_n\cos{n\omega_0t}+b_n\sin{n\omega_0t})} \]

其中

\[ \begin{aligned} \omega_0 &=\frac{2\pi}{T} \\ a_n &=\frac{2}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}{f(t)\cos{n\omega_0t\ dt}} \\ b_n &=\frac{2}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}{f(t)\sin{n\omega_0t\ dt}} \end{aligned} \]

在间断点处,有\(f(t)=\frac{f(t+0)+f(t-0)}{2}\)

傅里叶级数的物理含义

\(f(t)\)代表信号,

则由上可知,一个周期为\(T\)的信号可以分解为简谐波之和。

(非周期信号可以看作周期无穷大)

傅里叶变换

由于三角函数可以用自然底数的虚数幂表示,即\(e^{j\theta}=\cos{\theta}+j\sin{\theta}\)

所以可以将上面的傅里叶级数换一下表示方法,即有傅里叶变换与逆变换

\[ \begin{aligned} F(\omega) &=\int_{-\infty}^{+\infty}{f(t)e^{-j\omega t}\ dt} \\ F(\omega) &=\mathcal{F}[f(t)] \\ f(t) &=\mathcal{F}^{-1}[F(\omega)] \\ f(t) &=\frac{1}{2\pi}\int_{-\infty}^{+\infty}{F(\omega)e^{j\omega t}\ d\omega} \end{aligned} \]

\(F(\omega)\)称为像函数,\(f(t)\)称为原函数。

离散傅里叶变换 DFT

对于\(f(t)\)进行采样,得到离散信号\(x(n)\)。与上面连续函数的傅里叶变换类似,将积分改成求和,将周期\(T\)改为采样次数\(N\),将连续时间\(t\)变为采样的次数\(n\),得离散傅里叶变换和逆变换

\[ \begin{aligned} X[k] &=\sum_{n=0}^{N-1}{x[n]e^{-j\frac{2\pi nk}{N}}} \\ x[n] &=\frac{1}{N}\sum_{k=0}^{N-1}{X[k]e^{j\frac{2\pi nk}{N}}} \end{aligned} \]

其中\(k\in[0,N-1]\)

在ACM中常用的是对于多项式的离散傅里叶变换,但在讨论它之前,先讲一些补充知识。之后在下一节对于快速傅里叶变换的讨论主要就是着眼于对多项式的处理。

补充1 多项式的表示法

多项式有两种表示方法:系数表示法和点值表示法。

(下面讨论\(n-1\)次多项式)

系数表示法

\[ y=A(x)=\sum_{i=0}^{n-1}{a_i\cdot x^i} \]

点值表示法

将互不相同的\((x_0,x_1,\cdots,x_n)\)代入\(A(x)\)(插值),得\((y_0,y_1,\cdots,y_n)\).

\(n\)个不同点进行插值得到\(n\)组点值,即可唯一确定这个\(n-1\)次多项式。

多项式乘法

对于\(A(x)\)\(B(x)\)作节点\((x_0,x_1,\cdots,x_n)\)的插值,

分别得点值向量\((y_a0,y_a1,\cdots,y_an)\)\((y_b0,y_b1,\cdots,y_bn)\)

那么对于多项式\(C(x)=A(x)\cdot B(x)\),其点值也可由相乘得出,即\(y_{ci}=y_{ai}\cdot y_{bi}\)

容易想到,用点值形式进行多项式运算更快。

补充2 复数单位根的定义与性质

单位根

终点在复平面的单位圆上的复数向量。

即上面提到过的虚数幂表示的三角函数\(e^{i\theta}\)

n次单位向量

将单位圆均分成\(n\)份,原点为起点。

幅角为正且最小的向量为\(n\)次单位向量,记为\(\omega_n^1\)

\(n\)次单位向量的\(k\)次幂即\(\omega_n^k\)

性质一(折半引理)

\[ \omega_{2n}^{2k}=\omega_n^k \]

性质二(消去引理)

\[ \omega_n^{k+\frac{n}{2}}=-\omega_n^k \]

快速傅里叶变换 FFT

对于多项式\(A(x)\),其系数向量为\((a_0,a_1,\cdots,a_n)\)

\(n\)次单位向量的\(0\sim n-1\)次幂作为节点代入\(A(x)\)作插值得点值向量\((A(\omega_n^0),A(\omega_n^1),\cdots,A(\omega_n^n))\)

这个过程即多项式的离散傅里叶变换。

然而,这个过程是\(O(n^2)\)的,太慢了。快速傅里叶变换就是利用上面补充的性质,加快运算的方法。

快速傅里叶变换

对于

\[ A(x)=a_0+a_1x^1+a_2x^2+\cdots+a_{n-1}x^{n-1} \]

按奇偶性分组(假设\(n\)为偶数)

\[ \begin{aligned} A(x) &=(a_0+a_2x^2+a_4x^4+\cdots+a_{n-2}x^{n-2})+x\cdot(a_1+a_3x^2+a_5x^4\cdots+a_{n-1}x^{n-2}) \\ A(x) &=A1(x^2)+x\cdot A2(x^2) \end{aligned} \]

作插值\(x=\omega_n^k\),对上标分段讨论。

\(k\in[0,\frac{n}{2}-1]\)

\[ A(\omega_n^k)=A1(\omega_n^{2k})+\omega_n^k\cdot A2(\omega_n^{2k}) \]

根据折半引理,有

\[ A(\omega_n^k)=A1(\omega_{\frac{n}{2}}^k)+\omega_n^k\cdot A2(\omega_{\frac{n}{2}}^k) \]

对于另一半,即\((k+\frac{n}{2})\in[\frac{n}{2},n-1]\)

用消去原理,有

\[ A(\omega_n^{k+\frac{n}{2}})=A1(\omega_{\frac{n}{2}}^k)-\omega_n^k\cdot A2(\omega_{\frac{n}{2}}^k) \]

可见,只要知道\(A1(x),A2(x)\)\(\omega_{\frac{n}{2}}^0,\omega_\frac{n}{2}^1,\cdots,\omega_{\frac{n}{2}}^{\frac{n}{2}-1}\)的值,就可以\(O(n)\)时间求\(A(x)\)

每次都是一半的规模,整体的时间复杂度为\(O(n\log{n})\)

快速傅里叶逆变换

证明过程略。

\(A(x)\)由插值节点\((\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1})\)作离散傅里叶变换,得点值向量\((d_0,d_1,\cdots,d_{n-1})\)

将这个点值向量作为系数向量,由插值节点\((\omega_n^0,\omega_n^{-1},\cdots,\omega_n^{-(n-1)})\)再作一次变换,之后除以\(n\)

得到的\((\frac{c_0}{n},\frac{c_1}{n},\cdots,\frac{c_{n-1}}{n})\)就是多项式的系数向量\((a_0,a_1,\cdots,a_{n-1})\)

快速傅里叶变换的实现方式

一般用Cooley-Tukey算法,略,详见前言的参考链接。

做题用的模板如下

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
//FFT模板
//最高次为n、m次的多项式相乘
//系数为正整数,若有负数注意最后结果的取整方式
#include<bits/stdc++.h>
using namespace std;
inline int read(){
int f=1,x=0;char ch;
do{ch=getchar();if(ch=='-')f=-1;}while(ch<'0'||ch>'9');
do{x=x*10+ch-'0';ch=getchar();}while(ch>='0'&&ch<='9');
return f*x;
}
//_head
#define pi acos(-1)
#define N 3000005
typedef complex<double> cplx;
int n, m, l, r[N];
cplx a[N], b[N];
//_basic

void init(){
//注意按题目设置m
m += n;
for(n=1 ; n<=m ; n<<=1)
l++;
for(int i=0 ; i<n ; i++)
r[i]=(r[i>>1]>>1) | ((i&1)<<(l-1));
}

void fft(cplx *a, int flag){
for(int i=0 ; i<n ; i++)
if(i < r[i])
swap(a[i], a[r[i]]);
for(int i=1 ; i<n ; i<<=1){
cplx wn(cos(pi/i), flag*sin(pi/i));
for(int p=i<<1, j=0 ; j<n ; j+=p){
cplx w(1, 0);
for(int k=0 ; k<i ; k++, w*=wn){
cplx x=a[j+k], y=w*a[j+k+i];
a[j+k] = x+y;
a[j+k+i] = x-y;
}
}
}
//逆变换自己除n,范围0~m
if(flag == -1)
for(int i=0 ; i<=m ; i++)
a[i] /= n;
}

int main(){
//多项式最高n、m次
n = read();
m = read();
for(int i=0 ; i<=n ; i++)
a[i] = read();
for(int i=0 ; i<=m ; i++)
b[i] = read();

//初始化r[],并将n、m变为点值表示法、系数表示法的最高次数
init();

//DFT
fft(a, 1);
fft(b, 1);
//点值相乘0~n
for(int i=0 ; i<=n ; i++)
a[i] = a[i]*b[i];
//IDFT
fft(a, -1);

for(int i=0 ; i<=m ; i++){
printf("%d",(int)(a[i].real()+0.5));//向下取整
if(i == m)
printf("\n");
else
printf(" ");
}

}
/*
example
input>
1 2
1 2
1 2 1
output>
1 4 5 2
*/