赞
踩
例如
所以
那么我们该如何来储存一个多项式呢?最一般的方法就是储存多项式的系数,我们只需要把系数映射到一个表中,按顺序存储,这样的话,表中的第k个数字正好对应多项式的k阶项系数,我们称这种表示方法为多项式的系数表示法。
一般地,给定两个d阶多项式,
二者乘积应为2d阶多项式,并且这一算法,如果用naive的乘法分配律来算,时间复杂度为,因为多项式A的每一项都会跟B的每一项相乘。
问题来了,我们能改进这个算法吗?
这篇文章最奇妙的第一个想法,首先我们站在略微不同的角度来看待多项式,我们来看最简单的多项式,即一阶多项式(线性函数)
我们可以用两个系数来表示线性函数,即截距(0阶系数)和斜率(1阶系数),注意:一对系数唯一确定一条直线。那么,还有没有其他方法可以唯一确定一条直线呢?
这样的表示是很多的,我们这里用直线的两点表示,几何上我们学过:两点确定一条直线,事实上,高阶多项式也有类似的性质,即:任一d阶多项式都由d+1个点唯一确定,比如三个不同的点,只能唯一确定一个二次函数曲线。为了让大家信服这个性质,我们来证明一下它,假设我们有P(x)这个d阶多项式,通过了给定的d+1个点。
我们要证明,这个多项式的系数是唯一的,所以我们把这d+1个点代入多项式的方程,得到d+1个方程
将这个方程组用矩阵表示
矩阵
是范德蒙德矩阵,我们记作M。矩阵M可逆,只要这些点两两都不相同,故我们只需要证明矩阵M的行列式(范德蒙德矩阵)不为0,即可得出d阶行列式的系数是被d+1个点唯一确定的,从而多项式也是唯一的。
总结一下,我们现在有两种多项式的表示法,第一种是系数表示法,第二种是d+1个点表示法,称为值表示法。利用值表示法,多项式乘法变得非常简单,我们回到最开始的例子,
我们知道乘出来的多项式C(x)为4阶,所以需要5个点表示,我们现在只需要挑出5个点,将对应每个点的两个值相乘,得到C在每个点的函数值,根据我们前面证明的性质,这5个点确定了唯一的多项式,
可以看到,利用值表示法,我们计算多项式乘法的时间从一下子缩短到了
!
给定两个系数表示的d阶多项式
系数值 值
系数
我们已经知道了值表示法计算多项式乘法更快,所以我们可以先计算两个多项式在2d+1个点上的值,然后将函数值一对对地乘起来,从而得到乘出来的值表示,然后最后一步,想办法将值转换回系数表示。
大体思路有了,问题是具体怎么实现?也就是:
如何把系数表示转换成值表示?
如何把值表示转换成系数表示?
对这两个问题的解答也就是这篇文章的目标——FFT和IFFT。
我们先来解决第一个问题:从系数表示到值表示,我们称这个问题为求值,给定多项式
和n个点,我们想计算多项式在这n个点上的值,n大于等于d+1,最粗暴的做法就是在X轴上随便挑选n个点,一个一个计算函数值,
我们当然可以这样做,但是你如果仔细看看,就发现这样的话,反而又回到了最开始的地方,因为每个点的计算都是,加起来是
, 我们现在的问题是要寻找一个简单的方法,我们尝试着简化一下这个问题:举例比如
,在n=4个点进行求值,那我们来想一下,有没有这样一对点:当你知道一个点的函数值时,另一个点的也可以立刻知道? 答案是肯定的,有!
如果我们选了,
,那
的值我们立马知道了,也是1,同样地,如果我们选了
,那
的值我们也立马知道了,都是4。
由此可见,我们只要取互为相反数的数对,这其中的机理非常简单:因为是个偶函数。
OK,问题来了,如果我们取,这个技巧还使适用吗?事实上,这个技巧依旧正确,但是我们在计算的时候需要在-x的函数值前面加一个负号。所以在这两个偶/奇函数的情形中,我们只需要计算一半数量的点就好了,因为剩下一半从奇偶性中就可以得到了。
对于一般的多项式,我们怎样使用这个技巧呢?举个例子
我们把这个多项式分为奇/偶函数两部分,我们把奇函数部分的x提出来,那么括号里的正是两个偶多项式函数,我们将两个括号内的部分分别命名为和
,需要注意的是,这两个部分都是偶函数,这样的分解使得我们在计算一对相反数的函数值时,只需要计算其中一个,同时我们可以把两个部分都看做
的函数,这样的话
和
的阶数都降到了2,比原来的5阶降低了许多。
推广一下这个想法,如果我们有一个n-1阶多项式,我们想知道它在n个点的函数值,我们可以把多项式分为奇/偶两部分,每部分都只有n/2-1阶,所以对于这俩只有原来一半阶数的多项式,我们只需要按照之前的求值方法进行计算即可,不同的是这次需要计算我们所给的点的平方的函数值,但我们只需要取一对对地相反数,那么在平方以后我们就只剩n/个点了。(这里其实就可以开始使用递归了,埋个伏笔)
总结:我们想计算多项式P在n个点上的值,这n个点是一对对的相反数,我们把多项式拆分成两部分,
从而我们有两个阶为n/2-1的多项式,每个多项式也只需要求n/2个点的值,我们只需要递归地求这两个小多项式的值,利用下方这两个公式带回去,就可以得到原多项式的n个值。
最终我们将得到原多项式的值表达,如果这一切都走得通,那么我们地时间复杂度仅为,因为这两个子问题都是原问题规模的一半,而且我们只需要线性时间来计算原多项式的值,相比于原来的
是一个质的飞跃。
不过这里有一点问题,实际上在递归步骤:我们假设了每个多项式我们都使用相反数对来求值,但是对两个子问题而言,每个求值点都是平方数,所以都是正的,递归不成立。
那我们能不能把这些新的求值点也搞成相反数对呢?而这里,才是FFT最精妙的地方,只有一条路可走:我们把这些点都取成复数,我们还要专门挑一些复数,使得它们平方之后,依旧是正负成对出现的,那我们怎么取这些复数呢?我们来看一个例子逆向推导来找到答案。
需要
个点来求值
这些点应该正负成对出现,我们写作,我们的递归算法需要在
这两个点求子函数的值,那我们需要
应当也是一对正负数对,所以有
,再考虑一层递归,那么我们只有一个点:
。
我们假设,那么这四个数应该是1,-1,i,-i,第二行就是1和-1,最下方就是1。从另一个角度来看,我们实际上是解了方程
,而最上方四个点必须满足他们的四次方都等于1,事实上,我们知道方程
有四个根,称为1的四个四次方根。
至少需要
个点来求值
我们的递归算法每次递归后,点的个数要折半,所以我们取一个2的次方数,,那我们需要寻找8个数,形成4个正负对,并保正每个数字的8次方都是1,显然这8个数是1的八个八次方根。
需要
个点去求值,我们选
,这些点就相当于1的n个n次方根。关于为什么这些数可以,我做一些复数方面知识的补充:
补充:1的n个n次方根可以被解释为复平面上沿着单位圆等距排布的一系列点,任一相邻两点间的夹角为,从而这些点的最便捷的表达方法就是用欧拉公式表示成复指数。
我们定义,那么不同的i取值就代表了不同的单位复n次方根,所以当我们想在1的n次方根求值时,我们实际上是在
上求值。
所以回到原来的问题:为啥n次方根就可以用于递归过程呢?有两个主要原因:
首先,我们要求的点是正负成对的,这里正好是这么一对;
另外,当我们考虑对应的递归低阶奇/偶函数时,n次方根有一个好性质,平方以后我们得到的正好是n/2个n/2次方根,他们也是正负配对的,完美适合递归;我们再平方,又得到n/4个n/4次方根,适用于下一层的递归,最终直到我们只剩一个点。
下面终于是激动人心的FFT算法了!
FFT算法的输入是多项式的n个系数其中n是2的幂次,我们取
为1的n次方根,递归的底层情况是n=1,此时我们只在一个点上求多项式的值。
递归的主要命令就是把多项式拆分成奇/偶函数两部分,两部分就是调用函数两次,此时这两个子多项式函数都是n/2阶的,所以对应的求值点将是1的n/2次方根。
我们假设递归调用的两部分函数值为和
,然后我们把这两部分的函数值结合一下,计算出原多项式函数在对应的点的函数值,结合起来的核心思想就是利用正负对,不过这里要为n次方根略微修改一下我们的表达式,(为了方便写程序,这里的指标都是从0开始的。)此时
应该写为
,所以公式改写如下:
同时我们观察到,这里的正负点对是,所以我们可以把第二个式子再改写如下:
还有一个不可忽略的小细节:和
上的第j个元素,正好对应了
和
那么我们可以把我们的公式右侧全部写作和
,这样就方便写程序了。
最终FFT输出的是原多项式在n个n次方根的值。下面我们用程序来实现上面的过程。
- def FFT(P):
- # P- [po,P1,...,Pn-1]系数表示
- n=len(P)#n是2的幂次
- if n== l:
- return P
- W=e^(2πi/n)
- Pe,Po=P[::2],P[1::2]
- Ye,Yo= FFT(Pe), FFT(Po)
- y=[0]*n
- for j in range(n/2):
- y[j] =ye[j] + wjyo[j]
- y[j + n/2] = ye[j]一wjyo[j]
- return y
把值表示转换为系数表示,这个过程叫做插值。
参考前面将求值表示为矩阵-向量乘积,我们用系数向量乘以点对应的(范德蒙德)矩阵,从而得到对应的值,因为FFT中是对应的第k个n次方根,所以我们把这个矩阵乘法改写如下,改写后的矩阵称为离散傅里叶变换(DFT)矩阵,而我们前面所说的FFT就是快速计算DFT和向量乘法的算法,在n次方根处求值其实就是DFT矩阵和向量乘法的特例。
插值实际上就是给定位置的函数值,计算函数的系数,写作矩阵就是DFT矩阵的逆和值表示向量相乘。
上下两矩阵结构一致,我们只需要把上矩阵的换成
,就可以变为下矩阵了,这表明我们几乎可以在插值问题上套用FFT。
我们来对比一下区别,在求值过程中我们通过多项式系数和FFT算法,计算n次单位根处的函数值,这实际上是DFT矩阵和系数向量的乘法,其中,而在插值中,我们实际上是在使用逆FFT算法,逆FFT算法的输入是多项式在n次单位根处的函数值,输出多项式的系数,即反向操作FFT算法。
前面我们已经知道我们需要做一个DFT矩阵的逆和向量的乘法,求DFT矩阵的逆,我们只需要把原来矩阵的换成
,这意味着逆FFT实际上和FFT的操作一模一样,只是原本的系数向量变成了值表示向量,而且
换成了
,就这样我们就得到了计算插值的逆FFT算法。超级简单!
下面展示一下IFFT算法的程序:
- def IFFT(P):
- # P- [po,P1,...,Pn-1] 值表示
- n=len(P)#n是2的幂次
- if n== l:
- return P
- W=(1/n)*e^(-2πi/n)
- Pe,Po=P[::2],P[1::2]
- Ye,Yo= IFFT(Pe), IFFT(Po)
- y=[0]*n
- for j in range(n/2):
- y[j] =ye[j] + wjyo[j]
- y[j + n/2] = ye[j]一wjyo[j]
- return y
可以看出IFFT与FFT高度相同,只是调用函数从FFT变为IFFT,然后还有第6行略有区别,其他完全一样。
把多项式用对应点的值表示
利用相反数来保证点成对出现
用1的n次单位复根,这样每次平方完,它们还是正负成对出现,
IFFT算法和FFT算法除了第6行以外一模一样
附视频链接:https://www.youtube.com/watch?v=h7apO7q16V0&feature=emb_logo
感兴趣的同学可以自行看看原视频,真的很通透,看完之后立马对傅里叶变换的理解上一个大台阶
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。