A Continued Fraction Algorithm for Approximating All Real Polynomial Roots

Abstract

It is almost presumptuous to write another paper on approximating the real roots of a polynomial with real coefficients because so much is known, theoretically and practically, about the subject. However, the advent of computers makes it worthwhile to look again at some long forgotten theorems that can now be efficiently implemented on a computer to find all the real roots (and their multiplicities) of a real polynomial. Since a real number can be accommodated explicitly on the computer only by rational approximations, the generality of our method is restricted by the computer to polynomials with rational coefficients. Indeed, all the examples we tested had integral coefficients because a polynomial with rational coefficients is easily made equivalent to a polynomial with integral coefficients. For all practical purposes, therefore, our polynomials have integral coefficients. The theorem that lies at the heart of our algorithm was proved in 1836 by Vincent ([6], [7]) and uses a continued fraction algorithm to separate the roots. Once a root is separated out, we continue with a continued fraction algorithm that approximates the root to any desired accuracy. An interesting feature of the program is that there is no problem with round-off error because not only are the coefficients integers, but integer arithmetic is used and all computations are exact. We programmed in APL because it handles polynomials so efficiently. Our algorithm for finding all real roots first tests for rational roots, using the rational root theorem, and then factors out the linear factors, leaving a polynomial whose real roots are all irrational. The multiplicity of each rational root is also noted. The algorithm then tests the new polynomial for multiple irrational roots using an elegant procedure described in Uspensky [6, p. 65]. Let X,.denote the product of all linear factors corresponding to simple roots; let X2 denote the product of all linear factors corresponding to double roots, etc. If there are no roots of multiplicity k, we set Xk = 1. If, for example, P(x) = (x 1)(x 2)(x 3)2(x 4)2(x 5)3, then P(x) can be written in the form P(x) = X2X2X3, where X, = (x 1)(x 2), X2 = (x 3)(x 4), and XA = (x 5). In general a polynomial P(x) with r multiple roots can be written in the form P(x) = a,,X,X2 ... X'where a,, is a constant and X, = (x b,) ... (x bj) is a polynomial whose simple roots are all the roots of multiplicity i of P. The X, are obtained as follows: Let D, = gcd{P, P'} (P, P'), where P' is the derivative of P. Generate a sequence of polynomials D2= (D, D ), D3 =(D2, D2), etc. An easy calculation shows that