FFT&NTT学习笔记
参考资料
互动视频学习&实现FFT (建议配合视频阅读本章)
[TOC]
概述
什么是傅里叶变化？
通俗来说傅里叶变换用于将时域函数变换为频域函数
原来的时域函数可以用多个不同频率的sin or cos 函数叠加拟合出来
上图形象的描述了傅里叶变换在做什么
不过此处不多深究傅里叶变换数学上的推导及其含义(主要是👴忘了) 详细可以看这里
这里重点讨论 :1.FFT具体实现,2.它解决什么问题 3.它与信号处理有何关系
前置知识
key words : 单位根 , 多项式
单位根
可以理解为在复平面上的向量
一个单位圆上的复数 $z=cosθ+isinθ$
用欧拉定理 $e^{iθ}=cosθ+isinθ$ 其中θ为幅角 因为在单位圆上所以模长为1
而单位根$w_n^d=(e^{iθ})^d$ 其中 $θ=2pi/n$
可以想象成一个单位圆被平分成n分从x轴的 $w_n^0$逆时针数第d个点为$w_n^d$
关于单位根有三个重要的引理在之后实现FFT时会用到：
上面是严格按照定义证明，下面按照我的理解证明
Lemma.1 消去引理：$w_{dn}^{dk}=w_n^k$
可以想象被分成dn分的单位圆从x轴逆时针数第dk个 和 被分成n个的单位元从x轴逆时针数第k个 是同一个点
而这个引理主要用到这个公式的转换$w_n^{n/2}=w_2^1=-1$
Lemma.2 折半引理：$(w_n^{k+n/2})^2=(w_n^k)^2$
这也显然某个点转过半圈它们的实部长度相等 即 $|cos(θ)|=|cos(θ+pi/2)|$
Lemma.3求和引理:$∑_{i=0}^{n-1}(w_n^k)^i=[k==n]*n$
这个引理用于求IDFT
若将这个式子考虑为向量相加,单位圆上被平分成n+1个点,
当k=r*n时(r = any int) 等效于单位圆上只有$w_0^0=1+0i$ 它被加了n次所以求和结果为n
当k!=r*n时，无论分成奇数个还是偶数个点这几个向量相加还是0 如下图
多项式
多项式F的次数degree(F) = 最大幂次项的幂级数， 次数界>degree(F)
多项式的系数表示法:
$F(x)=a_0+a_1x+a_2x^2+….+a_nx^n=∑_{i=0}^n a_ix^i$ (F(x)的degree(F)=n)
一般可以用一个一维矩阵[a0,a1,a2,…an]仅仅保存系数来表示一个多项式
多项式的点值表示法:
顾名思义用一系列点(n+1个点)表示一个多项式函数图像
{$(x_0,F(x_0)),(x_1,F(x_1)),(x_2,F(x_2)),…,(x_n,F(x_n))$}
点值表示法有什么优势呢
假如求$A(x)=F(x)×G(x)$ 其中F(x)的系数以ai表示共m项 G(x)的系数以bi表示共n项
A(x)的系数以ci表示 共n+m项即得到的$A(x)$中degree(A)=degree(F)+degree(G)
用系数表示法$c_i=∑_{j=0}^{i} a_j*b_{i-j}$ 需要$On^2$
而用点值表示法我我们要求$A(x_i)$只用算$F(x_i)×G(x_i)$ 对应点相乘就好了
点值表示法可以$On$求出两多项式的乘积
如果有什么黑科技能快速将多项式从系数表示法变为点值表示法 求得两个多项式的乘积后再变回系数表示法
那么就能加速多项式乘法了,FFT算法解决了这个痛点
解决什么问题
FFT主要用于在$Onlogn$的复杂度内快速地 将多项式从系数表示法变为点值表示法 or 逆变换操作
由此来加速多项式乘法 比如模板题
给定一个 n 次多项式 F(x)，和一个 m 次多项式 G(x)。请求出 F(x) 和 G(x) 的卷积。
这里求两个多项式卷积实际上指的是两个多项式相乘的结果用系数表示法输出各个项的系数
假如求$A(x)=F(x)×G(x)$ 其中F(x)的系数以ai表示共m项 G(x)的系数以bi表示共n项
A(x)的系数以ci表示 共n+m项即得到的$A(x)$中degree(A)=degree(F)+degree(G)
用系数表示法求A(x)的系数$ci$需要进行卷积运算即：$c_i=∑_{j=0}^{i} a_j*b_{i-j}$
而解决这个问题我们所需要做的就是求A(x)
具体的求： $A(x)=IDFT(FFT(F)*FFT(G))$
其中FFT(X)表示将某多项式X从系数表示变为点值表示 IDFT(X)表示将某多项式X从点值表示变为系数表示
意思是我们需要将F(x),G(x)转化为点值式 然后$On$的对应点相乘 再转为系数表达式整体复杂度$Onlogn$
插入语
这里稍微提一下实现的细节问题，下文会具体贴出实现代码
这里需要注意的是由于degree(A)=degree(F)+degree(G)最后A(x)的系数有n+m个
于是我们直接读入长度为n+m系数矩阵a和b(原来没有的系数为0) 这样FFT和IDFT处理时就不会出错了
另外,在实现FFT时原来的数组含义为存第i项的系数
变为点值表达式后数组的含义为当多项式取$x=w_n^i$时对应的$G(w_n^i)$ 值
IDFT后数组又变成了第i项的系数
##FFT算法实现
key words:实现原理,优化加速
本章主要为了了解如何实现&使用FFT 因此推导过程略显潦草
实现原理
最终目的：使用FFT将$A(x)=∑_{i=0}^{n-1}a_ix^i$ 用点值法表示出来
而快速傅里叶变换使用单位根$w_n^i,i∈[0,n-1]$做为点代入求得点值表示法中的{$(x_0,y_0),..,(x_i,y_i)$} 集合
也就是说我们要求出$y_i=A(w_n^i)=∑_{j=0}^{n-1}a_jw_n^{ij}$ $i∈[0,n-1]$
对于多项式A(x) 它的系数构成的一维矩阵$a=[a_0,a_1,a_2….a_{n-1}]$
将这个系数矩阵按奇偶分组 对应系数下标表示此项为${a_ix^i}$
$a^{[0]}=[a_0,a_2,…a_{n-2}]$ 对应多项式为$A^{[0]}(x)$
$a^{[1]}=[a_1,a_3,…a_{n-1}]$对应多项式为$A^{[1]}(x)$
尝试寻找这三个矩阵(多项式)的关系
得到$A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)$ (中间的详细过请看鹤翔万里的互动视频)
分别带入单位根$w_n^k$ 和$w_n^{k+n/2}$得到下图两个式子
利用折半引理和消去引理得到
于是求解$A(x)$的问题被分成两个规模更小(只有A(x)一半系数)的子问题 我们只要求$A^{[0]}(x),A^{[1]}(x)$ 就能得到$A(x)$
同理利用递归分治将求解$A^{[0]}(x),A^{[1]}(x)$ 于是如何将多项式从系数表达式变为点值表达式的问题解决了
接下来考虑如何进行逆变换
以下截取自快速傅里叶变换详解(另外加了些注释讲解)
$(y_0,y_1,y_2,\dots,y_{n-1})$为多项式A(x)系数矩阵 $(a_0,a_1,a_2,\dots,a_{n-1})$的傅里叶变换（即点值表示）
设有另一个向量$(c_0,c_1,c_2,\dots,c_{n-1})$满足$c_k=∑_{i=0}^{n−1}y_i(ω_n^{−k})^i $
即多项式$B(x)=y_0,y_1x,y_2x^2,\dots,y_{n-1}x^{n-1}$在$\omega_n^{0},\omega_n^{-1},\omega_n^{-2},\dots,\omega_{n-1}^{-(n-1)}$处的点值表示
$(c_0,c_1,c_2,\dots,c_{n-1})$满足$c_k=∑_{i=0}^{n-1}y_i(w_n^{-k})^i$
这样就将点值表示法的矩阵$y$重新视为一个系数表示法的多项式系数
将$y$做一个替换
$c_k=∑_{i=0}^{n-1}(w_n^{-k})^i(∑_{j=0}^{n-1}a_j(w_n^i)^j)$
因为本来$yi$就是将$w_n^i$带入多项式A(x)的系数表达式$∑_{j=0}^{n-1}a_j(w_n^i)^j$求得的值
其中$aj$是我们最终将A(x)从点值表示变为系数表示所要求的系数
接着运算化简得到： $c_k=∑_{j=0}^{n−1}a_j(∑_{i=0}^{n−1}(ω_n^{j−k})^i) $
这里两个$∑$累加变量i 和 j里外互换是比较常见的操作
实际上可以理解成二维矩阵是一列一列加和还是一行一行加和 结果都是矩阵里所有数的和
设$S(x)=∑_{i=0}^{n-1}x^i$ 带入$w_n^k$ 通过Lemma.3求和引理可得当且仅当k=n时$S(w_n^k)$有值
应此$c_k=∑_{j=0}^{n−1}a_j(∑_{i=0}^{n−1}(ω_n^{j−k})^i) $被化简为$c_k=na_k$
其中$c_k$为点值表示法中矩阵中第k个元素 求得$a_k$为系数表示法中矩阵第k个元素
回头看IDFT的操作实际上是将A(x)的点值表示法视为B(x)系数表示法再进行FFT变换
需要修改的地方就是带入的点为{$(w_n^{-0},B(w_n^{-0})),(w_n^{-1},B(w_n^{-1})),…,(w_n^{-k},B(w_n^{-k})),…$}
并且得到B(x)的点值表示法矩阵ci需要 /n才等于要求的A(x)的系数表示法矩阵ai
至此解决了FFT和IDFT的推导&实现，现在能实现递归版的FFT&IDFT
1 | ///来自洛谷的模板 |
然而递归版的FFT仍然不够优秀
接下来要讲如何将递归变为迭代 提高效率
优化加速
主要有两种技巧优化加速
1.逆序位置换，2.蝴蝶操作
1.逆序位置换
不断递归时我们发现，最初的系数矩阵会被按照奇偶分组
递归到最后一层时当前系数下标与原数组的下标有一种规律
上一行是变换前的系数下标，下面为递归到最后的系数对应下标
经过观察得到规律：
于是我们可以舍弃递归直接使用位逆序置换操作预处理得到和递归一样的系数下标信息
根据公式
$$
A(w_n^k)=A^{[0]}(w_{n/2}^k)+ w_n^k A^{[1]}(w_{n/2}^k)
$$
$$
A(w_n^{k+n/2})=A^{[0]}(w_{n/2}^k)-w_n^kA^{[1]}(w_{n/2}^k)
$$
得到如下迭代运算过程
另外公式中$w_n^kA^{[1]}(w_{n/2}^k)$ 复数乘法被计算两次 可以将它计算一次后存入一个变量避免重复操作
这个优化被成为蝴蝶操作
于是迭代版的FFT也就出来了
1 | struct Complex{ |
还是之前的模板题模板题
我们要添加上逆序位置换 然后实现$A(x)=IDFT(FFT(F)*FFT(G))$
1 |
|
这里还有个小技巧三次FFT可以被优化成两次FFT同样得到结果
我们可以把$G(x)$放到$F(x)$的虚部上求出$F^2(x)$ ,然后取出$F(x)$的虚部全部/2就是答案了
正确性证明：$(a+bi)^2=(a^2-b^2)+(2abi)$
不过上面的代码中还有一些疑问比如为什么要+0.5上取整调整精度这是哪步产生的
到此为止FFT的实现优化已经结束。
接下来将要讨论NTT是什么以及FFT与信号处理的关系。
NTT
key words:原根,具体实现
由于FFT需要用到复数 使用double精度较差
而NTT使用原根替代FFT中的单位根进行点值转换 其余原理相同
原根
首先介绍一个概念”阶“
若a,p互素且p>1 对于$a^n≡1 (modp)$最小的n 我们称之为a模p的阶,记作$ δ_p(a)=n $
例如(似乎必须大于0)：
$δ_7(2)=3$，
$2^1≡2(mod7)$
$2^2≡4(mod7)$
$2^3≡1(mod7)$
然后就是要讨论的原根了
原根的定义：设p是正整数，a是整数，若$δ_p(a)$等于$ϕ(p)$，则称a为模p的一个原根
$δ_7(3)=6=ϕ(7)$，也就是说$3^6mod7=1$ 其中6是3模p的阶, 因此3是模7的一个原根
($ϕ(x)$是欧拉函数表示小于x且与x互素的数有几个不包括1,默认$ϕ(1)=1$)
注意原根的个数是不唯一的 如果模数p有原根，那么它一定有$ϕ(ϕ(p))$个原根
根据原根的定义我们得到 对于质数p,假若g是p的原根 则有$g^{ϕ(p)}=1(modp)$
因为$δ_p(g)=ϕ(p)$ ,进一步可以得到$g^imodP(1<g<P,0<i<P)$的结果两两不同
对于质数$p=k×2^N+1$,设其原根为g,我们令$g_n=g^{(p-1)/n} (modp)$ ,注意这里的$g_n$是定义出来的
在数值上等于$g^{(p-1)/n}$ 与g无关，在NTT中原根扮演了单位根的角色
$g_n$同样也满足单位根的三个性质 所以它能够替代单位根
证明参考此处
####细节
为什么质数p一定要取$p=k×2^N+1$
定义中$g_n=g^{(p-1)/n} (modp)$, $(p-1)/n$必须为整数
考虑FFT的过程是不断将系数数组二分的过程(递归时n>>1)
且在使用FFT前就要确保n(n就是系数数组的长度)为2的次幂
因而p-1要有足够大的2^N 才能保证(p-1)能被n整除 保证递归顺利进行
- 在逆变换时需要求 $g_n^{-1}$，在模意义下，也就是 $g_n$ 的逆元，预处理即可
p建议取998244353，它的原根为3。(注意读题有的题里会出现993244353这个数它的原根不是3)
如何求任意一个质数的原根呢？
对于质数p，质因子分解p−1，若$g^{(p−1)/pi}≠1(modp)$恒成立，g为p的原根
具体实现
1 |
|
FFT与信号处理的关系
参考资料
在上文讨论的FFT中 丝毫没有看见傅里叶变换最原始的用途：时域函数与频域函数的转换的影子
实际上多项式的点值表示法对应时域描述,而系数表示法则对应频域描述
回忆多项式$F(x)=∑a_kx^k$当带入点{$w_n^0, w_n^1, ..w_n^k…$}时得到了对应点{$(w_n^k ,F(w_n^k))$}
对于多项式而言$y_k=F(w_n^k)=∑_{j=0}^{n-1}a_jw_n^{kj}$ ，其中$w_n^k=(e^{2pi/n})^k$
而某函数$f(x)$的傅里叶级数展开$f(x)=∑_{n=-inf}^{inf}c_ke^{ikx}$
多项式带入一个自变量$x=w_n^k$ 就相当于函数的傅里叶级数表达式带入自变量$x=k*2pi/n$
当我们将{$w_n^0, w_n^1, ..w_n^k…$}带入多项式 时我们就得到了下图左侧傅里叶级数对应的函数值
也就是时域上一系列点值{$(w_n^k ,F(w_n^k))$}
傅里叶级数是将函数用不同的频率分量($cosθ+isinθ$ )×权值($c_i$)求和得到的一个时域表达式
(时域函数图像就是不同时刻(x)对应的信号强度(y) )
而傅里叶变换是通过 时域上的点值{$(w_n^k ,F(w_n^k))$} 求得 频域上不同频率分量的权值 $c_i$
逆变换是通过频域上的权值$c_i$求得 时域上的一系列点值
(频域函数图像就是不同频率分量(频率为x)的幅值(权值为y) )
相当于通过系数表达式 求得 点值表达式
只不过我们将多项式中的自变量替换成了$w_n^k$
这就相当于将多项式中的x当成不同频率分量而系数则成为权值