不忘历史
八百年前学了 FFT \text{FFT} FFT,因vicky过于垃圾,遂放弃。
七百年前重拾 FFT \text{FFT} FFT,勉强搞懂了它的递归写法,因vicky再一次懒癌附体,遂连板题都没写就弃疗了。
历史的今天 (是今天才怪),vicky重学了
FFT
\text{FFT}
FFT但仍然没写板题。
今天的昨天,vicky和蔼的教练给她讲了
FWT
\text{FWT}
FWT,vicky当然是继续弃疗。
于是vicky决定复习一下 FFT \text{FFT} FFT。
这时报应来了,vicky发现她的FFT学得很杂。(就vicky这有兴致就看看没兴致就不学了的学法她学得不杂才怪
于是vicky开了这个坑但是选择继续摆烂。
让我们采访一下vicky: QwQ
-
Q:请问您什么时候更完?
-
A:vicky:别指望我更完了,关于vicky,她死了。
u p d 23.4.27 upd 23.4.27 upd23.4.27:what should i do what should i do what should i do what should i do我学不会FWT呜呜呜呜呜
有些不得不说在前面的话
因为这些都是我边口胡边学的,其中一部分内容甚至是我在晚上睡前冥想的时候脑子里面突然冒出来、第二天上竞赛课时加以草率的验证才得到的。
所以这篇blog难免会出现很多很多的锅,各位发现了我在一些内容上理解/描述有误的话请私信or直接在评论区下面给我指出来,真的很感谢www QwQ \text{QwQ} QwQ。
因为这篇文章主要是让我用来复习 FFT&NTT&FWT \text{FFT\&NTT\&FWT} FFT&NTT&FWT 用的,所以一些我觉得无关紧要的概念名词在这篇blog中有一定概率不会被我说起。(即提到相关概念时我也不一定会专门告诉你这个概念的专有名词叫啥)
Reference meterial \text{Reference meterial} Reference meterial
-
小学生都能看懂的FFT!!!—— Written by \text{Written by } Written by 胡小兔
-
教练PPT(不知道教练给不给放所以不放了)
O(n log 2 n) \text{O(n log}_2\text{n)} O(n log2n) 的多项式乘法: FFT \text{FFT} FFT
Part 0. \text{Part 0.} Part 0. 点值表示&系数表示
0.1 \text{0.1} 0.1 点值表示
-
对于次数为 n n n 的多项式,我们很容易发现,只要我们代入 n + 1 n+1 n+1 个不同的 x x x 值求得 n + 1 n+1 n+1 个函数值,这 n + 1 n+1 n+1 个函数值完全可以代替原来多项式的系数表示,让我们能够通过拉格朗日插值法 O ( n ) O(n) O(n) 求出任意一个函数值。
-
也就是说想要通过函数值表示一个次数为 n n n 的多项式我们只需要 n + 1 n+1 n+1 个 x x x 的取值不同的函数值就可以表示出这个多项式。
-
此时我们用 n + 1 n+1 n+1 个点值表示一个次数为 n n n 的多项式,故称其为点值表示。
PS:拉插不是
F
F
T
FFT
FFT 的重点,所以这里为了日后vicky复习
F
F
T
FFT
FFT 时不犯迷糊着想,就不对拉插进行过多阐述。想学的话自己搜blog吧,这玩意还是蛮好学的。
0.2 \text{0.2} 0.2 系数表示
我写由题意得应该不会被骂吧。- 具体例子就是: A ( x ) = a 0 + a 1 x + x 2 x 2 + ⋯ + a n x n A(x)=a_0+a_1x+x_2x^2+\cdots+a_nx^n A(x)=a0+a1x+x2x2+⋯+anxn,这就是系数表示。
Part 1.FFT \text{Part 1.FFT} Part 1.FFT 思路概述
现在我们有三个多项式:
A ( x ) = a 0 + a 1 x + a 2 x 2 + ⋯ + a n − 1 x n − 1 , B ( x ) = b 0 + ⋯ + b n − 1 x n − 1 A(x)=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1},B(x)=b_0+\cdots+b_{n-1}x^{n-1} A(x)=a0+a1x+a2x2+⋯+an−1xn−1,B(x)=b0+⋯+bn−1xn−1。
C ( x ) = A ( x ) ∗ B ( x ) = ∑ i = 0 2 n − 2 c i x i C(x)=A(x)*B(x)=\sum\limits_{i=0}^{2n-2}c_ix^i C(x)=A(x)∗B(x)=i=0∑2n−2cixi 。
对于多项式 C ( x ) = A ( x ) ∗ B ( x ) C(x)=A(x)*B(x) C(x)=A(x)∗B(x) 而言,我们如果用点值表示这个多项式,因为它的次数为 2 n − 2 2n-2 2n−2 ,所以我们只需要求出 2 n − 1 2n-1 2n−1 个 C ( x ) C(x) C(x) 的函数值,就可以在 O ( n ) O(n) O(n) 的时间复杂度内求出多项式 A ( x ) A(x) A(x) 和 B ( x ) B(x) B(x) 的乘积。
那么此时求两个多项式相乘的思路其实就很明确了:
- 将多项式 A ( x ) A(x) A(x) 和 B ( x ) B(x) B(x) 由系数表示转化为点值表示,然后然后将它们俩相乘。【每求解一个 C ( x ) C(x) C(x) 所需要的时间均为 O ( n ) O(n) O(n) , C ( x ) C(x) C(x) 若要使用点值表示则需要求 2 n + 1 2n+1 2n+1 个 C ( x ) C(x) C(x) 的点值,所以这里如果使用朴素算法转化点值表示的话时间复杂度就是 O ( n 2 ) O(n^2) O(n2) 】
- 将 C ( x ) C(x) C(x) 的点值表示转化为系数表示,然后我们就可以非常自然地求得任意 C ( x ) C(x) C(x) 了。(这里用拉插的话时间复杂度又是 O ( n 2 ) O(n^2) O(n2) )
- 那么, F F T FFT FFT 要做的其实就是用一些神奇技巧来将两个转化过程简化。
Part 2. \text{Part 2.} Part 2. 关于复数
这里坚持住的话后面的 F F T FFT FFT 理解起来其实就很简单了。
首先我们了解一些概念。
对于一个复数
a
+
b
i
a+bi
a+bi,我们有两种理解方式:(因为这些都是vicky边口胡边自学的,所以有错记得跟我说一声
- 一个起点为原点的向量。
- 平面直角坐标系上的点 ( a , b ) (a,b) (a,b)。
两个复数相乘的规则如下:模长相乘,幅角相加。
用代数式进行表示的话就是:
( a + b i ) ( c + d i ) = a c + b d i 2 + a d i + c b i = ( a c − b d ) + ( a d + c b ) i (a+bi)(c+di)=ac+bdi^2+adi+cbi=(ac-bd)+(ad+cb)i (a+bi)(c+di)=ac+bdi2+adi+cbi=(ac−bd)+(ad+cb)i
模长就是点 ( a , b ) (a,b) (a,b) 到原点的距离,这里也可以理解为向量的模长。幅角为该复数所对应的向量与 x x x 轴的逆时针夹角。
如下图:
图中 P 6 P_6 P6 的幅角为 270 270 270 度, ( P 3 ) 2 = P 6 (P_3)^2=P_6 (P3)2=P6。
幅角相加的概念要理解到位。
介绍了基本虚数的概念之后,我们此时再引入一个新的概念: ω n i \omega _n^i ωni。
我们以原点为圆心作一个圆,取该圆与 x x x 轴正半轴的交点所表示的那个复数作为第一个点,我们将其命名为 ω n 0 \omega_n^0 ωn0 。
我们将这个圆平均分为 n n n 份,每一份都代表着一个在圆上的复数,即 ω n i \omega_n^i ωni ,我们将复数 ω n 0 \omega_n^0 ωn0 看作一个向量,对它进行若干次旋转,每次旋转都转 360 n \frac{360}{n} n360度,那么当我们进行了 i i i 次旋转之后我们就会得到新的复数 ω n i % n \omega_n^{i\%n} ωni%n ,其幅角为 360 n ⋅ i \frac{360}{n}\cdot i n360⋅i 度。
如上图, P 0 P_0 P0 为 ω 8 0 \omega_8^0 ω80 , P 5 P_5 P5 为 ω 8 5 \omega_8^5 ω85 ,即 ω 8 0 \omega_8^0 ω80 经过 5 5 5 次旋转之后得到的幅角为 360 8 ∗ 5 = 225 \frac{360}{8}*5=225 8360∗5=225 度的复数。
总结一下, ω n i \omega_n^i ωni 的意义即为:
- 将一个以原点为圆心的圆等分为 n n n 份,每一份对应圆上的一个代表复数的点。
- 定义该圆与x轴正半轴的交点为 ω n 0 \omega_n^0 ωn0。
- 该复数的幅角度数为 360 n ⋅ i \frac{360}{n}\cdot i n360⋅i 度,也可以理解为它在 ω n 0 \omega_n^0 ωn0 幅角的基础上加了 i i i 份度数为 360 n \frac{360}{n} n360 度的幅角(即 ω n 0 \omega_n^0 ωn0 等角度地旋转 i i i 次得到 ω n i \omega_n^i ωni )。
那么,这个 ω n k \omega_n^k ωnk 该怎么计算呢?
套公式即可:
ω
n
k
=
cos
(
2
π
⋅
k
n
)
+
i
⋅
sin
(
2
π
⋅
k
n
)
\omega_n^k=\cos(2\pi\cdot \frac{k}{n})+i\cdot \sin(2\pi\cdot \frac{k}{n})
ωnk=cos(2π⋅nk)+i⋅sin(2π⋅nk)
具体为啥是这个公式,只要稍微了解一下三角函数的弧度制即可理解该公式。
且根据复数相乘幅角相加的规则以及 ω n i \omega_n^i ωni 的定义,我们很容易得到不同复数间的一些关系:
-
ω n 2 i = ( ω n i ) 2 = ω n 2 i \omega_n^{2i}=(\omega_n^i)^2=\omega_{\frac{n}{2}}^i ωn2i=(ωni)2=ω2ni
-
ω n i = ω n i % n \omega_n^i=\omega_n^{i\%n} ωni=ωni%n
-
ω n i + n 2 = − ω n i \omega_n^{i+\frac{n}{2}}=-\omega_n^i ωni+2n=−ωni
垃圾markdown逼我用三级标题是吧。
其中最后一条关系还牵及到一个概念:共轭复数。
我们称
ω
n
i
+
n
2
\omega_n^{i+\frac{n}{2}}
ωni+2n 和
ω
n
i
\omega_n^i
ωni 互为共轭复数。正式一点地描述,即实部相等,虚部互为相反数的两复数互为共轭复数。从几何意义上来说,两个共轭复数所对应的向量是共线的。(其实共不共线什么的也不重要就是了
上面的两个式子直接套 ω n i \omega_n^i ωni 上下标的定义就可以很容易证得了吧。QwQ
离散傅里叶变换
我知道你很慌但你先别慌。
离散傅里叶变换其实只是一个概念而已,没啥难理解的地方。
上面说过了,对于一个次数为 n n n 的多项式 f ( x ) f(x) f(x) ,若我们要对其用点值表示,则需要取 n + 1 n+1 n+1 个不同的 x x x 值代入得到 n + 1 n+1 n+1 个函数值来表示这个函数 f ( x ) f(x) f(x) 。
那么这 n + 1 n+1 n+1 个不同的 x x x 值我们为何不去取 ( ω n + 1 0 , ⋯ , ω n + 1 n ) (\omega_{n+1}^0,\cdots,\omega_{n+1}^{n}) (ωn+10,⋯,ωn+1n) 上的值呢?
因为正经人谁会想到这种奇怪的取法啊喂! QAQ
离散傅里叶变换的定义其实就是对于一个次数为 n n n 的函数 f ( x ) f(x) f(x) ,我们取函数 f ( x ) f(x) f(x) 在 ( ω n + 1 0 , ⋯ , ω n + 1 n ) (\omega_{n+1}^0,\cdots,\omega_{n+1}^{n}) (ωn+10,⋯,ωn+1n) 上的取值作为点值表示。
简而论之,一个函数在 n + 1 n+1 n+1 个特殊位置上取值并将其作为该函数的点值表示,这个点值表示就是离散傅里叶变换。
再说一遍!一个函数的离散傅里叶变换是该函数的点值表示而不是单独一个点值!(写给我自己看的,没有任何侮辱各位智商的意思在。
Part 3. \text{Part 3.} Part 3. 复数在 FFT \text{FFT} FFT 的两个变换过程中的运用
搞定了复数!我们其实离搞懂 FFT \text{FFT} FFT 就差推几个式子的功夫啦!QwQ
下面的式子推导对我这样的数学弱智来说都还是很友好的,所以说其他人完全不用怕搞不懂 FFT \text{FFT} FFT 的式子推导就是了。
实在推不出来你把结论记下来感觉其实也没啥问题。
Part 3.0 \text{Part 3.0} Part 3.0 回顾
咱们先回顾一下 FFT \text{FFT} FFT 的两个转化过程:
- 系数表示 -> 点值表示。(具体来说是利用多项式 A ( x ) A(x) A(x) 和 B ( x ) B(x) B(x) 的系数表示求得 C ( x ) C(x) C(x) 的点值表示。)
- 点值表示 -> 系数表示。(将 C ( x ) C(x) C(x) 的点值表示转化为系数表示。)
Part 3.1 \text{Part 3.1} Part 3.1 系数表示 -> 点值表示
我们现有多项式 C ( x ) = A ( x ) ∗ B ( x ) = c 0 + ⋯ + c 2 n − 2 ⋅ x 2 n − 2 C(x)=A(x)*B(x)=c_0+\cdots+c_{2n-2}\cdot x^{2n-2} C(x)=A(x)∗B(x)=c0+⋯+c2n−2⋅x2n−2。
现在我们再弄两个多项式出来:
C
1
(
x
)
=
c
0
+
c
2
⋅
x
+
⋯
+
c
2
n
−
2
⋅
x
n
−
1
C_1(x)=c_0+c_2\cdot x+\cdots+c_{2n-2}\cdot x^{n-1}
C1(x)=c0+c2⋅x+⋯+c2n−2⋅xn−1
C
2
(
x
)
=
c
1
+
c
3
⋅
x
+
⋯
+
c
2
n
−
1
⋅
x
n
−
1
C_2(x)=c_1+c_3\cdot x+\cdots+c_{2n-1}\cdot x^{n-1}
C2(x)=c1+c3⋅x+⋯+c2n−1⋅xn−1
此时我们就可以将 C ( x ) C(x) C(x) 表示如下:
C ( x ) = C 1 ( x 2 ) + x ⋅ C 2 ( x 2 ) C(x)=C_1(x^2)+x\cdot C_2(x^2) C(x)=C1(x2)+x⋅C2(x2)
那么若我们将 C ( x ) C(x) C(x) 的系数表示转化为 C ( x ) C(x) C(x) 对应的离散傅里叶变换。
即系数表示 -> 点值表示。
为了下面的式子推导看上去更简洁,我们将 C ( x ) C(x) C(x) 的次数设为 n − 1 n-1 n−1 ,即函数 C ( x ) C(x) C(x) 的离散傅里叶变换在 ( ω n 0 , ⋯ , ω n n ) (\omega_n^0,\cdots,\omega_n^n) (ωn0,⋯,ωnn) 处取值。(不然满屏的 ω n + 1 i \omega_{n+1}^i ωn+1i谁受得了啊喂!QAQ
那么对于任意 C ( x ) C(x) C(x) 在 ω n i \omega_n^i ωni 处的点值求值,我们可以将点值求值过程分类讨论后进行转化:
温馨提示: ω n 0 = ω n n = 1 , ω n n 2 = − 1 \omega_n^0=\omega_n^n=1,\omega_n^\frac{n}{2}=-1 ωn0=ωnn=1,ωn2n=−1
- 0 ≤ i < n 2 0\leq i<\frac{n}{2} 0≤i<2n:
C
(
ω
n
i
)
=
C
1
(
(
ω
n
i
)
2
)
+
ω
n
i
⋅
C
2
(
(
ω
n
i
)
2
)
C(\omega_n^i)=C_1((\omega_n^i)^2)+\omega_n^i\cdot C_2((\omega_n^i)^2)
C(ωni)=C1((ωni)2)+ωni⋅C2((ωni)2)
=
C
1
(
ω
n
2
i
)
+
ω
n
i
⋅
C
2
(
ω
n
2
i
)
=C_1(\omega_n^{2i})+\omega_n^i\cdot C_2(\omega_n^{2i})
=C1(ωn2i)+ωni⋅C2(ωn2i)
=
C
1
(
ω
n
2
i
)
+
ω
n
i
⋅
C
2
(
ω
n
2
i
)
=C_1(\omega_\frac{n}{2}^{i})+\omega_n^i\cdot C_2(\omega_\frac{n}{2}^{i})
=C1(ω2ni)+ωni⋅C2(ω2ni)
- n 2 ≤ i ≤ n \frac{n}{2}\leq i\leq n 2n≤i≤n,设 i = k + n 2 ( 0 ≤ k ≤ n 2 ) i=k+\frac{n}{2}(0\leq k\leq \frac{n}{2}) i=k+2n(0≤k≤2n) :
C
(
ω
n
i
)
=
C
(
ω
n
k
+
n
2
)
C(\omega_n^i)=C(\omega_n^{k+\frac{n}{2}})
C(ωni)=C(ωnk+2n)
=
C
1
(
ω
n
2
k
+
n
)
+
ω
n
k
+
n
2
⋅
C
2
(
ω
n
2
k
+
n
)
=C_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}\cdot C_2(\omega_n^{2k+n})
=C1(ωn2k+n)+ωnk+2n⋅C2(ωn2k+n)
=
C
1
(
ω
n
2
k
+
n
2
)
+
ω
n
k
+
n
2
⋅
C
2
(
ω
n
2
k
+
n
2
)
=C_1(\omega_{\frac{n}{2}}^{k+{\frac{n}{2}}})+\omega_n^{k+\frac{n}{2}}\cdot C_2(\omega_{\frac{n}{2}}^{k+{\frac{n}{2}}})
=C1(ω2nk+2n)+ωnk+2n⋅C2(ω2nk+2n)
=
C
1
(
ω
n
2
k
⋅
ω
n
n
)
+
ω
n
k
+
n
2
⋅
C
2
(
ω
n
2
k
⋅
ω
n
n
)
=C_1(\omega_{\frac{n}{2}}^{k}\cdot \omega_n^n)+\omega_n^{k+\frac{n}{2}}\cdot C_2(\omega_{\frac{n}{2}}^{k}\cdot \omega_n^n)
=C1(ω2nk⋅ωnn)+ωnk+2n⋅C2(ω2nk⋅ωnn)
=
C
1
(
ω
n
2
k
)
+
(
ω
n
k
⋅
ω
n
n
2
)
⋅
C
2
(
ω
n
2
k
)
=C_1(\omega_{\frac{n}{2}}^{k})+(\omega_n^{k}\cdot \omega_n^{\frac{n}{2}})\cdot C_2(\omega_{\frac{n}{2}}^{k})
=C1(ω2nk)+(ωnk⋅ωn2n)⋅C2(ω2nk)
=
C
1
(
ω
n
2
k
)
+
(
ω
n
k
⋅
−
1
)
⋅
C
2
(
ω
n
2
k
)
=C_1(\omega_{\frac{n}{2}}^{k})+(\omega_n^{k}\cdot -1)\cdot C_2(\omega_{\frac{n}{2}}^{k})
=C1(ω2nk)+(ωnk⋅−1)⋅C2(ω2nk)
=
C
1
(
ω
n
2
k
)
−
ω
n
k
⋅
C
2
(
ω
n
2
k
)
=C_1(\omega_{\frac{n}{2}}^{k})- \omega_n^{k}\cdot C_2(\omega_{\frac{n}{2}}^{k})
=C1(ω2nk)−ωnk⋅C2(ω2nk)
我们对比一下两种情况下 C ( ω n i ) C(\omega_n^i) C(ωni) 最终化出来的式子:
- C ( ω n i ) = C 1 ( ω n 2 i ) + ω n i ⋅ C 2 ( ω n 2 i ) C(\omega_n^i)=C_1(\omega_\frac{n}{2}^{i})+\omega_n^i\cdot C_2(\omega_\frac{n}{2}^{i}) C(ωni)=C1(ω2ni)+ωni⋅C2(ω2ni)
- C ( ω n i ) = C 1 ( ω n 2 k ) − ω n k ⋅ C 2 ( ω n 2 k ) C(\omega_n^i)=C_1(\omega_{\frac{n}{2}}^{k})- \omega_n^{k}\cdot C_2(\omega_{\frac{n}{2}}^{k}) C(ωni)=C1(ω2nk)−ωnk⋅C2(ω2nk)
发现没有,我们通过式子推导,使得我们只需要求函数 C 1 , , C 2 C_1,,C_2 C1,,C2 在 ( ω n 2 0 , ⋯ , ω n 2 n 2 − 1 ) (\omega_{\frac{n}{2}}^0,\cdots,\omega_{\frac{n}{2}}^{\frac{n}{2}-1}) (ω2n0,⋯,ω2n2n−1)处的取值就可以求得函数 C C C 的点值了。
现在问题就变得很简单了,我们只需递推实现即可,碰到递归边界 n = 1 n=1 n=1 时我们直接套 C ( ω n i ) = A ( ω n i ) ∗ B ( ω n i ) C(\omega_n^i)=A(\omega_n^i)*B(\omega_n^i) C(ωni)=A(ωni)∗B(ωni) 求值返回即可。
现在我们看回时间复杂度。
对于每一个 C ( ω n i ) ( 0 ≤ i < n ) C(\omega_n^i)(0\leq i<n) C(ωni)(0≤i<n) ,如果用递归实现,递归边界为显然为 C ( ω n i ) C(\omega_n^i) C(ωni) 中的 n n n 的值为 1 1 1 。
每次 C ( ω n i ) C(\omega_n^i) C(ωni) 中的 n n n 都会缩小一半,不难得到我们要进行 l o g 2 n log_2n log2n 次递归,即每次求 C ( ω n i ) C(\omega_n^i) C(ωni) 的点值的时间复杂度为 O(log 2 n) \text{O(log}_2\text{n)} O(log2n),我们一共要求的点值总数与 n n n 成正比,所以求 C ( x ) C(x) C(x) 的点值表示所需的时间复杂度为 O(n log n) \text{O(n log n)} O(n log n) 。
但是,我们不能忽视的是, O(nlog 2 n) \text{O(nlog}_2\text{n)} O(nlog2n) 仅仅只是递归实现 FFT \text{FFT} FFT 的渐进时间复杂度而已,也就是说,递归实现的 FFT \text{FFT} FFT 在某些凉心畜题人的数据下还是很容易被卡成 T L E TLE TLE 的。
事实上我们有更快的非递归
FFT
\text{FFT}
FFT 实现方式,且vicky个人认为非递归版本并不比递归版本难写多少。
都学到这里了不把 FFT \text{FFT} FFT 的优化学完不就前功尽弃了嘛!QAQ
但是为了日后vicky看这篇帖子不会被冲昏头脑,所以我们先把 FFT \text{FFT} FFT 的基本思想讲完再进一步对 FFT \text{FFT} FFT 的代码实现进行讲解。QwQ
Part 3.2 \text{Part 3.2} Part 3.2 点值表示 -> 系数表示
这部分主要也是推式子,而且相对上一部分来说这部分的式子推导还是友好很多的。
但是上一部分的其实也不难,直接套
ω
\omega
ω 的定义就能推出来了对吧。
说实话其实这部分直接记结论也完全没有问题的。
在上一部分中我们求的是 C ( x ) C(x) C(x) 在 ( ω n 0 , ⋯ , ω n n − 1 ) (\omega_n^0,\cdots,\omega_n^{n-1}) (ωn0,⋯,ωnn−1) 处取值的点值表示,也就是它的离散傅里叶变换,那么我们是否可以这些在点值上进行一些瞎搞去求得 C ( x ) C(x) C(x) 的系数表示呢?
答案是可以的。
为了下面的式子推导相对简洁,我们先将 C ( x ) C(x) C(x) 的离散傅里叶变换 ( C ( ω n 0 ) , ⋯ , C ( ω n n − 1 ) ) (C(\omega_n^0),\cdots,C(\omega_n^{n-1})) (C(ωn0),⋯,C(ωnn−1)) 分别设为数组 ( d 0 , d 1 , ⋯ , d n − 1 ) (d_0,d_1,\cdots,d_{n-1}) (d0,d1,⋯,dn−1) 。
我们将 C ( x ) C(x) C(x) 的离散傅里叶变换 ( d 0 , ⋯ , d n − 1 ) (d_0,\cdots,d_{n-1}) (d0,⋯,dn−1) 作为另一个次数为 n − 1 n-1 n−1 的多项式 D ( x ) D(x) D(x) 的系数表示,即:
D
(
x
)
=
C
(
ω
n
0
)
+
C
(
ω
n
1
)
⋅
x
+
⋯
+
C
(
ω
n
n
−
1
)
⋅
x
n
−
1
D(x)=C(\omega_n^0)+C(\omega_n^1)\cdot x+\cdots+C(\omega_n^{n-1})\cdot x^{n-1}
D(x)=C(ωn0)+C(ωn1)⋅x+⋯+C(ωnn−1)⋅xn−1
=
d
0
+
d
1
x
+
⋯
+
d
n
−
1
x
n
−
1
=d_0+d_1x+\cdots+d_{n-1}x^{n-1}
=d0+d1x+⋯+dn−1xn−1
根据 d i d_i di 的定义,显然有 d i = C ( ω n i ) = ∑ j = 0 n − 1 c j ⋅ ( ω n i ) j d_i=C(\omega_n^i)=\sum\limits_{j=0}^{n-1}c_j\cdot (\omega_n^i)^j di=C(ωni)=j=0∑n−1cj⋅(ωni)j。
我们再将 D ( x ) D(x) D(x) 的离散傅里叶变换表示成 ( e 0 , ⋯ , e n − 1 ) (e_0,\cdots,e_{n-1}) (e0,⋯,en−1),但是这里需要注意的是,我们不再是在 ( ω n 0 , ⋯ , ω n n − 1 ) (\omega_n^0,\cdots,\omega_n^{n-1}) (ωn0,⋯,ωnn−1) 处取值,而是在 ( ω n 0 , ω n − 1 , ⋯ , ω n 1 − n ) (\omega_n^{0},\omega_n^{-1},\cdots,\omega_n^{1-n}) (ωn0,ωn−1,⋯,ωn1−n) 处进行取值,即 e k = D ( ω n − k ) = ∑ i = 0 n − 1 d i ⋅ ( ω n − k ) i e_k=D(\omega_n^{-k})=\sum\limits_{i=0}^{n-1}d_i\cdot (\omega_n^{-k})^i ek=D(ωn−k)=i=0∑n−1di⋅(ωn−k)i 。
那么对于每一个 e i e_i ei ,我们可以尝试一下对它进行式子推导,看看能否找出一种通过 e i e_i ei 快速求出 c i c_i ci 的办法:
e
k
=
D
(
ω
n
−
k
)
=
∑
i
=
0
n
−
1
d
i
⋅
(
ω
n
−
k
)
i
=
∑
i
=
0
n
−
1
D
(
ω
n
i
)
⋅
(
ω
n
−
k
)
i
e_k=D(\omega_n^{-k})=\sum\limits_{i=0}^{n-1}d_i\cdot (\omega_n^{-k})^i=\sum\limits_{i=0}^{n-1}D(\omega_n^{i})\cdot (\omega_n^{-k})^i
ek=D(ωn−k)=i=0∑n−1di⋅(ωn−k)i=i=0∑n−1D(ωni)⋅(ωn−k)i
=
∑
i
=
0
n
−
1
(
∑
j
=
0
n
−
1
c
j
⋅
(
ω
n
i
)
j
)
⋅
(
ω
n
−
k
)
i
=\sum\limits_{i=0}^{n-1}\left(\sum\limits_{j=0}^{n-1}c_j\cdot (\omega_n^i)^j\right)\cdot (\omega_n^{-k})^i
=i=0∑n−1(j=0∑n−1cj⋅(ωni)j)⋅(ωn−k)i
=
∑
i
=
0
n
−
1
∑
j
=
0
n
−
1
c
j
⋅
(
ω
n
i
)
j
⋅
(
ω
n
−
k
)
i
=
∑
i
=
0
n
−
1
∑
j
=
0
n
−
1
c
j
⋅
ω
n
i
⋅
j
⋅
ω
n
−
k
⋅
i
=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}c_j\cdot (\omega_n^i)^j\cdot (\omega_n^{-k})^i=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}c_j\cdot \omega_n^{i\cdot j}\cdot \omega_n^{-k\cdot i}
=i=0∑n−1j=0∑n−1cj⋅(ωni)j⋅(ωn−k)i=i=0∑n−1j=0∑n−1cj⋅ωni⋅j⋅ωn−k⋅i
=
∑
i
=
0
n
−
1
∑
j
=
0
n
−
1
c
j
⋅
ω
n
i
⋅
(
j
−
k
)
=
∑
i
=
0
n
−
1
∑
j
=
0
n
−
1
c
j
⋅
(
ω
n
j
−
k
)
i
=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}c_j\cdot \omega_n^{i\cdot (j-k)}=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}c_j\cdot (\omega_n^{j-k})^i
=i=0∑n−1j=0∑n−1cj⋅ωni⋅(j−k)=i=0∑n−1j=0∑n−1cj⋅(ωnj−k)i
=
∑
j
=
0
n
−
1
c
j
⋅
∑
i
=
0
n
−
1
(
ω
n
j
−
k
)
i
=\sum\limits_{j=0}^{n-1}c_j\cdot \sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i
=j=0∑n−1cj⋅i=0∑n−1(ωnj−k)i
k k k 显然可以理解成一个定值,我们枚举到每一个 j j j 的时候, j j j 其实也算作是一个定值,此时不难发现 ∑ i = 0 n − 1 ( ω n j − k ) i \sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i i=0∑n−1(ωnj−k)i 其实就是个等比数列求和。
因为 j j j 的值是在 [ 0 , n − 1 ] [0,n-1] [0,n−1] 范围内进行枚举的,因此我们可以对 j − k j-k j−k 进行分类讨论,然后再求 ∑ i = 0 n − 1 ( ω n j − k ) i \sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i i=0∑n−1(ωnj−k)i 的值:
- j − k = 0 j-k=0 j−k=0,即 j = k j=k j=k :
此时显然有:
∑ i = 0 n − 1 ( ω n j − k ) i = ∑ i = 0 n − 1 ( ω n 0 ) i = ∑ i = 0 n − 1 1 i = n \sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i=\sum\limits_{i=0}^{n-1}(\omega_n^0)^i=\sum\limits_{i=0}^{n-1}1^i=n i=0∑n−1(ωnj−k)i=i=0∑n−1(ωn0)i=i=0∑n−11i=n
- j − k ≠ 0 j-k\neq 0 j−k=0:
根据等比数列求和公式 s u m = a 1 ⋅ ( 1 − q n ) 1 − q sum=\dfrac{a_1\cdot (1-q^n)}{1-q} sum=1−qa1⋅(1−qn) ,我们易得:
PS:公式中的 a 1 a_1 a1 为数列首项, q q q 为公比, n n n 为项数。
为了下面的式子推导看起来更简洁一些,我们设 W = j − k W=j-k W=j−k 。
注意,下面式子中的
j
j
j 与
k
k
k 我们均看作定值。(不知道为啥的应该好好反思一下自己到底看没看懂上面的内容了……QwQ
∑
i
=
0
n
−
1
(
ω
n
j
−
k
)
i
=
∑
i
=
0
n
−
1
(
ω
n
W
)
i
\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i=\sum\limits_{i=0}^{n-1}(\omega_n^W)^i
i=0∑n−1(ωnj−k)i=i=0∑n−1(ωnW)i
=
1
⋅
(
1
−
(
ω
n
W
)
n
)
1
−
ω
n
W
=
1
−
(
ω
n
n
)
W
1
−
ω
n
W
=\dfrac{1\cdot (1-(\omega_n^W)^n)}{1-\omega_n^W}=\dfrac{1-(\omega_n^n)^W}{1-\omega_n^W}
=1−ωnW1⋅(1−(ωnW)n)=1−ωnW1−(ωnn)W
=
1
−
1
W
1
−
ω
n
W
=
1
−
1
1
−
ω
n
W
=
0
=\dfrac{1-1^W}{1-\omega_n^W}=\dfrac{1-1}{1-\omega_n^W}=0
=1−ωnW1−1W=1−ωnW1−1=0
也就是说,当且仅当 ∑ j = 0 n − 1 \sum\limits_{j=0}^{n-1} j=0∑n−1 枚举到 k k k 时, ∑ i = 0 n − 1 ( ω n j − k ) i \sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i i=0∑n−1(ωnj−k)i 才会有值,且此时 ∑ i = 0 n − 1 ( ω n j − k ) i \sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i i=0∑n−1(ωnj−k)i 的值为 n n n。
那么我们就可以进一步转化 e k e_k ek :
e k = ∑ j = 0 n − 1 c j ⋅ ∑ i = 0 n − 1 ( ω n j − k ) i e_k=\sum\limits_{j=0}^{n-1}c_j\cdot \sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i ek=j=0∑n−1cj⋅i=0∑n−1(ωnj−k)i
我们按上面分类讨论的两种情况进行求值:
e
k
=
(
∑
j
=
0
,
j
≠
k
n
−
1
c
j
⋅
∑
i
=
0
n
−
1
(
ω
n
j
−
k
)
i
)
+
c
k
⋅
∑
i
=
0
n
−
1
(
ω
n
0
)
i
e_k=\left(\sum\limits_{j=0,j\neq k}^{n-1}c_j \cdot \sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i\right)+c_k\cdot \sum\limits_{i=0}^{n-1}(\omega_n^0)^i
ek=
j=0,j=k∑n−1cj⋅i=0∑n−1(ωnj−k)i
+ck⋅i=0∑n−1(ωn0)i
=
(
∑
j
=
0
,
j
≠
k
n
−
1
c
j
⋅
0
)
+
c
k
⋅
n
=\left(\sum\limits_{j=0,j\neq k}^{n-1}c_j \cdot 0\right)+c_k\cdot n
=
j=0,j=k∑n−1cj⋅0
+ck⋅n
=
c
k
⋅
n
=c_k\cdot n
=ck⋅n
这个时候要怎么求 c k c_k ck 就不用我说了吧。 QwQ
c k = e k n c_k=\dfrac{e_k}{n} ck=nek
此时我们运用上一部分中的系数表示->点值表示的方法即可在 O(log 2 n) \text{O(log}_2\text{n)} O(log2n) 的时间复杂度下求每个 e k e_k ek 也就是 D ( ω n − k ) D(\omega_n^{-k}) D(ωn−k) 的值,每求出一个 e k e_k ek ,我们就可以运用 c k = e k n c_k=\dfrac{e_k}{n} ck=nek 这一关系 O(1) \text{O(1)} O(1) 求出 C ( x ) C(x) C(x) 的第 k k k 项系数。
而 C ( x ) C(x) C(x) 的项数与 n n n 成正比,所以我们运用这一方法实现点值表示->系数表示的渐进时间复杂度即为 O(nlog 2 n) \text{O(nlog}_2\text{n)} O(nlog2n)。
Part 3.3 \text{Part 3.3} Part 3.3 回顾 FFT \text{FFT} FFT 的全过程
FFT \text{FFT} FFT 是一种可以在 O(n log 2 n) \text{O(n log}_2\text{n)} O(n log2n) 的渐进时间复杂度内求次数都为 n n n 的两个多项式 A ( x ) A(x) A(x) 和 B ( x ) B(x) B(x) 的卷积 C ( x ) C(x) C(x) 的值的算法。
-
Part 3.1 \text{Part 3.1} Part 3.1
- 设 C ( x ) = A ( x ) ∗ B ( x ) C(x)=A(x)*B(x) C(x)=A(x)∗B(x),分治求 A ( x ) , B ( x ) A(x),B(x) A(x),B(x) 的离散傅里叶变换,然后求对应的 C ( x ) C(x) C(x) 的离散傅里叶变换,渐进时间复杂度为 O(nlog 2 n) \text{O(nlog}_2\text{n)} O(nlog2n)。
-
Part 3.2 \text{Part 3.2} Part 3.2
- 以 C ( x ) C(x) C(x) 的离散傅里叶变换作为另一个与 C ( x ) C(x) C(x) 同次的多项式 D ( x ) D(x) D(x) 的系数,每次用 Part 3.1 \text{Part 3.1} Part 3.1 中的办法 O(log 2 n) \text{O(log}_2\text{n)} O(log2n) 求每个 D ( x ) D(x) D(x) 在 ( ω 2 n − 1 0 , ω 2 n − 1 − 1 , ⋯ , ω 2 n − 1 − ( 2 n − 2 ) ) (\omega_{2n-1}^{0},\omega_{2n-1}^{-1},\cdots,\omega_{2n-1}^{-(2n-2)}) (ω2n−10,ω2n−1−1,⋯,ω2n−1−(2n−2)) 处取值的得到点值 e i e_i ei,根据 c i = e i n c_i=\dfrac{e_i}{n} ci=nei 的关系求 C ( x ) C(x) C(x) 的每一项系数,求 C ( x ) C(x) C(x) 的 2 n − 2 2n-2 2n−2 项系数的渐进时间复杂度为 O(n log 2 n) \text{O(n log}_2\text{n)} O(n log2n) 。
- 求出 C ( x ) C(x) C(x) 的每一项系数之后,因为 C ( x ) = A ( x ) ∗ B ( x ) C(x)=A(x)*B(x) C(x)=A(x)∗B(x) ,所以对于多项式 A ( x ) A(x) A(x) 和 B ( x ) B(x) B(x) 的乘积,我们直接将 x x x 代入系数表示下 C ( x ) C(x) C(x) 在 O(n) \text{O(n)} O(n) 渐进时间复杂度下进行求解即可。
所以 FFT \text{FFT} FFT 的渐进时间复杂度为 O(nlog 2 n) \text{O(nlog}_2\text{n)} O(nlog2n) 。
放一下板题(P3803 【模板】多项式乘法(FFT))递归写法的代码:
#include<bits/stdc++.h>
#define N 4000005
#define P acos(-1.0)
using namespace std;
int aa,bb;
struct cp{
double x,y;
cp (double xx=0,double yy=0){ x=xx,y=yy; }
}a[N],b[N];
cp operator +(cp a,cp b){ return cp(a.x+b.x,a.y+b.y);}
cp operator -(cp a,cp b){ return cp(a.x-b.x,a.y-b.y);}
cp operator *(cp a,cp b){ return cp(a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y);}
void fft(cp*,int,int);
int main(){
scanf("%d%d",&aa,&bb);
for(int i=0;i<=aa;i++)
scanf("%lf",&a[i].x);
for(int i=0;i<=bb;i++)
scanf("%lf",&b[i].x);
int n=1;
while(n<=aa+bb) n<<=1;
fft(a,n,1),fft(b,n,1);
for(int i=0;i<n;i++)
a[i]=a[i]*b[i];
fft(a,n,-1);
for(int i=0;i<=aa+bb;i++)
printf("%d ",(int)(a[i].x/n+0.5));
return 0;
}
void fft(cp *a,int n,int inv){
if(n==1) return ;
cp a1[(n>>1)+5],a2[(n>>1)+5];
for(int i=0;i<n;i+=2)//这里是i+=2!!!
a1[i>>1]=a[i],a2[i>>1]=a[i+1];
fft(a1,n>>1,inv),fft(a2,n>>1,inv);
for(int i=0;i<(n>>1);i++){
cp x(cos(2.0*P*i/n),inv*sin(2.0*P*i/n)),ax=x*a2[i];
a[i]=a1[i]+ax,a[i+(n>>1)]=a1[i]-ax;
}
}
Part 3.4 \text{Part 3.4} Part 3.4 FFT \text{FFT} FFT 优化
优化一:递归->迭代
迭代即用循环模拟递归的过程。
我们观察一下在做FFT的时候多项式 C ( x ) C(x) C(x) 的每一项系数的位置变化。
0 1 2 3 4 5 6 7
0 2 4 6|1 3 5 7
0 4|2 6|1 5|3 7
0|4|2|6|1|5|3|7
观察出来没有,(当然是没有啊!) 对于每一个系数,它的最终下标为它原来下标的二进制翻转。
如原来处在位置 6 ( 110 ) 6(110) 6(110) 的系数最终到了位置 3 ( 011 ) 3(011) 3(011) 。
那么此时我们其实就可以把要转换的多项式 C ( x ) C(x) C(x) 的系数一开始就调换到它的最终位置,然后再通过循环模拟递归的过程。
代码也很简单:
int sum_len = len_a+len_b, len = 0;
while (n <= sum_len)
++len, n <<= 1;
for (int i = 1; i < n; i++)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (len - 1));
for (int i = 0; i < n; i++)
if (i < r[i]) swap(a[i], a[r[i]]);
for (int len = 2; len <= n; len <<= 1) {
cp wn(cos(P * 2 / len), inv * sin(P * 2 / len));
for (int st = 0; st + len <= n; st += len) {
cp ax(1, 0);
for (int j = st; j < st + (len >> 1); j++, ax = ax * wn) {
cp a1 = a[j], a2 = a[j + (len >> 1)] * ax;
a[j] = a1 + a2, a[j + (len >> 1)] = a1 - a2;
}
}
}
优化二:减少复数乘法运算次数
而且实现FFT时还有一点十分重要,就是转化过程中的复数乘法运算。
有图有真相:D
让我们来一起看看两份代码有什么区别吧:D
5.42 s 5.42s 5.42s:
for(int len=2;len<=n;len<<=1){
for(int st=0;st+len-1<n;st+=len){
for(int j=st;j<=st+(len>>1)-1;j++){
cp x(cos(2.0*P*(j-st)/len),inv*sin(2.0*P*(j-st)/len)),ax=x*a[j+(len>>1)],a1=a[j];
a[j]=a1+ax,a[j+(len>>1)]=a1-ax;
}
}
}
2.83 s 2.83s 2.83s:
for(int len=2;len<=n;len<<=1){
complex wn(cos(2*P/len),inv*sin(2*P/len));
for(int st=0;st+len<=n;st+=len){
complex a1,a2,ax(1,0);
for(int j=st;j<st+(len>>1);j++,ax=ax*wn)
a1=a[j],a2=a[j+(len>>1)]*ax,
a[j]=a1+a2,a[j+(len>>1)]=a1-a2;
}
}
两份代码的其余部分均相同。
两份代码仅仅差了几次的复数运算,最终的运行结果却相差甚远。
所以说,复数乘法这种东西,咱们还是能少算几次就少算几次。
在第二份代码中我们就少算了很多次 (cos(2*P/len),inv*sin(2*P/len))
的值,故常数优化了不少。
因为这里是我感性理解的,并没有找理性证明的博客进行学习QwQ,所以如果我说错了的话请在评论区指出我的错误。
优化三:三次FFT变为两次FFT
在我们求 A ( x ) A(x) A(x) 和 B ( x ) B(x) B(x) 的卷积 C ( x ) C(x) C(x) 时,我们会对它们分别做 FFT \text{FFT} FFT 转化为点值表示,然后对应每项乘起来得到 C ( x ) C(x) C(x) 的点值表示,最后再做一次 FFT \text{FFT} FFT 将 C ( x ) C(x) C(x) 由点值表示转化为系数表示。
这样的话实际上我们一共进行了三次 FFT \text{FFT} FFT ,于是我们考虑能否减少做 FFT \text{FFT} FFT 的次数。
当然可以啊,不然我怎么可能写这个板块啊。
我们弄两个系数为复数的多项式,按正常 FFT \text{FFT} FFT 的做法,我们应该分别把多项式 A ( x ) A(x) A(x) 和 B ( x ) B(x) B(x) 的系数放到两个多项式系数的实部上面去。
但是在三次变两次 FFT \text{FFT} FFT 优化中,我们只需要弄一个系数为复数的多项式,并将多项式 A ( x ) A(x) A(x) 和 B ( x ) B(x) B(x) 的系数分别放到该多项式系数的实部和虚部中,然后对这个多项式跑一边 FFT \text{FFT} FFT ,然后对它做平方,再做一遍 FFT \text{FFT} FFT 逆变换回来,最后取它每一位系数虚部上的数字再 / 2 n /2n /2n 即可(该多项式次数为 n n n)。
证明也很简单,知道完全平方公式&复数概念即可证明。
( a + b i ) 2 = a 2 + ( b 2 i 2 ) + 2 a b i = a 2 − b 2 + 2 a b i (a+bi)^2=a^2+(b^2i^2)+2abi=a^2-b^2+2abi (a+bi)2=a2+(b2i2)+2abi=a2−b2+2abi
此时我们的常数就可以优化到约为原来的 2 3 \frac{2}{3} 32 了。
这是三个优化都加上了的代码跑出来的:D
最后放个加了优化的代码QwQ
#include<bits/stdc++.h>
#define P acos(-1.0)
#define N 1000005
using namespace std;
int a1,b1,n=1,r[N<<2],ans[N<<2],aa1;
struct cp{
double x,y;
cp(double xx=0,double yy=0){ x=xx,y=yy;}
}a[N<<2],b[N<<2];//这里的复数我用的是结构体表示,实际上复数运算也可以用STL
cp operator +(cp a,cp b){ return cp(a.x+b.x,a.y+b.y);}
cp operator -(cp a,cp b){ return cp(a.x-b.x,a.y-b.y);}
cp operator *(cp a,cp b){ return cp(a.x*b.x-b.y*a.y,a.y*b.x+b.y*a.x);}
void pre(),fft(cp*,int);
int main(){
scanf("%d%d",&a1,&b1);
for(int i=0;i<=a1;i++) scanf("%lf",&a[i].x);
for(int i=0;i<=b1;i++) scanf("%lf",&a[i].y);
pre(),fft(a,1);
for(int i=0;i<n;i++) a[i]=a[i]*a[i];
fft(a,-1);
for(int i=0;i<=a1+b1;i++) printf("%d ",(int)(a[i].y/n/2+0.5));//取虚部求最终的系数表示
return 0;
}
void pre(){
aa1=max(a1,b1)<<1;
while(n<=aa1) n<<=1;
int len=1;
while((1<<len)<=aa1) ++len;
for(int i=1;i<n;i++) r[i]=((r[i>>1]>>1|(i&1)<<(len-1)));//预处理每一个系数的最终位置(二进制翻转)
}
void fft(cp *a,int inv){
for(int i=0;i<n;i++)
if(i<r[i]) swap(a[r[i]],a[i]);
//用循环模拟递归的过程
for(int len=2;len<=n;len<<=1){
cp wn(cos(P*2/len),inv*sin(P*2/len));//提前算好单位复数根,减少重复运算
for(int st=0;st+len-1<n;st+=len){
cp ax(1,0);
for(int j=st;j<st+(len>>1);j++,ax=ax*wn){
cp a1=a[j],a2=ax*a[j+(len>>1)];
a[j]=a1+a2,a[j+(len>>1)]=a1-a2;
}
}
}
}