Download Polynomial GCDs using Linear Algebra: A Modern Approach and more Study notes Introduction to Business Management in PDF only on Docsity!
Polynomial GCDs by Linear Algebra
Barry Dayton
Northeastern Illinois University
www.neiu.edu/∼bhdayton
March 2004
This is a brief exposition on how to find the greatest common divisor of two univariate polynomials
using linear algebra rather than the Euclidean Algorithm. The motivation for this writeup is the
work of my colleague Z. Zeng. Here coefficients will be from the field F which will typcally be
the rationals, Q, the reals, R, or possibly the complex numbers, C.
Convolution Matrices The set PPP (^) n = {a 0 + a 1 x + a 2 x
2
n |aj ∈ F} of polynomials of
degree n or less forms a vector space. For a given polynomial p(x) = p 0 + p 1 x + · · · + pmx
m
of degree m, or less, multiplication of polynomials in PPP (^) n by p(x) gives a linear transformation
PPP (^) n → PPP (^) m+n because of the distributive law p(f + g) = pf + pg and the commutative law
p(rf ) = r(pf ) for constants r.
If we take as an ordered basis for PPP (^) n the list 1 , x, x
2 ,... , x
n , and likewise for PPP (^) m+n the convolution
matrix Cn(p) is the m + n + 1 × n + 1 matrix of the transformation f 7 → pf with respect to these
bases. For example, if p(x) = p 0 + p 1 x + p 2 x
2 and n = 3 then Cn(p) is given by
C 3 (p) =
p 0 0 0 0
p 1 p 0 0 0
p 2 p 1 p 0 0
0 p 2 p 1 p 0
0 0 p 2 p 1
0 0 0 p 2
Note that the rows correspond to the monomials 1 , x, x
2 , x
3 , x
4 , x
5 respecively and the first column
thus represents p, the second xp, the third x
2 p and the fourth x
3 p. If f = a 0 + a 1 x + a 2 x
2
3
then multiplying pf is just adding a 0 p + a 1 xp + a 2 x
2 p + a 3 x
3 p which means taking the appropriate
linear combination of columns, that is forming the matrix product
p 0 0 0 0
p 1 p 0 0 0
p 2 p 1 p 0 0
0 p 2 p 1 p 0
0 0 p 2 p 1
0 0 0 p 2
a 0
a 1
a 2
a 3
p 0 a 0
p 1 a 0 + p 0 a 1
p 2 a 0 + p 1 a 1 + p 0 a 2
p 2 a 1 + p 1 a 2 + p 0 a 3
p 2 a 2 + p 1 a 3
p 2 a 3
the result of which is the vector Coeff(f, m + n) of the coefficients of pf with respect to the basis
of PPP (^) m+n.
A variation on this is the division matrix where, in the example above, we augment with two
additional columns at the front with just a 1 in the diagonal spot. This makes a square upper
triangular matrix which will be non-singular assuming p 2 6 = 0.
Div 3 (p) =
1 0 p 0 0 0 0
0 1 p 1 p 0 0 0
0 0 p 2 p 1 p 0 0
0 0 0 p 2 p 1 p 0
0 0 0 0 p 2 p 1
0 0 0 0 0 p 2
We see that
Div 3 (p)
r 0
r 1
a 0
a 1
a 2
a 3
= Coeff(r 0 + r 1 x, m + n) + Coeff(pf, m + n))
If, instead, we replace the right hand side by Coeff(g, m + n) where g ∈ PPP (^) m+n and solve for the
vector on the left we see that we have simply performed the division algorithm to write the 5th
degree polynomial g as g = pf + r for some polynomial f of degree 3 and remainder r of degree
1 or less.
Generalizing this actually gives a proof of the division algorithm, since Divn(p) is non-singular if p
has exactly degree m and hence there is a unique solution. Moreover, this triangular system is easy
to solve by substitution, the work being equivalent to the usual division algorithm but more user
friendly to a computer. It should be remarked that if one is working over R or C and one knows
that the remainder should be zero then one can use the convolution matrix Cn(p) and solve as a
least square problem for a very accurate approximate solution. This was one of Zeng’s motivations
for replacing the division algorithm with this matrix formulation.
The Sylvester Matrix and the Resultant In a certain sense the use of the convolution matrix idea
to study the GCD goes back at least a century, but because of difficulty in actually performing
linear algebra computations without a computer, mathematicians only went half way, to decide if
gcd(f, g) = 1 without actually calculating the GCD.
The idea is now simple for us: given a polynomial f of degree exactly m and a polynomial g of
degree exactly n one forms a square m + n × m + n matrix Syl(f, g) = [Cn− 1 (f ) | Cm− 1 (g)]
formed by adjoining the two convolution matrices Cn− 1 (f ) and Cm− 1 (g). If we form a vector by
putting a coefficient vector Coeff(v, n − 1 ) of a, possibly unknown, polynomial v of degree n − 1
on top of a vector Coeff(w, m − 1 ) of coeffecients of a, also possibly unknown, polynomial w
of degree m − 1 then multiplying by the Sylvester matrix gives the coefficient vector of the sum
f v + gw.
The first follows simply by definition of nullspace, the second because if gcd(v, w) 6 = 1 then
dividing this GCD out would give v, w of smaller degree making the first equation true which
would have made an earlier Sj (f, g) rank deficient.
But the two properties above imply that v divides f and w divides g, moreover, dividing both sides
by −vw shows the quotients are equal, call this equal quotient u. u must be the GCD, it is certainly
a common divisor, and a larger common divisor would again contradict the minimality of degrees
of v, w. Using matrices, u can be found using linear algebra by the division matrices above.
If no such rank deficient Sj (f, g) is found by the time one gets to Sm− 1 (f, g) then, since Sm− 1 (f, g)
is square it is essentially the classical Sylvester matrix and the arguments above show gcd(f, g) =
Example We consider the same example as above: f (x) = x
3 +2x
2 − 4 x+1 and g(x) = x
2 − 3 x+2.
A trick is to shuffle the columns, that is to alternate columns from f and g rather then put all of
those from f first. We then have, writing, for example, S 0 (f, g) = S 0 ,
S 0 =
, S 1 =
, S 2 =
Using a computer we can just check the ranks of the matrices in order. By hand we could use
Gaussian elimination row operations directly on S 2 relying on the band matrix structure to avoid
doing actual work on the later columns unless we need to. We are looking for the first column
which is a linear combination of the earlier ones. If it is the second we note that we already
had this information in S 0 , the third or 4th was already in S 1 In this particular example we were
somewhat unlucky as it is the last column. Then completing Gauss-Jordan we can find a vector
[2, − 1 , − 1 , 3 , 0 , 1]
T in the nullspace. Since we shuffled the columns the coefficient vector for w is
[2, − 1 , 0]
T while that for v is [− 1 , 3 , 1] so −w = 2 − x and v = −1 + 3x + x
2
. Dividing by long
division or division matrices gives gcd(f, g) = u = x − 1.
The Extedned GCD. One advantange in this linear algebra approach is in finding the extended
GCD, that is find polynomials a(x), b(x) so that gcd(f, g) = a(x)f (x) + b(x)g(x). This is not
pleasant using the Euclidean Algorithm. Try it with the example above.
But using the linear algebra method it is fairly easy. The polynomials v, w obtained by the method
were seen to have gcd(v, w) = 1. Thus one can find the coefficients of v, w simply by using the
classical Sylvester matrix and solving
Syl(f, g)
a 0
. . .
an
b 0
. . .
bm
But multiplying av + bw = 1 by u gives af + bg = u as desired.
Example In the example above we saw v := −1 + 3x + x
2 and w = x − 2. So we need only solve
the 3 × 3 system
a 0
b 0
b 1
getting a = 1/ 9 , b = − 5 / 9 − 1 / 9 x and hence
f − (
x)g = x − 1
A final Word of Warning! When using linear algebra it is tempting to use decimal approximations
of numbers rather than exact fractions. Or you may wish to use a more sophisticated rank finding
method such as singular values. This latter method is the one used by Z. Zeng in his work. But
the result is then the approximate GCD which may not be the same as the exact GCD. One must
check to see if the computed GCD has the right kind of coefficients. For example, the fact the
GCD can be found using the Euclidean Algorithm or the method above using exact arithmetic says
that the GCD of polynomials with rational (fraction) coefficients also has rational coefficients. It
can be shown using number theory that the GCD of integer polynomials will always be an integer
polynomial, even if one is working in a larger coefficient field. The following example shows what
can happen.
Example Let
f = 225075x
5 − 20295
g = 88555x
4 − 12920
Applying Maple’s GCD algorithm gcd(f,g) one gets the constant polynomial 5, howver using
the command gcd(evalf(f),evalf(g)), which does the same problem using floating point
numbers, the answer is x − 0. 618034. But the number 0. 618034 is a known approximation to a
very irrational number, the reciprocal of the golden mean. So no multiple of this polynomial could
be an integer polynomial. Both answers are correct, however. The first answer is the exact GCD
as used in algebra and number theory, the second is the approximate GCD which might be used to
show that f, g have a common approximate real root of 0.618034.
