Iowa Hills Software
P51 Polynomial Root Finder and
Quartic and Cubic Root Finder
This page gives the source code for two root finders. The first is the P51 Root Finder which is for polynomials with real coefficients, up to 100th order. The second is a Quartic and Cubic root finder, which solves for roots algebraically. The C++ source code for both is here. P51 Source Code.
The P51 Root Finder is a top to bottom rewrite of the Jenkins Traub Root Finder. It uses some of the same mathematics used by Jenkins and Traub, but it is not a translation of their algorithm into C. The P51 algorithm is quite different in many respects.
The Jenkins Traub root finder was written in the early 70's, making it more than 40 years old. In 1970 the Jenkins Traub root finder was state of the art technology, and to its credit, many authors claim that it is still the root finder by which all others are judged. The problem however, is that this archaic piece of FORTRAN code is virtually indecipherable. Its liberal use of goto statements and global variables make it incredibly difficult to follow.
The development of P51 started as an effort to do a straight conversion of Jenkins Traub from FORTRAN into C. The original intent was to get rid of the global variables and goto statements, but otherwise make as few changes as possible. For numerous reasons however, this simple goal was unattainable, and we ultimately made numerous substantive changes.
For example, this algorithm no longer scales the polynomial coefficients at the beginning of the algorithm. Jenkins and Traub did this as a way to avoid overflow and underflow, but this scaling led to a convergence test that is indecipherable, which may be the ultimate reason for the algorithm's inaccuracies.
The P51 algorithm does not use special forms of recursion, or combine fixed shift iterations with variable shift iterations. Another important difference is the degree of, and the calculation of, the K polynomial used to find 2nd order factors (sometimes referred to as an H polynomial). We use an N-2 degree polynomial where Jenkins Traub used an N-1 degree polynomial. We initialize K with the 2nd derivative of P, while J-T initialized K to the 1st derivative.
Ultimately, these changes led to a more linear algorithm, making it much easier to follow and understand. The underlying mathematics used in P51 originated with Jenkins and Traub in this paper, (eq. 9.8 is of particular interest), but we did not implement their math as they did in their FORTRAN algorithm, and made significant changes to the equations used to update the 2nd order coefficients u and v.
In a word, accuracy is the biggest difference between this algorithm and Jenkins Traub. In virtually all of our tests, this algorithm gave more accurate results than J-T (unless the error was insignificant). The tests were run by building polynomials from known roots, and then comparing the calculated roots with the originals.
The error was calculated as follows:
The errors from P51 are typically 1% of those from Jenkins Traub, but as with all root finder tests, this result depends on the type of polynomial tested (degree, root types, and root locations). In some cases, the P51's results were only 10X more accurate, while in other cases the difference would be more than 1000X.
The accuracy tests were done by taking the average error result from 250 tests, across a wide range of polynomial orders and root locations. The J-T algorithm gave slightly better accuracy in some cases, but this was not typical.
Unfortunately, an accuracy test such as this is fraught with problems. For example, we can't know with certainty that the J-T code we were testing was done properly. We downloaded two different algorithms, which were converted to C in dramatically different ways, but both gave similar results. Unfortunately, both were done using 64 bit double variables, instead of 80 bit long double's as is used in P51. We converted the J-T code to long doubles, but this had little effect. We may have also done this conversion improperly as the convergence criteria used in J-T is quite complicated.
We would urge the reader to test this algorithm, or any other root finder, before using it in your application. No root finder is absolutely robust, so you need to test them with the types of polynomials you typically work with.
The author of this algorithm is an electrical engineer, so the testing was biased toward polynomials with roots primarily located between the two circles with radii of 0.01 and 100. For a given polynomial, the variation in polynomial coefficients was limited to about 50 orders of magnitude.
To facilitate your testing, we included two test functions in the source code files. The first function will build a polynomial from a set of roots. This makes it easy to generate a test polynomial from known roots, then the found roots can be compared to the originals.
The second function sorts roots. The problem with comparing roots is that the roots returned by the root finder will not be arranged in the same order as the originals. This simple function will align the original roots with the found roots in their respective arrays, making it easy to compare the two.
Speed is only important as an indicator of a properly functioning algorithm. Otherwise, speed isn't an issue, unless it's absurdly slow. We needed better accuracy from a root finder, so that's what we emphasized. Nevertheless, P51 can be slightly faster, or slightly slower than J-T, depending on root locations.
Links to Jenkins Traub in C.
These two sources are close representations of the original J-T algorithm written in C.
Quartic and Cubic Root Finder
The source code file for P51 also contains a 3rd and 4th order algebraic root finder. This code originated at TOMS, Algorithm 326 which is based on the paper "Roots of Low Order Polynomials" by Terence R.F.Nonweiler CACM (April 1968).
We made extensive changes to the TOMS code however. Among other things, we use scaling and polynomial reversals as a way to greatly improve the algorithm's accuracy and range. We also made extensive changes to the test conditions used in the algorithm. This code should be useful if a relatively light and fast algorithm is needed for low order polynomials. Since this code doesn't use iteration, it will find 4th order roots about seven times faster than P51.
P51 was named after the American P-51 Mustang used to hunt Zeros in WWII.
A Utility Program with P51 Implemented
You can try P51 with this Filter Polynomial utility program. This program's primary purpose is to generate the 2nd order coefficients needed for analog filters, but the P51 root finder can be accessed directly as shown in these two screen shots. The polynomial coefficients are entered in the left pane.
Copyright 2015 Iowa Hills Software