Suppose one has a polynomial such that two (or more) of its roots
have the same real part (this happens for the complex roots of a
real polynomial, for example).  I assert that for any polynomial,
there exists a rotation of the x,y coordinates (z=x+iy) such that
no two roots (except in the case of a multiple root) have the same
real part.
     A rotation of the coordinates is done as follows.  Suppose
the polynomial is P(z)=a(n) z^n, n=N to 0, where N is the order of
the polynomial.  Then
                     z=exp(i A) w                               (1)
where A is the rotation angle, transforms P(z) to
              a(n) exp(i nA) w^n = P(w) = b(n) w^n
P(w) has, in general, complex coefficients b(n) even when the a(n)
are real.  One chooses A so that no two roots of P(w) have the
same real part (except for multiple roots).
     It turns out that the scheme used in cpolyd.h for solving
complex polynomials produces one root when a root of P(w) is
of odd multiplicity (including 1), but does not produce a root
when a root of P(w) is of even multiplicity.  Therefore, the
program runs until it cannot find another root, and if the
proper number (N) of roots has not been found, the square root
of the polynomial is taken.  And so on in an obvious way.
     After the w roots are found, the z roots are calculated
from (1) above.
     Now, about the scheme used in cpolyd.h.
     A polynomial equation
            P(z)=a(n) z^n=a(n) (x+iy)^n =0                      (1)
where the a(n) are complex, has real and imaginary parts
            P(x,i) y^i = 0   and Q(x,i) y^i = 0.                (2)
where P(x,i) and Q(x,i) are polynomials in x.  These polynomials
are easily calculated according to the following rule:
            P(x,0) = (real part of a(n)) x^n
            Q(x,0) = (imaginary part of a(n)) x^n
            P(x,i) =  - dQ(x,i-1)/dx / i
            Q(x,i) =    dP(x,i-1)/dx / i
For a particular x, (2) is two simultaneous polynomial equations, each
of degree N, in y.  By linearly combining them, a polynomial equation
equation of degree N-1 in y can be found.  Multiplying this by y and
linearly combining the result with either of eqns (2) (say the first)
another polynomial equation of degree N-1 can be found.  Then the process
is repeated until two polynomials equations of degree 1 in y is found.
And finally a polynomial equation of degree 0 (call it lhs(x)) in y is
produced.  Unless this last equation is 0=0 (that is, lhs(x)=0), x is not
the real part of a root of (1). If lhs(x)=0, the two polynomials equations
of degree 1 are consistent and the imaginary part of z, namely y, can be
calculated from either of them.
     One varies x until lhs(x)=0.  When a zero of lhs(x) is found,
one computes y and, of course, checks that z=x+iy is a root.  One divides
out the root (P/(z-root)), resulting in a polynomial degree one less),
and repeats until N roots have been found.  If not enough roots
have been found, it is because all remaining roots are of even
multiplicity, so one takes the square root of the polynomial as
often as necessary to find N roots in all (allowing for the multiplicity
of all roots, that is, allowing for the number of times the square root
has been taken).
     In the variation of x, an upper bound can be placed on what the
magnitude of x.  One finds this bound by solving the polynomial equation
         abs(a(n))x^n=abs(a(n-1)x^(n-1)+...+abs(a(0).
One then varies x from -bound to +bound over as fine a mesh as one
cares to use.
     One learns a lot from plotting lhs(x) vs x for particular
polynomials (with no multiple roots, roots of odd multiplicity, and
roots with even multiplicity).
     The real polynomial solver polyd.h works in pretty much the same
way with some differences in detail.  Real roots are searched for
first, then complex roots of multiplicity 1, then real roots of
even multiplicity, and finally complex roots of even multiplicity.
The square root trick is not used in polyd.h.
     Polynomial solving is a difficult game.  The author would like
to hear from anyone who can devise a polynomial (of degree less than
40 or so) which this scheme will not solve.  The author has solved
some pretty high order simple polynomials (no multiple roots) and
the Wilkinson polynomial (please don't suggest (Wilkinson polynomial)^2
as an insoluble polynomial--the scheme works but the precision and time
required would be astronomical no doubt).  And, of course, there is
no use solving polynomials whose coefficients are not known to
sufficient precision.
     Comments or questions?  e-mail gammeljl@slu.edu or
gammeljl@att.net, who, if still able (maybe up to 2003--who knows?),
will reply (if he is unable to reply, you are on your own!).