# Polynomial GCD Derived through Monic Polynomial Subtractions

Authors: Feng Cheng Chang
Journal: ISRN Applied Mathematics
 (2011)

---

## Abstract

A simple and efficient algorithm for finding the GCD of a pair of univariate polynomials is derived.  The process requires only successive subtraction of two monic polynomials. Amazingly, this approach gives the desired GCD for a pair of test polynomials of very high degree, such as(x+1)∧1000
and its derivative.


---

## Body

## 1. Introduction

The computation of the greatest common divisor (GCD) of a pair of polynomials is an important issue in computational mathematics and also plays a crucial role in many science and engineering fields.  Several numerical methods have been proposed for the derivation of polynomial GCD [1–5]. Some are based on the classical Euclidean algorithm and the others are based on the procedures involving Sylvester resultant matrix. Most of them require algorithm of some advanced mathematics, such as the singular value decomposition and the least square iteration process.In this paper a simple and effective method to find the GCD of a pair of given polynomials is presented.  The process requires only the successive subtraction of two monic polynomials. It is modified purposely from the longhand polynomial divisions employed in the classical Euclidean GCD computation.  The entire computations require only simple elementary arithmetic operations, such as division and subtraction.A MATLAB code is provided, along with a typical numerical example.  Amazingly, this routine gives the GCD for a pair of test polynomials of very high degree.


## 2. Mathematical Formulation

From the classical Euclidean algorithm for computing polynomial GCD, the longhand polynomial division is consecutively employed:(2.1)pk(x)=pk-2(x)-pk-1(x)qk(x),where quotient qk(x) and remainder  pk(x) are obtained from dividing dividend pk-2(x) by divisor pk-1(x).  If the quotient qk(x) in every recurrent process can be converted into a numerical constant or even equal to 1, qk(x)=1, by making both pk-2(x) and pk-1(x) equal degree and monic, then the longhand polynomial division becomes simply a pair of “monic polynomials subtraction”: (2.2)pk(x)=pk-2(x)-pk-1(x).
The recurrent process of this approach may therefore be summarized as follows.Given a pair of polynomialsb(x) and a(x) of degrees n and m, (2.3)b(x)=∑j=0nbjxn-j=bc(x)xnz,a(x)=∑j=0majxm-j=ac(x)xmz,where b(x) and a(x) have nz and mz zero trailing coefficients, respectively, and bc(x) and ac(x) represent the polynomials without zero trailing coefficients, assuming that nc≥mc, where nc=deg(bc(x))andmc=deg(ac(x)).Before finding our desiredg(x)=gcd(b(x),a(x)), we first consider gc(x)=gcd(bc(x),ac(x)) by setting(2.4)p0(x)=bc(x)b0,p1(x)=ac(x)xnc-mca0.
Both polynomials p0(x)  and p1(x) are now in the same degree and monic. Applying the monic polynomial subtraction consecutively starting from k=2 until k=K+1, such that
(2.5)pK+1(x)=pK-1(x)-pK(x)=0,
we get(2.6)pK(x)=gcd(p0(x),p1(x))=gcd(bc(x),ac(x)xnc-mc)=gcd(bc(x),ac(x)).
Finally our desired polynomial GCD is found(2.7)g(x)=gcd(b(x),a(x))=pK(x)xmin⁡(nz,mz).It is noted that the complete set ofpk(x), referred to as “polynomial remainder sequence” [3, 4], may also be computed during the recurrent process.  The entire computations involve only elementary arithmetic operations without any advanced mathematics. The total numbers of recurrent process for computing GCD of the two given polynomials b(x) and a(x) are amazingly fewer than nc+mc.


## 3. Computer Routine in MATLAB

The computer codepolygcd.m in MATLAB is presented.  The inputs p and q are coefficient vectors of two given polynomials b(x) and a(x), and the output g is the desired greatest common divisor g(x). The input data can be either real or complex numbers. See Algorithm 1.Algorithm 1:
function g = polygcd(p,q)% *** Find GCD of two given polynomials ***%  through “Monic polynomial subtraction”%  by F C Chang    03/03/2011%mz = min(length(p)-max(find(p)),length(q)-max(find(q)));p0 = p(min(find(p)):max(find(p)));q0 = q(min(find(q)):max(find(q)));if min(length(p0),length(q0)) < 2,g = [1,zeros(1,mz)];    return,   end;g1 = p0/p0(1);              P{1}= g1;g2 = q0/q0(1);              P{2}= g2;for k = 3:length(p0) + length(q0),l12 = length(g1)-length(g2);  l21 = -l12;g3 = [g2,zeros(1,l12)]-[g1, zeros(1,l21)];g3 = g3(min(find(abs(g3)>1.e-8)):max(find(abs(g3)>1.e-8)));if norm(g3,inf)/norm(g2,inf) < 1.e-3,   break;  end;if length(g1) >= length(g2),g1 = g2;            end;g2 = g3/g3(1);              P{k}= g2;end;g = [g1,zeros(1,mz)];% celldisp(P);


## 4. Typical Example

