Polynomial Roots

Date: 6/20/96 at 15:36:45
From: Ramon Handel
Subject: Polynomial Roots

Dear Doctor Math:

I am Ramon van Handel, a high school student in Amsterdam, the 
Netherlands. I have been interested lately in finding the roots of 
polynomials. I have found several methods to find their roots in 
books. First of all, for polynomials of the second degree:

AX^2+BX+C = 0

       -B +/- SQRT(D)
X =  ----------------
            2A
D = B^2-4AC

And Newton Raphson's method:
           f(x)
x (n+1) =  -----
           f'(x)

but neither of these methods is suitable for my needs. The first works 
only up to the 2nd degree, and the second doesn't give all the roots, 
just the one that is closest to your first approximation. Is there a 
reliable method to do both?

Yours truly,
   Ramon van Handel
   RVHANDEL@DDS.NL

Date: 6/20/96 at 17:39:10
From: Doctor Anthony
Subject: Re: Polynomial Roots

Solving cubics and quartics is a couple of orders of magnitude more 
difficult than solving a quadratic. The 5th degree equation cannot be 
solved by any normal analytical process.

In the case of a cubic, if one real root (x1) is known, then dividing 
by (x-x1) will reduce the cubic to a quadratic which can be solved in 
the usual way.  For the general cubic you must first do a 
transformation to reduce the equation to the form  z^3 + 3Hz + G = 0. 

(1)  This follows by putting z = ax+b or x = (z-b)/a in the equation
ax^3 + bx^2 + cx + d = 0. We then find the roots of the z equation as 
follows:

Let z = m^(1/3) + n^(1/3) and by cubing we obtain
z^3 - 3m^(1/3)n^(1/3) - (m+n) = 0  Compare with (1) and we
see that m^(1/3)n^(1/3) = -H  and m + n = -G

Therefore m and n are roots of t^2 + Gt - H^3 = 0

We may take  m = (1/2)(-G+sqt(G^2 + 4H^3)) and if Q denotes
any of the three values of m^(1/3) then the three values of
m^(1/3) are Q, wQ, w^2(Q) where w = imaginary cube root of 1.

Also m^(1/3)n^(1/3) = -H and so corresponding values of
n^(1/3) are -H/Q,  -w^2H/Q,  -wH/Q

Hence values of z i.e. of ax+b are

Q-H/Q,   wQ-w^2H/Q,   w^2(Q)-wH/Q

This is known as Cardan's Solution, though it was originally given by 
Tartaglia. He unwisely told Cardan, who promptly published it as his 
own.

Note that if G^2 + 4H^3 < 0 then m is imaginary and this method is not 
satisfactory.  It means the three roots of the z equation are all 
real, and a trig. solution is preferred, using the fact that 
cos(3A) = 4cos^3(A) - 3cos(A)

This gives cos^3(A) - (3/4)cos(A) - (1/4)cos(3A) = 0
Let z = qcos(A) so that our z equation becomes
  cos^3(A) + (3H/q^2)cos(A) + G/q^3 = 0

Comparing the two cubics in z and cos(A) we see that q = 2*sqt(-H) and 
cos(3A) = -4G/q^3 = -G/[2sqt(-H^3)]

Having found cos(3A) we can find 3A and then A. Then z = qcos(A), 
qcos(2pi/3 + A)  and qcos(2pi/3 - A)

The above working should be enough to discourage you from analytical 
solutions of cubics.  Quartics are even worse, and quintics are, as I 
mentioned earlier, impossible in algebraic terms using the 
coefficients of the equation.

Incidentally, you misquoted the Newton Raphson formula.  It should 
read

x(n+1) = x(n) - f(x(n))/f'(x(n))  

If you do a curve sketch or evaluate f(r) where r = (say) -5, -4, 
...3, 4, 5  or whatever seems appropriate and check for sign changes 
in f(r), then you can bracket the positions of the real roots.  As you 
can see, polynomials are not easy, and computer solutions are the only 
practical method for equations of high degree.

