初步要求 学有余力 / 集训后继续消化 Learn for fun :) 记 使用 Iverson 括号约定:设 多项式的规模定义为多项式的次数加一.特别的,零多项式的规模为 给定两个至多 代入任意 点值意义下的多项式乘法是 若计算至多 朴素计算任意指定 Lagrange[1] 插值给出了 能否选取 离散傅里叶变换(Discrete Fourier Transform, DFT)接受一个至多 得益于单位根的特殊运算性质,二者均有被称为快速傅里叶变换(Fast Fourier Transform, FFT)的快速算法. 复数域上的 所有单位根模长均为 Remark (Euler’s formula for nerds). Euler 公式 Theorem 1 (消去引理) Theorem 2 (折半引理) 消去 / 折半引理将在 FFT 的推导中使用. Theorem 3 (求和引理) 求和引理的证明使用了等比数列求和公式.将在 IDFT 的推导中用到. 考虑将至多 记 由单位根的消去引理可证,DFT 矩阵 故快速计算 IDFT 的方法与 FFT 几乎一致,只需将计算 DFT 时使用的本原单位根 怎么计算卷积? 怎么快速求值? 怎么快速插值? 方便起见,我们只处理 Drawbacks? 迭代地实现 FFT 不仅在常数上更加优秀,亦更便于使用 C++ 的容器进行封装.这并不困难,只需自底向上模拟 FFT 递归过程即可. 唯一的问题——最底层的顺序? 来观察一轮 你发现了什么? 在计算点值前, 我们有 请自行实现更易用的容器封装版本. FFT 的缺点?浮点数带来的大常数与精度问题. 我们指出,在系数和点值模 我们指出,对于满足 只需修改单位根定义,把复数运算改为整数取模,就得到了能算 NTT 的 FFT 的实现. 关于 DFT 关于 NTT 十进制数可拆解为多项式表示,计算卷积后处理进位即可.由于数字最大只是 两个背包的合并就是多项式卷积. 反转多项式的系数数组再做卷积,可以快速得到两个多项式滑动窗口式的内积. 有些位运算可以写成卷积的形式.模 通过巧妙设定字符串距离函数,FFT 可解决更广泛的字符串匹配问题. 当值域较小时,将待计算的值放在多项式次数上统计贡献次数,可以绕开某些极难求解的数值问题.例如 Vandermonde 行列式的快速计算. 对多个长度相同的多项式的卷积,分治地卷起来可降低时间复杂度.长度不一时,挑小的先卷也可减小常数(用堆维护). 另有一种 CDQ 风格的分治 FFT.CDQ 长于处理带偏序的二元点对贡献,在处理形如 当卷积的前后项存在依赖关系时,也可使用此法保证处理顺序恰当.然而此类依赖问题往往也可通过解生成函数方程的方法求得封闭形式. 时间复杂度均为 生成函数是一种对数列的操作技巧.通过将数列表示为多项式或形式幂级数,数列间复杂的和式操作可用简单函数的乘法、复合等运算进行表示,从而大大降低了数列变换技巧的使用门槛. 生成函数在组合数学中应用广泛,且生成函数的部分操作在组合意义下也有较为直观的理解.本节将带大家初窥其中的奥妙.限于篇幅和主讲人能力,我们仅以题带点地讲解,期冀为大家建立构造生成函数的直觉.请感兴趣的同学下来做进一步研究. 熟悉 Taylor 展开的同学或能较快上手此部分内容. 序列 Exercise 1 写出下列数列的 OGF.下标从 OGF 相乘,是背包,是卷积,是两块无标号组合对象的有序拼接. Exercise 2 写出下列计数问题的 OGF,均以 序列 Exercise 3 写出下列序列的 EGF,下标从 观察两个 EGF 的乘积 Exercise 4 写出下列计数问题的 EGF,均以 本节介绍 NTT 的原理. FFT 加速卷积算法的核心,一是多项式的求值插值原理,二是单位根带来的分治快速算法.我们将在本节中证明,模 NTT 原理涉及原根等数论内容.本讲的目标是建立 DFT 变换和 FFT 算法的通用数学框架,而非具体研究其某一特例.故我们只讲解 NTT 所需的基础数论知识,无关的细节则略过处理.对数论感兴趣的同学可前往 OI Wiki 学习. 模 模 模 之后记系数均在 Lemma 1 ( Theorem 4 (Lagrange 定理) 设 Corollary 1 设 推论告诉我们,欲确定 回顾 DFT 中复数域 事实上所有的 我们把这一类重要的单位根称为本原单位根.抽象的来说, 如何定义 称在模 Theorem 5 (Fermat 小定理) 若 Theorem 6 (Euler 定理) 若 简化剩余系的大小即 设 补充讨论 当 当 需要注意的是,Euler 定理只给出了 称 Theorem 7 (原根存在定理) 下面再介绍原根的两个定理. Theorem 8 (原根判定定理) 若 必要性显然.充分性,反证出所有 Theorem 9 (原根个数定理) 若 设 使用数论中 Bézout 定理, 遗憾的是,一般的 求和引理是保障 NTT 逆变换对应矩阵确为 为保证 综上,逆变换矩阵确为 为满足 FFT 对 结合前述关于 当 我们找到一篇有关该问题的参考文献[7, Secs. 3 and appendix B],但尚不确定其证明的正确性.友情提示读者: 后续数学内容导读 NTT 原理虽已非常”数学”,但也只是 DFT 在有限域上的一个实例.本节往后,我们要尝试为多项式系数位于复数域 后续数学内容不再要求掌握.望同学们在纷繁的定理定义中抓住要旨,窥见抽象数学背后蕴藏的规律.熟悉高等代数和抽象代数的同学或会对某些内容感到熟悉.抽象地讨论 FFT 的资料并不多见,后续内容多为主讲人的新进探索,或有谬误,望不吝指正. 本节将重新审视已经熟知的多项式,把抽象的、代数的多项式和具体的、分析的多项式函数区分开来.我们指出,多项式和多项式函数不同但关联紧密,形式幂级数与幂级数亦有此类联系.这些抽象的讨论将帮助我们剖析多项式求值插值的基本原理. 除常见代数书目(如[8]),也推荐参考 OI Wiki 的多项式基础[9]和 Wikipedia 的形式幂级数[10]. 设(无穷)序列 多项式环上的加法、乘法的定义已经为大家所熟知.系数所处的环保证了多项式加法和乘法的良定义,而在这两种运算下, Remark (群、环、域). 群、环、域都是常见的代数结构,其中的元素在给定运算下封闭,并满足特定的运算性质.简单来说,环13上定义了加法和可能不可逆的、不一定交换的乘法,域上定义了加减乘除所有四则运算.交换环中的乘法满足交换律.除环中的所有元素有乘法逆元.域是交换除环. 习惯上也会将多项式 多项式 下面额外为多项式定义一种新的运算.多项式 形式幂级数定义与多项式的唯一区别是其不要求 由于涉及无限次运算,形式幂级数的复合运算需考虑环 Remark. DFT, NTT 与多项式环 DFT/FFT 加速的多项式卷积在复数域 整环是无零因子的交换幺环.所谓无零因子,即环中任意元素 定义有带余除法的环被称为 Euclid 整环.域上的多项式环都是 Euclid 整环[8, 第 7 章第 2 节定理 3, p. 11].值得注意的是, 如在带余除法中保证除式是首一多项式17,则带余除法的良定义和进行过程也均可在整环上实现. 刚刚强调,多项式 / 形式幂级数只是定义了加法和乘法的序列.现在介绍多项式函数和幂级数.它们不是序列,而是映射18. 多项式 与多项式 / 形式幂级数不同,这里的 多项式函数的加法和乘法定义为函数的加法和乘法,即 多项式和多项式函数似乎在许多情况下有着平行的关系.下面介绍一个较直观的结论. 若环 多项式和多项式函数的这一关系为多项式在任意点的求值操作提供了理论基础.在多项式函数的 需要注意的是,前述结论的逆命题不一定成立,即环 该逆命题的本质是通过多项式函数的所有函数值反过来确定多项式(系数)的过程.而如果在确定时只使用一部分函数值,就是所谓的多点插值过程.很多时候待求多项式的次数是已知的,这在相当程度上缩小了待定多项式的范围.我们指出,只要 证明的关键是用带余除法24讨论多项式根与其一次因式的关系,即多项式余式定理或小 Bézout 定理[8, 第 7 章第 6 节定理 6, p. 35]. 需要强调,上述结论只能说明,只有那些确实可通过多项式 到这里,多项式的求值与插值的线性表示已经呼之欲出了.若将交换环 若将对 本节的核心是多项式和多项式函数的区别与联系,两个方向的”确定”分别给出了多项式多点求值和多点插值的理论基础. 由于多项式和多项式函数这种若即若离的关系,往往在记号上也有意无意地混淆了它们,某些情况下加大了区分的难度.本篇使用的记号体系将尽量用单个字母 前面讨论了在多项式任意点处求值插值的基本原理,但 DFT/FFT 的运行只需在单位根处求值和插值.本节将进一步放宽对多项式环的限制,介绍定义在有主要单位根的环上的一般的 DFT 及其快速算法 FFT. 本节内容主要参考了 [12, Secs. 2, pp. 983–984] 和 [13]. 定义环 定义环 本原单位根在许多情况下与主要单位根等价,但亦非完全相同. Proposition 1 若环 Proof. 对任意正整数 Proposition 2 若环 Proof. 反证.若存在一正整数 分别根据定义和数论中的 Bézout 定理,容易证明主要单位根和本原单位根的消去引理:设 注意到 环上的 设 除无法在任意环上使用 设 考虑对模 注意到 因此,令 上一节中,我们建立了在有主要单位根的环上的 DFT 及其快速算法 FFT 的相关理论,但由于放宽了对多项式环 在求值与插值部分已经介绍求值变换 设 设 令 前述讨论已经证明29 我们指出,当要求所作变换可逆时,卷积定理反过来也要求所作变换是一类似 DFT 映射的变换30. Exercise 5 尝试将任意序列 DFT 两次,观察结果.证明你的结论. Exercise 6 求 DFT 矩阵的行列式.尽可能缩小可行解范围. 有哪些? 咋推的? 有啥用? 怎么讲? 给定一多项式 多项式逆元存在的充分必要条件是其常数项非零(这是因为边界条件 以下简记 两边平方得 给定一多项式 次数为 推导是容易的.方程两侧同时求导得 感谢 感谢队友 主讲人练题少,仅供参考. SPOJ-TSUM Triple Sums BZOJ3513-MUTC2013 Idiots 上面两道题都是 OGF 消序,较 EGF 消序困难.一般的方法是使用 Polya 计数原理. ABC291G OR Sum 也是滤波器的应用. 百度之星 2023 初赛第二场 T8 容斥后需要计算若干一次多项式乘积,分治 NTT 即可. 求 洛谷 P4721 【模板】分治 FFT CDQ 偏序化处理前后项依赖.也可解生成函数方程再多项式求逆. 一种常见的证法是使用 Vandermonde 行列式证明矩阵可逆.后面会介绍多项式环风格的证明.↩︎ 有的文献定义 这也表明适当归一化后的 DFT 矩阵是一个酉矩阵.↩︎ NTT 原理需较多笔墨,稍后介绍.↩︎ 否则只有 多项式求逆等多项式进阶操作,我们后续讲解.↩︎ 对质数 这是后文所述定理“整环上的本原单位根也是主要单位根”在一般环上的一个反例.↩︎ 网传此模数由 UOJ 站长 vfleaking 提出并推广.在所有需要取模的题目中使用该模数,可使选手无法通过模数判断题目的做法.↩︎ 本篇中环的定义包含乘法单位元,即幺环.↩︎ 无零因子的交换幺环,稍后介绍.↩︎ 可用前述多项式乘积次数公式证明.↩︎ 试试在 首项为 函数和映射几乎是等价名词.有时函数特指值域包含于复数域 为良定义 这里再次涉及环 这种保持结构不变的映射被称为同态(homomorphism).↩︎ 其证明了域上一元多项式环的通用性质.仿照该证明应可证明环上的版本,从而证明这一同态关系.↩︎ 该定理是下方高亮定理的一个自然的推论.↩︎ 由于一次因式均为首一多项式,可以在整环上对其使用带余除法.↩︎ 模是定义在环上的“线性空间”.↩︎ Lagrange 插值的构造用到了除法,且行列式非零推出矩阵可逆仅在域上的线性空间中适用,因此必须要求 使得 证明使用主要单位根的定义(求和引理)即可.↩︎ 虽然前面用到了 具体来说,该变换只能是 DFT 矩阵的某个行置换.在 1 Forewords
卷积,但不止卷积 - FFT 漫谈
推荐食用方法
记号说明
2 FFT/NTT in a nutshell
2.1 FFT
多项式卷积
系数 - 点值 - 系数
系数 - 点值 - 系数 - 快速转换?
Discrete Fourier Transform
复数域单位根
复数域单位根 - 三个重要性质3
复数域单位根 - 三个重要性质
Fast Fourier Transform
DFT 的矩阵表示
IDFT
梳理
FFT 递归实现 - DFT 部分
#include<bits/stdc++.h>
#include<complex>
using namespace std;
typedef long long ll; typedef complex<double> CP;
const ll MXN=4E6+5; const double PI=3.14159265358979323846l;
CP tmp[MXN];
void _DFT(CP A[],ll n,ll typ){
n/=2; if(n==0) return;
for(ll k=0;k<n;k++) tmp[k]=A[2*k],tmp[n+k]=A[2*k+1];
for(ll k=0;k<2*n;k++) A[k]=tmp[k];
_DFT(A,n,typ); _DFT(A+n,n,typ);
CP w(cos(2*PI/(2*n)),typ*sin(2*PI/(2*n))), wk=1;
for(ll k=0;k<n;k++){
tmp[ k]=A[k]+wk*A[n+k];
tmp[n+k]=A[k]-wk*A[n+k];
wk*=w;
} for(ll k=0;k<2*n;k++) A[k]=tmp[k];
} void DFT(CP A[],ll n,ll typ){
_DFT(A,n,typ); if(typ==-1) for(ll i=0;i<n;i++) A[i]*=1.0l/n;
}
FFT 递归实现 - 卷积部分
// alternatively, use std::__lg() in GCC
ll log2ceil(ll n){ll cnt=0; for(ll t=1;t<n;t<<=1) cnt++; return cnt;}
CP A[MXN],B[MXN],C[MXN]; ll outC[MXN];
ll* conv(ll inA[],ll aN,ll inB[],ll bN){
ll n=1LL<<log2ceil(aN+bN-1);
for(ll i=0;i<aN;i++) A[i]=inA[i];
for(ll i=0;i<bN;i++) B[i]=inB[i];
DFT(A,n,1); DFT(B,n,1);
for(ll i=0;i<n;i++) C[i]=A[i]*B[i];
DFT(C,n,-1); for(ll i=0;i<n;i++) outC[i]=round(C[i].real());
return outC;
}
FFT 迭代
蝶形运算
void spawnrev(ll n){
rev[0]=0;
for(ll i=1;i<(1<<n);i++)
rev[i]=(rev[i>>1]>>1)+((i&1)<<(n-1));
}
FFT 迭代实现
void DFT(CP A[],ll n,ll typ){ // rev[] should be spawned in advance
for(ll i=0;i<n;i++) if(i<rev[i]) swap(A[i],A[rev[i]]); // a one-to-one permutation
for(ll hf=1;hf<n;hf*=2){
CP w(cos(2*PI/(2*n)),typ*sin(2*PI/(2*n))), wk=1;
for(ll i=0;i<n;i+=hf*2){
CP wk=1;
for(ll k=0;k<hf;k++){
CP x=A[i+k],y=wk*A[i+hf+k];
A[i+k]=x+y; A[i+hf+k]=x-y;
wk=wk*w;
}
}
}
if(typ==-1) for(ll i=0;i<n;i++) A[i]*=1.0l/n;
}
2.2 NTT
NTT 速成5
const PR=3,MOD=998244353;
ll w=qpow(PR,(MOD-1)/(hf*2)); if(typ==-1) w=inv(w);
FFT/NTT in a nutshell - 小结:概念区分
3 Applications
3.1 基本应用
基本应用
基本应用 - 分治 FFT
3.2 生成函数初步
生成函数初步 - 导言
Ordinary Generating Function
OGF 组合意义
Exponential Generating Function
EGF 组合意义
EGF 组合意义
4 Mathematics behind
4.1 NTT 原理
NTT 原理 - 导言
模
本原单位根
Euler 定理的证明 - 简化剩余系
原根10
求阶和原根
int
型中的单次加减操作不会溢出,是 OI/XCPC 计数题中不可多得的优秀模数12.从分析到代数
4.2 求值与插值
求值与插值 - 导言
多项式
多项式
形式幂级数
带余除法
多项式函数 / 幂级数
求值
插值
求值与插值的线性表示
求值与插值 - 小结
4.3 环上的 DFT
环上的 DFT - 导言
环上的单位根 - 两个定义
环上的单位根 - 区别与联系
环上的单位根 - 区别与联系
环上的单位根 - 其它性质
环上的 DFT
环上的 FFT
环上的 FFT
4.4 循环卷积与卷积定理
循环卷积与卷积定理 - 导言
循环卷积
循环卷积
卷积定理
5 Advanced Operations
多项式全家桶 - 序言
多项式求逆
多项式求逆
多项式求逆 - 实现
Poly inv(const Poly &A){
ll n=A.len(); Poly B(1); B[0]=inv(A[0]);
for(ll hf=1;hf<n;hf<<=1){
B=B*2-B*B*A.subpoly(0,hf*2); B.resize(hf*2);
} B.resize(n);
return B;
}
多项式
多项式
Poly drv(Poly A){ // derivative
for(ll i=0;i<A.len();i++) A[i]=(i+1)*A[i+1]_;
A.pop_back(); return A;
}
Poly itg(Poly A,ll c){ // integral
A.push_back(0); for(ll i=A.len();i>=1;i--) A[i]=A[i-1]*inv(i)_;
A[0]=c; return A;
}
Poly ln(const Poly &A){
return itg((drv(A)*inv(A)).subpoly(0,A.len()-1),0/*log(A[0])*/);
}
Acknowledgements
keke_046
学长教授 FFT、集合幂级数与生成函数.微言大义,博大精深,至今仍在消化.ItzDesert
提供位运算典题一道并提供内容编排建议.题单
References
[1] OI-Wiki, “拉格朗日插值.” https://oi-wiki.org/math/numerical/lagrange/.
[2] 张筑生, “数学分析新讲(重排本)(第二册),” 2nd ed.北京: 北京大学出版社, 2021, pp. 256–262.
[3] 杨树森, “三角函数的严格定义.” https://zhuanlan.zhihu.com/p/58814328/, 2023.
[4] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, “算法导论(原书第三版),” 北京: 机械工业出版社, 2013.
[5] OI-Wiki, “离散对数.” https://oi-wiki.org/math/number-theory/discrete-logarithm/.
[6] OI-Wiki, “原根.” https://oi-wiki.org/math/number-theory/primitive-root/.
[7] R. Agarwal and C. Burrus, “Fast convolution using fermat number transforms with applications to digital filtering,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 22, no. 2, pp. 87–97, 1974, doi: 10.1109/TASSP.1974.1162555.
[8] 丘维声, “高等代数 下册,” 3rd ed.北京: 高等教育出版社, 2015.
[9] OI-Wiki, “多项式与生成函数简介.” https://oi-wiki.org/math/poly/intro/.
[10] Wikipedia, “Formal power series.” https://en.wikipedia.org/wiki/Formal_power_series.
[11] Wikipedia, “Lagrange polynomial.” https://en.wikipedia.org/wiki/Lagrange_polynomial.
[12] M. Fürer, “Faster integer multiplication,” SIAM Journal on Computing, vol. 39, no. 3, pp. 979–1005, 2009, doi: 10.1137/070711761.
[13] Wikipedia, “Discrete fourier transform over a ring.” https://en.wikipedia.org/wiki/Discrete_Fourier_transform_over_a_ring.
[14] I. Baraquin and N. Ratier, “Uniqueness of the discrete fourier transform,” Signal Processing, vol. 209, p. 109041, 2023.
[15] H. S. Wilf, “Hadamard determinants möbius functions, and the chromatic number of a graph,” 1968.
Footnotes