题解看不到的神奇优化,好写好记(确信)的代码和公式,尽在此篇!
update 2024.3.22 添加了几种之前没有提到的 FFT 优化方式。
update 2024.3.25 添加了关于 DIF-FFT 的说明和使用。
update 2024.3.28 去除了错误的描述,添加了分裂基和效率对比。
update 2024.11.6 更正了 FFT 的英文名。
KARATSUBA 分治乘法
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
(注:这一部分并不实用,只是前置知识少,门槛低。)
设两数的位数为 nnn,则正常情况下进行乘法需要把 xxx 的每一位分别乘以 yyy 的每一位,故时间复杂度为 O(n2)O(n^2)O(n2),运行慢,我们可以利用分治思想优化(下面的所有形如 n/2n/2n/2 的数均向下取整):
设 x=a×10n/2+bx = a \times 10^{n/2} + bx=a×10n/2+b,y=c×10n/2+dy = c \times 10^{n/2} + dy=c×10n/2+d,则:
x×y=ac×10n+(ad+bc)×10n/2+bdx \times y = ac \times 10^{n} + (ad+bc) \times 10^{n/2} + bdx×y=ac×10n+(ad+bc)×10n/2+bd。
所以我们只需要处理 acacac、bdbdbd、(a+b)(c+d)(a+b)(c+d)(a+b)(c+d) 三组长度为 n/2n/2n/2 的数相乘,总计算量为 3T(n/2)+O(n)3T(n/2) + O(n)3T(n/2)+O(n),继续用同样的方法计算 T(n/2)T(n/2)T(n/2),就能优化更多,设数据规模为 nnn 时运算量为 T(n)T(n)T(n),则
T(n)=3T(n/2)+O(n)T(n) = 3T(n/2) + O(n) T(n)=3T(n/2)+O(n)
由主定理分析可知复杂度为 O(nlog23)O(n^{\log_2 3})O(nlog2 3)。
参考代码(丑陋):
优化
如果你试过运行上述代码,就会发现它的效率实际很低,这是因为在实际实现中,大量的拆分使这份代码带上了巨大的常数,一个简单有效的优化就是当 x,yx,yx,y 的位数小于一定值时就直接改用暴力算法,而不是直接拆到不能再拆。只需将
改为
即可。
还没结束!我们还可以采取一个优化:long long 型的储存范围为 [−263,263)[-2^{63},2^{63})[−263,263),只存 0∼90 \sim 90∼9 位太可惜了,我们可以用一个元素存 10810^8108 位!只需将上述代码中 basebasebase 的值改为 10810^8108 即可。
至此,我们的代码已经可以轻松通过 P1919:
快速傅里叶变换(FAST FOURIER TRANSFORM,FFT)
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
前置知识
* nnn:一个常数,具体值与实现该算法没有直接联系。
* 复数:指形如 a+bia+bia+bi 的数,其中 i=−1i=\sqrt{-1}i=−1 ,a,ba,ba,b 为实数。
* 单位根:如果有一个数 ω\omegaω,满足 ωn=1\omega^n=1ωn=1,称 ω\omegaω 为 nnn 次单位根。
* e2πine^{\frac{2\pi i}{n}}en2πi 是一个 nnn 次单位根,而且这个单位根称为主 nnn 次单位根或本原单位根,下文大部分的 ωn\omega_nωn 都指主 nnn 次单位根。
* 主 nnn 次单位根的性质:
* ωnn=1\omega_n^n=1ωnn =1(定义)
* ωnk+n/2=−ωnk\omega_n^{k+n/2}=-\omega_n^kωnk+n/2 =−ωnk (上两条的推论)
* 多项式的点值表示:对于一个一元 nnn 次多项式 f(x)f(x)f(x),如果能知道它在 xxx 为 n+1n+1n+1 个不同的值时的取值,就可以确定唯一一个多项式。
* 显然,在点值表示法下,求两个多项式乘积只需要把相同 xxx 个 xxx 值时多项式的值分别相乘即可。
正文
(为了可读性,本篇大部分代码均按照公式书写,而不进行过度卡常或依据编程语言的规则进行优化,a=a+b 等类似写法也会频繁出现。)
快速傅里叶变换可以快速(O(nlogn)O(n\log n)O(nlogn))计算多项式在 xxx 为 nnn 个单位根的幂时的值,以此做到快速将多项式转化为点值表示来快速(O(n)O(n)O(n))相乘。
令 f(x)=a0+a1x+a2x2+⋯+an−1xn−1f(x)=a_0+a_1x+a_2x^2+\dots+a_{n-1}x^{n-1}f(x)=a0 +a1 x+a2 x2+⋯+an−1 xn−1,其中 nnn 为 222 的整数次幂(原因等下会讲),则:
f(ωnk)=a0+a1ωnk+a2ωn2k+⋯+an−1ωn(n−1)kf(\omega_n^k)=a_0+a_1\omega_n^k+a_2\omega_n^{2k}+\dots+a_{n-1}\omega_n^{(n-1)k} f(ωnk )=a0 +a1 ωnk +a2 ωn2k +⋯+an−1 ωn(n−1)k
令:
f[0](x)=a0+a2x+a4x2+⋯+an−2xn/2−1f^{[0]}(x)=a_0+a_2x+a_4x^2+\dots+a_{n-2}x^{n/2-1} f[0](x)=a0 +a2 x+a4 x2+⋯+an−2 xn/2−1
f[1](x)=a1+a3x+a5x2+⋯+an−1xn/2−1f^{[1]}(x)=a_1+a_3x+a_5x^2+\dots+a_{n-1}x^{n/2-1} f[1](x)=a1 +a3 x+a5 x2+⋯+an−1 xn/2−1
则:
f(ωnk)=f[0](ωn/2k)+ωnkf[1](ωn/2k)f(\omega_n^k)=f^{[0]}(\omega_{n/2}^k)+\omega_n^k f^{[1]}(\omega_{n/2}^k) f(ωnk )=f[0](ωn/2k )+ωnk f[1](ωn/2k )
f(ωnk+n/2)=f[0](ωn/2k)−ωnkf[1](ωn/2k)f(\omega_n^{k+n/2})=f^{[0]}(\omega_{n/2}^k)-\omega_n^k f^{[1]}(\omega_{n/2}^k) f(ωnk+n/2 )=f[0](ωn/2k )−ωnk f[1](ωn/2k )
这样我们只需要代入一半的单位根的幂就可以顺带求出另一半的点值啦!对于 f[0]f^{[0]}f[0] 和 f[1]f^{[1]}f[1] 显然可以递归地继续用同样的方法,因为通过一半求出另一半需要每次都严格的将多项式分成相等长度的两部分,所以 nnn 必须为 222 的整次幂。
蝶形变换
实现中,反复地分多项式是很费时间的,我们可以在算法开始前就把每个系数放到最后到达的位置,再递归计算,这样就避免了反复复制。
下面这段解释取自 OI wiki。
> 规律:其实就是原来的那个序列,每个数用二进制表示,然后把二进制翻转对称一下,就是最终那个位置的下标,我们称这个变换为位逆序置换(bit-reversal permutation)。
很多题解在说完这一点后就直接得出了迭代 FFT,但实际上用这种方法优化完后递归版一样跑的很快。
代码
优化
* 众所周知:(a+bi)2=a2−b2+2abi(a+bi)^2=a^2-b^2+2abi(a+bi)2=a2−b2+2abi
故我们可以将两个多项式分别放在同一多项式的实部和虚部,FFT 后再计算结果的平方,虚部的一半就是答案。
* STL 的 complex 在不开 O2 的情况下常数略大于手写。
* 模板递归!将多项式长度作为模板参数传入,使编译器可以在编译阶段确定其值并内联展开。
参考代码:
频率抽取(DIF)
刚才的 FFT 是按照下标奇偶性分的,其实也可以按前后分,令 k′=n−kk'=n-kk′=n−k 为偶数:
f(ωnk)=f[0](ωn/2k)+ωnkf[1](ωn/2k)f(\omega_n^k)=f^{[0]}(\omega_{n/2}^k)+\omega_n^k f^{[1]}(\omega_{n/2}^k) f(ωnk )=f[0](ωn/2k )+ωnk f[1](ωn/2k )
f(ωnk+n/2)=f[0](ωn/2k)−ωnkf[1](ωn/2k)f(\omega_n^{k+n/2})=f^{[0]}(\omega_{n/2}^k)-\omega_n^k f^{[1]}(\omega_{n/2}^k) f(ωnk+n/2 )=f[0](ωn/2k )−ωnk f[1](ωn/2k )
刚才的 FFT(时间抽取,DIT)是先递归计算,再进行合并,而现在,我们可以先直接计算每组 f[0](ωn/2k)f^{[0]}(\omega_{n/2}^k)f[0](ωn/2k ) 和 f[1](ωn/2k)f^{[1]}(\omega_{n/2}^k)f[1](ωn/2k ) 的值,再递归地计算两部分,稍微画蝶形图分析一下,我们这样计算可以直接传入原数组,而得到的结果是二进制逆序的,参考代码:
这一方法的使用将在快速傅里叶逆变换中提到。
分裂基-DIT
以上的 FFT 是单次将多项式分成 222 份的,我们可以称之为基 222 FFT,相应的,还有基 444,基 888 甚至基 161616 等多种实现方法。
而分裂基 FFT 就是其中一种,它是将基 222 FFT 中下标为奇数的部分再次拆分(似乎可以算是不太严格的基 444 FFT?),即在:
f(ωnk)=f[0](ωn/2k)+ωnkf[1](ωn/2k)f(\omega_n^k)=f^{[0]}(\omega_{n/2}^k)+\omega_n^k f^{[1]}(\omega_{n/2}^k) f(ωnk )=f[0](ωn/2k )+ωnk f[1](ωn/2k )
的基础上,再用同种方法使:
f[1](x)=f[2](x2)+xf[3](x2)f^{[1]}(x)=f^{[2]}(x^2)+x f^{[3]}(x^2) f[1](x)=f[2](x2)+xf[3](x2)
代入原式,得到:
f(ωnk)=f[0](ωn/2k)+ωnkf[2](ωn/4k)+ωn3kf[3](ωn/4k)f(\omega_n^k)=f^{[0]}(\omega_{n/2}^k)+\omega_n^k f^{[2]}(\omega_{n/4}^k)+\omega_n^{3k} f^{[3]}(\omega_{n/4}^k) f(ωnk )=f[0](ωn/2k )+ωnk f[2](ωn/4k )+ωn3k f[3](ωn/4k )
这样分,同时兼具基 222 FFT 的灵活性和单次拆分多部分的高效率。
参考代码:
分裂基-DIF
(为了公式,代码与介绍一致且便于叙述,下面称原下标偶数部分为左半,奇数部分为右半。)
思想同上,还是把 f(ωnk)f(\omega_n^k)f(ωnk ) 写成一个基 222 FFT 加两个基 444 FFT 的形式即可。
这里需要注意:分裂基 DIT 在递归完后,因为左右都处理完了,而右半采取的是基 444 FFT 的拆分,所以用类似基 222 FFT 的方式合并,而分裂基 DIF 在拆分之前,因为左半和右半都需要继续递归,所以左半按基 222 拆,右半按基 444 拆。
参考代码:
快速数论变换(NTT)
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
前置知识
* FFT 的思想。
* OI wiki-数论基础中关于同余的部分。
* 对于质数 ppp 的原根 ggg,gp−1ng^{\frac{p-1}{n}}gnp−1 在模意义下与 nnn 次单位根有相同的性质。
正文
对于质数 ppp 的原根 ggg,gp−1ng^{\frac{p-1}{n}}gnp−1 在模意义下与 nnn 次单位根有相同的性质,所以可以用原根的幂代替复数单位根代入进行模意义下的 FFT,以此规避浮点数计算并加快速度。
为了 p−1n\frac{p-1}{n}np−1 为整数,我们的 nnn 应取形如 2k2^k2k 的数,ppp 应取形如 c⋅2k+1c\cdot2^k+1c⋅2k+1 的质数,比如 998244353=119⋅223+1998244353=119\cdot2^{23}+1998244353=119⋅223+1,1004535809=479⋅221+11004535809=479\cdot2^{21}+11004535809=479⋅221+1 等,这两个的原根都是 333。
如果选择 998244353998244353998244353 作模数,则多项式长度不能超过 2232^{23}223。
如果原多项式的各项的系数不超过 10910^9109,长度为 10510^5105,则相乘后各项系数不超过 102310^{23}1023,注意开 long long 和快速乘,同时选择合适的模数。
优化
* 没有了浮点数误差,我们可以选一个大模数然后压位,实测选择 998244353998244353998244353 时可以压 999 位。
* 两个不超过 p−1p-1p−1 的数进行加减,结果在区间 [0,2p)[0,2p)[0,2p) 中,我们可以用判断和加减代替部分模运算。
* 预处理原根的幂并用快速幂优化。
* 同时对两个多项式 NTT/FFT。
代码
快速傅里叶逆变换/快速数论逆变换(IFFT/INTT)
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
笔者才疏学浅,这里直接给出结论:用相乘后的各个点值作系数构成新多项式,用单位根的倒数(或原根的逆元)代替单位根(或原根)进行 FFT/NTT,得到的各项点值再除以多项式的长度(或乘以长度的逆元)就是各项系数。
DIF 和 DIT 混合使用加速
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
DIF 需要最后进行二进制逆序置换,DIT 需要开始时二进制逆序置换,我们用 DIF 的 FFT,DIT 的 IFFT 就相当于两个都不需要二进制逆序置换了,代码又少写一点。
极端优化
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
前排提醒,这一部分的内容实用性远低于前面各部分,且不给出任何形式的实现,仅供大佬们参考。
* (OI 禁止使用)使用 AVX2,SIMD 等指令集。
* 在多项式长度较短时改用迭代 FFT。
* 人工内联展开多项式较短的情况(多写几百行且效率提升低)。
* 预处理甚至打表单位根和其幂。
效率对比
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
(以 洛谷 P3803 为测试题目,均采用快读和三次变两次优化。)
省流:纯递归最慢,迭代次之,模板递归最快,DIF 和 DIT 混用优化和分裂基都很有较大提升,总计时间较最慢的写法(DIT,递归,基 222)减少了约 40%40\%40%。
参考文献
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
[1] 从多项式乘法到快速傅里叶变换
[2] 快速傅里叶变换(FFT)详解
[3] Split-radix FFT
[4] OI wiki - 快速傅里叶变换
[5] NTT 与多项式乘法
特别鸣谢
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
所有参考文献作者及网站维护人员。
我的同学,他提供给了我部分数学知识。
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
> 此页面来源洛谷!!
> 原页面:https://www.luogu.com.cn/article/rj58c2ebfor 洛谷
> 作者:zhangbo1000
> 本文作者未含有权限
> 此文为AI-Deepseek转换markdown格式