Polynomial GCD Derived through Monic Polynomial Subtractions

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 (𝑥+1)^∧1000 and its derivative.

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(𝑥)𝑞_𝑘(𝑥),(2.1)where quotient 𝑞_𝑘(𝑥) and remainder 𝑝_𝑘(𝑥) are obtained from dividing dividend 𝑝_𝑘−2(𝑥) by divisor 𝑝_𝑘−1(𝑥). If the quotient 𝑞_𝑘(𝑥) in every recurrent process can be converted into a numerical constant or even equal to 1, 𝑞_𝑘(𝑥)=1, by making both 𝑝_𝑘−2(𝑥) and 𝑝_𝑘−1(𝑥) equal degree and monic, then the longhand polynomial division becomes simply a pair of “monic polynomials subtraction”: 𝑝_𝑘(𝑥)=𝑝_𝑘−2(𝑥)−𝑝_𝑘−1(𝑥).(2.2) The recurrent process of this approach may therefore be summarized as follows.

Given a pair of polynomials 𝑏(𝑥) and 𝑎(𝑥) of degrees 𝑛 and 𝑚, 𝑏(𝑥)=^𝑛_𝑗=0𝑎(𝑥)=_𝑚𝑏_𝑗𝑥^𝑛−𝑗=𝑏_𝑐(𝑥)𝑥^𝑛𝑧,_𝑗=0𝑎_𝑗𝑥^𝑚−𝑗=𝑎_𝑐(𝑥)𝑥^𝑚𝑧,(2.3)where 𝑏(𝑥) and 𝑎(𝑥) have 𝑛_𝑧 and 𝑚_𝑧 zero trailing coefficients, respectively, and 𝑏_𝑐(𝑥) and 𝑎_𝑐(𝑥) represent the polynomials without zero trailing coefficients, assuming that 𝑛_𝑐≥𝑚_𝑐, where 𝑛_𝑐=deg(𝑏_𝑐(𝑥))and𝑚_𝑐=deg(𝑎_𝑐(𝑥)).

Before finding our desired 𝑔(𝑥)=gcd(𝑏(𝑥),𝑎(𝑥)), we first consider 𝑔_𝑐(𝑥)=gcd(𝑏_𝑐(𝑥),𝑎_𝑐(𝑥)) by setting𝑝₀(𝑥)=𝑏^𝑐(𝑥)𝑏₀𝑝₁,(𝑥)=𝑎^𝑐(𝑥)𝑥^{𝑛𝑐−𝑚𝑐}𝑎₀.(2.4) Both polynomials 𝑝₀(𝑥) and 𝑝₁(𝑥) are now in the same degree and monic. Applying the monic polynomial subtraction consecutively starting from 𝑘=2 until 𝑘=𝐾+1, such that 𝑝_𝐾+1(𝑥)=𝑝_𝐾−1(𝑥)−𝑝_𝐾(𝑥)=0,(2.5) we get𝑝_𝐾(𝑥)==ggccdd𝑝𝑏_0𝑐(𝑥),𝑝₁(𝑥)=gcd𝑏_𝑐(𝑥),𝑎_𝑐(𝑥)𝑥^{𝑛𝑐−𝑚𝑐}(𝑥),𝑎_𝑐(𝑥).(2.6) Finally our desired polynomial GCD is found𝑔(𝑥)=gcd(𝑏(𝑥),𝑎(𝑥))=𝑝_𝐾(𝑥)𝑥^{min(𝑛𝑧,𝑚𝑧)}.(2.7)

It is noted that the complete set of 𝑝_𝑘(𝑥), 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 𝑏(𝑥) and 𝑎(𝑥) are amazingly fewer than 𝑛_𝑐+𝑚_𝑐.

3. Computer Routine in MATLAB

The computer code polygcd.m in MATLAB is presented. The inputs p and q are coefficient vectors of two given polynomials 𝑏(𝑥) and 𝑎(𝑥), and the output g is the desired greatest common divisor 𝑔(𝑥). The input data can be either real or complex numbers. See 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);

Algorithm 1  

4. Typical Example

