多项式的点值表示(Point Value Representation)
设多项式的系数表示(Coefficient Representation):
Pa(x)=a0+a1x+a2x2+⋯+an−1xn−1=n−1∑i=0aixiPa(x)=a0+a1x+a2x2+⋯+an−1xn−1=n−1∑i=0aixi
则我们对上面的式子可以代入不同的 nn 个 xx 的值,构成一个 nn 维向量:
[Pa(x0)Pa(x1)Pa(x2)⋮Pa(xn−1)]=[1x0x20⋯xn−101x1x21⋯xn−111x2x22⋯xn−12⋮⋮⋮⋱⋮1xn−1x2n−1⋯xn−1x−1][a0a1a2⋮an−1]⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣Pa(x0)Pa(x1)Pa(x2)⋮Pa(xn−1)⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦=⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣1x0x20⋯xn−101x1x21⋯xn−111x2x22⋯xn−12⋮⋮⋮⋱⋮1xn−1x2n−1⋯xn−1x−1⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣a0a1a2⋮an−1⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦
更简洁的写法:
Pa=XαPa=Xα
对上式观察后发现,XX 是所谓的范德蒙德矩阵(Vandermonde's Matrix),在 nn 个 xx 的值不同的情况下,其行列式的值为:
det(X)=∏0⩽i<j⩽n−1(xj−xi)det(X)=∏0⩽i<j⩽n−1(xj−xi)
很明显,当所有 nn 个 xx 取值不同时,其行列式不为零,因此 XX 可逆。
所以我们可以唯一确定多项式系数构成的向量 αα:
α=X−1Paα=X−1Pa
也就是说,多项式 Pa(x)Pa(x) 还可以由 nn 个 xx 代入得到的 nn 个点值来唯一表示:
{[x0,P(x0)],[x1,P(x1)],[x2,P(x2)],⋯,[xn−1,P(xn−1)]}{[x0,P(x0)],[x1,P(x1)],[x2,P(x2)],⋯,[xn−1,P(xn−1)]}
这就是多项式的点值表示。
多项式的点值表示是指,对于 nn 次多项式,可以用 nn 个不同的 xx 和与之对应的多项式的值 P(x)P(x) 构成一个长度为 nn 的序列,这个序列唯一确定多项式,并且能够与系数表示相互转化。
nn 次单位根
了解了多项式的点值表示,一个很自然的问题是:如何选择 xx 的值,来防止其指数大小爆炸型增长呢?这里可以借用复数的单位根。
简单回顾一下,复数有两种表示方法:迪卡尔积坐标表示和极坐标表示,这里我们用到的是后者:
z=reiθz=reiθ
ii 是虚数单位,rr 表示模长,θθ 表示相角。
复数的 nn 次单位根 ωω 需要满足条件:
ωnn=1ωnn=1
了解复数乘法及其几何意义的同学知道,复数相乘则相角相加,相当于复数点逆时针转动;nn 个复数相乘则说明有 nn 个相角相加,nn 次逆时针转动。因此:
ωnn=1=ei⋅2πωnn=1=ei⋅2π
则 nn 次单位根为:
n√ωn=ei⋅2πnn√ωn=ei⋅2πn
很容易想到,nn 等分 2π2π 相当于 nn 等分圆。下图是 n=16n=16 的例子:

引入了 nn 次单位根,我们就可以找到任意 nn 个不同的点值 xx,并且不用担心指数增长带来的大小爆炸性增长的问题。
DFT
设长度为 nn 的离散序列 αT=[a0,a1,⋯,an−1]αT=[a0,a1,⋯,an−1],构建多项式:
Pa(x)=n−1∑i=0aixi=xαTPa(x)=n−1∑i=0aixi=xαT
其中,xT=[1,x,x2,⋯,xn−1]xT=[1,x,x2,⋯,xn−1]
用 nn 次单位根生成 nn 个不同的值:ω0n, ω1n, ω2n, ⋯, ωn−1nω0n, ω1n, ω2n, ⋯, ωn−1n
对应的点值表示可以用下面的矩阵运算完成:
[Pa(ω0n)Pa(ω1n)Pa(ω2n)⋮Pa(ωn−1n)]=[111⋯11ω1nω2n⋯ωn−1n1ω2nω4n⋯ω2(n−1)n⋮⋮⋮⋱⋮1ωn−1nω2(n−1)n⋯ω(n−1)(n−1)n][a0a1a2⋮an−1]⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣Pa(ω0n)Pa(ω1n)Pa(ω2n)⋮Pa(ωn−1n)⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦=⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣111⋯11ω1nω2n⋯ωn−1n1ω2nω4n⋯ω2(n−1)n⋮⋮⋮⋱⋮1ωn−1nω2(n−1)n⋯ω(n−1)(n−1)n⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣a0a1a2⋮an−1⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦
更简洁的写法:
Pa=ΩαPa=Ωα
序列 PaPa 被称为 αα 的离散傅里叶变换(DFT),也称为频域序列。
很明显,DFT 的计算时间复杂度是 O(n2)O(n2)
IDFT
离散傅里叶反变换,就是根据 DFT 得到的频域序列算出多项式的系数(也称时域序列)。这可以根据单位根矩阵 ΩΩ 的逆矩阵 Ω−1Ω−1 得到
α=Ω−1Paα=Ω−1Pa
一般来说,计算矩阵的逆的时间复杂度高达 O(n3)O(n3)。所幸,单位根矩阵的逆 Ω−1Ω−1 有个优良的性质可以省去庞大的计算量:
[111⋯11ω1nω2n⋯ωn−1n1ω2nω4n⋯ω2(n−1)n⋮⋮⋮⋱⋮1ωn−1nω2(n−1)n⋯ω(n−1)(n−1)n]−1=1n[111⋯11ω−1nω−2n⋯ω−(n−1)n1ω−2nω−4n⋯ω−2(n−1)n⋮⋮⋮⋱⋮1ω−(n−1)nω−2(n−1)n⋯ω−(n−1)(n−1)n]⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣111⋯11ω1nω2n⋯ωn−1n1ω2nω4n⋯ω2(n−1)n⋮⋮⋮⋱⋮1ωn−1nω2(n−1)n⋯ω(n−1)(n−1)n⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦−1=1n⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣111⋯11ω−1nω−2n⋯ω−(n−1)n1ω−2nω−4n⋯ω−2(n−1)n⋮⋮⋮⋱⋮1ω−(n−1)nω−2(n−1)n⋯ω−(n−1)(n−1)n⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦
也就是说,求 Ω−1Ω−1 就是将里面元素的指数改成其相反数,再对所有元素除以 nn。
有了这个性质,我们就可以得到:
[a0a1a2⋮an−1]=1n[111⋯11ω−1nω−2n⋯ω−(n−1)n1ω−2nω−4n⋯ω−2(n−1)n⋮⋮⋮⋱⋮1ω−(n−1)nω−2(n−1)n⋯ω−(n−1)(n−1)n][Pa(ω0n)Pa(ω1n)Pa(ω2n)⋮Pa(ωn−1n)]⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣a0a1a2⋮an−1⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦=1n⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣111⋯11ω−1nω−2n⋯ω−(n−1)n1ω−2nω−4n⋯ω−2(n−1)n⋮⋮⋮⋱⋮1ω−(n−1)nω−2(n−1)n⋯ω−(n−1)(n−1)n⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦⎡⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣Pa(ω0n)Pa(ω1n)Pa(ω2n)⋮Pa(ωn−1n)⎤⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦
很明显,DFT 和 IDFT 的计算时间复杂度都是 O(n2)O(n2),效率并不高。
FFT
FFT 是用分治法(Divide & Conquer)的思想,用来优化 DFT 计算矩阵相乘的时间复杂度过高这一问题的算法。
设 nn 次多项式 P(x)P(x):
P(x)=a0+a1x+a2x2+a3x3+⋯+an−1xn−1P(x)=a0+a1x+a2x2+a3x3+⋯+an−1xn−1
我们把多项式 P(x)P(x) 按照系数下标的奇偶性分成两部分
Pe(x2)=a0+a2x2+a4(x2)2+⋯+an−2xn−2xPo(x2)=x[a1+a3x2+a5(x2)2+⋯+an−1xn−2]Pe(x2)=a0+a2x2+a4(x2)2+⋯+an−2xn−2xPo(x2)=x[a1+a3x2+a5(x2)2+⋯+an−1xn−2]
则多项式 P(x)P(x) 就是奇偶两部分之和
P(x)=Pe(x2)+xPo(x2)P(x)=Pe(x2)+xPo(x2)
从上式中可以看出,多项式 P(x)P(x) 可以由两个系数个数为 n/2n/2 的多项式 Pe(x2)Pe(x2) 和 Po(x2)Po(x2) 计算得到。
对于 Pe(x2)Pe(x2) 和 Po(x2)Po(x2),我们令 x=x2x=x2,就会发现这很明显是一个递归的过程:
Pe(x)=a0+a2x+a4x2+⋯+an−2xn−22Pee(x2)=a0+a4x2+a8x4+⋯+an−4xn−22−1xPeo(x2)=x[a2+a6x2+a10x4+⋯+an−2xn−22−1]Pe(x)=a0+a2x+a4x2+⋯+an−2xn−22Pee(x2)=a0+a4x2+a8x4+⋯+an−4xn−22−1xPeo(x2)=x[a2+a6x2+a10x4+⋯+an−2xn−22−1]
就可以求出
Pe(x)=Pee(x2)+xPeo(x2)Pe(x)=Pee(x2)+xPeo(x2)
同理
Po(x)=a1+a3x+a5x2+⋯+an−1xn−22Poe(x2)=a1+a5x2+a9x4+⋯+an−3xn−22−1xPoo(x2)=x[a3+a7x2+a11x4+⋯+an−1xn−22−1]Po(x)=a1+a3x+a5x2+⋯+an−1xn−22Poe(x2)=a1+a5x2+a9x4+⋯+an−3xn−22−1xPoo(x2)=x[a3+a7x2+a11x4+⋯+an−1xn−22−1]
同样可以求出
Po(x)=Poe(x2)+xPoo(x2)Po(x)=Poe(x2)+xPoo(x2)
递归的终止条件很明显,就是当遇到多项式中只有一个系数时返回该系数。
因此我们将 nn 个单位根 ω0n, ω1n, ω2n, ⋯, ωn−1nω0n, ω1n, ω2n, ⋯, ωn−1n 带入多项式 P(x)P(x):
P(ωkn)=Pe(ω2kn)+ωknPo(ω2kn)(k=0,1,⋯,n−1)P(ωkn)=Pe(ω2kn)+ωknPo(ω2kn)(k=0,1,⋯,n−1)
刚刚提到,多项式 Pe(x2)Pe(x2) 和 Po(x2)Po(x2) 的系数个数为 n/2n/2 ,恰好 nn 个单位根平方后也只剩下 n/2n/2 个不同的单位根,简单证明如下:
首先证明:
(ωkn)2=ei2πkn⋅2=ei2πkn/2=ωkn/2(ωkn)2=ei2πkn⋅2=ei2πkn/2=ωkn/2
因此相当于 nn 等分圆变成了 n/2n/2 等分圆。下面证明 k=0, 1, ⋯, n2−1k=0, 1, ⋯, n2−1:
设 k1=mk1=m, k2=m+n/2k2=m+n/2,m=0,1,⋯,n2−1m=0,1,⋯,n2−1。则
ωk1n/2=ωmn/2=ei2πmn/2ωk2n/2=ωm+n/2n/2=ei2π(m+n/2)n/2=ei(2πmn/2+2π)=ei2πmn/2ωk1n/2=ωmn/2=ei2πmn/2ωk2n/2=ωm+n/2n/2=ei2π(m+n/2)n/2=ei(2πmn/2+2π)=ei2πmn/2
即:
ωmn/2=ωm+n/2n/2(m=0,1,⋯,n2−1)ωmn/2=ωm+n/2n/2(m=0,1,⋯,n2−1)
由此可以说明,ωknωkn 是周期序列,其最小正周期为 n/2n/2。因此 k=0,1,⋯,n2−1k=0,1,⋯,n2−1
因此多项式 P(ωkn)P(ωkn) 可改写为
P(ωkn)=Pe(ωkn/2)+ωknPo(ωkn/2)(k=0,1,⋯,n2−1)P(ωkn)=Pe(ωkn/2)+ωknPo(ωkn/2)(k=0,1,⋯,n2−1)
从上式中可以看出,计算 Pe(ωkn/2)Pe(ωkn/2) 和 Po(ωkn/2)Po(ωkn/2) 各需要 n/2n/2 次乘法运算,即一共 nn 次乘法运算;而计算Pee(ωkn/4)Pee(ωkn/4) 、 Peo(ωkn/4)Peo(ωkn/4) 、Poe(ωkn/4)Poe(ωkn/4) 和 Poo(ωkn/4)Poo(ωkn/4) 各需要 n/4n/4 次乘法运算,即一共 nn 次乘法运算……以此类推,这种每层递归规模减半的递归深度很明显是 lognlogn,因此 FFT 算法的时间复杂度就是 O(nlogn)O(nlogn)
IFFT
快速傅里叶变换反变换同样是优化 IDFT 计算矩阵相乘的时间复杂度。由于 DFT 和 IDFT 核心操作一样,都是矩阵相乘,因此 FFT 和 IDFT 的核心操作就是利用分治的思想减少矩阵相乘的运算量。可以想到,FFT 的过程可以直接移植到 IFFT 中来,需要修改的两个地方是
- 单位根指数部分改成其相反数:ω−knω−kn
- 得到的结果均除以 nn
因此,IFFT 的时间复杂度也是O(nlogn)O(nlogn)
Refer: COMP3121/9101, CSE UNSW
Written with StackEdit.
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 一个奇形怪状的面试题:Bean中的CHM要不要加volatile?
· [.NET]调用本地 Deepseek 模型
· 一个费力不讨好的项目,让我损失了近一半的绩效!
· .NET Core 托管堆内存泄露/CPU异常的常见思路
· PostgreSQL 和 SQL Server 在统计信息维护中的关键差异
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 如何使用 Uni-app 实现视频聊天(源码,支持安卓、iOS)
· C# 集成 DeepSeek 模型实现 AI 私有化(本地部署与 API 调用教程)