1、FFT前言
FFT全称(Fast Fourier Transformation)快速傅里叶变换
NTT全称(Fast Number Theoretical Transformation)快速数论变换
1.1、引入问题
如何计算多项式 f ( x ) f(x) f(x)与 g ( x ) g(x) g(x)的乘积?
例如:计算 f ( x ) = 1 + x 2 + x 3 f(x)=1+x^2+x^3 f(x)=1+x2+x3与 g ( x ) = 1 + x 5 + x 6 + x 7 g(x)=1+x^5+x^6+x^7 g(x)=1+x5+x6+x7的乘积
通常而言,我们思考暴力的算法:
f
(
x
)
×
g
(
x
)
=
1
×
(
1
+
x
5
+
x
6
+
x
7
)
+
x
2
×
(
1
+
x
5
+
x
6
+
x
7
)
+
x
3
×
(
1
+
x
5
+
x
6
+
x
7
)
=
1
+
x
5
+
x
6
+
x
7
+
x
2
+
x
7
+
x
8
+
x
9
+
x
3
+
x
8
+
x
9
+
x
10
=
1
+
x
2
+
x
3
+
x
5
+
x
6
+
2
x
7
+
2
x
8
+
2
x
9
+
x
10
\begin{aligned} f(x)\times g(x)&=1\times(1+x^5+x^6+x^7)+x^2\times(1+x^5+x^6+x^7)+x^3\times(1+x^5+x^6+x^7) \\& =1+x^5+x^6+x^7+x^2+x^7+x^8+x^9+x^3+x^8+x^9+x^{10} \\& =1+x^2+x^3+x^5+x^6+2x^7+2x^8+2x^9+x^{10} \end{aligned}
f(x)×g(x)=1×(1+x5+x6+x7)+x2×(1+x5+x6+x7)+x3×(1+x5+x6+x7)=1+x5+x6+x7+x2+x7+x8+x9+x3+x8+x9+x10=1+x2+x3+x5+x6+2x7+2x8+2x9+x10
显而易见,这个算法的时间复杂度是
O
(
n
2
)
O(n^2)
O(n2)
那么有快速的计算方法吗?
当然是有的,于1965年由J.W.库利和T.W.图基提出FFT算法,可以在 O ( n l o g n ) O(nlogn) O(nlogn)的时间复杂度下计算出答案。
1.2、定义一些记号:
# f ( x ) f(x) f(x)表示一个多项式
**#**对于一个 n n n次多项式 f ( x ) f(x) f(x)的第 m m m项,记为 f [ m ] f[m] f[m];换言之, f ( x ) = ∑ i = 0 n f [ i ] x i f(x)=\sum_{i=0}^{n}f[i]x^i f(x)=∑i=0nf[i]xi
1.3、多项式乘法理解
例如:现有两个 n n n次多项式 f ( x ) f(x) f(x)与 g ( x ) g(x) g(x),需要计算这两个多项式的乘积.
不妨设 h ( x ) = f ( x ) × g ( x ) h(x)=f(x)\times g(x) h(x)=f(x)×g(x)
那么通过1.1中的例子,我们可以发现,其实 h [ k ] = ∑ i = 0 k f [ i ] g [ k − i ] h[k]=\sum_{i=0}^{k}f[i]g[k-i] h[k]=∑i=0kf[i]g[k−i]
换言之,我们可以理解为:要想乘出 x k x^k xk,需要 f ( x ) f(x) f(x)的 x i x^i xi项和 g ( x ) g(x) g(x)的 x k − i x^{k-i} xk−i项相乘,然后将他们的相乘的贡献累加便是 h [ k ] h[k] h[k]
1.4、卷积
将形如: h [ k ] = ∑ i ⊕ j f [ i ] g [ j ] h[k]=\sum_{i⊕j}f[i]g[j] h[k]=∑i⊕jf[i]g[j]的式子称为卷积,其中⊕在数学中称为假设符号(假设什么运算)
例如,在多项式乘法中,为加法卷积: h [ k ] = ∑ i + j f [ i ] g [ j ] h[k]=\sum_{i+j}f[i]g[j] h[k]=∑i+jf[i]g[j]
因此,对于卷积可以利用FFT构造多项式进行计算!
2、FFT思想
2.1、多项式的点值表达式
我们可以从
得出推论:在平面直角坐标系中, n + 1 n+1 n+1个点就能唯一确定一个 n n n次多项式。
因此,对于一个 n n n次多项式 f ( x ) f(x) f(x)而言,我们可以用 n + 1 n+1 n+1个点来描述一个多项式;换言之,两个多项式相乘,等价于两个多项式的对应点值相乘。(点值的乘法对应着整个多项式的乘法,也就是浓缩了整个多项式)
将能唯一表达一个 n n n次多项式 f ( x ) f(x) f(x)的 n + 1 n+1 n+1个点,称为 f ( x ) f(x) f(x)的点值表达式。
$$
例如:n次多项式f(x)的点值表达式为:\
\left{(x_1,f(x_1)),(x_2,f(x_2)),(x_3,f(x_3))…(x_{n+1},f(x_{n+1}))\right}
$$
由于
n
n
n次多项式
f
(
x
)
f(x)
f(x)与
n
n
n次多项式
g
(
x
)
g(x)
g(x),他们的乘积一定是一个
2
n
+
1
2n+1
2n+1次的多项式,所以我们只需要对
f
(
x
)
f(x)
f(x)和
g
(
x
)
g(x)
g(x)找到
2
n
+
1
2n+1
2n+1个对应的点,然后将他们对应点的纵坐标相乘,便可以确定多项式
f
(
x
)
×
g
(
x
)
f(x)\times g(x)
f(x)×g(x)的点值表达式;
例如:
2
n
+
1
次多项式
f
(
x
)
×
g
(
x
)
的点值表达式为:
{
(
x
1
,
f
(
x
1
)
g
(
x
1
)
)
,
(
x
2
,
f
(
x
2
)
g
(
x
2
)
)
,
(
x
3
,
f
(
x
3
)
g
(
x
3
)
)
.
.
.
.
.
.
(
x
2
n
+
1
,
f
(
x
2
n
+
1
)
g
(
x
2
n
+
1
)
)
}
例如:2n+1次多项式f(x)\times g(x)的点值表达式为: \\ \left\{(x_1,f(x_1)g(x_1)),(x_2,f(x_2)g(x_2)),(x_3,f(x_3)g(x_3))......(x_{2n+1},f(x_{2n+1})g(x_{2n+1}))\right\}
例如:2n+1次多项式f(x)×g(x)的点值表达式为:{(x1,f(x1)g(x1)),(x2,f(x2)g(x2)),(x3,f(x3)g(x3))......(x2n+1,f(x2n+1)g(x2n+1))}
最后我们将多项式
f
(
x
)
×
g
(
x
)
f(x)\times g(x)
f(x)×g(x)的点值表达式转化为多项式的一般表达式(系数表达式)便得到了我们想要的多项式乘积了。
而这也是FFT算法的流程
2.2、DFT&IDFT
DFT:把多项式的系数表达式转化为点值表达式的变换过程
IDFT:把多项式的点值表达式转化为系数表达式的变换过程
下面这张图便对比了FFT算法与暴力算法的区别
3、FFT的实现
3.1、DFT
由于我们想要多项式的点值表达式,所以我们需要利用一些方法获得多项式的上的点;但是细想一下会发现,
如果对于一个
n
次多项式
f
(
x
)
,我们任选
n
+
1
个点的横坐标
x
1
,
x
2
,
x
3
.
.
.
.
x
n
,
x
n
+
1
,
再分别计算其对应的纵坐标
f
(
x
1
)
,
f
(
x
2
)
,
f
(
x
3
)
.
.
.
f
(
x
n
)
,
f
(
x
n
+
1
)
的时间复杂度是
O
(
n
2
)
(
由于计算每个横坐标的时间复杂度是
O
(
n
)
,
但是有
n
个点
)
;
如果对于一个n次多项式f(x),我们任选n+1个点的横坐标x_1,x_2,x_3....x_n,x_{n+1},\\再分别计算其对应的纵坐标f(x_1),f(x_2),f(x_3)...f(x_n),f(x_{n+1})的时间复杂度是O(n^2)\\(由于计算每个横坐标的时间复杂度是O(n),但是有n个点);
如果对于一个n次多项式f(x),我们任选n+1个点的横坐标x1,x2,x3....xn,xn+1,再分别计算其对应的纵坐标f(x1),f(x2),f(x3)...f(xn),f(xn+1)的时间复杂度是O(n2)(由于计算每个横坐标的时间复杂度是O(n),但是有n个点);
这样一来,我们通过随便找
n
+
1
n+1
n+1个点的方式,会让我们的时间复杂度和暴力没什么区别!
那么有快速找到 n + 1 n+1 n+1个点的方法么?答案当然是有的!
我们可以这样考虑问题:
3.1.1、Good idea 1:函数的奇偶性
首先,我们不考虑一般多项式,只考虑一个简单的多项式 f ( x ) = x 2 f(x)=x^2 f(x)=x2
那么,如果我们要从 f ( x ) f(x) f(x)中选择8个点,我们如何选择?
这里我们需要注意 f ( x ) f(x) f(x)的一个优美的性质———— f ( x ) f(x) f(x)是偶函数,所以我们计算出一个点 ( x , f ( x ) ) (x,f(x)) (x,f(x))后,便可以立马知道点 ( − x , f ( x ) ) (-x,f(x)) (−x,f(x))也一定是 y ( x ) y(x) y(x)上的点。
那么,如果 f ( x ) = x 3 f(x)=x^3 f(x)=x3,应该如何处理喃?
显而易见,我们会发现 f ( x ) = x 3 f(x)=x^3 f(x)=x3是一个奇函数,它与偶函数的区别在于,奇函数上的点是关于原点对称的,换言之,如果我们计算出点 ( x , f ( x ) ) (x,f(x)) (x,f(x)),那就可以知道点 ( − x , − f ( x ) ) (-x,-f(x)) (−x,−f(x))也是多项式上的点!
既然这样,对于多项式 f ( x ) = 1 + 2 x 2 + 3 x 3 + 4 x 4 + 5 x 5 f(x)=1+2x^2+3x^3+4x^4+5x^5 f(x)=1+2x2+3x3+4x4+5x5
那么我可以将
f
(
x
)
f(x)
f(x)中按幂数分离,因此可得
f
(
x
)
=
(
1
+
2
x
2
+
4
x
4
)
+
(
3
x
3
+
5
x
5
)
f(x)=(1+2x^2+4x^4)+(3x^3+5x^5)
f(x)=(1+2x2+4x4)+(3x3+5x5)
不妨设
f
o
(
x
)
=
1
+
2
x
+
4
x
2
,
f
e
(
x
)
=
3
x
+
5
x
2
那么有:
f
(
x
)
=
(
1
+
2
x
2
+
4
x
4
)
+
x
(
3
x
2
+
5
x
4
)
=
f
o
(
x
2
)
+
x
f
e
(
x
2
)
不妨设f_o(x)=1+2x+4x^2,f_e(x)=3x+5x^2\\ 那么有: \begin{aligned} f(x)&=(1+2x^2+4x^4)+x(3x^2+5x^4) \\&=f_o(x^2)+xf_e(x^2) \end{aligned}
不妨设fo(x)=1+2x+4x2,fe(x)=3x+5x2那么有:f(x)=(1+2x2+4x4)+x(3x2+5x4)=fo(x2)+xfe(x2)
这样,对于
f
o
(
x
2
)
,
f
e
(
x
2
)
f_o(x^2),f_e(x^2)
fo(x2),fe(x2)都是关于
x
2
x^2
x2的多项式,我们便可以利用函数的奇偶性快速找到6个点,然后利用点值表达式的计算方式计算出该函数的点值表达式。
然后我们可以推广到一般多项式
f
(
x
)
=
a
0
+
a
1
x
+
a
2
x
2
+
.
.
.
.
.
a
n
x
n
f(x)=a_0+a_1x+a_2x^2+.....a_nx^{n}
f(x)=a0+a1x+a2x2+.....anxn
可以设
f
(
x
)
按照幂数分离,将偶幂数部分设为
f
e
(
x
2
)
,
将奇幂数部分提出一个
x
后的部分设为
f
o
(
x
2
)
所以:
f
(
x
)
=
f
e
(
x
2
)
+
x
f
o
(
x
2
)
然而对于,
f
e
(
x
)
与
f
o
(
x
)
而言,最高次都变成了
n
2
−
1
同理可以推出重要结论
{
f
(
x
i
)
=
f
e
(
x
i
2
)
+
x
i
f
o
(
x
i
2
)
f
(
−
x
i
)
=
f
e
(
x
i
2
)
−
x
i
f
o
(
x
i
2
)
这样一来,问题得到了简化:
分别求
f
e
(
x
2
)
,
f
o
(
x
2
)
在
x
1
2
,
x
2
2
.
.
.
.
x
n
2
2
的值
,
因为这样我们便可以得到
±
x
1
,
±
x
2
.
.
.
±
x
n
2
这
n
个点
可以设f(x)按照幂数分离,将偶幂数部分设为f_e(x^2),将奇幂数部分提出一个x后的部分设为f_o(x^2)\\ 所以:f(x)=f_e(x^2)+xf_o(x^2)\\ 然而对于,f_e(x)与f_o(x)而言,最高次都变成了\frac{n}{2}-1\\ 同理可以推出重要结论 \begin{cases} f(x_i)=f_e(x_i^2)+x_if_o(x_i^2)\\ f(-x_i)=f_e(x_i^2)-x_if_o(x_i^2) \end{cases} \\ 这样一来,问题得到了简化: \\分别求f_e(x^2),f_o(x^2)在x_1^2,x_2^2....x_{\frac{n}{2}}^2的值,因为这样我们便可以得到\pm x_1,\pm x_2...\pm x_\frac{n}{2}这n个点
可以设f(x)按照幂数分离,将偶幂数部分设为fe(x2),将奇幂数部分提出一个x后的部分设为fo(x2)所以:f(x)=fe(x2)+xfo(x2)然而对于,fe(x)与fo(x)而言,最高次都变成了2n−1同理可以推出重要结论{f(xi)=fe(xi2)+xifo(xi2)f(−xi)=fe(xi2)−xifo(xi2)这样一来,问题得到了简化:分别求fe(x2),fo(x2)在x12,x22....x2n2的值,因为这样我们便可以得到±x1,±x2...±x2n这n个点
而对于
f
e
(
x
)
与
f
o
(
x
)
f_e(x)与f_o(x)
fe(x)与fo(x),我们可以重复
f
(
x
)
f(x)
f(x)的奇偶分离的方式;这样我们可以把对于
n
n
n次多项式转化为求
n
2
\frac{n}{2}
2n次多项式(类似分治的思想),所以这个思路下,求多项式点值表达式的时间复杂度是
O
(
n
l
o
g
n
)
O(nlogn)
O(nlogn)
因此,只要我们最后找到的点为合法的点,那么算法就是合法的!
但是这里出现了一个问题:对于 f ( x ) f(x) f(x)中的 [ ± x 1 , ± x 2 . . . ± x n 2 ] [\pm x_1,\pm x_2...\pm x\frac{n}{2}] [±x1,±x2...±x2n]而言,每一个 ± x i \pm x_i ±xi都是 ± \pm ±配对的;但是对于 f o ( x 2 ) 与 f e ( x 2 ) f_o(x^2)与f_e(x^2) fo(x2)与fe(x2)中的 [ x 1 2 , x 2 2 . . x n 2 ] [x_1^2,x_2^2..x_n^2] [x12,x22..xn2]而言,没有一个 x i 2 x_i^2 xi2是配对的;所以在实数域下,是找不到满足条件的点的!
不过解决问题的的方法总比问题多,所以我们考虑把数域扩展到复数域!
3.1.2、Good idea 2:复数域下的单位根
3.1.2.1解决疑惑:为何复数域下,满足3.1.1中的分治找点
一个简单例子:求多项式
f
(
x
)
=
1
+
2
x
+
3
x
3
f(x)=1+2x+3x^3
f(x)=1+2x+3x3的点值表达式
在这里我们不妨设这四个点的横坐标为
:
x
1
,
x
2
,
x
3
,
x
4
那么根据分治的规则
:
那么必须保证
x
1
=
−
x
2
,
x
3
=
−
x
4
∴
x
1
2
=
x
2
…
…
2
,
x
3
2
=
x
4
2
∵
x
1
2
=
−
x
3
2
,
那么可令
x
1
=
1
∴
x
3
2
=
−
1
∵
在复数域下有
i
2
=
−
1
∴
x
3
=
i
,
x
4
=
−
i
∴
对于横坐标
x
1
=
1
,
x
2
=
−
1
,
x
3
=
i
,
x
4
=
−
i
是满足条件的
在这里我们不妨设这四个点的横坐标为:x_1,x_2,x_3,x_4\\ 那么根据分治的规则:那么必须保证x_1=-x_2,x_3=-x_4\\ \therefore x_1^2=x_2……2,x_3^2=x_4^2\\ \because x_1^2=-x_3^2,那么可令x_1=1\\ \therefore x_3^2=-1\\ \because 在复数域下有i^2=-1\\ \therefore x_3=i,x_4=-i\\ \therefore 对于横坐标x_1=1,x_2=-1,x_3=i,x_4=-i是满足条件的
在这里我们不妨设这四个点的横坐标为:x1,x2,x3,x4那么根据分治的规则:那么必须保证x1=−x2,x3=−x4∴x12=x2……2,x32=x42∵x12=−x32,那么可令x1=1∴x32=−1∵在复数域下有i2=−1∴x3=i,x4=−i∴对于横坐标x1=1,x2=−1,x3=i,x4=−i是满足条件的
3.1.2.2引入单位根
单位根:在复平面上,以单位圆点为起点,单位圆的 n n n 等分点为终点,在单位圆上可以得到 n n n 个复数,设幅角为正且最小的复数为 n n n次单位根,记作 ω n \omega_n ωn
**#**推导
ω
n
\omega_n
ωn
由于在复平面当中,故俯角最小为
2
π
n
,
所以可知在复平面中,
ω
n
=
c
o
s
2
π
n
+
i
s
i
n
2
π
n
又
∵
欧拉公式为
e
i
θ
=
c
o
s
θ
+
i
s
i
n
θ
∴
ω
n
=
e
i
2
π
n
由于在复平面当中,故俯角最小为\frac{2\pi}{n},所以可知在复平面中,\omega_n=cos\frac{2\pi}{n}+isin\frac{2\pi}{n}\\ 又\because 欧拉公式为 e^{i\theta}=cos\theta+isin\theta\\ \therefore\omega_n=e^{i\frac{2\pi}{n}}
由于在复平面当中,故俯角最小为n2π,所以可知在复平面中,ωn=cosn2π+isinn2π又∵欧拉公式为eiθ=cosθ+isinθ∴ωn=ein2π
**#**推导单位根的一些性质:
(1)
w
n
a
×
w
n
b
=
w
n
i
+
j
w_n^{a}\times w_n^{b}=w_n^{i+j}
wna×wnb=wni+j
w
n
a
×
w
n
b
=
e
a
i
2
π
n
×
e
b
i
2
π
n
=
e
(
a
+
b
)
i
2
π
n
=
(
e
i
2
π
n
)
a
+
b
=
w
n
a
+
b
\begin{aligned} w_n^{a}\times w_n^{b}&=e^{ai\frac{2\pi}{n}}\times e^{bi\frac{2\pi}{n}}\\ &=e^{(a+b)i\frac{2\pi}{n}}\\ &=(e^{i\frac{2\pi}{n}})^{a+b}\\ &=w_n^{a+b} \end{aligned}
wna×wnb=eain2π×ebin2π=e(a+b)in2π=(ein2π)a+b=wna+b
(2)
w
r
n
r
k
=
w
n
k
w_{rn}^{rk}=w_n^{k}
wrnrk=wnk
w
r
n
r
k
=
(
e
i
2
π
r
n
)
r
k
=
(
e
i
2
π
r
n
r
)
k
=
(
e
i
2
π
n
)
k
=
w
n
k
w_{rn}^{rk}=(e^{i\frac{2\pi}{rn}})^{rk}=(e^{i\frac{2\pi}{rn}r})^k=(e^{i\frac{2\pi}{n}})^k=w_n^{k}
wrnrk=(eirn2π)rk=(eirn2πr)k=(ein2π)k=wnk
(3)
w
n
2
k
=
w
n
2
k
w_n^{2k}=w_{\frac{n}{2}}^{k}
wn2k=w2nk
w
n
2
k
=
e
2
k
i
2
π
n
=
e
k
i
2
π
n
2
=
w
n
2
k
w_n^{2k}=e^{2ki\frac{2\pi}{n}}=e^{ki\frac{2\pi}{\frac{n}{2}}}=w_{\frac{n}{2}}^{k}
wn2k=e2kin2π=eki2n2π=w2nk
(4)如果
n
n
n是偶数,
w
n
k
+
n
2
=
−
w
n
k
w_n^{k+\frac{n}{2}}=-w_n^{k}
wnk+2n=−wnk
w
n
k
+
n
2
=
(
e
i
2
π
n
)
k
+
n
2
=
e
k
i
2
π
n
+
i
π
=
e
i
π
e
k
i
2
π
n
=
−
1
×
e
k
i
2
π
n
=
−
w
n
k
\begin{aligned} w_n^{k+\frac{n}{2}}&=(e^{i\frac{2\pi}{n}})^{k+\frac{n}{2}}\\ &=e^{ki\frac{2\pi}{n}+i\pi}\\ &=e^{i\pi}e^{ki\frac{2\pi}{n}}\\ &=-1\times e^{ki\frac{2\pi}{n}}\\ &=-w_n^{k} \end{aligned}
wnk+2n=(ein2π)k+2n=ekin2π+iπ=eiπekin2π=−1×ekin2π=−wnk
3.1.2.3代入单位根
现已知 f ( x ) = f e ( x 2 ) + x f o ( x 2 ) f(x)=f_e(x^2)+xf_o(x^2) f(x)=fe(x2)+xfo(x2)
(1)代入
x
=
w
n
k
,
k
<
n
2
x=w_n^{k},k<\frac{n}{2}
x=wnk,k<2n
f
(
w
n
k
)
=
f
e
(
w
n
2
k
)
+
w
n
k
f
o
(
w
n
2
k
)
=
f
e
(
w
n
2
k
)
+
w
n
k
f
o
(
w
n
2
k
)
\begin{aligned} f(w_n^{k})&=f_e(w_n^{2k})+w_n^{k}f_o(w_n^{2k})\\ &=f_e(w_{\frac{n}{2}}^{k})+w_n^{k}f_o(w_{\frac{n}{2}}^{k}) \end{aligned}
f(wnk)=fe(wn2k)+wnkfo(wn2k)=fe(w2nk)+wnkfo(w2nk)
(2)代入
x
=
w
n
k
+
n
2
,
k
<
n
2
x=w_n^{k+\frac{n}{2}},k<\frac{n}{2}
x=wnk+2n,k<2n
f
(
w
n
k
+
n
2
)
=
f
e
(
w
n
2
(
k
+
n
2
)
)
+
w
n
k
+
n
2
f
o
(
w
n
2
(
k
+
n
2
)
)
=
f
e
(
w
n
2
k
+
n
2
)
+
w
n
k
+
n
2
f
o
(
w
n
2
k
+
n
2
)
=
f
e
(
w
n
2
k
+
n
2
)
−
w
n
k
f
o
(
w
n
2
k
+
n
2
)
=
f
e
(
w
n
2
k
+
n
2
)
−
w
n
k
f
o
(
w
n
2
k
+
n
2
)
=
f
e
(
w
n
2
k
+
n
)
−
w
n
k
f
o
(
w
n
2
k
+
n
)
=
f
e
(
w
n
2
k
)
−
w
n
k
f
o
(
w
n
2
k
)
=
f
e
(
w
n
2
k
)
−
w
n
k
f
o
(
w
n
2
k
)
\begin{aligned} f(w_n^{k+\frac{n}{2}})&=f_e(w_n^{2(k+\frac{n}{2})})+w_n^{k+\frac{n}{2}}f_o(w_n^{2(k+\frac{n}{2})})\\ &=f_e(w_{\frac{n}{2}}^{k+\frac{n}{2}})+w_n^{k+\frac{n}{2}}f_o(w_{\frac{n}{2}}^{k+\frac{n}{2}})\\ &=f_e(w_{\frac{n}{2}}^{k+\frac{n}{2}})-w_n^kf_o(w_{\frac{n}{2}}^{k+\frac{n}{2}})\\ &=f_e(w_{\frac{n}{2}}^{k+\frac{n}{2}})-w_n^kf_o(w_{\frac{n}{2}}^{k+\frac{n}{2}})\\ &=f_e(w_n^{2k+n})-w_n^kf_o(w_n^{2k+n})\\ &=f_e(w_n^{2k})-w_n^kf_o(w_n^{2k})\\ &=f_e(w_{\frac{n}{2}}^{k})-w_n^kf_o(w_{\frac{n}{2}}^{k}) \end{aligned}
f(wnk+2n)=fe(wn2(k+2n))+wnk+2nfo(wn2(k+2n))=fe(w2nk+2n)+wnk+2nfo(w2nk+2n)=fe(w2nk+2n)−wnkfo(w2nk+2n)=fe(w2nk+2n)−wnkfo(w2nk+2n)=fe(wn2k+n)−wnkfo(wn2k+n)=fe(wn2k)−wnkfo(wn2k)=fe(w2nk)−wnkfo(w2nk)
那么得到( 1 )( 2 )两个结论,我们可以: 如果知道多项式 f e ( x ) , f o ( x ) 在 w n 2 0 , w n 2 1 , w n 2 2 . . . . w n 2 n 2 − 1 的点值表达式, 利用( 1 )可以 O ( n ) 求出 f ( x ) 在 w n 2 0 , w n 2 1 , w n 2 2 . . . . w n 2 n 2 − 1 的点值表达式 ; 利用( 2 )可以 O ( n ) 求出 f ( x ) 在 w n 2 n 2 , w n 2 n 2 + 1 , w n 2 n 2 + 2 . . . . w n 2 n − 1 的点值表达式 ; 那么得到(1)(2)两个结论,我们可以:\\ 如果知道多项式f_e(x),f_o(x)在w_{\frac{n}{2}}^0,w_{\frac{n}{2}}^1,w_{\frac{n}{2}}^2....w_{\frac{n}{2}}^{\frac{n}{2}-1}的点值表达式,\\ 利用(1)可以O(n)求出f(x)在w_{\frac{n}{2}}^0,w_{\frac{n}{2}}^1,w_{\frac{n}{2}}^2....w_{\frac{n}{2}}^{\frac{n}{2}-1}的点值表达式;\\ 利用(2)可以O(n)求出f(x)在w_{\frac{n}{2}}^{\frac{n}{2}},w_{\frac{n}{2}}^{\frac{n}{2}+1},w_{\frac{n}{2}}^{\frac{n}{2}+2}....w_{\frac{n}{2}}^{n-1}的点值表达式; 那么得到(1)(2)两个结论,我们可以:如果知道多项式fe(x),fo(x)在w2n0,w2n1,w2n2....w2n2n−1的点值表达式,利用(1)可以O(n)求出f(x)在w2n0,w2n1,w2n2....w2n2n−1的点值表达式;利用(2)可以O(n)求出f(x)在w2n2n,w2n2n+1,w2n2n+2....w2nn−1的点值表达式;
既然这样,那么我们再3.1.1中 O ( n l o g n ) O(nlogn) O(nlogn)求多项式的思路是可行的!!!!
3.1.3:优化分治(蝴蝶变换)
我们用下面一张图来表示3.1.1中的分治过程:
显然,对于上述过程是可以使用递归的写法,但是递归的写法常数过大,
但是我们需要观察每个数最终的位置,最终发现对于幂数
i
i
i而言,其实
i
i
i最终在序列中的位置其实是
i
i
i在二进制下翻转(蝴蝶变换)得到的数
设
i
通过蝴蝶变换后的位置为
r
e
v
e
r
s
e
[
i
]
,
i
在二进制下的位数为
l
e
n
(
i
)
那么有
r
e
v
e
r
s
e
[
i
]
=
⌊
r
e
v
e
r
s
e
⌊
x
2
⌋
2
⌋
+
(
i
×
(
m
o
d
2
)
)
×
(
l
e
n
(
i
)
2
)
设i通过蝴蝶变换后的位置为reverse[i],i在二进制下的位数为len(i) \\ 那么有reverse[i]=⌊\frac{reverse⌊\frac{x}{2}⌋}{2}⌋+(i\times(mod2))\times(\frac{len(i)}{2})
设i通过蝴蝶变换后的位置为reverse[i],i在二进制下的位数为len(i)那么有reverse[i]=⌊2reverse⌊2x⌋⌋+(i×(mod2))×(2len(i))
因此我们可以把每个数都放在其最终位置上,不断往上迭代即可!
**#**附一个蝴蝶变换的小结论:
蝴蝶变换过后的序列,任取三个数都不构成等差数列。
3.2、IDFT
现在我们已经知道 k k k次多项式 f ( x ) f(x) f(x)的点值表达式: { ( w n 0 , f ( w n 0 ) ) , ( w n 1 , f ( w n 1 ) ) , ( w n 2 , f ( w n 2 ) ) , . . . . , ( w n n , f ( w n n ) ) } ( k ≤ n , n = 2 m − 1 ) \{(w_n^0,f(w_n^0)),(w_n^1,f(w_n^1)),(w_n^2,f(w_n^2)),....,(w_n^n,f(w_n^n))\} (k\leq n,n=2^m-1) {(wn0,f(wn0)),(wn1,f(wn1)),(wn2,f(wn2)),....,(wnn,f(wnn))}(k≤n,n=2m−1)
那么可一把点值表达式改写成与之对应的矩阵的形式:
[
1
(
w
n
0
)
1
(
w
n
0
)
2
…
(
w
n
0
)
n
1
(
w
n
1
)
1
(
w
n
1
)
2
…
(
w
n
1
)
n
1
(
w
n
2
)
1
(
w
n
2
)
2
…
(
w
n
2
)
n
⋮
⋮
⋮
⋱
⋮
1
(
w
n
n
)
1
(
w
n
n
)
2
…
(
w
n
n
)
n
]
×
[
a
0
a
1
a
2
⋮
a
n
]
=
[
f
(
w
n
0
)
f
(
w
n
1
)
f
(
w
n
2
)
⋮
f
(
w
n
n
)
]
\begin{bmatrix} &1 &(w_n^0)^1 &(w_n^0)^2 &\dots&(w_n^0)^n\\ &1 &(w_n^1)^1 &(w_n^1)^2 &\dots&(w_n^1)^n\\ &1 &(w_n^2)^1 &(w_n^2)^2 &\dots&(w_n^2)^n\\ &\vdots &\vdots &\vdots &\ddots &\vdots\\ &1 &(w_n^n)^1 &(w_n^n)^2 &\dots&(w_n^n)^n\\ \end{bmatrix} \times \begin{bmatrix} &a_0\\ &a_1\\ &a_2\\ &\vdots\\ &a_n \end{bmatrix} =\begin{bmatrix} &f(w_n^0)\\ &f(w_n^1)\\ &f(w_n^2)\\ &\vdots\\ &f(w_n^n) \end{bmatrix}
111⋮1(wn0)1(wn1)1(wn2)1⋮(wnn)1(wn0)2(wn1)2(wn2)2⋮(wnn)2………⋱…(wn0)n(wn1)n(wn2)n⋮(wnn)n
×
a0a1a2⋮an
=
f(wn0)f(wn1)f(wn2)⋮f(wnn)
那么在这个矩阵乘法中:
[
a
0
a
1
a
2
⋮
a
n
]
即为多项式的系数表达式
\begin{bmatrix} &a_0\\ &a_1\\ &a_2\\ &\vdots\\ &a_n \end{bmatrix}即为多项式的系数表达式
a0a1a2⋮an
即为多项式的系数表达式
那么我们只需要两边同时乘以 [ 1 ( w n 0 ) 1 ( w n 0 ) 2 … ( w n 0 ) n 1 ( w n 1 ) 1 ( w n 1 ) 2 … ( w n 1 ) n 1 ( w n 2 ) 1 ( w n 2 ) 2 … ( w n 2 ) n ⋮ ⋮ ⋮ ⋱ ⋮ 1 ( w n n ) 1 ( w n n ) 2 … ( w n n ) n ] 的逆矩阵,便可以得到多项式的系数表达式 那么我们只需要两边同时乘以\begin{bmatrix} &1 &(w_n^0)^1 &(w_n^0)^2 &\dots&(w_n^0)^n\\ &1 &(w_n^1)^1 &(w_n^1)^2 &\dots&(w_n^1)^n\\ &1 &(w_n^2)^1 &(w_n^2)^2 &\dots&(w_n^2)^n\\ &\vdots &\vdots &\vdots &\ddots &\vdots\\ &1 &(w_n^n)^1 &(w_n^n)^2 &\dots&(w_n^n)^n\\ \end{bmatrix}的逆矩阵,便可以得到多项式的系数表达式 那么我们只需要两边同时乘以 111⋮1(wn0)1(wn1)1(wn2)1⋮(wnn)1(wn0)2(wn1)2(wn2)2⋮(wnn)2………⋱…(wn0)n(wn1)n(wn2)n⋮(wnn)n 的逆矩阵,便可以得到多项式的系数表达式
因为矩阵 [ f ( w n 0 ) f ( w n 1 ) f ( w n 2 ) ⋮ f ( w n n ) ] 是一个关于 f ( w n i ) 已知的矩阵 因为矩阵\begin{bmatrix} &f(w_n^0)\\ &f(w_n^1)\\ &f(w_n^2)\\ &\vdots\\ &f(w_n^n) \end{bmatrix}是一个关于f(w_n^i)已知的矩阵 因为矩阵 f(wn0)f(wn1)f(wn2)⋮f(wnn) 是一个关于f(wni)已知的矩阵
所以最后 I D F T 这一块的问题等价于求矩阵 [ 1 ( w n 0 ) 1 ( w n 0 ) 2 … ( w n 0 ) n 1 ( w n 1 ) 1 ( w n 1 ) 2 … ( w n 1 ) n 1 ( w n 2 ) 1 ( w n 2 ) 2 … ( w n 2 ) n ⋮ ⋮ ⋮ ⋱ ⋮ 1 ( w n n ) 1 ( w n n ) 2 … ( w n n ) n ] 的逆矩阵 在这里我们不妨考虑,该矩阵的每一行: 对于每 j 行有 : [ 1 ( w n j ) 1 ( w n j ) 2 … ( w n j ) n ] 肯定有逆矩阵的第 j 列与之相乘 [ 1 ( w n j ) 1 ( w n j ) 2 … ( w n j ) p … ( w n j ) n ] × [ x 0 x 1 x 2 ⋮ x q ⋮ x n ] ∴ 对于每一项有 w n j p × x q 这个我们考虑单位根性质: w n j × w n n − j = w n n = 1 ∴ 可以构造 x q = w n ( n − j ) p ∴ w n j p × x q = w n n p = 1 ∴ [ x 0 x 1 x 2 ⋮ x q ⋮ x n ] = [ w n ( n − j ) 0 w n ( n − j ) 1 w n ( n − j ) 2 ⋮ w n ( n − j ) p ⋮ w n ( n − j ) n ] 所以最后IDFT这一块的问题等价于求矩阵\begin{bmatrix} &1 &(w_n^0)^1 &(w_n^0)^2 &\dots&(w_n^0)^n\\ &1 &(w_n^1)^1 &(w_n^1)^2 &\dots&(w_n^1)^n\\ &1 &(w_n^2)^1 &(w_n^2)^2 &\dots&(w_n^2)^n\\ &\vdots &\vdots &\vdots &\ddots &\vdots\\ &1 &(w_n^n)^1 &(w_n^n)^2 &\dots&(w_n^n)^n\\ \end{bmatrix}的逆矩阵 \\ 在这里我们不妨考虑,该矩阵的每一行:\\ 对于每j行有:[1\space\space\space(w_n^j)^1\space\space\space(w_n^j)^2 \dots(w_n^j)^n] 肯定有逆矩阵的第j列与之相乘\\ \begin{bmatrix} &1 &(w_n^j)^1 &(w_n^j)^2 &\dots&(w_n^j)^p&\dots&(w_n^j)^n\\ \end{bmatrix} \times \begin{bmatrix} &x_0\\ &x_1\\ &x_2\\ &\vdots\\ &x_q\\ &\vdots\\ &x_n\\ \end{bmatrix} \\ \therefore对于每一项有w_n^{jp}\times x_q \\ 这个我们考虑单位根性质:w_n^j\times w_n^{n-j}=w_n^n=1\\ \therefore 可以构造x_q=w_n^{(n-j)p}\\ \therefore w_n^{jp}\times x_q=w_n^{np}=1 \\ \therefore \begin{bmatrix} &x_0\\ &x_1\\ &x_2\\ &\vdots\\ &x_q\\ &\vdots\\ &x_n\\ \end{bmatrix} =\begin{bmatrix} &w_n^{(n-j)0}\\ &w_n^{(n-j)1}\\ &w_n^{(n-j)2}\\ &\vdots\\ &w_n^{(n-j)p}\\ &\vdots\\ &w_n^{(n-j)n}\\ \end{bmatrix} 所以最后IDFT这一块的问题等价于求矩阵 111⋮1(wn0)1(wn1)1(wn2)1⋮(wnn)1(wn0)2(wn1)2(wn2)2⋮(wnn)2………⋱…(wn0)n(wn1)n(wn2)n⋮(wnn)n 的逆矩阵在这里我们不妨考虑,该矩阵的每一行:对于每j行有:[1 (wnj)1 (wnj)2…(wnj)n]肯定有逆矩阵的第j列与之相乘[1(wnj)1(wnj)2…(wnj)p…(wnj)n]× x0x1x2⋮xq⋮xn ∴对于每一项有wnjp×xq这个我们考虑单位根性质:wnj×wnn−j=wnn=1∴可以构造xq=wn(n−j)p∴wnjp×xq=wnnp=1∴ x0x1x2⋮xq⋮xn = wn(n−j)0wn(n−j)1wn(n−j)2⋮wn(n−j)p⋮wn(n−j)n
∴ [ 1 ( w n j ) 1 ( w n j ) 2 … ( w n j ) p … ( w n j ) n ] × [ w n ( n − j ) 0 w n ( n − j ) 1 w n ( n − j ) 2 ⋮ w n ( n − j ) p ⋮ w n ( n − j ) n ] = n 又 ∵ ∀ d ∈ Z , d ≠ 0 且 n − j + d > 0 有 [ 1 ( w n j ) 1 ( w n j ) 2 … ( w n j ) p … ( w n j ) n ] × [ w n ( n − j + d ) 0 w n ( n − j + d ) 1 w n ( n − j + d ) 2 ⋮ w n ( n − j + d ) p ⋮ w n ( n − j + d ) n ] = 0 \therefore\begin{bmatrix} &1 &(w_n^j)^1 &(w_n^j)^2 &\dots&(w_n^j)^p&\dots&(w_n^j)^n\\ \end{bmatrix} \times\begin{bmatrix} &w_n^{(n-j)0}\\ &w_n^{(n-j)1}\\ &w_n^{(n-j)2}\\ &\vdots\\ &w_n^{(n-j)p}\\ &\vdots\\ &w_n^{(n-j)n}\\ \end{bmatrix}=n\\ 又\because \forall d\in Z,d\neq 0且n-j+d>0有 \begin{bmatrix} &1 &(w_n^j)^1 &(w_n^j)^2 &\dots&(w_n^j)^p&\dots&(w_n^j)^n\\ \end{bmatrix} \times\begin{bmatrix} &w_n^{(n-j+d)0}\\ &w_n^{(n-j+d)1}\\ &w_n^{(n-j+d)2}\\ &\vdots\\ &w_n^{(n-j+d)p}\\ &\vdots\\ &w_n^{(n-j+d)n}\\ \end{bmatrix}=0 ∴[1(wnj)1(wnj)2…(wnj)p…(wnj)n]× wn(n−j)0wn(n−j)1wn(n−j)2⋮wn(n−j)p⋮wn(n−j)n =n又∵∀d∈Z,d=0且n−j+d>0有[1(wnj)1(wnj)2…(wnj)p…(wnj)n]× wn(n−j+d)0wn(n−j+d)1wn(n−j+d)2⋮wn(n−j+d)p⋮wn(n−j+d)n =0
∴ [ 1 ( w n 0 ) 1 ( w n 0 ) 2 … ( w n 0 ) n 1 ( w n 1 ) 1 ( w n 1 ) 2 … ( w n 1 ) n 1 ( w n 2 ) 1 ( w n 2 ) 2 … ( w n 2 ) n ⋮ ⋮ ⋮ ⋱ ⋮ 1 ( w n n ) 1 ( w n n ) 2 … ( w n n ) n ] × [ 1 ( w n n − 0 ) 1 ( w n n − 0 ) 2 … ( w n n − 0 ) n 1 ( w n n − 1 ) 1 ( w n n − 1 ) 2 … ( w n n − 1 ) n 1 ( w n n − 2 ) 1 ( w n n − 2 ) 2 … ( w n n − 2 ) n ⋮ ⋮ ⋮ ⋱ ⋮ 1 ( w n n − n ) 1 ( w n n − n ) 2 … ( w n n − n ) n ] = [ n 0 0 … 0 0 n 0 … 0 0 0 n … 0 ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 0 … n ] \therefore\begin{bmatrix} &1 &(w_n^0)^1 &(w_n^0)^2 &\dots&(w_n^0)^n\\ &1 &(w_n^1)^1 &(w_n^1)^2 &\dots&(w_n^1)^n\\ &1 &(w_n^2)^1 &(w_n^2)^2 &\dots&(w_n^2)^n\\ &\vdots &\vdots &\vdots &\ddots &\vdots\\ &1 &(w_n^n)^1 &(w_n^n)^2 &\dots&(w_n^n)^n\\ \end{bmatrix} \times \begin{bmatrix} &1 &(w_n^{n-0})^1 &(w_n^{n-0})^2 &\dots&(w_n^{n-0})^n\\ &1 &(w_n^{n-1})^1 &(w_n^{n-1})^2 &\dots&(w_n^{n-1})^n\\ &1 &(w_n^{n-2})^1 &(w_n^{n-2})^2 &\dots&(w_n^{n-2})^n\\ &\vdots &\vdots &\vdots &\ddots &\vdots\\ &1 &(w_n^{n-n})^1 &(w_n^{n-n})^2 &\dots&(w_n^{n-n})^n\\ \end{bmatrix}= \begin{bmatrix} &n &0 &0 &\dots&0\\ &0 &n &0 &\dots&0\\ &0 &0 &n &\dots&0\\ &\vdots &\vdots &\vdots &\ddots &\vdots\\ &0 &0 &0 &\dots&n\\ \end{bmatrix} ∴ 111⋮1(wn0)1(wn1)1(wn2)1⋮(wnn)1(wn0)2(wn1)2(wn2)2⋮(wnn)2………⋱…(wn0)n(wn1)n(wn2)n⋮(wnn)n × 111⋮1(wnn−0)1(wnn−1)1(wnn−2)1⋮(wnn−n)1(wnn−0)2(wnn−1)2(wnn−2)2⋮(wnn−n)2………⋱…(wnn−0)n(wnn−1)n(wnn−2)n⋮(wnn−n)n = n00⋮00n0⋮000n⋮0………⋱…000⋮n
∴ 矩阵的逆矩阵为: 1 n [ 1 ( w n n − 0 ) 1 ( w n n − 0 ) 2 … ( w n n − 0 ) n 1 ( w n n − 1 ) 1 ( w n n − 1 ) 2 … ( w n n − 1 ) n 1 ( w n n − 2 ) 1 ( w n n − 2 ) 2 … ( w n n − 2 ) n ⋮ ⋮ ⋮ ⋱ ⋮ 1 ( w n n − n ) 1 ( w n n − n ) 2 … ( w n n − n ) n ] \therefore 矩阵的逆矩阵为:\frac{1}{n}\begin{bmatrix} &1 &(w_n^{n-0})^1 &(w_n^{n-0})^2 &\dots&(w_n^{n-0})^n\\ &1 &(w_n^{n-1})^1 &(w_n^{n-1})^2 &\dots&(w_n^{n-1})^n\\ &1 &(w_n^{n-2})^1 &(w_n^{n-2})^2 &\dots&(w_n^{n-2})^n\\ &\vdots &\vdots &\vdots &\ddots &\vdots\\ &1 &(w_n^{n-n})^1 &(w_n^{n-n})^2 &\dots&(w_n^{n-n})^n\\ \end{bmatrix} ∴矩阵的逆矩阵为:n1 111⋮1(wnn−0)1(wnn−1)1(wnn−2)1⋮(wnn−n)1(wnn−0)2(wnn−1)2(wnn−2)2⋮(wnn−n)2………⋱…(wnn−0)n(wnn−1)n(wnn−2)n⋮(wnn−n)n
∴ 系数 [ a 0 a 1 a 2 ⋮ a n ] = 1 n [ 1 ( w n n − 0 ) 1 ( w n n − 0 ) 2 … ( w n n − 0 ) n 1 ( w n n − 1 ) 1 ( w n n − 1 ) 2 … ( w n n − 1 ) n 1 ( w n n − 2 ) 1 ( w n n − 2 ) 2 … ( w n n − 2 ) n ⋮ ⋮ ⋮ ⋱ ⋮ 1 ( w n n − n ) 1 ( w n n − n ) 2 … ( w n n − n ) n ] × [ f ( w n 0 ) f ( w n 1 ) f ( w n 2 ) ⋮ f ( w n n ) ] \therefore 系数\begin{bmatrix} &a_0\\ &a_1\\ &a_2\\ &\vdots\\ &a_n \end{bmatrix}= \frac{1}{n}\begin{bmatrix} &1 &(w_n^{n-0})^1 &(w_n^{n-0})^2 &\dots&(w_n^{n-0})^n\\ &1 &(w_n^{n-1})^1 &(w_n^{n-1})^2 &\dots&(w_n^{n-1})^n\\ &1 &(w_n^{n-2})^1 &(w_n^{n-2})^2 &\dots&(w_n^{n-2})^n\\ &\vdots &\vdots &\vdots &\ddots &\vdots\\ &1 &(w_n^{n-n})^1 &(w_n^{n-n})^2 &\dots&(w_n^{n-n})^n\\ \end{bmatrix}\times\begin{bmatrix} &f(w_n^0)\\ &f(w_n^1)\\ &f(w_n^2)\\ &\vdots\\ &f(w_n^n) \end{bmatrix} ∴系数 a0a1a2⋮an =n1 111⋮1(wnn−0)1(wnn−1)1(wnn−2)1⋮(wnn−n)1(wnn−0)2(wnn−1)2(wnn−2)2⋮(wnn−n)2………⋱…(wnn−0)n(wnn−1)n(wnn−2)n⋮(wnn−n)n × f(wn0)f(wn1)f(wn2)⋮f(wnn)
最后我们注意到对于 w n n − j w_n^{n-j} wnn−j,其实它为 w n j w_n^j wnj的共轭复数;因此在IDFT中,只需将DFT中的乘以单位根变为乘以其共轭复数即可,但是分治结束后,需要除以 n n n
其实上述需要求逆的矩阵是一个范德蒙德矩阵,也可以利用范德蒙德矩阵求逆的过程
好了!耗时三天终于敲完FFT的公式推导了!!!下面的代码实现就比较容易了!
4、FFT代码实现
直接上模版:
double PI = acos(-1.0);
//复数类
struct Complex
{
double r, i;
Complex() {}
Complex(double r, double i) : r(r), i(i) {}
inline void real(const double &x) { r = x; }
inline double real() { return r; }
inline Complex operator+(const Complex &rhs) const
{
return Complex(r + rhs.r, i + rhs.i);
}
inline Complex operator-(const Complex &rhs) const
{
return Complex(r - rhs.r, i - rhs.i);
}
inline Complex operator*(const Complex &rhs) const
{
return Complex(r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r);
}
inline Complex operator*(const double &rhs) const
{
return Complex(r * rhs, i * rhs);
}
inline void operator/=(const double &x)
{
r /= x, i /= x;
}
inline void operator*=(const Complex &rhs)
{
*this = Complex(r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r);
}
inline void operator+=(const Complex &rhs)
{
r += rhs.r, i += rhs.i;
}
inline Complex conj()
{
return Complex(r, -i);
}
};
// i的二进制翻转数组
int bitrev[maxn];
// FFT预处理
void fft_prepare(int n)
{
int bits = 0;
//由于FFT需要项数为2的整数次方倍,所以多项式(形式)次数为第一个大于n的二的正整数次方,bits表示多少位
while ((1 << bits) < n)
{
bits++;
}
for (int i = 0; i < n; i++)
{
//递推式:二进制翻转
// 0 1 2 3 4 5 6 7
// 0|4|2|6|1|5|3|7
bitrev[i] = bitrev[i >> 1] >> 1 | ((i & 1) << (bits - 1));
}
}
//快速傅里叶变换系数表达式转换为点值表达式
void fft(Complex *a, int n, int type = 1)
{
//求出要迭代的序列
for (int i = 0; i < n; i++)
{
if (i < bitrev[i])
{
swap(a[i], a[bitrev[i]]);
}
}
//待合并区间的中点
for (int mid = 1; mid < n; mid <<= 1)
{
//单位根
Complex Wn{cos(PI / mid), type * sin(PI / mid)};
// R是区间的右端点,j表示已经到哪个位置
for (int R = mid << 1, j = 0; j < n; j += R)
{
Complex w{1, 0}; //幂
//枚举左半部分
for (int k = 0; k < mid; k++, w *= Wn)
{
//蝴蝶效应
Complex x = a[j + k], y = w * a[j + mid + k];
a[j + k] = x + y;
a[j + k + mid] = x - y;
}
}
}
// type为-1则进行逆变换,即IDFT
if (type == -1)
{
for (int i = 0; i < n; i++)
{
a[i].r /= n;
}
}
}
代码有注释,就不细说了!
5、FFT的瓶颈
虽然FFT的在时间复杂度上非常优秀,但是计算机在计算时仍然会遇到一些问题!
例如:由于复数类中使用的是double类型,而且在计算单位根存在除法,除此之外计算机计算下 π \pi π始终是一个逼近值,因此必然也必然会存在误差!所以在特殊情况下计算机计算会出现一定的误差!
但是有没有解决办法喃?
可以说在一定情况下是可以避免误差的!这个时候就要引出原根
6、To NTT
原根:
设
m
∈
N
∗
,
a
∈
Z
,
若
g
c
d
(
a
,
m
)
=
1
,
且
δ
m
(
a
)
=
φ
(
m
)
,
则称
a
为模
m
的原根
注:
δ
m
(
a
)
为
a
模
m
的阶,
φ
为
m
的欧拉函数
设m\in N^*,a\in Z,若gcd(a,m)=1,且\delta_m(a)=\varphi(m),则称a为模m的原根\\ 注:\delta_m(a)为a模m的阶,\varphi为m的欧拉函数
设m∈N∗,a∈Z,若gcd(a,m)=1,且δm(a)=φ(m),则称a为模m的原根注:δm(a)为a模m的阶,φ为m的欧拉函数
通过定义我们可以知道,原根是一种定义在模意义下的定义,因此对于FFT遇到的的精度问题,我们也只能保证在模意义下的答案正确
我们通过推导可以发现原根具有和单位根几乎一样的性质,因此我们可以通过原根推出基本相同的结论
这样一来我们就不必考虑复数域了!!
也可以理解为NTT其实就是把FFT中的单位根换成了原根
6.1 NTT代码实现
const int g = 3, gi = 332748118, mod = 998244353;
// 998244353的一个原根为3且998244353-1=2^23*119,3在模998244353意义下的逆元为332748118
int qp(int a, int n)
{
int ans = 1;
while (n)
{
if (n & 1)
{
ans = ans * a % mod;
}
a = a * a % mod;
n >>= 1;
}
return ans;
}
// i的二进制翻转数组
int bitrev[maxn];
// NTT预处理
void ntt_prepare(int n)
{
int bits = 0;
//由于NTT需要项数为2的整数次方倍,所以多项式(形式)次数为第一个大于n的二的正整数次方,bits表示多少位
while ((1 << bits) < n)
{
bits++;
}
for (int i = 0; i < n; i++)
{
//递推式:二进制翻转
// 0 1 2 3 4 5 6 7
// 0|4|2|6|1|5|3|7
bitrev[i] = bitrev[i >> 1] >> 1 | ((i & 1) << (bits - 1));
}
}
// NTT,type=1时系数表示法转点值表示法,否则点值表示法转系数表示法
void NTT(int *a, int n, int type = 1)
{
for (int i = 0; i < n; i++)
{
if (i < bitrev[i])
{
swap(a[i], a[bitrev[i]]);
}
}
//在这之前NTT与FFT无异
for (int mid = 1; mid < n; mid <<= 1)
{
//单位原根g_n
int gn = qp(type == 1 ? g : gi, (mod - 1) / (mid << 1));
//枚举具体区间,j也就是区间右端点
for (int R = mid << 1, j = 0; j < n; j += R)
{
int g0 = 1;
//合并,记得要取模
for (int k = 0; k < mid; k++, g0 = gn * g0 % mod)
{
int x = a[j + k], y = g0 * a[mid + j + k] % mod;
a[j + k] = (x + y) % mod;
a[j + k + mid] = (x - y + mod) % mod;
}
}
}
}
可以发现,NTT和FFT代码大体逻辑是相同的!!!只是从单位根换成了原根
6.2 tips
在抽象代数中,原根就是循环群的生成元。这个概念只在模 m m m缩剩余系关于乘法形成的群中有“原根”这个名字,在一般的循环群中都称作“生成元”。
并非每个模 m m m缩剩余系关于乘法形成的群都是循环群,存在原根就表明它同构于循环群,如果不存在原根就表明不同构。