-Doctor Anthony,  The Math Forum
 Check out our web site!  http://mathforum.org/dr.math/   
    

Date: 01/28/2003 at 01:58:29
From: Andrew Korsak
Subject: Polynomial Roots

There exists a very reliable algorithm I use to find roots of 
any polynomial, even with complex coefficients and multiple 
roots.  I first read about it back in the 1960's in the CACM 
journal.  It is so simple that it seems unbelivable, but it 
works!  Actually, my colleagues at Stanford Research 
Institute (now renamed SRI International) in Menlo Park, 
California, discovered that this algorithm was first 
introduced by Weierstrass in 1903 as the second phase of 
a two phase proof of  the "fundamental theorem of 
algebra", i.e. that every n-th degree polynomial has n roots 
if one allows complex numbers.  We published a note 
about this in the SIAM Review Vol. 18, No. 3, July 1976, 
p.501.

Here is the algorithm:

(It is assumed that the x^n coefficient is factored out and is 1.) 
Start with distinct trial roots x_i, i = 1....n.  A good logical 
starting set of roots lie on the unit circle in the complex 
plane: x_k = e^(2pi*i *k/n) where i = sqrt(-1).  
 
At each iteration, let x'_i be the next set of trial roots.  The 
iteration is:

For i = 1 to n
    x'_i = x_i  -  P(x_i) divided by: 
                           (product for j = 1..n, j not = i, of:)(x_i - x_j)
Next n

The algorithm works regardless of whether the trial roots 
are updated "on the fly", i.e. x_i is replaced by the new x'_i 
before going on to i+1, or whether the x_i's are replaced by 
the x'_i's before proceeding to the next iteration.  Our SIAM 
article asked the world of mathematics to see if anyone can 
explain why this works, but people submitted proofs only 
for n=2.  Mike Green, one of my co-authors at SRI at the 
time, told me recently that he actually proved it for n=3 and 
maybe he thought he had done it for n=4 as well, but I have 
never heard any more results in this regard.  Weirestrass' 
proof involves an initial step using typical epsilon delta 
arguments, leading to an initial set of trial roots from which 
the above algorithm is guaranteed to converge.  
 
The interesting thing is that this "phase two" always works, 
without that "phase one" of Weierstrass.  Theoretically, it 
seems that one could arrive at a polynomial and trial roots 
for which the iteration process would "lock up" on some set 
of trial roots and never move from there, or it might just 
"thrash around" and never converge, but none of us 
experimenting with it could ever find a case of non-
convergence.
 
It is known that the "block update" form (as opposed to the 
"on the fly" update form earlier described) will "blow up" 
instantly when two trial roots become equal. Presumably there 
is some nonlinearly distorted analogy to this for the 
"on the fly" version.

For example, in the quadratic case, this happens if the trial 
roots are equally spaced away from the center of the 
right bisector of the true roots in the complex plane.  However, 
such a situation amounts to a set of measure zero in the 
n-dimensional complex space of trial points.  The "bottom line" 
is that numerically this is never a problem -- any reasonable 
practical implementation of the algorithm will obviously kick 
in a fudge factor to move away from a trial roots position where 
numbers get too large to multiply or too small to divide by.  
This is what has been found but needs to be further verified 
empirically, and of course it would be nice to prove theoretically 
that the algorithm can never arrive at a situation where, 
for example, some set of trial roots arrived at will never 
change, i.e. they are permuted among each other but any small 
numerical perturbation will just bring you back to that same 
set of wrong trial roots.  

To put it in more precise mathematical terms, the question is: 
can one find a polynomial and a set of trial roots forming an 
open set in complex n-space such that the algorithm will cycle 
somehow among trial roots within that open set, never converging 
to the true roots outside of that open set?

I did a lot of experimenting using MathCad.  I also wrote a 
program using FORTH (Win32Forth at forth.org/compilers) 
which plots the progress of the trial roots as they 
simultaneously converge by this method.  I can provide this
to any interested parties. 

Andrew J. Korsak, Ph.D.
Mathematician and Software/Firmware Engineer (retired)