FFT,即快速傅里叶变换,是离散傅里叶变换的快速方法,可以在很低复杂度内解决多项式乘积的问题(两个序列的卷积)
卷积通俗来说就一个公式(本人觉得卷积不重要)
Ck=∑i+j=kAi∗Bi
那么这个表达式是啥意思了:
有两个序列A和B,其中A={A1,A2,...},B={B1,B2,...}
A序列有a个元素,B序列有b个元素。于是,由这两个序列可以推出另一个序列C,C序列如何确定了?确定方法就按照卷积的公式得来的,即C=⎧⎨⎩∑i+k=0Ai∗Bj,∑i+k=1Ai∗Bj,...,∑i+k=2a+2bAi∗Bj⎫⎬⎭
这里只介绍一下卷积,下面来探究卷积和多项式之间的关系。
形如
f(x)=a0+a1x+a2x2+...+anxn
的函数关系式叫做关于x的多项式,其中多项式系数可以构成一个序列,即
A={a0,a1,a2,...,an}
假如有两个多项式,其中
f(x)=a0+a1x+a2x2+...+anxn
g(x)=b0+b1x+b2x2+...+bmxm
现在要求f(x)*g(x)的序列,很明显
f(x)∗g(x)=∑i+j=0ai∗bj+∑i+j=1ai∗bjx+∑i+j=2ai∗bjx2+...+∑i+j=m+nai∗bjxm+n
现在把f(x),g(x)以及f(x)∗g(x)三个多项式的系数拿出来写成三个序列,可得:
A={a0,a1,a2,...an}
B={b0,b1,b2,...bm}
C={∑i+j=0ai∗bj,∑i+j=1ai∗bj,∑i+j=2ai∗bj,...∑i+j=m+nai∗bj}
于是惊讶的发现,两个序列的卷积相当于两个多项式乘积后得到的多项式的系数序列。而FFT算法的目的,就是求两个多项式乘积后得到的多项式的系数序列。
多项式有两种表示方法,即系数表示法和点值表示法。
1.系数表示法是我们常用的表示方法,即
f(x)=a0+a1x+a2x2+...+anxn
的形式
2.点值表示法,顾名思义,就是在这个多项式上任意选取不重复的n+1个点,即
f(x)={(x0,y0),(x1,y1),...,(x2,y2)}
可以证明:任何n+1个点可以确定唯一一个最高次项为n次方的多项式(下面是证明,可以看看,也可以忽略)
将一个多项式的系数写成一个系数矩阵(当然,这些系数我们是不知道的)
⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣a0a1a2⋮an⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦
然后将刚刚选取的n+1个点写成下面两个矩阵
⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣1x0x20x30⋯xn01x1x21x31⋯xn11x2x22x32⋯xn2⋮⋮⋮⋮⋱⋮1xnx2nx3n⋯xnn⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣y0y1y2⋮yn⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦
可以得出以下三个矩阵的关系
⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣1x0x20x30⋯xn01x1x21x31⋯xn11x2x22x32⋯xn2⋮⋮⋮⋮⋱⋮1xnx2nx3n⋯xnn⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣a0a1a2⋮an⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦=⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣y0y1y2⋮yn⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦
解出系数矩阵
⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣a0a1a2⋮an⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦=⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣1x0x20x30⋯xn01x1x21x31⋯xn11x2x22x32⋯xn2⋮⋮⋮⋮⋱⋮1xnx2nx3n⋯xnn⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦−1⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣y0y1y2⋮yn⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦
只能唯一确定一个系数矩阵,即只能唯一确定一组系数,即只能唯一确定一个多项式,证明完毕
为了实现多项式乘法并得到系数序列,我们可以考虑实现的方法,如果直接暴力(通过定义直接算系数)是O(n2)复杂度,肯定会超时,于是我们这样考虑
首先可以选取n+1个点,把两个多项式转化为点值表示法。
然后将两个点值表示法的多项式相乘(复杂度为O(n)),然后将新得到的多项式的点值表示法转化为系数表示法。
注意:假设f(x)最高有n次方,g(x)最高有m次方,所以f(x)*g(x)最高有m+n次方,但是m+n可能不是2的幂次方,如果不是,则需要选取点的数量应该是大于m+n的2的幂次方个,假设这个值为1<<k,所以从系数表示法到点值表示法的转化过程中,我们要在f(x)和g(x)内选择1<<k的点,才能保证点值表示法的f(x)和g(x)之后,至少仍有m+n个点,才能确定唯一一个f(x)*g(x)的多项式。所以在选取点的数量的时候,要保证点的个数能确定f(x)*g(x)后的多项式。
如果任意的选取n+m个点然后转化,复杂度还是太高,所以我们需要巧妙的选取1<<k个点,从而使得系数表示法与点值表示法之间转化的复杂度降为O(nlogn),这就是FFT算法。
一个虚数,是由实部和虚部组成,表示为(a+bi),虚数的几何表示可以在坐标系中体现,如下图所示

