文章目录
FFT ( 快速傅里叶变换 ) 学习笔记
本文有错误之处请诸位大佬多多指正!
F F T FFT FFT :快速傅里叶变换的英文缩写,快速傅里叶变换是对离散傅里叶变换 D F T DFT DFT 的优化,用来解决多项式上的操作如 卷积 等问题。
参考文章:
https://www.luogu.com.cn/blog/command-block/fft-xue-xi-bi-ji
https://blog.csdn.net/weixin_43346722/article/details/116353343
https://www.cnblogs.com/zwfymqz/p/8244902.html
多项式
系数表示法
一般在 初中 数学上,表示一个多项式我们用的是系数表示法
设
A
(
x
)
A(x)
A(x) 表示一个关于
x
x
x 的
n
n
n 次的多项式,那么:
A
(
x
)
=
∑
i
=
0
n
a
i
∗
x
i
A(x)=\sum_{i=0}^{n}a_i*x^i
A(x)=i=0∑nai∗xi
展开后就是:
A
(
x
)
=
a
0
+
a
1
∗
x
+
a
2
∗
x
2
+
⋅
⋅
⋅
⋅
⋅
⋅
+
a
n
∗
x
n
A(x)=a_0+a_1*x+a_2*x^2+······+a_n*x^n
A(x)=a0+a1∗x+a2∗x2+⋅⋅⋅⋅⋅⋅+an∗xn
然后两个多项式相乘,每一项都要和另一个多项式的每一项相乘,时间复杂度
O
(
n
2
)
O(n^2)
O(n2)
点值表示法
我们找
n
+
1
n+1
n+1 组不同的
x
x
x 代入多项式,就可以得出
n
+
1
n+1
n+1 个
y
y
y 值,
我们把它们对应:
(
x
0
,
y
0
)
,
(
x
1
,
y
1
)
⋅
⋅
⋅
⋅
⋅
⋅
(
x
n
,
y
n
)
(x_0,y_0),(x_1,y_1)······(x_n,y_n)
(x0,y0),(x1,y1)⋅⋅⋅⋅⋅⋅(xn,yn)
其中
y
i
=
∑
j
=
0
n
−
1
a
j
∗
x
i
j
y_i=\sum\limits_{j=0}^{n-1}a_j*x_i^j
yi=j=0∑n−1aj∗xij
听说这玩意叫做
D
F
T
DFT
DFT,怎么和我的名字缩写一样
当然,你会发现时间复杂度还是没有变化,选点
O
(
n
)
O(n)
O(n),计算
O
(
n
)
O(n)
O(n),总共
O
(
n
2
)
O(n^2)
O(n2)
然后好像没有啥卵用
最后你可悲的发现,好像时间复杂度在
O
(
n
2
)
O(n^2)
O(n2)处达到了下限,没有什么优化的空间。
但是数学家们的智慧不是你能揣测的。
前方第一波高能!!!
复数
前置芝士
向量
具有大小又有方向的量,在物理上又称矢量,比如说力,加速度等都是向量。
在图上一般用一根箭头来表示。
弧度制
1
∘
=
π
180
r
a
d
1~^∘=\frac {\pi}{180}~~rad
1 ∘=180π rad
180
∘
=
1
π
r
a
d
180~^∘=1~~\pi~~rad
180 ∘=1 π rad
定义
我们设
a
,
b
∈
R
a,b\in R
a,b∈R,
i
2
=
−
1
i^2=-1
i2=−1,所有形如
a
+
b
i
a+bi
a+bi 的数通称为复数,其中
i
i
i 被称为虚数单位,
b
i
bi
bi 被称为复数的虚部,
a
a
a 被称为实部,如果一个复数只有
b
i
bi
bi 部分,那么这个复数被称为纯虚数。
复数域
C
C
C 是目前最大的域,实数集是复数集的真子集。
如何表示复数?
我们建立一个平面直角坐标系,称之为复平面, x x x 轴代表实轴, y y y 轴代表虚轴 (除原点外),那么 x x x 轴上的点表示的就是实数(类比数轴), y y y 轴除原点外表示的都是纯虚数,原点表示实数 0 0 0 。那么对于每一个复数 a + b i a+bi a+bi,我们都可以用复平面上的点 ( a , b ) (a,b) (a,b) 来表示这个复数,或者我们可以看做是起点为 ( 0 , 0 ) (0,0) (0,0) 终点为 ( a , b ) (a,b) (a,b) 的向量。
复数的模长:就是向量的大小,具体来说就是
a
2
+
b
2
\sqrt {a^2+b^2}
a2+b2
幅角:假以逆时针为正方向,从
x
x
x 轴正半轴到已知向量的转角的有向角叫做幅角.
运算法则
加法遵循向量的平行四边形法则,这个玩意长这样:(盗的别人的图)
A
B
⃗
+
B
C
⃗
=
A
D
⃗
+
D
C
⃗
=
A
C
⃗
\vec{AB}+\vec{BC}=\vec{AD}+\vec{DC}=\vec{AC}
AB+BC=AD+DC=AC
减法:加法的时候将向量变成反过来的(终点变起点,起点变终点)
乘法:
几何定义:模长相乘,幅角相加(虽然我也不懂)
证明:
对于每一个复数,都可以表示成
r
(
cos
x
+
i
∗
sin
x
)
r(\cos x+i*\sin x )
r(cosx+i∗sinx) 的形式,其中
r
r
r 为复数的模长,
x
x
x 是幅角 (不理解的建议学学三角函数 勾股定理 )
对于两个复数
r
1
(
cos
x
+
i
∗
sin
x
)
,
r
2
(
cos
y
+
i
∗
sin
y
)
r_1(\cos x+i*\sin x )~,~r_2(\cos y+i*\sin y )
r1(cosx+i∗sinx) , r2(cosy+i∗siny)
相乘,得到的答案: (
i
2
=
−
1
i^2=-1
i2=−1 所以第二项变成了减号 )
r
1
∗
r
2
∗
[
(
cos
x
cos
y
−
sin
x
sin
y
+
i
∗
(
cos
x
sin
y
+
cos
y
sin
x
)
]
r_1*r_2*[(\cos x\cos y-\sin x\sin y+i*(\cos x\sin y+\cos y\sin x)]
r1∗r2∗[(cosxcosy−sinxsiny+i∗(cosxsiny+cosysinx)]
根据三角恒等变换式 (不会的翻翻高中数学必修第一册最后一章)
cos
(
x
+
y
)
=
cos
x
cos
y
−
sin
x
sin
y
,
sin
(
x
+
y
)
=
sin
x
cos
y
+
sin
y
cos
x
\cos(x+y)=\cos x\cos y-\sin x\sin y~~,~~\sin(x+y)=\sin x\cos y+\sin y\cos x
cos(x+y)=cosxcosy−sinxsiny , sin(x+y)=sinxcosy+sinycosx
再看看上面的式子,我们发现可以化简:
r
1
r
2
∗
(
cos
(
x
+
y
)
+
i
∗
sin
(
x
+
y
)
)
r_1r_2*(\cos(x+y)+i*\sin(x+y))
r1r2∗(cos(x+y)+i∗sin(x+y))
所以就有了模长相乘 (显然),幅角相加 (还是显然)
代数定义:
(
a
+
b
i
)
∗
(
c
+
d
i
)
=
a
c
+
a
d
i
+
b
c
i
+
b
d
i
2
(a+bi)*(c+di)=ac+adi+bci+bdi^2
(a+bi)∗(c+di)=ac+adi+bci+bdi2
因为
i
2
=
−
1
i^2=-1
i2=−1 ,就有
(
a
c
−
b
d
)
+
(
a
d
+
b
c
)
i
(ac-bd)+(ad+bc)i
(ac−bd)+(ad+bc)i
搞定。
然后你会一头雾水的发现,这和多项式乘法有什么关系啊?
别急,真正的好东西来了。
单位根
注:以下提到的复数都是在单位圆上的复数。
在复平面上,以原点为圆心,
1
1
1 为半径作圆,所得的圆叫单位圆
。
又去盗了一张图
以圆点为起点,圆的
n
n
n 等分点为终点,做
n
n
n个向量,设幅角为正且最小的向量对应的复数为
ω
n
1
ω_n^1
ωn1,称为
n
n
n 次单位根。
就是
x
n
=
1
x^n=1
xn=1 在复数域上的
n
n
n 个解
根据复数乘法的运算,其余
n
−
1
n−1
n−1 个复数为
ω
n
2
,
ω
n
3
,
…
,
ω
n
n
ω^2_n,ω^3_n,…,ω^n_n
ωn2,ωn3,…,ωnn
其中
ω
n
0
=
ω
n
n
=
1
\omega^0_n=\omega_n^n=1
ωn0=ωnn=1(对应在
x
x
x 轴上长度为
1
1
1 的向量)
然后我们就有
ω
n
k
=
cos
k
2
π
n
+
i
∗
sin
k
2
π
n
\omega^k_n=\cos k~\frac{2\pi}{n}+i*\sin k~\frac{2\pi}{n}
ωnk=cosk n2π+i∗sink n2π
单位根的几条性质:
(
1
)
:
ω
n
k
=
(
ω
n
1
)
k
(1):~~~~\omega_n^k=(\omega_n^1)^k
(1): ωnk=(ωn1)k
证明:两个向量相乘,模长相乘还是 1 ,幅角相加,总共加了
k
k
k 次 (也可以看做向逆时针方向旋转了
k
n
\frac{k}{n}
nk 个圆周),所以就是
ω
n
k
\omega_n^k
ωnk (显而易见)
(2):
ω
n
i
∗
ω
n
j
=
ω
n
i
+
j
\omega_n^i*\omega_n^j=\omega_n^{i+j}
ωni∗ωnj=ωni+j
证明:
ω
n
i
=
(
ω
n
1
)
i
,
ω
n
j
=
(
ω
n
1
)
j
,
ω
n
i
∗
ω
n
j
=
(
ω
n
1
)
i
∗
(
ω
n
1
)
j
=
(
ω
n
1
)
i
+
j
=
ω
n
i
+
j
\omega_n^i=(\omega_n^1)^i~,~\omega_n^j=(\omega_n^1)^j~,~\omega_n^i*\omega_n^j=(\omega_n^1)^i*(\omega_n^1)^j=(\omega_n^1)^{i+j}=\omega_n^{i+j}
ωni=(ωn1)i , ωnj=(ωn1)j , ωni∗ωnj=(ωn1)i∗(ωn1)j=(ωn1)i+j=ωni+j
(3):
ω
n
k
=
ω
p
n
p
k
\omega_n^k=\omega_{pn}^{pk}
ωnk=ωpnpk
证明:把
ω
p
n
p
k
\omega_{pn}^{pk}
ωpnpk 这玩意代入到单位根的表示的式子里,发现分子分母同时乘了一个
p
p
p ,复数的值不变
(4):
ω
n
k
+
n
2
=
−
ω
n
k
\omega_n^{k+\frac{n}{2}}=-\omega_n^k
ωnk+2n=−ωnk
证明:
ω
n
n
2
=
cos
n
2
⋅
2
π
n
+
i
∗
sin
n
2
⋅
2
π
n
=
cos
π
+
i
∗
sin
π
\omega_n^\frac{n}{2}=\cos \frac{n}{2}·\frac{2\pi}{n}+i*\sin \frac{n}{2}·\frac{2\pi}{n}=\cos\pi+i*\sin\pi
ωn2n=cos2n⋅n2π+i∗sin2n⋅n2π=cosπ+i∗sinπ
众所周知 根据三角函数表,其中
cos
π
=
−
1
,
sin
π
=
0
\cos\pi=-1~,~\sin\pi=0
cosπ=−1 , sinπ=0
所以 ω n n 2 = − 1 , ω n k + n 2 = ω n k ∗ ω n n 2 = − ω n k \omega_n^\frac{n}{2}=-1~,~\omega_n^{k+\frac{n}{2}}=\omega_n^k*\omega_n^{\frac{n}{2}}=-\omega_n^k ωn2n=−1 , ωnk+2n=ωnk∗ωn2n=−ωnk
当然你可以看做把这个复数向逆时针方向旋转了 n 2 \frac{n}{2} 2n 个圆周 (就是一个半圆),然后这个新的复数和之前的复数就关于原点对称了 (建议画图理解)
(5):
ω
n
k
(
k
≥
n
)
=
ω
n
k
m
o
d
n
\omega_n^k~(k\geq n)=\omega_n^{k\bmod n}
ωnk (k≥n)=ωnkmodn
证明:很明显这个复数围绕单位圆转了
n
n
=
1
\frac{n}{n}=1
nn=1 圈后回到了之前的位置,对答案没有影响。
用三角函数的角度来理解就是终边相同的角三角函数值相等
正文真正的开始
快速傅里叶变换:
前文说到,对于一个
n
n
n 次多项式,我们可以用
n
+
1
n+1
n+1 个点来确定
对于一个
n
−
1
n-1
n−1 次多项式
A
(
x
)
A(x)
A(x),我们的超级神仙傅里叶 降下了一道救世之光,拯救了那些被多项式梦魇诅咒的善良之人 反手将这个多项式拍成了渣渣,然后按照奇偶性分成了两半,就有了这个东西:
A
(
x
)
=
(
a
0
+
a
2
x
2
+
⋅
⋅
⋅
+
a
n
−
2
x
n
−
2
)
+
(
a
1
x
+
a
3
x
3
+
⋅
⋅
⋅
+
a
n
−
1
x
n
−
1
)
A(x)=(a_0+a_2x^2+···+a_{n-2}x^{n-2})+(a_1x+a_3x^3+···+a_{n-1}x^{n-1})
A(x)=(a0+a2x2+⋅⋅⋅+an−2xn−2)+(a1x+a3x3+⋅⋅⋅+an−1xn−1)
这时候我们设
A
1
(
x
)
=
a
0
+
a
2
x
+
⋅
⋅
⋅
+
a
n
−
2
x
n
2
−
1
A_1(x)=a_0+a_2x+···+a_{n-2}x^{\frac{n}{2}-1}
A1(x)=a0+a2x+⋅⋅⋅+an−2x2n−1
A
2
(
x
)
=
a
1
+
a
3
x
+
⋅
⋅
⋅
+
a
n
−
1
x
n
2
−
1
A_2(x)=a_1+a_3x+···+a_{n-1}x^{\frac{n}{2}-1}
A2(x)=a1+a3x+⋅⋅⋅+an−1x2n−1
对应上面把
A
(
x
)
A(x)
A(x)切成两半的操作,我们 惊悚 惊喜地发现:
A
(
x
)
=
A
1
(
x
2
)
+
x
A
2
(
x
2
)
A(x)=A_1(x^2)+xA_2(x^2)
A(x)=A1(x2)+xA2(x2)
这时候我们把
ω
n
k
\omega_n^k
ωnk (
k
≤
n
2
k\leq \frac{n}{2}
k≤2n ) 代入这个式子:
A
(
ω
n
k
)
=
A
1
(
ω
n
2
k
)
+
ω
n
k
A
2
(
ω
n
2
k
)
A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})
A(ωnk)=A1(ωn2k)+ωnkA2(ωn2k)
然后再把
ω
n
k
+
n
2
\omega_n^{k+\frac{n}{2}}
ωnk+2n 代入,得到:
A
(
ω
n
k
+
n
2
)
=
A
1
(
ω
n
2
k
+
n
)
+
ω
n
k
+
n
2
A
2
(
ω
n
2
k
+
n
)
=
A
1
(
ω
n
2
k
)
+
ω
n
k
+
n
2
A
2
(
ω
n
2
k
)
A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k+n})=A_1(\omega_n^{2k})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k})
A(ωnk+2n)=A1(ωn2k+n)+ωnk+2nA2(ωn2k+n)=A1(ωn2k)+ωnk+2nA2(ωn2k)
再根据
ω
n
k
+
n
2
=
−
ω
n
k
\omega_n^{k+\frac{n}{2}}=-\omega_n^k
ωnk+2n=−ωnk,得到
A
(
ω
n
k
+
n
2
)
=
A
1
(
ω
n
2
k
)
−
ω
n
k
A
2
(
ω
n
2
k
)
A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k})-\omega_n^{k}A_2(\omega_n^{2k})
A(ωnk+2n)=A1(ωn2k)−ωnkA2(ωn2k)
我们惊喜的发现,两个式子只有一个常数项不同,那么你每一次枚举一个式子的时候,你可以在
O
(
1
)
O(1)
O(1) 时间复杂度之内求出第二个式子的值,那么你的问题就缩小了一半!这时候你还会发现,对于第一个式子,继续拆分,它的问题规模又可以减半,那就可以递归实现了 (是不是很玄学)
不要以为这样就结束了,逼事还有一大堆呢
快速傅里叶逆变换 (IFFT):
当你以为你搞定了前面的东西的时候,你发现你得到的是一堆点值表达式,好像并没有卵用
你悲哀的发现你只会系数表示法, 点值表示法你不会······
所以你就必须将这个点值表示法的多项式转换为系数表示法来表示。
所以神仙数学家们又搞了一个叫做快速傅里叶逆变换 (
I
F
F
T
IFFT
IFFT ) 的东西。
内容:
我们设
(
y
0
,
y
1
,
y
2
⋅
⋅
⋅
y
n
−
1
)
(y_0,y_1,y_2···y_{n-1})
(y0,y1,y2⋅⋅⋅yn−1) 为多项式
(
a
0
,
a
1
,
a
2
⋅
⋅
⋅
a
n
−
1
)
(a_0,a_1,a_2···a_{n-1})
(a0,a1,a2⋅⋅⋅an−1) 的傅里叶变换
这时候我们再设
(
c
0
,
c
1
,
c
2
⋅
⋅
⋅
c
n
−
1
)
(c_0,c_1,c_2···c_{n-1})
(c0,c1,c2⋅⋅⋅cn−1) 是
(
y
0
,
y
1
,
y
2
⋅
⋅
⋅
y
n
−
1
)
(y_0,y_1,y_2···y_{n-1})
(y0,y1,y2⋅⋅⋅yn−1) 在
ω
n
0
,
ω
n
−
1
⋅
⋅
⋅
,
ω
n
−
(
n
−
1
)
\omega_n^0,\omega_n^{-1}···,\omega_n^{-(n-1)}
ωn0,ωn−1⋅⋅⋅,ωn−(n−1) 处的傅里叶变换
那么就有
c
k
=
∑
i
=
0
n
−
1
y
i
∗
(
ω
n
−
k
)
i
c_k=\sum\limits_{i=0}^{n-1}y_i*(\omega_n^{-k})^i
ck=i=0∑n−1yi∗(ωn−k)i
推式子高能警告!!!
展开式子:
c
k
=
∑
i
=
0
n
−
1
(
∑
j
=
0
n
−
1
a
j
(
ω
n
i
)
j
)
∗
(
ω
n
−
k
)
i
)
c_k=\sum\limits_{i=0}^{n-1}(\sum \limits_{j=0}^{n-1}a_j(\omega_n^i)^j)*(\omega_n^{-k})^i)
ck=i=0∑n−1(j=0∑n−1aj(ωni)j)∗(ωn−k)i)
= ∑ i = 0 n − 1 ( ∑ j = 0 n − 1 a j ( ω n j ) i ) ∗ ( ω n − k ) i ) =\sum\limits_{i=0}^{n-1}(\sum \limits_{j=0}^{n-1}a_j(\omega_n^j)^i)*(\omega_n^{-k})^i) =i=0∑n−1(j=0∑n−1aj(ωnj)i)∗(ωn−k)i)
= ∑ i = 0 n − 1 ( ∑ j = 0 n − 1 a j ( ω n j ) i ∗ ( ω n − k ) i ) =\sum\limits_{i=0}^{n-1}(\sum\limits _{j=0}^{n-1}a_j(\omega_n^j)^i*(\omega_n^{-k})^i) =i=0∑n−1(j=0∑n−1aj(ωnj)i∗(ωn−k)i)
= ∑ i = 0 n − 1 ( ∑ j = 0 n − 1 a j ( ω n j − k ) i ) =\sum\limits_{i=0}^{n-1}(\sum \limits_{j=0}^{n-1}a_j(\omega_n^{j-k})^i) =i=0∑n−1(j=0∑n−1aj(ωnj−k)i)
= ∑ i = 0 n − 1 ( ∑ j = 0 n − 1 a j ( ω n j − k ) i ) =\sum\limits_{i=0}^{n-1}(\sum\limits _{j=0}^{n-1}a_j(\omega_n^{j-k})^i) =i=0∑n−1(j=0∑n−1aj(ωnj−k)i)
= ∑ j = 0 n − 1 a j ( ∑ i = 0 n − 1 ( ω n j − k ) i ) =\sum\limits_{j=0}^{n-1}a_j(\sum\limits _{i=0}^{n-1}(\omega_n^{j-k})^i) =j=0∑n−1aj(i=0∑n−1(ωnj−k)i) 注意 j , i j,i j,i 位置的调换
这时候我们设
S
(
x
)
=
∑
i
=
0
n
−
1
x
i
S(x)=\sum\limits_{i=0}^{n-1}x_i
S(x)=i=0∑n−1xi
代入
ω
n
k
\omega_n^k
ωnk ,有
S
(
ω
n
k
)
=
1
+
ω
n
k
+
(
ω
n
k
)
2
+
⋅
⋅
⋅
+
(
ω
n
k
)
n
−
1
S(\omega_n^k)=1+\omega_n^k+(\omega_n^k)^2+···+(\omega_n^k)^{n-1}
S(ωnk)=1+ωnk+(ωnk)2+⋅⋅⋅+(ωnk)n−1
当
k
!
=
0
k!=0
k!=0 时,这时候
ω
n
k
!
=
1
\omega_n^k~!=1
ωnk !=1
将这个式子乘上
ω
n
k
\omega_n^k
ωnk ,就有了:
ω
n
k
S
(
ω
n
k
)
=
ω
n
k
+
(
ω
n
k
)
2
+
(
ω
n
k
)
3
+
⋅
⋅
⋅
+
(
ω
n
k
)
n
\omega_n^kS(\omega_n^k)=\omega_n^k+(\omega_n^k)^2+(\omega_n^k)^3+···+(\omega_n^k)^{n}
ωnkS(ωnk)=ωnk+(ωnk)2+(ωnk)3+⋅⋅⋅+(ωnk)n
后式减去前式,得到:
ω n k S ( ω n k ) − S ( ω n k ) = ( ω n k ) n − 1 \omega_n^kS(\omega_n^k)-S(\omega_n^k)=(\omega_n^k)^n-1 ωnkS(ωnk)−S(ωnk)=(ωnk)n−1
S ( ω n k ) = ( ω n k ) n − 1 ω n k − 1 S(\omega_n^k)=\frac{( \omega_n^k )^n-1}{\omega_n^k-1} S(ωnk)=ωnk−1(ωnk)n−1
S ( ω n k ) = ( ω n n ) k − 1 ω n k − 1 S(\omega_n^k)=\frac{( \omega_n^n )^k-1}{\omega_n^k-1} S(ωnk)=ωnk−1(ωnn)k−1
因为
ω
n
n
=
1
\omega_n^n=1
ωnn=1,所以这个分数的分子为
0
0
0,分数值为
0
0
0
当
k
=
0
k=0
k=0 的时候,
S
(
ω
n
k
)
=
1
+
1
+
1
+
⋅
⋅
⋅
⋅
⋅
=
n
S(\omega_n^k)=1+1+1+·····=n
S(ωnk)=1+1+1+⋅⋅⋅⋅⋅=n
这时候我们看回刚刚这个式子:
c
k
=
∑
j
=
0
n
−
1
a
j
(
∑
i
=
0
n
−
1
(
ω
n
j
−
k
)
i
)
c_k=\sum\limits_{j=0}^{n-1}a_j(\sum\limits _{i=0}^{n-1}(\omega_n^{j-k})^i)
ck=j=0∑n−1aj(i=0∑n−1(ωnj−k)i)
=
∑
j
=
0
n
−
1
a
j
S
(
ω
n
j
−
k
)
~~~~~=\sum\limits_{j=0}^{n-1}a_jS(\omega_n^{j-k})
=j=0∑n−1ajS(ωnj−k)
当
j
−
k
!
=
0
j-k~!=0
j−k !=0 时,贡献为
0
0
0
当
j
=
k
,
j
−
k
=
0
j=k,j-k=0
j=k,j−k=0 时,贡献为
n
n
n
所以就有
c
k
=
n
∗
a
k
c_k=n*a_k
ck=n∗ak
最后
a
k
=
c
k
n
a_k=\frac{c_k}{n}
ak=nck
这样我们就得到了点值表示法和系数表示法之间的关系啦。
这个过程和上一个过程一样,都可以分治递归来实现
迭代优化
你慢悠悠的打出了一个递归,交上洛谷,然后只有
77
77
77 分
没错,递归常数太大,在
1
0
6
10^6
106 的常数面前,它没了······ (但这并不表示这个算法是不可行的)
所以我们把递归扔掉,使用迭代
去盗了一个大佬的图
有没有看出什么显而易见的性质?
没错,我们需要求的序列实际是原序列下标的二进制反转 (真的玄学)
因此没有必要对数列的奇偶性进行分类
然后我就不会证了
来看看这位 command_block 大佬的博客
原文地址:https://www.luogu.com.cn/blog/command-block/fft-xue-xi-bi-ji
最后迭代就是从下往上合并
悄悄话
累死我了,终于打完了······
(热泪盈眶)
代码:
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<queue>
#include<cstring>
#include<cmath>
#define r register
#define rep(i,x,y) for(r ll i=x;i<=y;++i)
#define per(i,x,y) for(r ll i=x;i>=y;--i)
using namespace std;
typedef long long ll;
const ll V=(1e7+10)/2;
const double pi=acos(-1.0);
ll in()
{
ll res=0,f=1;
char ch;
while((ch=getchar())<'0'||ch>'9')
if(ch=='-') f=-1;
res=res*10+ch-48;
while((ch=getchar())>='0'&&ch<='9')
res=res*10+ch-48;
return res*f;
}
ll n,m,f[V];
ll now,k;
struct complex
{
double x,y;
complex (double xx=0,double yy=0) { x=xx,y=yy; }
}a[V],b[V];
complex operator + (complex a,complex b) { return complex(a.x+b.x,a.y+b.y); }
complex operator - (complex a,complex b) { return complex(a.x-b.x,a.y-b.y); }
complex operator * (complex a,complex b) { return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x); }
void FFT(complex *a,ll v)
{
rep(i,0,now-1)
if(i<f[i]) swap(a[i],a[f[i]]);
for(r ll mid=1;mid<now;mid<<=1)
{
complex Wn(cos(pi/mid),v*sin(pi/mid));
ll R=mid<<1;
for(r ll j=0;j<now;j+=R)
{
complex val(1,0);
for(r ll k=0;k<mid;++k,val=val*Wn)
{
complex x=a[j+k],y=val*a[j+mid+k];
a[j+k]=x+y;
a[j+mid+k]=x-y;
}
}
}
}
int main()
{
scanf("%lld%lld",&n,&m);
rep(i,0,n) a[i].x=in();
rep(i,0,m) b[i].x=in();
now=1;//单位根中的k值,默认为2的自然数次幂
while(now<=n+m) now<<=1,++k;
rep(i,0,now-1)
f[i]=(f[i>>1]>>1)|((i&1)<<(k-1));
FFT(a,1);
FFT(b,1);
rep(i,0,now) a[i]=a[i]*b[i];
FFT(a,-1);
rep(i,0,m+n)
printf("%lld ",(ll)(a[i].x/now+0.49));
return 0;
}