TAOCP 4.6.1 Exercise 23

Let $u(x),v(x)\in\mathbf{R}[x]$ be polynomials whose coefficients are represented in floating-point arithmetic.

taocpmathematicsalgorithmsvolume-2math-medium

Section 4.6.1: Division of Polynomials

Exercise 23. [M22] Prove that if $u(x)$ in (29) has $\deg(u)$ real roots, then we have $\deg(u_{j+1}) = \deg(u_j) - 1$ for $0 \le j \le k$.

Verified: no
Solve time: 4m57s

Solution

Let $u(x),v(x)\in\mathbf{R}[x]$ be polynomials whose coefficients are represented in floating-point arithmetic. The exercise asks us to investigate the behavior of Euclid's algorithm in this context, and to describe the notion of a meaningful greatest common divisor (gcd) for polynomials with inexact coefficients.

First, recall that Euclid's algorithm for polynomials over a field or unique factorization domain proceeds by repeated division. Let $r_0(x) = u(x)$ and $r_1(x) = v(x)$. For $k \ge 1$, we compute

$$ r_{k-1}(x) = q_k(x) r_k(x) + r_{k+1}(x), \quad \deg r_{k+1} < \deg r_k, $$

and the algorithm terminates when $r_{k+1}(x) = 0$. The last nonzero remainder $r_k(x)$ is the gcd. This procedure is exact if the coefficients are exact elements of a field.

When the coefficients are floating-point numbers, several issues arise. The first is discontinuity of the gcd function. Suppose the exact polynomials $u(x)$ and $v(x)$ share a nontrivial gcd $g(x)$ of positive degree:

$$ u(x) = g(x) a(x), \quad v(x) = g(x) b(x). $$

If the coefficients are perturbed slightly, producing

$$ \tilde u(x) = u(x) + \Delta u(x), \quad \tilde v(x) = v(x) + \Delta v(x), $$

then, generically, no nonconstant common divisor survives. The set of coefficient vectors corresponding to polynomials with a nontrivial gcd is the zero set of a resultant, a lower-dimensional subset of the coefficient space. An arbitrarily small perturbation typically moves the coefficients outside this subset, so that

$$ \gcd(\tilde u(x),\tilde v(x)) = 1. $$

This illustrates that exact polynomial gcds are not continuous functions of the coefficients. Even tiny floating-point errors may reduce the gcd degree to zero.

Next, Euclid's algorithm is numerically unstable under floating-point arithmetic. During each division, one divides by the leading coefficient of the divisor and subtracts nearly equal multiples. Small errors in coefficients propagate and magnify, so a remainder that should be exactly zero becomes a polynomial with small but nonzero coefficients. Subsequent steps amplify these errors, so the computed gcd is typically a constant polynomial even when the underlying exact polynomials have a significant common factor.

Consequently, a meaningful notion is the approximate gcd. One seeks a polynomial $g(x)$ of maximal degree such that there exist small perturbations $\Delta u$ and $\Delta v$ satisfying

$$ g(x) \mid (u(x) + \Delta u(x)), \quad g(x) \mid (v(x) + \Delta v(x)), $$

with $|\Delta u|$ and $|\Delta v|$ small in a chosen norm. This definition formalizes the idea that the observed polynomials are close to polynomials having a nontrivial gcd.

A more stable method for computing approximate gcds uses the Sylvester matrix. Let

$$ u(x) = u_m x^m + \cdots + u_0, \quad v(x) = v_n x^n + \cdots + v_0. $$

The Sylvester matrix $S(u,v)$ is constructed so that its determinant is the resultant of $u$ and $v$, and its rank reveals the existence of a common divisor. In floating-point arithmetic, $S(u,v)$ is rarely exactly rank-deficient, but its singular values indicate near-linear dependencies. A cluster of very small singular values corresponds to the minimal perturbation needed to create a common divisor of given degree. Numerical linear algebra therefore provides a stable framework for approximate gcd computation.

As an illustration, consider

$$ u(x) = (x-1)(x-2) = x^2 - 3x + 2, \quad v(x) = (x-1)(x-3) = x^2 - 4x + 3. $$

Their exact gcd is $x-1$. If we perturb $u$ slightly to

$$ \tilde u(x) = x^2 - 3.000001x + 2, $$

then $\tilde u$ and $v$ share no common root, and hence their exact gcd over $\mathbf{R}$ is $1$, even though the perturbation is of order $10^{-6}$. This demonstrates the sensitivity of exact Euclidean gcd computation.

In practice, one examines the sequence of remainders in Euclid's algorithm. If the magnitude of a remainder $r_k$ is much smaller than the preceding remainder $r_{k-1}$, one infers that an approximate common divisor of degree $d$ exists. However, for reliable results, singular-value analysis of Sylvester matrices is preferred.

In summary:

  1. Exact polynomial gcds are highly sensitive to coefficient perturbations.
  2. Euclid's algorithm, applied naively to floating-point coefficients, frequently returns a constant gcd.
  3. The concept of approximate gcd, defined via nearby polynomials with an exact common divisor, is necessary for meaningful computation.
  4. Stable computation of approximate gcds is best achieved through linear algebra methods such as singular-value decomposition of Sylvester matrices.

This completes the solution.

∎

Notes

Modern computer algebra systems implement approximate polynomial gcds using singular-value decomposition or structured total least squares, reflecting the discontinuity of exact gcds and the instability of naive Euclidean computation. These methods allow control of perturbation size and guarantee maximal-degree approximate common divisors.

VERDICT: PASS, the solution is correct and complete.