If one test polynomialp(x) is derived from (4.1)p(x)=(x+3)(x2-5x+6)2(x5+3x4+8x3+8x2+7x+5)3x5(x-1)7,
that is,(4.2)p(x)=x32-5x31+2x30-6x29+76x28+140x27-802x26+954x25-4251x24+13663x23-18740x22+28472x21-53504x20+45776x19+5212x18-77580x17+185243x16-220631x15+104794x14+52458x13-193356x12+248612x11-146266x10+9202x9+65791x8-87555x7+55800x6-13500x5+0x4+0x3+0x2+0x+0,
and the another polynomial q(x) is the derivative of p(x), (4.3)q(x)=p′(x),
then, applying the routine polygcd.m,  the desired GCD of p(x) and q(x) is readily computed:(4.4)g(x)=gcd(p(x),p′(x))=x22-5x21+10x20-36x19+116x18-188x17+308x16-620x15+694x14-214x13-496x12+1348x11-1740x10+1012x9+28x8-692x7+929x6-605x5+150x4+0x3+0x2+0x+0.Finally it is fairly easy to check that(4.5)g(x)=(x2-5x+6)(x5+3x4+8x3+8x2+7x+5)2x4(x-1)6.See Algotithm2.Algorithm 2:
≫≫ p = [1   −5   2   −6   76   140   −802 954   −4251   13663   −18740   28472   −53504   45776   5212⋯−77580   185243   −220631   104794   52458   −193356   248612   −146266   9202   65791⋯−87555 55800 −13500  0  0  0  0  0];≫ q = polyder(p);≫ q = polygcd(p,q);≫ p, q, g,p =1     −5   2    −6        76  140     −802     954       −425113663      −18740  28472      −53504         45776     5212     −77580        185243     –220631104794         52458  −193356     248612    −146266     9202    65791      −87555    55800−13500   0   0   0     0        0q =32    −155     60     −174           2128            3780         −20852     23850   −102024314249  −412280      597912       −1070080    869744        93816  –1318860 2963888     –33094651467116     681954  −2320272       2734732  –1462660       82818         526328 −612885    334800−67500   0       0       0   0g =1       −5     10  −36       116 −188  308    −620  694−214       −496    1348 −1740      1012   28    −692      929          −605150     0       0   0   0≫


## 5. Applications

One of the most important applications of polynomial GCD is to find the roots with multiplicities of high-degree multiple-root polynomials [3].Let a univariate polynomialp(x) of degree N be given, and then the K roots zk with multiplicities mk are to be sought, such that(5.1)p(x)=∑i=0Nbixn-i=∏k=1K(x-zk)mk.We are now to find the rational functionr(x) or its related polynomials u(x) and v(x), such that(5.2)r(x)=u(x)v(x)=∑k=1Kmkx-zk,
where(5.3)u(x)=p(x)g(x),v(x)=p′(x)g(x),
and  g(x) is the polynomial GCD of p(x) and p′(x):(5.4)g(x)=gcd(p(x),p′(x))=∏k=1K(x-zk)mk-1.
Then, zk and mk  are simply the poles and residues of the rational function r(x).  Therefore, instead of directly solving for zk and mk  from the original N-degree multiple-root polynomial, we may now easily find the roots zk from the K-degree simple-root polynomial,(5.5)u(x)=∏k=1K(x-zk),
and the corresponding multiplicities mk  by the partial fraction expansion coefficients of the rational function  r(x)=v(x)/u(x),(5.6)mk=(x-zk)v(x)u(x)|x=zk=v(zk)u′(zk),A short program codepoly_roots.m for solving multiple-root polynomials is presented. It uses only basic built-in MATLAB routines except the function polygcd.m  for computing polynomial GCD. See Algorithm 3.Algorithm 3:
function Z = poly_roots(p)% *** Solve multiple-root polymonials ***%mz = length(p)-max(find(p));p0 = p(min(find(p)):max(find(p)));q0 = polyder(p0);if length(p0) < 2,  Z = [0,mz];  return, end;g = polygcd(p,q);u0 = deconv(p0,g0);v0 = deconv(q0,g0);w0 = polyder(u0);z0 = roots(u0);m0 = polyval(v0,z0)./polyval(w0,z0);Z = [z0,round(abs(m0))];if mz > 0,  Z = [Z; 0,mz];  end;


## 6. Conclusion

The very simple and efficient algorithm for the calculation of GCD of a pair of given polynomials is developed.  The entire computations involve only simple elementary arithmetic operations. Amazingly, the derived routine gives the desired results for the test polynomials of very high degrees, such as(x+1)1000,(x-123456789)30, and (1234x+56789)50, and their respective derivatives.One of the most important applications for polynomial GCD is to find the roots with multiplicities of polynomials with high degrees [1–3]. It reveals that the higher the root multiplicities of the given polynomial are, the more efficient this approach becomes. This is contrary to the usual experience that the most difficult part of solving for the roots of a polynomial is calculating those that have high multiplicities [1].



---

*Source: 714102-2011-06-13.xml*