UOJ Logo EntropyIncreaser的博客

博客

HALF-GCD 算法的阐述

2020-03-10 17:48:53 By EntropyIncreaser

翻了各课件,发现……

某个英文网站 是这么说的:

There is a way to calculate the GCD and resultants in O(nlog2n). To do this you should note that if you consider a(x)=a0+xka1 and b(x)=b0+xkb1 where k=min(dega,degb)/2 then basically the first few operations of Euclidean algorithm on a(x) and b(x) are defined by the Euclidean algorithm on a1(x) and b1(x) for which you may also calculate GCD recursively and then somehow memorize linear transforms you made with them and apply it to a(x) and b(x) to lower the degrees of polynomials. Implementation of this algorithm seems pretty tedious and technical thus it's not considered in this article yet.

你如果脑子一冲,傻乎乎地去直接分治一半算出线性变换,DEBUG 一下会让你冷静过来…… 你发现这个线性变换应用到整体多项式上没有任何效果……

而 Picks 的 Introduction to Polynomials 也写的比较晦涩,因此在这里稍微详细地谈谈,以正视听……


给多项式 a,bF[z],我们原本的问题是求 GCD 以及扩展欧几里得的解:x,yF[z],ax+by=gcd(a,b)

由于多项式取模最坏情况下每次只让一个多项式的度数减小,这种方法的复杂度最好也就是 Θ((dega)(degb)) 的。

我们考虑求解一个中间情况叫做 HALF-GCD,设 k=min(dega,degb),这个东西让我们对于 a,b 能得到一个矩阵 (x1y1x2y2),其中 degxk2+o(k),且 deg(ax1+by1),deg(ax2+by2)k2+o(k)(除非它们的 gcddeg 比这个大)。也就是说这个 HALF-GCD 可以理解成是将欧几里得过程执行到一半的产物。

睿智的你想必已经意识到了这个东西如何得出真正的欧几里得的解。我们只需要调用一次这个算法,然后就规约到了一个 n2 的子问题(可能需要额外支付 O(M(n)) 时间用于调整较大的多项式),那么我们就得到了一个 Θ(M(n))+T(n)+T(n/2)+T(n/4)+ 的方法。

那么接下来考虑通过分治来解决 HALF-GCD:设 k=min(dega,degb),我们首先算出 a/zk/2,b/zk/2 对应的矩阵,这样一个矩阵作用与原多项式有两部分:考虑将原来的表为 a=a0+a1zk/2,b=b0+b1zk/2,那么 ax+by 中,deg((a1x+b1y)zk/2)=34k+o(k),另外 deg(a0x+b0y)=34k+o(k) 也比较显然。接着我们将得到的两个 34k 次多项式除掉前 zk/4 项送进去,这样又再次保证了递归返回的又是一个 14k+o(k) 量级的结果,我们将两个矩阵相乘就得到了变换矩阵。而结果正好满足度数在 k2+o(k) 以内。时间复杂度 T(n)=2T(n/2)+M(n),即 M(n)logn。消耗了 M(n) 的一方面是进行矩阵乘法,另一方面是可能在一些情况减少两个多项式次数的差,同时也可能发生剪枝。

综上所述,我们在 M(n)logn 的时间内完成了两个多项式的欧几里得求算,高精度整数的做法也是相似的。

这里 是一份代码,这里处理的时候是令 k=max(dega,degb),理论上是差不多的,目前尚且不知道这个代码是否存在一些 edge case 上的问题,欢迎大家指正。

评论

sto EntropyIncreaser
sto EntropyIncreaser
sto EntropyIncreaser

发表评论

可以用@mike来提到mike这个用户,mike会被高亮显示。如果你真的想打“@”这个字符,请用“@@”。