Solving Polynomials

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)\), it is possible to find (with exact arithmetic) a set of ranges each containing exactly one root of \(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)\) can be trivially replaced with \(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 \(\DeclareMathOperator{\var}{var} \var(P(x))-2\,k\), where \(\var(P(x))\) is the number of sign changes between consecutive non-zero coefficients of \(P(x)\), and \(k\) is a non-negative integer.

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

If we don't get conclusive information from one of the tests (i.e. \(\var(A(x))>1\) or \(\var(B(x))>1\)), apply the same two transformations \((P(x) \mapsto P(x+1)\) and \(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)=\sum_{i=0}^n p_i\,x^i\), then

\(\displaystyle A(x)=\sum_{i=0}^n \left(\sum_{j=i}^n {j \choose i}\,p_j\right)\,x^i,\quad B(x)=\sum_{i=0}^n \left(\sum_{j=i}^n {j \choose 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+(\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.

Leave a Reply

Your email address will not be published. Required fields are marked *