TAOCP 4.6.1 Exercise 20

Section 4.6.1: Division of Polynomials

Exercise 20. [M40] Investigate approximate polynomial gcds and the accuracy of Euclid's algorithm: What can be said about calculation of the greatest common divisor of polynomials whose coefficients are floating point numbers?

Verified: no
Solve time: 3m34s

Setup

The exercise asks for an investigation rather than a theorem with a single conclusion. Let

$u(x),v(x)\in \mathbf{R}[x]$

have coefficients represented in floating-point arithmetic. We seek to understand:

Whether Euclid's algorithm remains reliable when the coefficients are only approximate.
How the notion of greatest common divisor should be interpreted in the presence of rounding errors.
What algorithms are appropriate for computing a meaningful gcd.

The exact Euclidean algorithm repeatedly forms remainders

$$ \deg r_{k+1}<\deg r_k, $$

and terminates when some remainder is zero. The last nonzero remainder is a gcd.

The difficulty is that floating-point coefficients are not exact elements of a unique factorization domain or field; they are approximations to real numbers.

Solution

The classical Euclidean algorithm is numerically unstable for polynomial gcd computation.

Assume that two exact polynomials possess a nontrivial common divisor,

$$ v(x)=g(x)b(x), $$

with $\deg g>0$. If the coefficients are perturbed slightly, producing

$$ \tilde v(x)=v(x)+\delta v(x), $$

then the common divisor generally disappears. The set of polynomial pairs having a nonconstant gcd is determined by the vanishing of certain polynomial expressions in the coefficients, namely the appropriate resultants. Such conditions define a lower-dimensional subset of coefficient space. An arbitrarily small generic perturbation therefore moves the coefficient vector away from that subset.

Consequently, two floating-point polynomials almost always have exact gcd equal to $1$, even when they arise from exact polynomials with a large common factor.

A simple example illustrates the phenomenon. Consider

$$ v(x)=(x-1)(x-3). $$

Their gcd is $x-1$. Replace $u$ by

$$$$

Then $\tilde u$ and $v$ have no common root, hence their exact gcd over the real numbers is $1$. Euclid's algorithm applied exactly to the perturbed coefficients returns a constant polynomial, although the perturbation is extremely small.

Therefore the meaningful problem is not exact gcd computation but rather computation of an approximate gcd. One seeks a polynomial $g(x)$ of largest possible degree for which there exist small perturbations

$$$$

such that

$$$$

while the coefficients of $\Delta u$ and $\Delta v$ remain small according to a chosen norm.

The behavior of Euclid's algorithm reflects this difficulty. In exact arithmetic a common divisor of degree $d$ causes the remainder sequence to terminate after the degree has dropped by $d$. Under floating-point perturbations the theoretically zero remainder is replaced by a polynomial whose coefficients are small but nonzero. Subsequent divisions magnify relative errors because each step divides by leading coefficients and subtracts nearly equal quantities. Tiny perturbations may therefore produce large changes in later remainders.

A practical criterion is obtained by examining the sizes of the remainders. Suppose the exact problem would yield

$$$$

In floating-point arithmetic one obtains instead

$$$$

whose coefficients are of magnitude comparable to the accumulated rounding error. If

$$$$

is much smaller than

$$$$

the computation suggests the existence of an approximate common divisor. The degree of that divisor is inferred from the location where the dramatic decrease occurs.

A more reliable approach uses the Sylvester matrix. For polynomials

$$ \qquad v(x)=v_nx^n+\cdots+v_0, $$

the existence of a common divisor of degree at least $d$ is equivalent to a rank deficiency in an appropriate Sylvester matrix. With floating-point data the matrix is seldom exactly rank deficient, but its singular values reveal near-dependencies. A cluster of very small singular values indicates that only a small perturbation of the coefficients is needed to create a common divisor. Numerical linear algebra therefore provides a stable framework for approximate gcd computation.

The principal conclusions are:

Exact polynomial gcds are not continuous functions of the coefficients.
Arbitrarily small perturbations usually destroy every nonconstant common divisor.
Euclid's algorithm, used naively with floating-point coefficients, frequently returns a constant gcd even when a substantial common factor is present in the underlying exact data.
Approximate gcds should be defined through nearby polynomials having an exact common divisor.
Stable algorithms are based on rank-revealing methods, particularly singular-value analysis of Sylvester matrices, rather than on direct interpretation of the final remainders of Euclid's algorithm.

Verification

Let

$$$$

The exact gcd is $x-1$.

Now replace $u$ by

$$$$

The roots of $\tilde u$ are $1.000001$ and $2$, while the roots of $v$ are $1$ and $3$. No root is common. Hence every common divisor of $\tilde u$ and $v$ has degree $0$, so

$$$$

The coefficient perturbation is of order $10^{-6}$, yet the gcd degree changes from $1$ to $0$. This demonstrates the discontinuity of exact gcd computation and explains the sensitivity of Euclid's algorithm.

Notes

Approximate polynomial gcds form an important topic in computer algebra and numerical analysis. Modern algorithms typically combine polynomial algebra with matrix computations. The central problem is often formulated as follows: Given floating-point polynomials $u$ and $v$, find the smallest perturbations that produce nearby polynomials having a common divisor of prescribed degree. The singular values of Sylvester matrices provide both a diagnostic tool and a computational method for solving this problem.