Solving Polynomials

2010-10-28 Tavian Barnes

A well known (if not by name) theorem is the Abel–Ruffini theorem, which states that there is no algebraic expression for the roots of polynomials with degree higher than 4.

A not-so-well-known fact is that for any polynomial P(x)P(x), it is possible to find (with exact arithmetic) a set of ranges each containing exactly one root of P(x)P(x). One such algorithm is due to James Victor Uspensky in 1948.

Uspensky's algorithm requires all roots to have multiplicity 1; P(x)P(x) can be trivially replaced with P(x)/gcd(P(x),P(x))P(x)/\gcd(P(x), P'(x)) to eliminate duplicate roots, if necessary. The algorithm relies on Descartes' rule of signs, which states that the number of positive solutions (counting multiplicities) to a polynomial is equal to var(P(x))2k\mathrm{var}(P(x)) - 2 k, where var(P(x))\mathrm{var}(P(x)) is the number of sign changes between consecutive non-zero coefficients of P(x)P(x), and kk is a non-negative integer.

This clearly implies that if var(P(x))\mathrm{var}(P(x)) is 0 or 1, then the polynomial must have exactly 0 or 1 root in (0,)(0, \infty), respectively. Otherwise, we can transform the polynomial and try again. Expanding A(x)=P(x+1)A(x) = P(x + 1), we can test var(A(x))\mathrm{var}(A(x)) to learn about the roots in (1,)(1, \infty). Expanding B(x)=(x+1)nP(1/(x+1))B(x) = (x + 1)^n P(1/(x + 1)), where nn is the degree of the polynomial P(x)P(x), we can similarly test var(B(x))\mathrm{var}(B(x)) to learn about the roots in (0,1)(0, 1). If A(0)=0A(0) = 0, then x=1x = 1 is a root of P(x)P(x).

If we don't get conclusive information from one of the tests (i.e. var(A(x))>1\mathrm{var}(A(x)) > 1 or var(B(x))>1\mathrm{var}(B(x)) > 1), apply the same two transformations (P(x)P(x+1)P(x) \mapsto P(x+1) and P(x)(x+1)nP(1/(x+1))P(x) \mapsto (x+1)^n P(1/(x+1)) to the offending polynomial(s). Uspensky proves that this recursive procedure always terminates for square-free polynomials. A reasonable bound on the recursion depth may also be shown.

From a CS perspective, the two transformations are easy to implement. If P(x)=i=0npixiP(x) = \sum_{i=0}^n p_i x_i, then

A(x)=i=0n(j=in(ji)pj)xi,B(x)=i=0n(j=in(ji)pnj)xi.A(x) = \sum_{i=0}^n \left( \sum_{j=i}^n \binom{j}{i} p_j \right) x^i, \quad B(x) = \sum_{i=0}^n \left( \sum_{j=i}{n} \binom{j}{i} p_{n-j} \right) x^i.

Another CS-related point is that one of the resultant intervals may have an infinite upper bound. Cauchy provides a general upper bound 1+(maxi=0n1pi)/pn1 + (\max_{i=0}^{n-1} |p_i|) / |p_n| which may be used to produce a finite upper bound.

A C implementation of this algorithm as a root-solver in my ray-tracer Dimension can be seen here.