-1、为什么学 FFT
退役(很早)之前听说 FFT 很神(e)奇(xin),Po姐来讲的时候也是膜(sha)了(ye)一(bu)发(dong),于是就放那里了。退役之后有(xian)了(de)时(mei)间(shi),并且在篮球赛之前立了赢一场的 flag 否则学FFT结果又双叒叕全输了,于是下面就是成果了……网上 FFT 讲解多得是不看也罢……
0、什么是 FFT?为什么 FFT?
0.1、为什么 FFT?
从一个简单问题说起:大整数乘法。在做 VijosP2000 的时候,看数据范围 O ( n 2 ) O(n^2) O(n2) 有点卡常?在做 HDOJ1402 的时候 O ( n 2 ) O(n^2) O(n2) 接着卡常?大整数乘法的时间复杂度只能是 O ( n 2 ) O(n^2) O(n2) 了?不是的。从算法优化上可以实现为 O ( n log 2 3 ) O(n^{\log_2 3}) O(nlog23)的方法( log 2 3 ≈ 1.585 \log_2 3\approx1.585 log23≈1.585,此方法为 Karatsuba 乘法,详解看这里),还是不太友好。那么就会用到 O ( n log 2 n ) O(n\log_2 n) O(nlog2n) 的 FFT 了。
0.2、什么是 FFT?
FFT(Fast Fourier Transform),全名快速傅里叶变换,这与傅里叶变换只有半毛钱的关系,傅里叶变换是分析波的成分的方法,通过推广傅里叶变换得到了优化的快速傅里叶变换,为了展示区别和联系,下面是傅里叶变换和傅里叶逆变换的公式:
傅里叶变换:
F
(
ω
)
=
F
[
f
(
t
)
]
=
∫
−
∞
∞
f
(
t
)
e
−
i
ω
t
d
t
F(\omega)=\mathcal{F}[f(t)]=\int\limits_{-\infty}^\infty f(t)e^{-i\omega t} dt
F(ω)=F[f(t)]=−∞∫∞f(t)e−iωtdt
傅里叶逆变换:
f ( t ) = F − 1 [ F ( ω ) ] = 1 2 π ∫ − ∞ ∞ F ( ω ) e i ω t d ω f(t)=\mathcal{F}^{-1}[F(\omega)]=\frac{1}{2\pi} \int\limits_{-\infty}^\infty F(\omega)e^{i\omega t} d\omega f(t)=F−1[F(ω)]=2π1−∞∫∞F(ω)eiωtdω
其中
f
(
t
)
f(t)
f(t) 为以
t
t
t 为周期的周期函数,
t
t
t 满足狄利赫里条件,
F
(
ω
)
F(\omega)
F(ω) 称为
f
(
t
)
f(t)
f(t) 的像函数,
f
(
t
)
f(t)
f(t) 称为
F
(
ω
)
F(\omega)
F(ω) 的原像函数。具体问题请参考《高等数学》,博主还是高中生……(高中高等数学233333)
然后,我们终于可以开始讲 FFT 了…
1 、(初中知识)多项式
1.1、定义
我们把形如
a
0
+
a
1
x
+
a
2
x
2
+
⋯
+
a
n
−
1
x
n
−
1
a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}
a0+a1x+a2x2+⋯+an−1xn−1 的式子叫做多项式,其中
a
0
,
a
1
,
⋯
,
a
n
−
1
a_0,a_1,\cdots,a_{n-1}
a0,a1,⋯,an−1 称为多项式的系数,多项式中单项式个数为多项式的项数。因为多项式由单项式组成,规定最高次项的次数叫做多项式的次数,这里
A
(
x
)
A(x)
A(x) 的次数为
n
−
1
n-1
n−1。
x
x
x 为未知数。
以上内容均为初二(我记得是)的内容,不明白的自觉面壁思过(初中数学老师:这个定义抄
10
10
10 遍给我)……
下面就是没学过的了:
为方便表示,我们记多项式
a
0
+
a
1
x
+
a
2
x
2
+
⋯
+
a
n
−
1
x
n
−
1
a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}
a0+a1x+a2x2+⋯+an−1xn−1 为
A
(
x
)
A(x)
A(x),即:
A ( x ) = ∑ i = 0 n − 1 a i x i A(x)=\sum\limits_{i=0}^{n-1}a_ix^i A(x)=i=0∑n−1aixi
我们规定:严格大于
n
n
n 的任意整数为这个多项式的次数界,即次数的范围。多项式中任意单项式次数均严格小于次数界。
哦,行……
下面无特殊说明,多项式均用
A
(
x
)
A(x)
A(x) 等类似形式表示。
1.2、多项式的两种表示法
1.2.1、系数表示法
我们知道,向量是个好东西(蛤?),矩阵也是个好东西(蛤?),于是我们把三者结合起来看(蛤玩意?)。
我们将
n
−
1
n-1
n−1 次多项式的系数看成
n
n
n 维向量,即
a
⃗
=
(
a
0
,
a
1
,
⋯
,
a
n
−
1
)
\vec a=(a_0,a_1,\cdots,a_{n-1})
a=(a0,a1,⋯,an−1),则称
a
⃗
\vec a
a 为多项式
A
(
x
)
A(x)
A(x) 的系数表示。
在选修 4-2 中,我们知道向量与矩阵有着密切联系,这就是为什么我们把多项式,向量和矩阵联系起来。(然而我们并没选 4-2 蛤蛤蛤)
于是这个
n
n
n 维向量可以表示为一个
n
n
n 维列向量,为:
[
a
0
a
1
⋮
a
n
−
1
]
\begin{bmatrix} a_0\\ a_1\\ \vdots \\ a_{n-1} \end{bmatrix}
⎣⎢⎢⎢⎡a0a1⋮an−1⎦⎥⎥⎥⎤
1.2.2、点值表示法
我们把多项式
A
(
x
)
A(x)
A(x) 看做一个函数,选取不同的
n
n
n 个值
x
0
,
x
1
,
⋯
,
x
n
−
1
x_0,x_1,\cdots,x_{n-1}
x0,x1,⋯,xn−1 代入
A
(
x
)
A(x)
A(x) 中,可以得到
A
(
x
)
A(x)
A(x) 的图像上
n
n
n 个不同的点,那么称点集
α
=
{
(
x
i
,
A
(
x
i
)
)
∣
i
∈
[
0
,
n
−
1
]
,
i
∈
Z
}
\alpha=\{ (x_i,A(x_i))\mid i\in[0,n-1],i\in \mathbb{Z}\}
α={(xi,A(xi))∣i∈[0,n−1],i∈Z} 为多项式
A
(
x
)
A(x)
A(x) 的点值表示。
到这里,是不是差不多忘了要干什么了……
1.3、多项式的运算
1.3.1、多项式求值
这个比较简单,拿题说。这是我还没出也不想出的一道题,拿出来娱乐一下……
题目描述 Description
都知道 FZ 酱数学不好,数学老师很着急,于是让她求个多项式的值。但是数学老师一不小心数据就出大了,他却并没有注意到这一点。老师共出了 T T T 道题,FZ 酱拿到了这 T T T 道题就崩溃了,请帮她尽快完成作业。
输入 Input
第一行一个整数 T T T,表示数据组数。
对于每组数据:
第一行有一个数 n n n 表示这个多项式的次数界;
第二行有 n n n 个整数 a 0 , a 1 , ⋯ , a n − 1 a_0,a_1,\cdots\ ,a_{n-1} a0,a1,⋯ ,an−1,表示多项式 ∑ i = 0 n − 1 a i x i \sum\limits_{i=0}^{n-1}a_ix^i i=0∑n−1aixi;
第三行,为未知数的值 x x x。
输出 Output
对于每组数据,输出一行,为答案。
样例输入 Sample Input
2
2
3 2
100
3
1 2 3
-9
样例输出 Sample Output
203
226
样例解释 Explanation
对于第一个多项式为 A ( x ) = 2 x + 3 A(x)=2x+3 A(x)=2x+3,代入后得 A ( 100 ) = 2 × 100 + 3 = 203 A(100)=2\times 100+3=203 A(100)=2×100+3=203;
对于第二个多项式为 A ( x ) = 3 x 2 + 2 x + 1 A(x)=3x^2+2x+1 A(x)=3x2+2x+1,代入后得 A ( − 9 ) = 226 A(-9)=226 A(−9)=226。
限制 Limits
对于 50 % 50\% 50% 的数据,答案在 [ 0 , 2 128 − 1 ] [0,2^{128}-1] [0,2128−1] 范围内;
对于 100 % 100\% 100% 的数据, T ≤ 2 T\le 2 T≤2, n ≤ 1 0 5 n\le 10^5 n≤105,保证输入均为绝对值不大于 2 31 − 1 2^{31}-1 231−1 的整数并且答案小于 3 × 1 0 6 3\times 10^6 3×106 位。
那么很好我们必须写高精度了……
2
128
−
1
≈
3.4
×
1
0
38
2^{128}-1\approx 3.4\times 10^{38}
2128−1≈3.4×1038,这里的
50
%
50\%
50% 数据是给不想写高精度的童鞋写 __int128
准备的(__int128
的范围是
[
0
,
2
128
−
1
]
[0,2^{128}-1]
[0,2128−1])。于是我们就有这样三种的方法。
记数字长度为
p
p
p。
1.3.1.1、方法一: O ( p n 2 ) → 50 pts O(pn^2)\to 50\text{pts} O(pn2)→50pts
暴力不会吗?
核心代码如下:
1.3.1.2、方法二: O ( p 2 n log n ) → [ 50 , 80 ] pts O(p^2n\log n)\to [50,80]\text{pts} O(p2nlogn)→[50,80]pts
这里挨个乘太慢了,快速幂处理会快一些。
(其实不一定会快多少)
核心代码如下:
1.3.1.3、方法三: O ( p 2 n ) → 100 pts O(p^2n)\to 100\text{pts} O(p2n)→100pts
学过必修三吗?
学过必修三还不会?
拿出数学必修三翻到 37 页你看到了什么?
秦九韶算法可以大量减少乘法和加法次数,并且把运算量简化为
O
(
n
)
O(n)
O(n) 的。
核心代码如下:
当然,还可以处理前缀积,也是
O
(
n
)
O(n)
O(n) 的。
核心代码如下:
不过没上面的优雅不是吗……(呵呵)
回到正题上,我们要谈的其实是……
1.3.2、多项式的加法和乘法
1.3.2.1、多项式加法
没啥难度,直接加即可:
C
(
x
)
=
A
(
x
)
+
B
(
x
)
=
∑
i
=
0
n
−
1
a
i
x
i
+
∑
j
=
0
n
−
1
b
j
x
j
=
∑
i
=
0
n
−
1
(
a
i
+
b
i
)
x
i
=
∑
i
=
0
n
−
1
c
i
x
i
\begin{aligned} C(x)&=A(x)+B(x)\\ &=\sum\limits_{i=0}^{n-1}a_ix^i+\sum\limits_{j=0}^{n-1}b_jx^j\\ &=\sum\limits_{i=0}^{n-1}(a_i+b_i)x^i\\ &=\sum\limits_{i=0}^{n-1}c_ix^i \end{aligned}
C(x)=A(x)+B(x)=i=0∑n−1aixi+j=0∑n−1bjxj=i=0∑n−1(ai+bi)xi=i=0∑n−1cixi
这个操作可以
O
(
n
)
O(n)
O(n) 完成。
1.3.2.2、多项式乘法
这就是个开括号的问题……
两个次数界为
n
a
,
n
b
n_a,n_b
na,nb 的多项式
A
(
x
)
,
B
(
x
)
A(x),B(x)
A(x),B(x) 相乘,得到的多项式
C
(
x
)
C(x)
C(x) 的次数界为
n
a
+
n
b
−
1
n_a+n_b-1
na+nb−1。
我们把不足的次数用系数
0
0
0 补充,变成两个齐次多项式。
考虑怎样的两项乘完后,结果的幂指数为同一个。即
a
i
x
i
×
b
j
x
j
=
c
k
x
k
a_ix^i\times b_jx^j=c_kx^k
aixi×bjxj=ckxk,
i
,
j
i,j
i,j 怎样取值才能使
k
k
k 为定值。显然
k
=
i
+
j
k=i+j
k=i+j。
那么乘完之后,
C
(
x
)
C(x)
C(x) 就可以知道了。
C
(
x
)
=
∑
i
=
0
2
n
−
2
(
∑
j
=
0
i
a
j
b
i
−
j
)
x
i
=
∑
i
=
0
2
n
−
2
c
i
x
i
\begin{aligned} C(x)&=\sum\limits_{i=0}^{2n-2}(\sum\limits_{j=0}^ia_jb_{i-j})x^i\\ &=\sum\limits_{i=0}^{2n-2}c_ix^i \end{aligned}
C(x)=i=0∑2n−2(j=0∑iajbi−j)xi=i=0∑2n−2cixi
我们把中间的那个加法,即
c
i
=
∑
j
=
0
i
a
j
b
i
−
j
c_i=\sum\limits_{j=0}^ia_jb_{i-j}
ci=j=0∑iajbi−j 叫做卷积。
这样,朴素算法时间复杂度为
O
(
n
2
)
O(n^2)
O(n2),显然不是太好。
但是点值表达式的表示就很好了,直接乘起来就好了嘛!即:
设
A
(
x
)
A(x)
A(x) 的点值表示为
α
=
{
(
x
i
,
A
(
x
i
)
)
∣
i
∈
[
0
,
n
−
1
]
,
i
∈
Z
}
\alpha=\{ (x_i,A(x_i))\mid i\in[0,n-1],i\in \mathbb{Z}\}
α={(xi,A(xi))∣i∈[0,n−1],i∈Z},设
B
(
x
)
B(x)
B(x) 的点值表示为
β
=
{
(
x
i
,
B
(
x
i
)
)
∣
i
∈
[
0
,
n
−
1
]
,
i
∈
Z
}
\beta=\{ (x_i,B(x_i))\mid i\in[0,n-1],i\in \mathbb{Z}\}
β={(xi,B(xi))∣i∈[0,n−1],i∈Z},然后就能得到
C
(
x
)
C(x)
C(x) 的点值表示为:
γ
=
{
(
x
i
,
A
(
x
i
)
B
(
x
i
)
)
∣
i
∈
[
0
,
n
−
1
]
,
i
∈
Z
}
\gamma=\{ (x_i,A(x_i)B(x_i))\mid i\in[0,n-1],i\in \mathbb{Z}\}
γ={(xi,A(xi)B(xi))∣i∈[0,n−1],i∈Z},但是这里只是两多项式的积。并不是真正的点值表示,不过也够用了。
1.4、插值-点值转换
就要接近重点了同志们加油!
先定义两个运算:
1、定义一种运算,可以将系数表示转化为点值表示,记为
DFT
(
A
(
x
)
)
=
α
\text{DFT}(A(x))=\alpha
DFT(A(x))=α;
2、定义其逆运算,即从一个多项式的点值表示确定其系数表示为插值。记为
DFT
−
1
(
α
)
=
IDFT
(
α
)
=
A
(
x
)
\text{DFT}^{-1}(\alpha)=\text{IDFT}(\alpha)=A(x)
DFT−1(α)=IDFT(α)=A(x)。
第一个运算在
x
i
x_i
xi 确定时唯一确定,那么第二个在
x
i
x_i
xi 确定时可不可以唯一确定
A
(
x
)
A(x)
A(x) 呢?
答案是肯定的。我们把这些点代入,可以得出一个
n
n
n 元一次方程组,当
x
i
x_i
xi 均不相同时,可以唯一确定
a
0
,
a
1
,
⋯
,
a
n
a_0,a_1,\cdots ,a_n
a0,a1,⋯,an。(如果有一个
x
i
x_i
xi 相同立刻变成不定方程)
当然可以拿矩阵证明网上多得是。
那么就很尴尬,
DFT
(
A
(
x
)
)
\text{DFT}(A(x))
DFT(A(x)) 的复杂度是
O
(
n
2
)
O(n^2)
O(n2) 的,它的逆运算是
O
(
n
3
)
O(n^3)
O(n3) 的(高斯消元)……还不如不优化了。
所以我们还需要一些姿势。
2、(高中知识)复数
一些奇奇怪怪的问题无法解决,可以把它从实数域扩展到复数域上解决,加入了复数这种类似向量的东西,问题就会变得简单。
噫,向量?点值表示不就是个向量吗!
2.1、定义
一堆定义就不说了,选修 2-2 内容,这里简单的强调一些运算:
i
2
=
−
1
i^2=-1
i2=−1
设
z
1
,
z
2
∈
C
z_1,z_2\in \mathbb{C}
z1,z2∈C,
z
1
=
a
+
b
i
,
z
2
=
c
+
d
i
z_1=a+bi,z_2=c+di
z1=a+bi,z2=c+di,则:
z
1
±
z
2
=
(
a
±
c
)
+
(
b
±
d
)
i
z_1\pm z_2=(a\pm c)+(b\pm d)i
z1±z2=(a±c)+(b±d)i
z
1
×
z
2
=
(
a
c
−
b
d
)
+
(
a
d
+
b
c
)
i
z_1\times z_2=(ac-bd)+(ad+bc)i
z1×z2=(ac−bd)+(ad+bc)i
z
1
z
2
=
a
c
+
b
d
c
2
+
d
2
+
b
c
−
a
d
c
2
+
d
2
i
\frac{z_1}{z_2}=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}i
z2z1=c2+d2ac+bd+c2+d2bc−adi
设
z
∈
C
z\in \mathbb{C}
z∈C,
z
=
a
+
b
i
z=a+bi
z=a+bi,则
z
z
z 的共轭复数为
z
‾
=
a
−
b
i
\overline z=a-bi
z=a−bi
这是学过的,下面就说说没学过的。
2.2、单位根
根?方程的根?什么方程的根?
定义:
x
n
=
1
x^n=1
xn=1 的所有根叫做
n
n
n 次单位根。
当我们把问题引入复数域,
x
x
x 就不能简单计算为
x
=
1
n
=
1
x=\sqrt[n] 1=1
x=n1=1了,而是在复平面内解决问题。
2.3、欧拉公式
一个很 6 的公式。内容是:
e
i
x
=
cos
x
+
i
sin
x
e^{ix}=\cos x+i\sin x
eix=cosx+isinx
式子很简单,证明如下:
由 Taylor 展开:
e
x
=
1
+
x
+
x
2
2
!
+
⋯
+
x
n
n
!
+
⋯
cos
x
=
1
−
x
2
2
!
+
x
4
4
!
−
x
6
6
!
+
⋯
sin
x
=
x
−
x
3
3
!
+
x
5
5
!
−
x
7
7
!
+
⋯
\begin{aligned} e^x&=1+x+\frac{x^2}{2!}+\cdots +\frac{x^n}{n!}+\cdots \\ \cos x&=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+\cdots \\ \sin x&=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\cdots \\ \end{aligned}
excosxsinx=1+x+2!x2+⋯+n!xn+⋯=1−2!x2+4!x4−6!x6+⋯=x−3!x3+5!x5−7!x7+⋯
将 Taylor 展开推广到复数域上,有复数运算
±
i
=
±
i
,
(
±
i
)
2
=
−
1
,
(
±
i
)
3
=
∓
i
,
(
±
i
)
4
=
1
\pm i=\pm i,(\pm i)^2=-1,(\pm i)^3=\mp i,(\pm i)^4=1
±i=±i,(±i)2=−1,(±i)3=∓i,(±i)4=1,可得:
e
i
x
=
1
+
i
x
−
x
2
2
!
−
i
x
3
3
!
+
x
4
4
!
+
i
x
5
5
!
+
⋯
=
(
1
−
x
2
2
!
+
x
4
4
!
−
⋯
)
+
i
(
x
−
x
3
3
!
+
x
5
5
!
−
⋯
)
=
cos
x
+
i
sin
x
\begin{aligned} e^{ix}&=1+ix-\frac{x^2}{2!}-i\frac{x^3}{3!}+\frac{x^4}{4!}+i\frac{x^5}{5!}+\cdots \\ &=(1-\frac{x^2}{2!}+\frac{x^4}{4!}-\cdots)+i(x-\frac{x^3}{3!}+\frac{x^5}{5!}-\cdots)\\ &=\cos x+i\sin x \end{aligned}
eix=1+ix−2!x2−i3!x3+4!x4+i5!x5+⋯=(1−2!x2+4!x4−⋯)+i(x−3!x3+5!x5−⋯)=cosx+isinx
证毕。
将
x
x
x 换成
π
\pi
π,就出现了
e
π
i
=
−
1
e^{\pi i}=-1
eπi=−1,即
e
π
i
+
1
=
0
e^{\pi i}+1=0
eπi+1=0,被誉为上帝创造的公式。
将
x
x
x 换成
2
k
π
(
k
∈
Z
)
2k\pi\ \ (k\in \mathbb{Z})
2kπ (k∈Z),就出现了
e
2
k
π
i
=
1
e^{2k\pi i}=1
e2kπi=1。
我们不是要解
x
n
=
1
x^n=1
xn=1 吗!
那么
x
=
e
2
k
π
i
n
(
k
∈
Z
)
x=e^{\frac{2k\pi i}{n}} \ \ (k \in \mathbb{Z})
x=en2kπi (k∈Z),记这些单位根为
ω
n
k
=
e
2
k
π
i
n
\omega _n^k=e^{\frac{2k\pi i}{n}}
ωnk=en2kπi。
下面就是一些数论问题了……
记
k
=
n
r
+
p
k=nr+p
k=nr+p,其中
p
∈
[
0
,
n
−
1
]
p\in[0,n-1]
p∈[0,n−1],则
ω
n
k
=
e
2
p
π
i
n
\omega _n^k=e^{\frac{2p\pi i}{n}}
ωnk=en2pπi
ω
n
k
=
e
2
k
π
i
n
=
e
2
(
n
r
+
p
)
π
i
n
=
e
2
n
r
π
i
n
+
2
p
π
i
n
=
e
2
n
r
π
i
n
×
e
2
p
π
i
n
=
e
2
π
i
r
×
e
2
p
π
i
n
=
1
r
×
e
2
p
π
i
n
=
e
2
p
π
i
n
\begin{aligned} \omega _n^k&=e^{\frac{2k\pi i}{n}}\\ &=e^{\frac{2(nr+p)\pi i}{n}}\\ &=e^{\frac{2nr\pi i}{n}+\frac{2p\pi i}{n}}\\ &=e^{\frac{2nr\pi i}{n}}\times e^{\frac{2p\pi i}{n}}\\ &=e^{2\pi ir}\times e^{\frac{2p\pi i}{n}}\\ &=1^r\times e^{\frac{2p\pi i}{n}}=e^{\frac{2p\pi i}{n}} \end{aligned}
ωnk=en2kπi=en2(nr+p)πi=en2nrπi+n2pπi=en2nrπi×en2pπi=e2πir×en2pπi=1r×en2pπi=en2pπi
这不一样吗……
所以可以证明
k
∈
[
0
,
n
−
1
]
k\in [0,n-1]
k∈[0,n−1]。也就是
n
n
n 次单位根有
n
n
n 个,不难发现,这
n
n
n 个单位根就是单位圆上的
n
n
n 个
n
n
n 等分点。
于是记这些根为
ω
n
0
,
ω
n
1
,
⋯
,
ω
n
n
−
1
\omega_n^0,\omega_n^1,\cdots ,\omega_n^{n-1}
ωn0,ωn1,⋯,ωnn−1,这些根各不相同,并且在乘法意义下构成一个群。即:
ω
n
i
ω
n
j
=
ω
n
i
+
j
=
ω
n
(
i
+
j
)
m
o
d
n
\omega_n^i\omega_n^j=\omega_n^{i+j}=\omega_n^{(i+j)\bmod n}
ωniωnj=ωni+j=ωn(i+j)modn,且
ω
n
−
1
=
ω
n
n
−
1
\omega_n^{-1}=\omega_n^{n-1}
ωn−1=ωnn−1,证明如下:
ω
n
−
1
=
1
×
ω
n
−
1
=
e
2
π
i
×
e
−
2
π
i
n
=
e
2
π
i
−
−
2
π
i
n
=
e
2
(
n
−
1
)
π
i
n
=
ω
n
n
−
1
\begin{aligned} \omega_n^{-1}&=1\times \omega_n^{-1}\\ &=e^{2\pi i}\times e^{\frac{-2\pi i}{n}}\\ &=e^{2\pi i-\frac{-2\pi i}{n}}\\ &=e^{\frac{2(n-1)\pi i}{n}}=\omega_n^{n-1} \end{aligned}
ωn−1=1×ωn−1=e2πi×en−2πi=e2πi−n−2πi=en2(n−1)πi=ωnn−1
问题得证。
2.4、三个引理
2.4.1、消去引理
内容:当
d
≠
0
d\not =0
d=0 时,
ω
d
n
d
k
=
ω
n
k
\omega_{dn}^{dk}=\omega_n^k
ωdndk=ωnk。
证明:
ω
d
n
d
k
=
e
2
d
k
π
i
d
n
=
e
2
k
π
i
n
=
ω
n
k
\begin{aligned} \omega_{dn}^{dk}&=e^{\frac{2dk\pi i}{dn}}\\ &=e^{\frac{2k\pi i}{n}}=\omega_n^k \end{aligned}
ωdndk=edn2dkπi=en2kπi=ωnk
很简单……也很智障
2.4.2、折半引理
内容:如果
n
>
0
n>0
n>0 且
n
n
n 为偶数,
(
ω
n
k
)
2
=
ω
n
2
k
(\omega_n^k)^2=\omega_{\frac{n}{2}}^k
(ωnk)2=ω2nk。
证明:
(
ω
n
k
)
2
=
(
e
2
k
π
i
n
)
2
=
e
2
k
π
i
n
2
=
ω
n
2
k
\begin{aligned} (\omega_{n}^{k})^2&=(e^{\frac{2k\pi i}{n}})^2\\ &=e^{\frac{2k\pi i}{\frac{n}{2}}}=\omega_{\frac{n}{2}}^k \end{aligned}
(ωnk)2=(en2kπi)2=e2n2kπi=ω2nk
还行吧……可以看出,
n
n
n 个
n
n
n 次单位根平方的集合就是
n
2
\frac{n}{2}
2n 个
n
2
\frac{n}{2}
2n 次单位根的集合,并且每个元素出现两次。可以通过单位根的对称性来理解。
更一般的结论:
ω
n
k
+
n
2
=
ω
n
n
2
ω
n
k
=
−
ω
n
k
\omega_n^{k+\frac{n}{2}}=\omega_n^{\frac{n}{2}}\omega_n^k=-\omega_n^k
ωnk+2n=ωn2nωnk=−ωnk
ω
n
n
2
=
e
π
i
=
−
1
\omega_n^{\frac{n}{2}}=e^{\pi i}=-1
ωn2n=eπi=−1
2.4.3、求和引理
内容:对于
∀
k
∈
(
0
,
n
)
\forall k\in(0,n)
∀k∈(0,n),有
∑
j
=
0
n
−
1
(
ω
n
k
)
j
=
0
\sum\limits_{j=0}^{n-1}(\omega_n^k)^j=0
j=0∑n−1(ωnk)j=0。
证明:
∑
j
=
0
n
−
1
(
ω
n
k
)
j
=
1
−
(
ω
n
k
)
n
1
−
ω
n
k
(
等
比
数
列
求
和
)
=
1
−
(
ω
n
n
)
k
1
−
ω
n
k
=
0
\begin{aligned} \sum\limits_{j=0}^{n-1}(\omega_n^k)^j &=\frac{1-(\omega_{n}^{k})^n}{1-\omega_{n}^{k}}(等比数列求和)\\ &=\frac{1-(\omega_{n}^{n})^k}{1-\omega_{n}^{k}}=0 \end{aligned}
j=0∑n−1(ωnk)j=1−ωnk1−(ωnk)n(等比数列求和)=1−ωnk1−(ωnn)k=0
注意:
ω
n
n
=
1
\omega_n^n=1
ωnn=1,化简省略了部分很多步骤。
现在全都就绪了,终于要到目标了……
3、离散傅里叶变换(DFT)
DFT,即 Discrete Fourier Transform,离散傅里叶变换,就是我们之前定义的那个,现在有了复数,我们就可以在
O
(
n
log
n
)
O(n\log n)
O(nlogn) 的时间内做完 DFT 运算。运用到了 Cooley-Tukey 算法。
我们将多项式按照指数的奇偶分类,记原多项式为
A
(
x
)
A(x)
A(x),那么构造两个多项式:
A
0
(
x
)
=
a
0
+
a
2
x
+
a
4
x
2
+
⋯
+
a
n
−
2
x
n
2
−
1
=
∑
i
=
0
n
2
−
1
a
2
i
x
i
A
1
(
x
)
=
a
1
+
a
3
x
+
a
5
x
2
+
⋯
+
a
n
−
1
x
n
2
−
1
=
∑
i
=
0
n
2
−
1
a
2
i
+
1
x
i
\begin{aligned} A_0(x)&=a_0+a_2x+a_4x^2+\cdots +a_{n-2}x^{\frac{n}{2}-1}\\ &=\sum\limits_{i=0}^{\frac{n}{2}-1}a_{2i}x^i\\ A_1(x)&=a_1+a_3x+a_5x^2+\cdots +a_{n-1}x^{\frac{n}{2}-1}\\ &=\sum\limits_{i=0}^{\frac{n}{2}-1}a_{2i+1}x^i\\ \end{aligned}
A0(x)A1(x)=a0+a2x+a4x2+⋯+an−2x2n−1=i=0∑2n−1a2ixi=a1+a3x+a5x2+⋯+an−1x2n−1=i=0∑2n−1a2i+1xi
通过观察可知,
A
(
x
)
=
A
0
(
x
2
)
+
x
A
1
(
x
2
)
A(x)=A_0(x^2)+xA_1(x^2)
A(x)=A0(x2)+xA1(x2),根据折半引理,我们考虑将
n
n
n 次单位根作为
x
x
x 代入计算。又因为折半引理,我们即可分治操作。
那么,对于
k
<
n
2
k<\frac{n}{2}
k<2n,有:
A
(
ω
n
k
)
=
A
0
(
(
ω
n
k
)
2
)
+
ω
n
k
A
1
(
(
ω
n
k
)
2
)
=
A
0
(
ω
n
2
k
)
+
ω
n
k
A
1
(
ω
n
2
k
)
=
∑
i
=
0
n
2
−
1
a
2
i
ω
n
2
k
i
+
ω
n
k
∑
i
=
0
n
2
−
1
a
2
i
+
1
ω
n
2
k
i
A
(
ω
n
k
+
n
2
)
=
A
0
(
(
ω
n
k
)
2
)
+
ω
n
k
+
n
2
A
1
(
(
ω
n
k
)
2
)
=
A
0
(
ω
n
2
k
)
−
ω
n
k
A
1
(
ω
n
2
k
)
=
∑
i
=
0
n
2
−
1
a
2
i
ω
n
2
k
i
−
ω
n
k
∑
i
=
0
n
2
−
1
a
2
i
+
1
ω
n
2
k
i
\begin{aligned} A(\omega_n^k)&=A_0((\omega_n^k)^2)+\omega_n^kA_1((\omega_n^k)^2)\\ &=A_0(\omega_{\frac{n}{2}}^k)+\omega_n^kA_1(\omega_{\frac{n}{2}}^k)\\ &=\sum\limits_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}+\omega_n^k\sum\limits_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki}\\ A(\omega_n^{k+\frac{n}{2}})&=A_0((\omega_n^k)^2)+\omega_n^{k+\frac{n}{2}}A_1((\omega_n^k)^2)\\ &=A_0(\omega_{\frac{n}{2}}^k)-\omega_n^kA_1(\omega_{\frac{n}{2}}^k)\\ &=\sum\limits_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}-\omega_n^k\sum\limits_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki} \end{aligned}
A(ωnk)A(ωnk+2n)=A0((ωnk)2)+ωnkA1((ωnk)2)=A0(ω2nk)+ωnkA1(ω2nk)=i=0∑2n−1a2iω2nki+ωnki=0∑2n−1a2i+1ω2nki=A0((ωnk)2)+ωnk+2nA1((ωnk)2)=A0(ω2nk)−ωnkA1(ω2nk)=i=0∑2n−1a2iω2nki−ωnki=0∑2n−1a2i+1ω2nki
这样,由奇偶性使得需要代入的值范围减半,递归做下去就行了,时间复杂度为:
T
(
n
)
=
2
T
(
n
2
)
+
O
(
n
)
=
O
(
n
log
n
)
T(n)=2T(\frac{n}{2})+O(n)=O(n\log n)
T(n)=2T(2n)+O(n)=O(nlogn)
因为永远是奇偶合并,如果想要更方便地合并必然
n
n
n 取
2
2
2 的整数次幂,否则合并到一定时候左右不一定都是相同次数(可以考虑一个完全二叉树)。所以做的时候取
2
2
2 的整数次幂作为
n
n
n。
4、逆离散傅里叶变换(IDFT)
IDFT,即 Inverse Discrete Fourier Transform,逆离散傅里叶变换,还是我们之前定义的那个。我们在
O
(
n
log
n
)
O(n\log n)
O(nlogn) 的时间里搞定了 DFT,可是 IDFT 还是个
O
(
n
3
)
O(n^3)
O(n3)的,所以我们继续观(tui)察(da)性(shi)质(zi)。
我们说过,IDFT 相当于解个
n
n
n 元一次方程组,这个方程组写成矩阵形式是这样的:
[
(
ω
n
0
)
0
(
ω
n
0
)
1
(
ω
n
0
)
2
⋯
(
ω
n
0
)
n
−
1
(
ω
n
1
)
0
(
ω
n
1
)
1
(
ω
n
1
)
2
⋯
(
ω
n
1
)
n
−
1
⋮
⋮
⋮
⋱
⋮
(
ω
n
n
−
1
)
0
(
ω
n
n
−
1
)
1
(
ω
n
n
−
1
)
2
⋯
(
ω
n
n
−
1
)
n
−
1
]
×
[
a
0
a
1
⋮
a
n
−
1
]
=
[
A
(
ω
n
0
)
A
(
ω
n
1
)
A
(
ω
n
2
)
⋮
A
(
ω
n
n
−
1
)
]
\begin{bmatrix} (\omega_n^0)^0& (\omega_n^0)^1&(\omega_n^0)^2 &\cdots &(\omega_n^0)^{n-1}\\ (\omega_n^1)^0& (\omega_n^1)^1&(\omega_n^1)^2&\cdots &(\omega_n^1)^{n-1}\\ \vdots&\vdots &\vdots &\ddots &\vdots \\ (\omega_n^{n-1})^0 & (\omega_n^{n-1})^1& (\omega_n^{n-1})^2&\cdots & (\omega_n^{n-1})^{n-1} \end{bmatrix} \times \begin{bmatrix} a_0\\a_1\\ \vdots \\a_{n-1} \end{bmatrix}= \begin{bmatrix} A(\omega_n^0)\\ A(\omega_n^1)\\ A(\omega_n^2)\\ \vdots \\ A(\omega_n^{n-1})\\ \end{bmatrix}
⎣⎢⎢⎢⎡(ωn0)0(ωn1)0⋮(ωnn−1)0(ωn0)1(ωn1)1⋮(ωnn−1)1(ωn0)2(ωn1)2⋮(ωnn−1)2⋯⋯⋱⋯(ωn0)n−1(ωn1)n−1⋮(ωnn−1)n−1⎦⎥⎥⎥⎤×⎣⎢⎢⎢⎡a0a1⋮an−1⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡A(ωn0)A(ωn1)A(ωn2)⋮A(ωnn−1)⎦⎥⎥⎥⎥⎥⎤
这同时也是做 DFT 时我们乘的矩阵。只不过这次
a
i
a_i
ai 变成未知数了。
记上面的系数矩阵为
D
D
D,再考虑以下矩阵:
[
(
ω
n
0
)
0
(
ω
n
0
)
1
(
ω
n
0
)
2
⋯
(
ω
n
0
)
n
−
1
(
ω
n
−
1
)
0
(
ω
n
−
1
)
1
(
ω
n
−
1
)
2
⋯
(
ω
n
−
1
)
n
−
1
⋮
⋮
⋮
⋱
⋮
(
ω
n
−
(
n
−
1
)
)
0
(
ω
n
−
(
n
−
1
)
)
1
(
ω
n
−
(
n
−
1
)
)
2
⋯
(
ω
n
−
(
n
−
1
)
)
n
−
1
]
\begin{bmatrix} (\omega_n^0)^0& (\omega_n^0)^1&(\omega_n^0)^2 &\cdots &(\omega_n^0)^{n-1}\\ (\omega_n^{-1})^0& (\omega_n^{-1})^1&(\omega_n^{-1})^2&\cdots &(\omega_n^{-1})^{n-1}\\ \vdots&\vdots &\vdots &\ddots &\vdots \\ (\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1& (\omega_n^{-(n-1)})^2&\cdots & (\omega_n^{-(n-1)})^{n-1} \end{bmatrix}
⎣⎢⎢⎢⎡(ωn0)0(ωn−1)0⋮(ωn−(n−1))0(ωn0)1(ωn−1)1⋮(ωn−(n−1))1(ωn0)2(ωn−1)2⋮(ωn−(n−1))2⋯⋯⋱⋯(ωn0)n−1(ωn−1)n−1⋮(ωn−(n−1))n−1⎦⎥⎥⎥⎤
记此矩阵为
W
W
W,则记
E
=
D
×
W
E=D\times W
E=D×W。可知:
E
i
j
=
∑
k
=
0
n
−
1
D
i
k
W
k
j
=
∑
k
=
0
n
−
1
ω
n
−
i
k
ω
n
k
j
=
∑
k
=
0
n
−
1
ω
n
k
(
j
−
i
)
\begin{aligned} E_{ij}&=\sum\limits_{k=0}^{n-1}D_{ik}W_{kj}\\ &=\sum\limits_{k=0}^{n-1}\omega_n^{-ik}\omega_n^{kj}\\ &=\sum\limits_{k=0}^{n-1}\omega_n^{k(j-i)} \end{aligned}
Eij=k=0∑n−1DikWkj=k=0∑n−1ωn−ikωnkj=k=0∑n−1ωnk(j−i)
当
i
=
j
i=j
i=j 时:
E
i
j
=
∑
k
=
0
n
−
1
ω
n
k
(
j
−
i
)
=
∑
k
=
0
n
−
1
ω
n
0
=
n
\begin{aligned} E_{ij}&=\sum\limits_{k=0}^{n-1}\omega_n^{k(j-i)}\\ &=\sum\limits_{k=0}^{n-1}\omega_n^0=n\\ \end{aligned}
Eij=k=0∑n−1ωnk(j−i)=k=0∑n−1ωn0=n
当
i
≠
j
i\not =j
i=j 时,由求和引理:
E
i
j
=
∑
k
=
0
n
−
1
ω
n
k
(
j
−
i
)
=
∑
k
=
0
n
−
1
(
ω
n
j
−
i
)
k
=
0
\begin{aligned} E_{ij}&=\sum\limits_{k=0}^{n-1}\omega_n^{k(j-i)}\\ &=\sum\limits_{k=0}^{n-1}(\omega_n^{j-i})^k=0\\ \end{aligned}
Eij=k=0∑n−1ωnk(j−i)=k=0∑n−1(ωnj−i)k=0
记单位矩阵为
I
I
I,所以
I
=
1
n
E
I=\frac{1}{n}E
I=n1E,即
1
n
D
=
W
−
1
\frac{1}{n}D=W^{-1}
n1D=W−1。
将由方程组推得那个矩阵左乘
1
n
D
\frac{1}{n}D
n1D,可得:
[
a
0
a
1
⋮
a
n
−
1
]
=
[
(
ω
n
0
)
0
(
ω
n
0
)
1
(
ω
n
0
)
2
⋯
(
ω
n
0
)
n
−
1
(
ω
n
−
1
)
0
(
ω
n
−
1
)
1
(
ω
n
−
1
)
2
⋯
(
ω
n
−
1
)
n
−
1
⋮
⋮
⋮
⋱
⋮
(
ω
n
−
(
n
−
1
)
)
0
(
ω
n
−
(
n
−
1
)
)
1
(
ω
n
−
(
n
−
1
)
)
2
⋯
(
ω
n
−
(
n
−
1
)
)
n
−
1
]
×
[
A
(
ω
n
0
)
A
(
ω
n
1
)
A
(
ω
n
2
)
⋮
A
(
ω
n
n
−
1
)
]
\begin{bmatrix} a_0\\a_1\\ \vdots \\a_{n-1} \end{bmatrix}= \begin{bmatrix} (\omega_n^0)^0& (\omega_n^0)^1&(\omega_n^0)^2 &\cdots &(\omega_n^0)^{n-1}\\ (\omega_n^{-1})^0& (\omega_n^{-1})^1&(\omega_n^{-1})^2&\cdots &(\omega_n^{-1})^{n-1}\\ \vdots&\vdots &\vdots &\ddots &\vdots \\ (\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1& (\omega_n^{-(n-1)})^2&\cdots & (\omega_n^{-(n-1)})^{n-1} \end{bmatrix} \times \begin{bmatrix} A(\omega_n^0)\\ A(\omega_n^1)\\ A(\omega_n^2)\\ \vdots \\ A(\omega_n^{n-1})\\ \end{bmatrix}
⎣⎢⎢⎢⎡a0a1⋮an−1⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡(ωn0)0(ωn−1)0⋮(ωn−(n−1))0(ωn0)1(ωn−1)1⋮(ωn−(n−1))1(ωn0)2(ωn−1)2⋮(ωn−(n−1))2⋯⋯⋱⋯(ωn0)n−1(ωn−1)n−1⋮(ωn−(n−1))n−1⎦⎥⎥⎥⎤×⎣⎢⎢⎢⎢⎢⎡A(ωn0)A(ωn1)A(ωn2)⋮A(ωnn−1)⎦⎥⎥⎥⎥⎥⎤
这和 DFT 的操作很像,所以用 DFT 的过程实现 IDFT 即可。
负号?把
ω
n
1
\omega_n^1
ωn1 改成
ω
n
−
1
\omega_n^{-1}
ωn−1 即可;
1
n
\frac{1}{n}
n1?做完之后再除个
n
n
n 即可。现在 IDFT 也是
O
(
n
log
n
)
O(n\log n)
O(nlogn) 的了。
至此,所有的时间都是
O
(
n
log
n
)
O(n\log n)
O(nlogn) 的了!即 FFT 的复杂度为
O
(
n
log
n
)
O(n\log n)
O(nlogn)!
5、实现问题
说是一回事,写是另一回事。
奇偶分组这个玩意就很坑,容易发现,设
r
i
r_i
ri 为将
i
i
i 的高位翻转到低位得到的数字(如 0011 就变成 1100),那么初始位置为
i
i
i 的数就会被移动到
r
i
r_i
ri,所以我们可以预先完成移动,然后合并计算。时间就会加快很多了。
不过 FFT 本身常数巨大,主要因为 double
类型的计算时间,还有使用
sin
,
cos
\sin,\cos
sin,cos 函数带来的计算时间问题,所以常数问题一定要注意。
我们看到,DFT 做完后,得到的是将
n
n
n 次单位根代入后表达式的点值表示,IDFT 做完后,得到的是两多项式相乘的系数表示。
由于 FFT 涉及复数和除法,容易出现精度问题,所以考虑一种整数在模的前提下的 FFT,就出现了 FNT,也叫 NTT,即快速数论变换(Fast Number-theory Transform 或 Number Theory Transform)。这里有学习笔记……
如果我们计算卷积的时候
a
,
b
a,b
a,b 的下标不是相加为定值,而是一或二元逻辑运算(与,或,异或等等)后为定值,则可应用快速沃尔什变换(Fast Walsh Hadamard Transform)。学习笔记在这里
FFT真是大毒瘤……
FFT 和 FWT 都在图像处理,音频处理方面有较大应用。
FFT 的板子在这里,以 HDOJ1402 为例……Code
VijosP2000 那道题除了写 FFT 还得压位……或者用一种神奇的
O
(
n
2
)
O(n^2)
O(n2) 乘法:Code
终于把我的 flag 填完了我好开心啊!!!