引言 快速计算两个多项式的乘法

1.多项式乘法

例如

所以

那么我们该如何来储存一个多项式呢？最一般的方法就是储存多项式的系数，我们只需要把系数映射到一个表中，按顺序存储，这样的话，表中的第k个数字正好对应多项式的k阶项系数,我们称这种表示方法为多项式的系数表示法。

一般地，给定两个d阶多项式，

二者乘积应为2d阶多项式，并且这一算法，如果用naive的乘法分配律来算，时间复杂度为,因为多项式A的每一项都会跟B的每一项相乘。

问题来了，我们能改进这个算法吗？

2.多项式乘法的改进

这篇文章最奇妙的第一个想法，首先我们站在略微不同的角度来看待多项式，我们来看最简单的多项式，即一阶多项式(线性函数)

我们可以用两个系数来表示线性函数，即截距（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个点确定了唯一的多项式，

可以看到，利用值表示法，我们计算多项式乘法的时间从一下子缩短到了！

3.快速多项式乘法

给定两个系数表示的d阶多项式

系数值 值系数

我们已经知道了值表示法计算多项式乘法更快，所以我们可以先计算两个多项式在2d+1个点上的值，然后将函数值一对对地乘起来，从而得到乘出来的值表示，然后最后一步，想办法将值转换回系数表示。

大体思路有了，问题是具体怎么实现？也就是：

1. 如何把系数表示转换成值表示？

1. 如何把值表示转换成系数表示？

对这两个问题的解答也就是这篇文章的目标——FFT和IFFT。

一 FFT（快速傅里叶变换）推导

1.求值

我们先来解决第一个问题：从系数表示到值表示，我们称这个问题为求值，给定多项式

和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，-1，i，-i，第二行就是1和-1，最下方就是1。从另一个角度来看，我们实际上是解了方程，而最上方四个点必须满足他们的四次方都等于1，事实上，我们知道方程有四个根，称为1的四个四次方根。

1. 至少需要个点来求值

我们的递归算法每次递归后，点的个数要折半，所以我们取一个2的次方数，，那我们需要寻找8个数，形成4个正负对，并保正每个数字的8次方都是1，显然这8个数是1的八个八次方根。

1. 需要个点去求值,我们选，这些点就相当于1的n个n次方根。关于为什么这些数可以，我做一些复数方面知识的补充:

补充：1的n个n次方根可以被解释为复平面上沿着单位圆等距排布的一系列点，任一相邻两点间的夹角为，从而这些点的最便捷的表达方法就是用欧拉公式表示成复指数。

我们定义,那么不同的i取值就代表了不同的单位复n次方根，所以当我们想在1的n次方根求值时，我们实际上是在上求值。

所以回到原来的问题：为啥n次方根就可以用于递归过程呢？有两个主要原因：

• 首先，我们要求的点是正负成对的，这里正好是这么一对；

• 另外，当我们考虑对应的递归低阶奇/偶函数时，n次方根有一个好性质，平方以后我们得到的正好是n/2个n/2次方根，他们也是正负配对的，完美适合递归；我们再平方，又得到n/4个n/4次方根，适用于下一层的递归，最终直到我们只剩一个点。

2.FFT算法及程序

下面终于是激动人心的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 PW=e^(2πi/n)Pe,Po=P[::2],P[1::2]Ye,Yo= FFT(Pe), FFT(Po)y=[0]*nfor j in range(n/2):y[j] =ye[j] + wjyo[j]y[j + n/2] = ye[j]一wjyo[j]return y

三 IFFT算法

1.插值

把值表示转换为系数表示，这个过程叫做插值。

参考前面将求值表示为矩阵-向量乘积，我们用系数向量乘以点对应的（范德蒙德）矩阵，从而得到对应的值，因为FFT中是对应的第k个n次方根，所以我们把这个矩阵乘法改写如下，改写后的矩阵称为离散傅里叶变换（DFT）矩阵，而我们前面所说的FFT就是快速计算DFT和向量乘法的算法，在n次方根处求值其实就是DFT矩阵和向量乘法的特例。

插值实际上就是给定位置的函数值，计算函数的系数，写作矩阵就是DFT矩阵的逆和值表示向量相乘。

2.IFFT算法及程序

上下两矩阵结构一致，我们只需要把上矩阵的换成，就可以变为下矩阵了，这表明我们几乎可以在插值问题上套用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 PW=(1/n)*e^(-2πi/n)Pe,Po=P[::2],P[1::2]Ye,Yo= IFFT(Pe), IFFT(Po)y=[0]*nfor 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. 把多项式用对应点的值表示

1. 利用相反数来保证点成对出现

1. 用1的n次单位复根，这样每次平方完，它们还是正负成对出现，

1. IFFT算法和FFT算法除了第6行以外一模一样

附视频链接：https://www.youtube.com/watch?v=h7apO7q16V0&feature=emb_logo

感兴趣的同学可以自行看看原视频，真的很通透，看完之后立马对傅里叶变换的理解上一个大台阶