# 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 `$\mathrm{var}(P(x)) - 2 k$`

, where `$\mathrm{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 `$\mathrm{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 `$\mathrm{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 `$\mathrm{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. `$\mathrm{var}(A(x)) > 1$`

or `$\mathrm{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

```
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 + (\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.