如图所示,r为此虚数的模长,θ为此虚数的幅角,于是虚数还能表示为r(cosθ,isinθ)的形式
若z1=r1(cosθ1,isinθ1),z2=r2(cosθ2,isinθ2)
于是z1∗z2=r1(cosθ1,isinθ1)∗r2(cosθ2,isinθ2)=r1r2(cos(θ1+θ2),isin(θ1+θ2))
可以得到虚数相乘的几何意义:模长相乘,幅角相加
n次单位根wn,即xn=1在虚数范围内的解,即wnn=1,故n次单位根wn也是一个复数,且模长为1,所以n次单位根在乘方的时候模长不变,幅角等差增大
如图所示,将一个单位圆分成n块,幅角为正且最小的那个虚数即为n次单位根(w1n,红色点),下图是n=8的情况

还要注意,w0n=wnn=1
由图及n次单位根的性质(n次单位根模长为1,幅角为2πn)可得
wkn=cos2kπn+isin2kπn
比较简单的性质:
1.wm+nn=wmn∗wnn=wmn
2.(wmn)n=(wnn)m=1,m∈[0,n−1]
比较复杂(重要)的性质:
1.w2k2n=wkn,如图所示

2.wkn=−wk+n2n
通过上面那个图也可以看出来,注意一点:虚数乘-1,那么虚数的实部与虚部都要乘上-1,所以−wkn在单位圆上的点和wkn在单位圆上的点关于原点对称
3.(wkn)2=wkn2
推导:(wkn)2=(−wk+n2n)2=w2∗kn=wkn2
更加复杂(重要)的性质
n−1∑k=0wikn={0imodn≠0nimodn=0
证明:
当imodn≠0时,相当如等比数列前n项和求和,故n−1∑k=0wikn=w0n(1−(win)n)1−win=1−(wnn)i1−win=0
当imodn=0时,wikn=1,故n−1∑k=0wikn=n∗1=n
至此,你应该知道,多项式从系数表示法转化为点值表示法时,选取的n个点的x值即为n次单位根序列{w0n,w1n,....,wn−1n}
注意,此处的n是f(x)*g(x)的多项式最高次项的次数,不是f(x)或者g(x)最高此项的系数,前面(多项式乘法与FFT)讲过要选取1<<k(再次强调这个值比f(x)*g(x)的多项式最高次的次数m+n还要大,且是2的幂次方)个点
假设现在多项式为
f(x)=a0+a1x+a2x2+...+ak−1xk−1,注意,多项式最高次项为k-1次方,共k项,k如果不等于2的n次幂,就凑成2的n次幂(加系数为0的项)
g(x)=b0+b1x+b2x2+...+bh−1xh−1,注意,多项式最高次项为h-1次方,共h项,h如果不等于2的n次幂,就凑成2的n次幂(加系数为0的项)
那么现在,我们从多项式中选取了这样一些点
(w0n,f(w0n)),(w1n,f(w1n)),(w2n,f(w2n)),...,(wn−1n,f(wn−1n)),共n个点
注意n是大于k+h的最小2次幂,而不是等于k+h!!!!!
此时选取得点数n大于k+h,故可以确定f(x)*g(x)这个多项式,现在,只需确定f(w0n),f(w1n),...,f(wn−1n)的值即可。
如何算值?我们可以观察一下f(win)的值(0≤i<n)
f(win)=a0+a1win+a2(win)2+a3(win)3+...+ak−1(win)k−1
=a0+a2(win)2+a4(win)4+...+ak−2(win)k−2+win(a1+a3(win)2+...+ak−1(win)k−2)
根据(wkn)2=wkn2可得
=a0+a2win2+a4(win2)2+...+ak−2(win2)k−22+win(a1+a3win2+...+ak−1(win2)k−22)
=f1(win)+winf2(win)
再算f1(win)的值和f2(win)的值的时候,我们又可以按照奇偶分开,像算f(win)一样变成两个问题
同时我们发现
f(wi+n2n)=f(−win)
=a0+a2(−win)2+a4(−win)4+...+ak−2(−win)k−2−win(a1+a3(−win)2+...+ak−1(−win)k−2)
=a0+a2win2+a4(win2)2+...+ak−2(win2)k−22−win(a1+a3win2+...+ak−1(win2)k−22)
=f1(win)−winf2(win)
又发现,如果我们要算f(win)的值,我们会先算f1(win)和f2(win)的值,但是算出f1(win)和f2(win)的值之后,不仅能算出f(win)的值,还能同时算出f(wi+n2n)的值
总的来说,如果n=8算出了f(w0n)的值,就算出了f(w4n)的值;算出了f(w1n),就算出来f(w5n)的值;...;算出了f(w3n)的值,就算出了f(w7n)的1值
同理,算f1(win)和f2(win)的值的时候,也是算出一半,得到另一半的值------每次解决问题,都会变成解决一半问题然后直接得到另一半的答案,在解决这一半问题的时候,也变成了解决一半的一半的问题,另一半的一半的问题又被直接得到答案。
这样,就可以把f(w0n),f(w1n),f(w2n),...,f(wn−1n)的值全部算出来,得到了n个点的值;同理算出g(w0n),g(w1n),f(w2n),...,g(wn−1n)的值。得到了f(x)和g(x)的点值表示法,点值表示法相乘,O(n)复杂度就得到了f(x)*g(x)的点值表示法。
代码如下:

FFT变换(递归版)
这里再来详细解释一下,上面代码中的n值
加入开始两个序列是n次多项式和m次多项式,那么开始n值就等于大于等于(m+n)的最小2次幂,如下代码所示
此时lim值就是n的初值
强调一遍wkn=cos2kπn+isin2kπn
FFT逆变换则是将点值表示法转化为系数表示法的步骤。在逆变换之前,我们已经对两个多项式系数序列进行了正变换,即
如上所示,只要将a序列进行逆变换,就可以得到新多项式的系数序列。在FFT正变换之中,我们一直在求f(win)的值。
假设所求多项式系数序列为k0,k1,...kn−1
现在对于新的多项式a(新得到的序列,不是原来的a序列)来求卷积ch=n−1∑i=0ai(w−hn)i
相当于对a序列又进行FFT变换,只不过,我们现在以新序列a为系数的多项式在x=w−in的一系列值
相当于求C(x)=a0+a1x+a2x2+...+an−1xn−1在x=w0n,w−n1,w−2n,...w−(n−1)n的点值表示法
而我们知道ax=k0+k1wxn+k2(wkn)2+...=∑n−1j=0kj(wxn)j
n−1∑i=0ai(w−hn)i
=n−1∑i=0(n−1∑j=0kj∗(win)j)(w−hn)i
=n−1∑j=0kj(n−1∑i=0(wj−hn)i)
=n−1∑j=0kj(n−1∑i=0w(j−h)∗in)
由于当且只当j=h时,(j−h)才为i的倍数,其余时候为0,根据n次单位根更加复杂的性质可得(不知道的看看n次单位根预备性质)
n−1∑j=0kj(n−1∑i=0w(j−h)∗in)
=n∗kh
所以
kh=∑n−1i=0ai(w−hn)ih=chn
说明:w−in与win实部相同,虚部相反。现在知道代码中inv的作用了吧
当inv=1,表示FFT变换;当inv=-1,表示FFT逆变换
为啥要加0.5捏:首先我们确定答案肯定是整数,但是在代码实现过程中,由于有π介入计算,导致答案可能变小了(误差很小),所以加一个0.5,然后强制类型转换切掉小数部位
比如说本来答案是1,但是代码实现的结果可能是0.999999999999,此时加上了0.5,然后转换为int,答案就是1了。
强调一遍wkn=cos2kπn+isin2kπn
FFT代码

FFT变换
在进行FFT之前,我们有两个多项式系数序列,用数组存系数的不要用普通的数组存,而要用复数数组,因为在计算的时候系数涉及复数计算,所以这样
这样,复数的实部存系数的值,虚部为0
一个算多项式相乘后系数序列的程序

多项式相乘
测试用例

测试
在讲FFT迭代版本之前,我打算先用一个例子来展示FFT变换,加深认识
假设
多项式f(x)=3+2x+x2
多项式g(x)=2+x+2x2
在草稿纸上算一下,这两个多项式相乘的结果为h(x)=6+7x+10x2+5x3+2x4
先得到f(x)与g(x)多项式的系数序列,注意应该写成复数形式的序列(复数形式为(real,image))
多项式f(x)的系数序列为A={(3,0),(2,0),(1,0)}
多项式g(x)的系数序列为B={(2,0),(1,0),(2,0)}
然后判断两个多项式相乘后有多少项,应该不超过3+3=6项。而大于6的最小2次幂应该为8,所以要把两个多项式序列补成8项。
故(用(0,0)来补项)
多项式f(x)的系数序列为A={(3,0),(2,0),(1,0),(0,0),(0,0),(0,0),(0,0),(0,0)}
多项式g(x)的系数序列为B={(2,0),(1,0),(2,0),(0,0),(0,0),(0,0),(0,0),(0,0)}
对于每一个f(wi8),可以寻得规律
f(wi8)=A0+A1wi8+A2(wi8)2+A3(wi8)3+A4(wi8)4+A5(wi8)5+A6(wi8)6+A7(wi8)7
=A0+A2(wi8)2+A4(wi8)4+A6(wi8)6+wi8(A1+A3(wi8)2+A5(wi8)4+A7(wi8)6)
=A0+A2wi4+A4(wi4)2+A6(wi4)3+wi8(A1+A3wi4+A5(wi4)2+A7(wi4)3)
=A0+A4(wi4)2+wi4(A2+A6(wi4)2)+wi8(A1+A5(wi4)2+wi4(A3+A7(wi4)2))
=A0+A4wi2+wi4(A2+A6wi2)+wi8(A1+A5wi2+wi4(A3+A7wi2))
然后
f(wi+48)=f(−wi8)=A0+A1(−wi8)+A2(−wi8)2+A3(−wi8)3+A4(−wi8)4+A5(−wi8)5+A6(−wi8)6+A7(−wi8)7
=A0+A2(wi8)2+A4(wi8)4+A6(wi8)6−wi8(A1+A3(wi8)2+A5(wi8)4+A7(wi8)6)
只要算出来A0+A2wi4+A4(wi4)2+A6(wi4)3和A1+A3wi4+A5(wi4)2+A7(wi4)3的值,就可以算出f(wi8)和f(wi+48)的值
同理,算出A0+A4(wi4)2和A2+A6wi2的值,不仅可以得到A0+A2wi4+A4(wi4)2+A6(wi4)3的值,还可以得到A0+A2wi+24+A4(wi+24)2+A6(wi+24)3
然后将每一个序列分成长度为2的小序列(下面是用1~8来编号的)

后面的是重点,对于每一个长度为2的小序列(设为(x,y)):
第一个元素x重新赋值为x+w01∗y
第二个元素y重新赋值为x−w01∗y
(这个w11这样来的:$w_{当前序列一半长度l}^{一半长度l的第i个元素,i从0到l-1}$)
重新赋值完之后就可以合并了,每两个长度为2的相邻序列合并为一个长度为4的序列
对于每一个长度为4的序列(设为(a,b,c,d)):
第一个元素a重新赋值为a+w02∗c
第三个元素c重新赋值为a−w02∗c
第二个元素b重新赋值为b+w12∗d
第四个元素d重新赋值为b−w12∗d
重新赋值完之后就可以合并了,每两个长度为4的相邻序列合并为一个长度为8的序列(此时全部小序列已经合并为一个序列了)
对于整个序列(长度为8),按照上面的赋值方法重新赋值(可以参考代码),于是FFT变换就ok了
两个序列就变成了
A={A′0,A′1,A′2,A′3,A′4,A′5,A′6,A′7}
B={B′0,B′1,B′2,B′3,B′4,B′5,B′6,B′7}
将两个序列对应元素相乘,得到新序列
C={A′0∗B′0,A′1∗B′1,A′2∗B′2,A′3∗B′3,A′4∗B′4,A′5∗B′5,A′6∗B′6,A′7∗B′7}
然后将C序列进行FFT逆变换,就是多项式相乘后的系数序列了。
ps:FFT逆变换就是把当前序列再来一次FFT变换,只不过在重新赋值环节的时候是乘上w−1mid,而不是wimid,最后得到的序列虚部为0,将实部全部除以8即可(这里不懂参考代码或前面逆变换的讲解)
综上所述,FFT的序列更新是有规律的,先把序列分解为序列长度为2的小段,每一段更新完毕后就合并小段,继续一组一组更新,再合并,最后合并成一个序列,然后再最后一次更新完毕。
如图所示:

最后序列{A′′′0,A′′′1,...,A′′′n−1}为原序列进过FFT变换后的序列。
不像递归版本,迭代版本是直接从每一个小区间开始更新,然后合并,然后更新...,但是我们需要找到原系数序列经过奇偶分离的一系列操作后每个值的位置发生了什么变化。
假设原系数序列为{A0,A1,A2,A3,A4,A5,A6,A7},下图是奇偶分离的示意图

没错,这不是巧合,序列经过奇偶分离之后,前后位置数的二进制是对称的,根据这一点我们就不需要用递归来实现,可以直接来得到奇偶分离最终序列,然后再进行FFT
FFT迭代版比递归版快的不是一点。。
关于序列奇偶分离可以用下面方法
于是就可以实现迭代版本FFT

FFT变换(迭代版)
FFT递归版本

FFT变换(迭代版)
FFT迭代版本

FFT变换(迭代版)
__EOF__
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· DeepSeek 解答了困扰我五年的技术问题
· 为什么说在企业级应用开发中,后端往往是效率杀手?
· 用 C# 插值字符串处理器写一个 sscanf
· Java 中堆内存和栈内存上的数据分布和特点
· 开发中对象命名的一点思考
· PPT革命!DeepSeek+Kimi=N小时工作5分钟完成?
· What?废柴, 还在本地部署DeepSeek吗?Are you kidding?
· DeepSeek企业级部署实战指南:从服务器选型到Dify私有化落地
· 程序员转型AI:行业分析
· 重磅发布!DeepSeek 微调秘籍揭秘,一键解锁升级版全家桶,AI 玩家必备神器!