深入探究多项式乘法的快速算法.docx
《深入探究多项式乘法的快速算法.docx》由会员分享，可在线阅读，更多相关《深入探究多项式乘法的快速算法.docx（20页珍藏版）》请在冰点文库上搜索。
深入探究多项式乘法的快速算法
深入探究多项式乘法的快速算法
焦作市第一中学闵梓轩
一、高精度、多项式与生成函数
1.1高精度
在OI中我们有时会碰到一些问题的必要数值超出64位整形的范围，这个时候我们就需要用到高精度方式存储。
而高精度数的思想是进制思想的一个具体体现，出于正常人类的习惯，我们所使用的高精度数都采用10进制，即每一位都表示十进制上的一个数，从0~9，更进一步，为了优化高精度数运算所花费的时间与空间，我们采用了万进制，即每一位存0~9999的数，这样同时优化了程序效率，同时在输出上也没有什么太大的问题（每一位不足1000补0即可）。
当然，我们也可以用三进制、五进制、450进制，8964进制的高精度数，虽然因为在输出时会变得非常麻烦而没有人去用，但是它们的可行性正对应了进制的一种思想，比如一个十进制数12450，它的算数含义是
二进制数10010，它的算数含义是
（把为0的位忽略），这样形如
的每一位上的数字在数值表示上都乘上了某个数的一个幂的数正是进制思想的基础。
在编程实现上这样的一个数我们通常用整形数组来表示，a[i]表示i次项的系数，如果数组长度为n，那么学过高精度的人都知道两个数相加的时间复杂度是θ（n），两个数相乘的时间复杂度是O（n^2），在信息学竞赛中，这样的时间复杂度足以满足大部分题目的需求，因为一般来说我们的数值都不会达到10^100000次方这么大。
1.2多项式
熟悉数学的我们能够发现上面这样的一个式子，如果忽略了括号中的内容的限制，那么我们可以发现这样的式子其实就是我们所学的n次多项式
，比如十进制数12450就是
当x=10的时候的数值嘛。
所以，当一个值b代入多项式A（x）时，这个式子也就变成了一个值A（b）。
但是要注意的是多项式的系数是没有限制的，所以多项式可以用浮点数组表示，而且我们可以惊奇地发现多项式的加法和乘法在代码上除了不需要进位之外和高精度是一样的。
所以说，我们所见的b进制数值，就是一个当x=b的多项式的取值而已。
但是在多项式中，x的意义仅仅是一个符号而已，ai*x^i你可以理解为ai在数组的第i个位置。
我们需要注意的是，n次多项式的数组表示需要用到n+1个数，为什么？
因为有n个含x的项和一个常数项，所以我们一般把多项式A（x）的最高次项的次数+1称作为这个多项式的次数界（次数界的真正意义是系数不为零的最高次项的次数+1，下文中提到的“次数界“为了更清晰的算法表述，将其定义为按需求扩展后的数组长度，此时最高次项可能是为零的。
）。
所以，存储一个次数界为n的多项式只需要开一个长度为n的数组就行了。
现在我们需要严格地定义多项式的运算，以下假定A、B的次数界分别为n、m：
注意D的次数界：
两个次数界为n和m的多项式相乘时，积的次数界为n+m-1，D被称为A和B的卷积。
这样一来求和就是从0到n-1枚举i，将ai和bi相加的结果作为ci即可，时间复杂度θ（max（n,m））。
求积就是从0到n-1两层枚举i和j，将ai和bj相乘的结果加到d（i+k）上，时间复杂度θ（nm）。
说白了就是不用进位的高精度乘法。
1.3生成函数
然而多项式在OI中有什么用呢？
好像除了高精度之外没有什么卵用。
但是熟悉组合数学的小伙伴们会了解一个叫做生成函数的东西，生成函数是什么呢？
是
，看到了木有，和多项式简直一模一样，而且如果把∞换成n，那这个式子不就真的变成多项式了吗？
在组合数学中，数列A的生成函数的i次项系数ai就是这个数列的第i项。
然而我们在生成函数解决问题时是用得到乘法的，有人可能会说这TM有无穷项你乘个卵？
但是我们真正关心的没有那么多项，比如我们可能会只关注这个数列的第n项，这个时候第n+1项以及往后的数值我们就没有必要再乘了（相当于对x^（n+1）取模，同余原理懂吧？
懂吧？
就像组合计数时要求你对一个大质数取模避免高精度一样），这个时候，生成函数的乘法就变成了多项式乘法！
O（n^2）的乘法或许对于高精度来说足够了，但是对于只需要算常数次乘法，并且只关心其中一项的生成函数来说，你需要计算的，可能有很多位。
比如我们需要计算一个一个数列的第n项，也就是这个数列生成函数的第n位，那么如果当n=300000的时候你的乘法就会瞬间爆炸，那么，有没有更快的算法呢？
二、分而治之
分治思想对于计算机科学而言就像心脏，而数组这样的形式是很容易引起我们的分治欲望的，我们不如来看看是否能够通过对数组进行分治而来实现优化。
2.1裂项相乘
对于数组来说，最基本的分治就是从中分开了，一般的数组分治就是取左端点右端点加起来除以二省略余数作为分治的界限，但是如果数组的长度正好是偶数的话思路分治算法会更加方便一些，直接从中间分开就行了，两边长度一样的。
我们注意到，生成函数如果只有有限项不为0的话，那么它就是多项式，如果我们反过来想，这个性质就变成了可以通过往多项式最高位后面补0来无限地增加数组的长度。
比如我们如果有249个项，可以往最高位补个0来让它变成250个项，这不就容易分治了吗？
那么现在关键是如何分治，数组分治就是把数组分成连续的多份嘛，而由于很多原因（比如基于二进制的位运算常数很小，二分比多分编程实现更容易而且渐进复杂度是一样的）我们通常采用二分，也就是把数组一分为二，此时由于上面的性质我们不妨假设我们所需要面对的多项式已经有偶数项了，我们设它为2n，然后把它从中间（当然有的基于数组的分治可能不是从正中间，比如下面我们要讨论的）分为两个有n项的数。
由于运算是两个多项式的事情，所以我们需要另一个数组，把第一个数组分成A和B，第二个数组分成C和D。
这里为了画图方便，n=5。
那么我们该如何分治呢？
我们可以注意到第一个多项式的数组表示是B数组右移n位和A数组拼接起来的，那么这在多项式中意味着什么呢？
第一个多项式可以表示成（A+B*x^n），不是吗？
那么第二个多项式就可以表示成（C+D*x^n）。
那么两个多项式相乘就变成了
，此处乘以x的n次方只是表示将这个多项式乘出来之后往右移n位。
这样的话我们需要计算AC、BD、AD、BC四组多项式的乘积并进行相加，需要注意的是如果要继续对这几组多项式乘法施展递归的方法，我们需要n也是偶数，而继续递归下去就会发现我们最初的数组长度必须是2的幂，但这没什么，来人，给我添0。
所以，在以下的叙述中我们都假定n是2的幂假定n是2的幂假定n是2的幂（因为很重要所以说三遍）算一下这个算法的时间复杂度T（n）=4*T（n/2）+O（n），主定理解递归式得T（n）=O（n^2）（主定理是什么自己查），嗯，好像并没有什么卵用，常数上还因为递归的缘故比原来大了，但这并不代表我们的思路没有卵用。
现在考虑是否可以通过将式子变形来少计算几次乘积呢？
这里AC和BD是必须要计算出来的乘积（思考为什么），所以现在考虑（AD+BC）的变形，由于我们已经计算出了AC和BD，那么是否可以直接把这个结果利用上呢？
先相加起来AC+BD+AD+BC=（A+B）（C+D），所以AD+BC=（A+B）（C+D）-AC-BD，那么整个式子就变成了：
这样一来，我们只需要计算AC、BD、（A+B）（C+D）三次乘法就行了，然后进行加加减减，时间复杂度T（n）=3*T（n/2）+O（n），解得：
这样的复杂度看起来可能有些玄学，那么它究竟是怎样的呢？
举个栗子，当n=100000（一般来说十万这个数量级是O（nlogn）与O（n^2）的分界线）时，这个式子约等于0.8亿，仍然可以满足，但是已经接近极限了，再加上递归的巨大常数，是很勉强的。
而当n=300000时，这个式子约等于4.5亿，已经爆掉了。
虽然我们将多项式乘法通过这样的分治算法进行了优化，但仍然无法满足我们的需求。
我们可以发现这样的分治没有达到对数级别，这意味着一个坏消息：
我们还没做到位，以及一个好消息：
我们可以做得更好。
2.2另一种分治
既然从中间分治不够完美，我们不妨退一步，从多项式的表达式中寻求规律。
从中间分治试过了，那如果我们按照奇偶性分开呢？
我们发现如果按照奇偶项分开，左右两边的形式都是一样的，而且都可以被表示成一个多项式的形式，也就是说，A（x）可以这样被表示：
那么这样的分治有什么卵用呢？
答：
取值的时候有用。
如果我们要将n个不同的x值代入到A（x）中，那么就会产生n个值，这样的操作叫做多项式的求值。
求一个值的时间复杂度是O（n）（利用秦九韶算法），那么求n个值的时间复杂度就是O（n^2），但我们可以从上面这个式子中发现这样一个性质：
这意味着什么？
我们可以选择n/2个x值，把n/2个x^2代入到AL和AR中，然后根据上面的性质，就能求出n个A（xi）值了，这n个A（x）的值分别是n/2个x和n/2个-x代入A（x）的值，按照这样的算法，我们可以得出这样分治求值的复杂度：
好的，现在我们成功地偏离了主题：
说好的求乘积呢？
不慌,不慌。
三、变换与反演
3.1另一种表示方式
上面我们讨论到了从多项式的系数表达求出了值，也就是说我们可以快速地找到对于已知系数ai的多项式A（x），n个不同的x值到n个取值的对应关系，即：
这像一个化学反应方程式，现在想：
如果我们知道了左边的和右边的，我们是否可以求出中间的呢？
也就是说如果我们知道有一个次数界为n的多项式，知道n个不同的x值以及它们代入所得的n个不同的取值，是否可以求出这个多项式的系数从而确定这个多项式呢？
这看起来像个数学题，然而我们好像做过这样的数学题：
已知二次函数f（x）经过三个点，求f（x）的表达式。
这就变得有意思了吗？
上面那个例子正是知道（x0,f（x0））,（x1,f（x1））,（x2,f（x3））三个点，求一个次数界为2+1=3的多项式的系数。
数学老师告诉我们：
这是可行的，只不过可能会有些难算。
但这没什么，毕竟计算机的一个优势就在于计算速度快，现在我们关心的是这样的一个性质是否对于任意一个n都成立呢？
要想严格地证明这个问题，我们就要先严格地描述这个问题。
如果我们用yi这样一个直观的东西来表示A（xi）的话，上面的求值可以表示为：
上面所叙述的求值操作就是在已知X和A的情况下求Y，那么我们要求的问题就是已知X和Y来求A，那这个问题就很简单了，我们只需要求出来X的逆矩阵就行了，这样
。
不幸的是，不是每一个矩阵都有逆矩阵的，但幸运的是，X是必然有逆矩阵的。
像X这样有个性的矩阵当然得有个自己的名字了，它叫做范德蒙矩阵，而它的行列式是这样的：
至于如何推导……我也不知道
可以看到，如果n个xi都不相同的话，范德蒙矩阵的行列式就不为零，而我们知道矩阵有逆当且仅当其行列式不为零。
现在，我们知道了，一个次数界为n的多项式的系数与其在n个不同x值处的取值存在一一对应的关系：
既然存在这样一个一一对应的关系，那么我们就可以用值而不是系数的形式表示一个多项式，也就是说现在多项式有两种表示方式，一种是箭头左边的用n个系数来表示的系数表示，另一种是箭头右边的用n个不同x值以及其取值的点值表示。
点值表示，是系数表示的一种变换，转化的操作叫做求值;而从点值表示到系数表示，即是反演，我们将其叫做插值。
点值表示同样需要长度与次数界相同的数组，不过是两个，一个是x[]一个是y[]，虽然经过后面的探究之后我们几乎没可能这么用，但这样会有助于我们的理解。
现在我们新发现了一种与系数表示截然不同的表示方法，现在我们不妨开始研究它的性质，既然本文的目的是研究快速算法，那么我们不妨可以看看点值表达的多项式是如何进行运算的，以下给定次数界为n的多项式A和B（即我们已知其系数表示或点值表示）之间的运算关系：
我们发现，如果有A和B的点值表示，那么我们就可以求出在Θ（n）的时间内求出它们的和的点值表示，甚至它们的积的点值表示，这不就是我们梦寐以求的快速乘法吗？
让我们梳理一下思路：
研究这么半天，我们发现了一条捷径，离成功只有一步之遥了，让我们来看看插值的运算吧：
插值有两步：
第一步矩阵求逆，普通的算法Θ（n^3），对于范德蒙矩阵，更快的算法有Θ（n^2）。
第二步逆求值，因为逆矩阵并不一定是范德蒙矩阵的形式，所以我们不能用上文分治的方法Θ（nlogn）的时间逆求值，此时进行插值的复杂度为Θ（n^2），这样一来插值的时间复杂度为：
Θ（n^3）。
陷入困境了呢。
3.2问题的发现
我们一步一步地后退，审视我们在思考的过程中是否漏掉了什么东西。
好消息是我们少考虑到一个“问题”，坏消息是它是一个“problem”：
在求值的过程中会涉及到一个数的i次方，而i是可能会达到100000的，这就造成了一个灾难：
数值的精度问题。
如果我们用整数代入的话，按照指数函数的增长速度来看，那么用不了几次方就会爆掉64位长整形，需要借助高精度，但问题是我们研究的就是特么的高精度的变种，这不就是个没用边界的递归，死循环了吗。
如果我们用浮点数的话，那么这个浮点数的精度……十万次方，看着办吧。
在OI中，每当我们所需要计算的结果会很大时，会让我们把结果对一个大质数取模，从而避免高精度而关注算法本身。
注意，这样的问题通常都是组合数学问题，这跟我们要解决的问题所属的范围恰好一致。
如果在求值的时候运用取模来把乘法降为O
（1）就可以解决这个灾难了，但是问题在于取模的话运算结果是否正确呢？
并不一定。
而且求逆的问题依然没有得到解决。
该怎么办呢？
实际上，到最后我们确实会将这个算法置于模运算的意义之上，但这并不是我们现在要讨论的问题。
模运算是在循环群上进行的，循环这个词给了我们一个解决精度的思路，但以普通高中生的知识层面是无法沿着这个思路解决问题的，不过不虚，我们不是普通高中生，我们是OIER。
所以我们现在需要做的是放下对算法的研究，补充基础知识。
这个过程是很生硬的知识传输，但我会尽量以一种亲民的方法来讲。
学习新知识是很痛苦的，但是学会了会有一种豁然开朗的欣喜和自豪，这在自然科学这条道路上行走时是无法避免的。
3.3复杂的数值问题要用复杂的数解决
为了解决这个复杂的问题（complexproblem）接下来我们要以毒攻毒（好像毒的都是自己），介绍一个超出常识的概念：
复数（complexnumber）。
复数的定义是这样的：
复数c（complex的首字母）=a+bi，其中a和b均属于实数集，而i也不是程序员们平时所说的循环变量，在这里i被赋予了一个超越常识的意义：
。
因为c是由两个数组成的，所以我们也可以把复数c=a+bi表示成（a,b）。
这样像不像一个坐标表示？
从某种意义上来说，是的。
既然知道了复数集这样一个集合，那么我们需要研究一下它们的运算关系：
加法和乘法，减法和除法实际上就是加（乘）上另一个数的加（乘）法逆元：
有了这样的运算法则，我们就可以在计算机上对复数进行运算了，一个复数可以定义为一个含有两个double的结构体，然后按上面的式子来定义它们的加减乘除函数（最直观的方法还是重载运算符）即可。
不过，在c++中内置了复数类complex，它定义好了加减乘除法，但具体要用时请参考别人的代码，使用时我们需要含入库，需要使用时声明complexc;即可（嫌麻烦的可以用typedef）。
注意，库里内置的东西在评测时不开-O2常数是很大的，最好还是尝试自己写一下。
（C++实现见附录代码）
我们发现复数的加法和向量的加法是一模一样的，这是否能启发到我们的乘法呢？
不慌，我们再来看以下的内容：
复数的指数形式定义是这样的：
，其中e为自然常数。
为了阅读方便我们定义g=e^i（在其它书籍或资料中没有这样的定义，g只是我为了方便大家理解胡乱起的），这样的话
。
那么我们就可以发现这样：
感觉好爽啊。
这样的话我们就可以把复数表示成向量了，而它们的相乘就是倾斜角相加！
或许你可能在哪里听到过世界上最美的数学公式：
，现在我们可以解释它了，而且可以按我们的思路发现
，这些性质来自于我们学过骂过的三角函数的各种性质。
下面我们再来讨论真正用得着的东西：
单位复数根。
同样超越常识的它的定义是这样的：
在超越常识的复数域内，存在有数，这个数的n次方等于1，而且这样的数不止一个，而是有n个，分别为
，其中
被称为主n次单位复数根，被简写为
。
按照定义，我们可以直接找到这n个数，看图：
（截自第四版的算导）
这张图已经很清晰了吧？
那么我就直接给公式了：
。
没错，就是将一个单位圆分成n份取其中的第k份。
接下来介绍一个重要的定理：
这个定理在算法导论上被称消去引理，由此我们可以联想到，如果n是2的幂的话，那么将
平方其实就相当于将n除以2。
回头想一想，我们到现在为止利用到了复数的哪些性质解决了哪些问题？
我们遇到的问题有二：
数值的精度问题，以及插值运算的问题。
现在看来，数值的精度问题已经得到了完美地解决，因为运用单位复数根的幂指数是循环的，这是因为
。
这样一来，只要我们取
代入到多项式中求值就可以了。
3.4快速傅里叶变换
快速傅里叶变换（fastFouriertransform,FFT），这就是本文所要介绍的算法真正的名字，它所蕴含的思想正是上文所介绍的运用分治，在Θ（nlogn）的时间内完成n个不同的值的求值运算，并采用n个n次单位复数根作为n个值代入多项式来解决精度问题。
而且如果给定n的话，n个单位复数根就确定了，那么这个多项式的点值表示就仅用一个数组就可以做到，所以我们在计算过程中直接把系数数组变换为点值数组就行了。
具体做法我用中文伪代码来解释：
（c++实现见附录代码）
我们来用归纳法来证明，对于次数界为n的复数数组（实数其实也是复数）a，以及它的长度n（n为2的幂），该程序在运行之后能够正确地计算出a的傅里叶变换，也就是将多项式的系数数组a输入并运行之后，a数组就从系数数组变为数值数组，其中
。
程序第2、3行给出了递归边界，当n=1时，相当于多项式只有一个常数项，那么不管代入什么值它的值都不会变，所以直接结束递归调用即可。
当n=2^q（q>1）时，程序第6~9行创建aL和aR数组并对它们进行递归调用求值。
第4~5行为初始化单位复数根，程序第10~13行是通过分治的方法，利用已得出的两个子数组的值，进行求值，其中第13行的乘法保证了在第k次循环执行时，
。
看不懂的，请自行回头看上文介绍的分治求值方法。
好了，现在我们已经成功地解决了利用单位复数根进行求值的问题，至于两个点值表示相乘，那都不会你就来看这个？
3.5快速傅里叶逆变换（我一般习惯把逆变换叫做反演）
现在我们就差一步了，那就是插值，我们需要做到对单位复数根的范德蒙矩阵求并逆求值，看起来一点也不简单。
不急，我们再来介绍一个性质：
这是一个很好理解的引理，因为当n和k互质时，n个j和k相乘会得到n个不同的值（这涉及到数论的理论），根据循环的性质，这n个值也就是n个不同的n次单位根，而上面的那个单位圆的图告诉了我们一个很明显的事实，n个n次单位根的和为0。
至于另一种更严格的证明，参见《算法导论》。
好了，矩阵乘法是啥？
，现在我们知道
（中括号为艾弗森约定，当括号内的表达式为真时取值为1，否则为0），现在要做的是求b，既然原矩阵的取值是n次单位复数根，而这n个在乘法意义下是循环的，那我们为啥不假设b也是n次单位复数根？
我们先假设
，此时
，这个形式和求和引理很相近了，我们可以发现
。
而我们要求的是
。
那么我们就可以发现，只需要令上面c’那个式子忠的i+j变为i-j，就能根据求和引理满足后半句了，此时前半句的式子依然为n，但这并不麻烦，我们只需要将其除以n就行了，所以我们得出
，此时我们再进行验证：
。
完美。
那么，逆矩阵求出来了，而且逆矩阵的形式和一个范德蒙矩阵接近，逆求值呢？
逆求值和求值多么相近，我们只需要对y数组运行快速傅里叶变换，并且把单位复数根取反，最后得出的结果再除以n，那么我们就能够得出a！
现在，我们对上面的程序进行改进，增加一个参数，使得f为1时对a执行傅里叶变换，当f为-1时对a执行傅里叶逆变换：
新的程序比原来的程序多出了14~16行，这是因为逆变换后数组元素要除以n。
并且请注意第四行，这么改的原理是这样的：
大功告成了，让我们来用伪代码梳理一下整个过程：
貌似真的很简单，就和上面的框图解释的一样：
算法第2~3行为求值，第4~5行为相乘，第6行为插值，求出成绩c并在第7行返回结果。
需要注意的是，长度n是不小于c的长度的2的幂。
我们需要枚举出最小的大于等于a,b长度和的2的幂，当然大点也无所谓，只要你不介意让自己的程序跑得更慢一些。
3.6优美的蝶代
作为oier，发现一个递归的算法之后，首先要想到能不能将其转化为迭代的形式，因为迭代的实现往往比迭代的要搞笑。
首先我们从递归树入手，首先来看当n=8时的情形：
在一层层的递归一直到底层之后，最后数组元素所在的位置变换为了它二进制翻转后的位置。
比如原来的下标为4（100）的元素，变到了第1（001）的位置，而0（000）、2（010）、5（101）、7（111）由于它们的二进制是回文串，所以递归后位置并没有发生改变。
该如何解释这一神奇的现象呢？
其实并不难，在第一次递归调用的时候，我们将原数组按照奇偶划分为了左右两个部分，偶数在左，奇数在右，而奇偶的判断依据是二进制末位是否为0，这也就相当于将原数组按照二进制末位为第一关键字重新排序。
而进行递归调用时将n除以2，也代表着在以后的划分中将最末尾抹去不参与了（因为都一样），所以，倒数第二位，也就相当于第二关键字。
在我们平时的排序中，是按照二进制最高位（第一位）为第一关键字，第二位为第二关键字以此类推进行排序的，而这样的递归划分是从最末尾倒着来的，这也就很好地解释了为什么递归进行完之后所有元素的下标都进行了二进制翻转。
我们既然知道了这样一个规律，那么就可以直接将数组元素进行二进制翻转，然后从底层向上进行求值即可。
将二进制翻转的操作，我们叫做蝶形变换，因为它真实的递归树像蝴蝶一样美（当然文科生无法理解这种数学的美）。
中文伪代码如下：
循环5~14行枚举数组长度从底向上迭代（从2开始是因为最底一行不用计算），第7~13行枚举本层各个区间的起点，第9~12行和递归版本一样进行计算。
现在该解决的问题都解决了，甚至我们还对这个算法进行了迭代的优化，到这里为止快速傅里叶变换就算基本结束了。
但是上文还是留下了一个问题。
3.6*快速数论变换
还记得本文开头所说的研究多项式快速乘法的目的是什么吗？
并不是解决高精度问题，而是要解决组合数学中生成函数的乘法问题，而这一类问题往往由于结果过大是需要对一个大质数取模的，这个大质数往往会是2的幂的倍数+1。
当遇到这类问题时，直接进行复数运算最后再进行取模的话，可能会承受无法估计的精度问题所带来的不确定性，那么，对此我们是否有解决的办法呢？
当然有，因为在我们提出利用复数来解决幂次方的精度问题之前，便想到了整数群的循环性质，而单位复数根的运算恰好运用了这样的性质。
现在我们研究完了基于复数的快速傅里叶变换，是否可以回头看看整数群也拥有这些性质呢？
快速傅里叶变换核心的核心在于复数群有n个n次单位复数根，而且它们之间的运算满足循环的性质。
那么整数群呢？
是否也有n个不同的n次方为1的数呢？
由这个性质，我们可以知道g的p个幂次方正好对应了p个不同的数，再看上面的那个式子，我们再来找n次方为1的数不就简单了吗？
我们只需令主n次单位根
。
由于n是2的幂，所以很显然n|（p-1）。
我们不妨再来验证一下整数群上n次单位根的性质：
很好，这样，看来我们只是正确地取了
，就完美地解决了一切问题，发现了一个新的算法。
《算法导论》上这种算法被称为基于模运算的FF
升级会员 