If one test polynomial 𝑝(𝑥) is derived from 𝑝(𝑥)=(𝑥+3)𝑥²−5𝑥+6²𝑥⁵+3𝑥⁴+8𝑥³+8𝑥²+7𝑥+5³𝑥⁵(𝑥−1)⁷,(4.1) that is,𝑝(𝑥)=𝑥³²−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,(4.2) and the another polynomial 𝑞(𝑥) is the derivative of 𝑝(𝑥), 𝑞(𝑥)=𝑝^(𝑥),(4.3) then, applying the routine polygcd.m, the desired GCD of 𝑝(𝑥) and 𝑞(𝑥) is readily computed:𝑔(𝑥)==g𝑥c₂₂d𝑝(𝑥),𝑝^(𝑥)−5𝑥²¹+10𝑥²⁰−36𝑥¹⁹+116𝑥¹⁸−188𝑥¹⁷+308𝑥¹⁶−620𝑥¹⁵+694𝑥¹⁴−214𝑥¹³−496𝑥¹²+1348𝑥¹¹−1740𝑥¹⁰+1012𝑥⁹+28𝑥⁸−692𝑥⁷+929𝑥⁶−605𝑥⁵+150𝑥⁴+0𝑥³+0𝑥²+0𝑥+0.(4.4)

Finally it is fairly easy to check that𝑔(𝑥)=𝑥²−5𝑥+6𝑥⁵+3𝑥⁴+8𝑥³+8𝑥²+7𝑥+5²𝑥⁴(𝑥−1)⁶.(4.5)

See Algotithm 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       −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 =      32    −155     60     −174           2128            3780         −20852     23850   −102024     314249  −412280      597912       −1070080    869744        93816  –1318860 2963888     –3309465      1467116     681954  −2320272       2734732  –1462660       82818         526328 −612885    334800     −67500   0       0       0   0    g =     1       −5     10  −36       116 −188  308    −620  694       −214       −496    1348 −1740      1012   28    −692      929          −605    150     0       0   0   0 ≫

Algorithm 2  

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 polynomial 𝑝(𝑥) of degree 𝑁 be given, and then the 𝐾 roots 𝑧_𝑘 with multiplicities 𝑚_𝑘 are to be sought, such that𝑝(𝑥)=^𝑁_𝑖=0𝑏_𝑖𝑥^𝑛−𝑖=^𝐾_𝑘=1𝑥−𝑧_𝑘^𝑚𝑘.(5.1)

We are now to find the rational function 𝑟(𝑥) or its related polynomials 𝑢(𝑥) and 𝑣(𝑥), such that𝑟(𝑥)=𝑢(𝑥)𝑣(𝑥)=^𝐾_𝑘=1𝑚_𝑘𝑥−𝑧_𝑘,(5.2) where𝑢(𝑥)=𝑝(𝑥)𝑔(𝑥),𝑣(𝑥)=𝑝^(𝑥)𝑔(𝑥),(5.3) and 𝑔(𝑥) is the polynomial GCD of 𝑝(𝑥) and 𝑝^(𝑥):𝑔(𝑥)=gcd𝑝(𝑥),𝑝^(𝑥)=^𝐾_𝑘=1𝑥−𝑧_𝑘^𝑚𝑘−1.(5.4) Then, 𝑧_𝑘 and 𝑚_𝑘 are simply the poles and residues of the rational function 𝑟(𝑥). Therefore, instead of directly solving for 𝑧_𝑘 and 𝑚_𝑘 from the original 𝑁-degree multiple-root polynomial, we may now easily find the roots 𝑧_𝑘 from the 𝐾-degree simple-root polynomial,𝑢(𝑥)=^𝐾_𝑘=1𝑥−𝑧_𝑘,(5.5) and the corresponding multiplicities 𝑚_𝑘 by the partial fraction expansion coefficients of the rational function 𝑟(𝑥)=𝑣(𝑥)/𝑢(𝑥),𝑚_𝑘=𝑥−𝑧_𝑘𝑣(𝑥)𝑢(𝑥)||||_{𝑥=𝑧𝑘}=𝑣𝑧^𝑘𝑢^𝑧_𝑘,(5.6)

A short program code poly_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.

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;

Algorithm 3  

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 (𝑥+1)¹⁰⁰⁰,(𝑥−123456789)³⁰, and (1234𝑥+56789)⁵⁰, 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].