深入探究多项式乘法的快速算法

深入探究多项式乘法的快速算法焦作市第一中学闵梓轩一、高精度、多项式与生成函数1.1高精度在OI中我们有时会碰到一些问题的必要数值超出64位整形的范围，这个时候我们就需要用到高精度方式存储。而高精度数的思想是进制思想的一个具体体现，出于正常人类的习惯，我们所使用的高精度数都采用10进制，即每一位都表示十进制上的一个数，从0~9，更进一步，为了优化高精度数运算所花费的时间与空间，我们采用了万进制，即每一位存0~9999的数，这样同时优化了程序效率，同时在输出上也没有什么太大的问题（每一位不足1000补0即可）。当然，我们也可以用三进制、五进制、450进制，8964进制的高精度数，虽然因为在输出时会变得非常麻烦而没有人去用，但是它们的可行性正对应了进制的一种思想，比如一个十进制数12450，它的算数含义是0123410*010*510*410*210*1二进制数10010，它的算数含义是142*12*1（把为0的位忽略），这样形如),0(*0Naxaxaiiniii的每一位上的数字在数值表示上都乘上了某个数的一个幂的数正是进制思想的基础。在编程实现上这样的一个数我们通常用整形数组来表示，a[i]表示i次项的系数，如果数组长度为n，那么学过高精度的人都知道两个数相加的时间复杂度是θ(n)，两个数相乘的时间复杂度是O(n^2)，在信息学竞赛中，这样的时间复杂度足以满足大部分题目的需求，因为一般来说我们的数值都不会达到10^100000次方这么大。1.2多项式熟悉数学的我们能够发现上面这样的一个式子，如果忽略了括号中的内容的限制，那么我们可以发现这样的式子其实就是我们所学的n次多项式0*)(iiixaxA，比如十进制数12450就是05421234xxxx当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：220022010)(，那么)()()(D定义)()(，那么)()()(定义nkkikikinkkkjijiniiiixbaxbaxDxBxAxxbaxCxBxAxC注意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中有什么用呢？好像除了高精度之外没有什么卵用。但是熟悉组合数学的小伙伴们会了解一个叫做生成函数的东西，生成函数是什么呢？是0*iiixa，看到了木有，和多项式简直一模一样，而且如果把∞换成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）。那么两个多项式相乘就变成了nnnnxBCADxBDACxDCxBA*)(*)*)(*(2，此处乘以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，那么整个式子就变成了：nnnnxBDACDCBAxBDACxDCxBA*)))(((*)*)(*(2这样一来，我们只需要计算AC、BD、(A+B)(C+D)三次乘法就行了，然后进行加加减减，时间复杂度T(n)=3*T(n/2)+O(n)，解得：)()()3()(58.13loglog22nOnOOnTn这样的复杂度看起来可能有些玄学，那么它究竟是怎样的呢？举个栗子，当n=100000（一般来说十万这个数量级是O(nlogn)与O(n^2)的分界线）时，这个式子约等于0.8亿，仍然可以满足，但是已经接近极限了，再加上递归的巨大常数，是很勉强的。而当n=300000时，这个式子约等于4.5亿，已经爆掉了。虽然我们将多项式乘法通过这样的分治算法进行了优化，但仍然无法满足我们的需求。我们可以发现这样的分治没有达到对数级别，这意味着一个坏消息：我们还没做到位，以及一个好消息：我们可以做得更好。2.2另一种分治既然从中间分治不够完美，我们不妨退一步，从多项式的表达式中寻求规律。11332211011......)(nniniixaxaxaxaaxaxA从中间分治试过了，那如果我们按照奇偶性分开呢？))(...)()()(())(...)()()(()...()...()(12/2122512302112/222241221202145230122442200nnnnnnnnxaxaxaxaxxaxaxaxaxaxaxaxaxxaxaxaxaxA我们发现如果按照奇偶项分开，左右两边的形式都是一样的，而且都可以被表示成一个多项式的形式，也就是说，A(x)可以这样被表示：niiiRniiiLRLxxAxxAxxAxAxAi是奇数，i是偶数，22a)(,a)(，其中)()()(那么这样的分治有什么卵用呢？答：取值的时候有用。如果我们要将n个不同的x值代入到A(x)中，那么就会产生n个值，这样的操作叫做多项式的求值。求一个值的时间复杂度是O(n)（利用秦九韶算法），那么求n个值的时间复杂度就是O(n^2)，但我们可以从上面这个式子中发现这样一个性质：)()()(),()()(2222iRiiLiiRiiLixAxxAxAxAxxAxA这意味着什么？我们可以选择n/2个x值，把n/2个x^2代入到AL和AR中，然后根据上面的性质，就能求出n个A(xi)值了，这n个A(x)的值分别是n/2个x和n/2个-x代入A(x)的值，按照这样的算法，我们可以得出这样分治求值的复杂度：)log()()2/(*2)(nnOnOnTnT好的，现在我们成功地偏离了主题：说好的求乘积呢？不慌,不慌。三、变换与反演3.1另一种表示方式上面我们讨论到了从多项式的系数表达求出了值，也就是说我们可以快速地找到对于已知系数ai的多项式A(x)，n个不同的x值到n个取值的对应关系，即：))(),...,(),(),((),...,,,(),...,,,(121012101210nnnxAxAxAxAaaaaxxxx这像一个化学反应方程式，现在想：如果我们知道了左边的和右边的，我们是否可以求出中间的呢？也就是说如果我们知道有一个次数界为n的多项式，知道n个不同的x值以及它们代入所得的n个不同的取值，是否可以求出这个多项式的系数从而确定这个多项式呢？)),...,,,()))(,()),...,(,()),(,()),(,((1210???11221100nnnaaaaxAxxAxxAxxAx这看起来像个数学题，然而我们好像做过这样的数学题：已知二次函数f(x)经过三个点，求f(x)的表达式。这就变得有意思了吗？上面那个例子正是知道(x0,f(x0)),(x1,f(x1)),(x2,f(x3))三个点，求一个次数界为2+1=3的多项式的系数。数学老师告诉我们：这是可行的，只不过可能会有些难算。但这没什么，毕竟计算机的一个优势就在于计算速度快，现在我们关心的是这样的一个性质是否对于任意一个n都成立呢？要想严格地证明这个问题，我们就要先严格地描述这个问题。如果我们用yi这样一个直观的东西来表示A(xi)的话，上面的求值可以表示为：为矩阵）、、（，即Y110110112111011121110110201000AXYXAyyyaaaxxxxxxxxxxxxnnnnnnnnn