快速傅里叶变换FFT
作为我5.1出来培训这几天的原因,学会了之后当然要记录一下awa
虽然这是我唯一学会的内容qwq
首先,FFT是干什么的:
如果有这么两个东西:
A(x)=n−1∑i=0aixi
B(x)=n−1∑i=0bixi
让你求出 A(x)∗B(x)
啊,我会我会!直接 n2 枚举,进行多项式乘法就好了awa
是的,但是如果 n<=1e6 呢?
你会发现 n2 就解决不了了,而且,枚举的方法显然已经无法继续优化了,所以我们需要另辟蹊径。
首先我们设 C(x)=A(x)∗B(x)
然后如果我们能快速的求出 A(x) 和 B(x) 的点值,那么我们就可以把他们乘在一起,之后就可以得到 C(x) 的点值,然后就可以逆向唯一确定 C(x) 的表达式了
(原理:n 次多项式唯一的由 n + 1 个不同点的值确定。)
(求点值就是随便带入一个 x ,然后求出它在函数图像上的 y)
所以FFT就是干这个的:
(分为以下三步)
1.快速求出 A(x) 和 B(x) 的点值。(DFT)
2.将求出的 A(x) 和 B(x) 相乘
3.通过得到的 C(x) 的 n 个点值求出 C(x) 的 n 个系数 (IDFT)
(可以发现第三步求出的数就是答案)
开始求解:
我们可以记单位根 ωn=e2πin
( i 就是虚数,一般为 √−1 )
(选择 e2πin 的原因是 e2πi=1 ,一会方便计算awa)
那么我们就可以将选的 n 个点值设定为:
x=ω0n,ω1n,ω2n……ωn−1n
之后我们将 A(x)=n−1∑i=0aixi 拆成 A1(x2)+xA2(x2)
其中:
A1(x)=a0+a2x+a4x2+……+an−2xn2−1
A2(x)=a1+a3x+a5x2+……+an−1xn2−1
那我们再把我们设的 x=ωkn 带入:
A(ωkn)=A1((ωkn)2)+ωknA2((ωkn)2)
=A1((e2πi∗2n)k)+ωknA2((e2πi∗2n)k)
=A1((e2πin2)k)+ωknA2((e2πin2)k)
=A1(ωkn2)+ωknA2(ωkn2)
因为我们设的是x=ω0n,ω1n,ω2n……ωn−1n,所以当 ωi 是 n2 阶的时候, k 一定是小于 n2 的,所以这只是前一半。
然后对于另一半,我们可以带入带入 x=ωk+n2n ,可以得出:
A(ωk+n2n)=A1((ωkn)2∗ωnn)+ωk+n2nA2((ωkn)2∗ωnn)
=A1((e2πi∗2n)k∗1)+ωn2n∗ωknA2((e2πi∗2n)k∗1)
=A1((e2πin2)k)−ωknA2((e2πin2)k)
=A1(ωkn2)−ωknA2(ωkn2)
也就是说,如果我们已经求出了 A1,A2 在阶为 n2 时,某个单位根的点值,那么我们可以通过 O(n) 计算出 A(x) 在阶为 n 时,相应单位根的点值。
然后一直递归下去,当 n=1 时,ω1=e2πi1=e2πi=1 。答案就是 A(x) 的系数,然后我们再一层一层的返回,就可以求出最终答案。
这样就会发现,复杂度就是 T(n)=2T(n/2)+O(n)=O(nlogn) 的 awa
然后我们发现当上式可做时,最好需要 n=2k , 那我们就不妨让 n=2k ,然后多出来的位数就让它的系数为 0 就好了。
所以我们就通过分治的方法,在 O(logn) 的时间内,解决了第一步(DFT)awa
之后我们就分别对 A(x) 和 B(x) 求一遍,再将它们对应的点值相乘,就求出了 C(x) 的点值( 复杂度: O(n) )
第二步解决awa
然后是第三步。
乍一看第三步十分的困难 qwq ,但是我们从另一个角度去看第一步(DFT)
可以发现:我们实际上是把一个单位根矩阵 M=[ωij] 左乘在了系数向量上, 那么实现 IDFT 实际上我们需要找到这个单位根矩阵的逆。
(上面一句一听就不是人话,让我们来形象的翻译一下)
首先,我们可以发现, DFT 其实干了下面这个事情:
⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣(ω0n)0(ω0n)1…(ω0n)n−2(ω0n)n−1(ω1n)0(ω1n)1…(ω1n)n−2(ω1n)n−1⋮⋱⋮(ωn−2n)0(ωn−2n)1…(ωn−2n)n−2(ωn−2n)n−1(ωn−1n)0(ωn−1n)1…(ωn−1n)n−2(ωn−1n)n−1⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦×⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣a0a1⋮an−2an−1⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦=⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣A(ω0n)A(ω1n)⋮A(ωn−2n)A(ωn−1n)⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦
对应上面的:“我们实际上是把一个单位根矩阵 M=[ωij] 左乘在了系数向量上”
只不过我们用了一些比较魔幻的方法将它加速到了 O(nlogn) awa。
为了方便表示,我们设:
Mω=⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣(ω0n)0(ω0n)1…(ω0n)n−2(ω0n)n−1(ω1n)0(ω1n)1…(ω1n)n−2(ω1n)n−1⋮⋱⋮(ωn−2n)0(ωn−2n)1…(ωn−2n)n−2(ωn−2n)n−1(ωn−1n)0(ωn−1n)1…(ωn−1n)n−2(ωn−1n)n−1⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦
Ma=⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣a0a1⋮an−2an−1⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦
MA=⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣A(ω0n)A(ω1n)⋮A(ωn−2n)A(ωn−1n)⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦
相应的,Mc 和 MC 同理
然后现在的情况是我们通过某些奇妙的方法求出了 MC ,并且我们已经知道了 Mω 然后现在我们就是要求出 Mc 。
怎么求呢 qwq?
可以发现其实就是在 MC 上乘上 Mω 的逆矩阵(不知道为什么的可以先去学一下矩阵相关的芝士)
我们先设 M′ω 为 Mω 的逆矩阵。
那么 M′ω 是什么呢?
我们先来浅猜一手:
M′ω=[ω−ij]
(就是把 Mω 的每一项的次数都 ×(−1) )
然后我们发现 Mω×M′ω=n×I0
( I0 就是单位矩阵,即除了对角线是 1 ,其他都是 0 该不会有人不知道吧awa )
所以我们就找到了 Mω 的逆矩阵 M′ω=[ω−ijn]
第三步解……
Q:哎!不对!!还没解决呢!!!
Q:你这样进行矩阵乘法,不就是 O(n2) 的了吗?
Q:你这复杂的有问题啊。
不不不,我们可以发现,其实 IDFT 就是把 DFT 里面的 Mω 换成了 M′ω , Ma 换成了 MC
其他的没有任何的改变。
所以可以如法炮制,直接向上面一样 O(nlogn) 解决 awa
所以第三步(IDFT)解决.
所以,我宣布,FFT,搞定 awa !
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】博客园携手 AI 驱动开发工具商 Chat2DB 推出联合终身会员
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 一个适用于 .NET 的开源整洁架构项目模板
· API 风格选对了,文档写好了,项目就成功了一半!
· 【开源】C#上位机必备高效数据转换助手
· .NET 9.0 使用 Vulkan API 编写跨平台图形应用
· MyBatis中的 10 个宝藏技巧!