Algorithm for finding polynomial roots
In mathematics, Graeffe's method or Dandelin–Lobachesky–Graeffe method is an algorithm for finding all of the roots of a polynomial. It was developed independently by Germinal Pierre Dandelin in 1826 and Lobachevsky in 1834. In 1837 Karl Heinrich Gräffe also discovered the principal idea of the method. The method separates the roots of a polynomial by squaring them repeatedly. This squaring of the roots is done implicitly, that is, only working on the coefficients of the polynomial. Finally, Viète's formulas are used in order to approximate the roots.
Dandelin–Graeffe iteration.
Let "p"("x") be a polynomial of degree n
formula_1
Then 
formula_2
Let "q"("x") be the polynomial which has the squares formula_3 as its roots,
formula_4
Then we can write:
formula_5
"q"("x") can now be computed by algebraic operations on the coefficients of the polynomial "p"("x") alone. Let: 
formula_6
then the coefficients are related by
formula_7
Graeffe observed that if one separates "p"("x") into its odd and even parts:
formula_8
then one obtains a simplified algebraic expression for "q"("x"):
formula_9
This expression involves the squaring of two polynomials of only half the degree, and is therefore used in most implementations of the method. 
Iterating this procedure several times separates the roots with respect to their magnitudes. Repeating "k" times gives a polynomial of degree n:
formula_10
with roots 
formula_11 
If the magnitudes of the roots of the original polynomial were separated by some factor formula_12, that is, formula_13, then the roots of the "k"-th iterate are separated by a fast growing factor 
formula_14.
Classical Graeffe's method.
Next the Vieta relations are used
formula_15
If the roots formula_16 are sufficiently separated, say by a factor formula_12, formula_18, then the iterated powers formula_19 of the roots are separated by the factor formula_20, which quickly becomes very big.
The coefficients of the iterated polynomial can then be approximated by their leading term,
formula_21
formula_22 and so on,
implying
formula_23
Finally, logarithms are used in order to find the absolute values of the roots of the original polynomial. These magnitudes alone are already useful to generate meaningful starting points for other root-finding methods.
To also obtain the angle of these roots, a multitude of methods has been proposed, the most simple one being to successively compute the square root of a (possibly complex) root of formula_24, "m" ranging from "k" to 1, and testing which of the two sign variants is a root of formula_25. Before continuing to the roots of formula_26, it might be necessary to numerically improve the accuracy of the root approximations for formula_25, for instance by Newton's method.
Graeffe's method works best for polynomials with simple real roots, though it can be adapted for polynomials with complex roots and coefficients, and roots with higher multiplicity. For instance, it has been observed that for a root formula_28 with multiplicity "d", 
the fractions 
formula_29 tend to formula_30
for formula_31. This allows to estimate the multiplicity structure of the set of roots.
From a numerical point of view, this method is problematic since the coefficients of the iterated polynomials span very quickly many orders of magnitude, which implies serious numerical errors. One second, but minor concern is that many different polynomials lead to the same Graeffe iterates.
Tangential Graeffe method.
This method replaces the numbers by truncated power series of degree 1, also known as dual numbers. Symbolically, this is achieved by introducing an "algebraic infinitesimal" formula_32 with the defining property formula_33. Then the polynomial
formula_34 has roots formula_35, with powers
formula_36
Thus the value of formula_37 is easily obtained as fraction formula_38
This kind of computation with infinitesimals is easy to implement analogous to the computation with complex numbers. If one assumes complex coordinates or an initial shift by some randomly chosen complex number, then all roots of the polynomial will be distinct and consequently recoverable with the iteration.
Renormalization.
Every polynomial can be scaled in domain and range such that in the resulting polynomial the first and the last coefficient have size one. If the size of the inner coefficients is bounded by "M", then the size of the inner coefficients after one stage of the Graeffe iteration is bounded by formula_39. After "k" stages one gets the bound formula_40 for the inner coefficients.
To overcome the limit posed by the growth of the powers, Malajovich–Zubelli propose to represent coefficients and intermediate results in the "k"th stage of the algorithm by a scaled polar form
formula_41
where formula_42 is a complex number of unit length and formula_43 is a positive real. Splitting off the power formula_44 in the exponent reduces the absolute value of "c" to the corresponding dyadic root. Since this preserves the magnitude of the (representation of the) initial coefficients, this process was named renormalization.
Multiplication of two numbers of this type is straightforward, whereas addition is performed following the factorization formula_45, where formula_46 is chosen as the larger of both numbers, that is, formula_47. Thus
formula_48 and formula_49 with formula_50
The coefficients formula_51 of the final stage "k" of the Graeffe iteration, for some reasonably large value of "k", are represented by pairs formula_52, formula_53. By identifying the corners of the convex envelope of the point set formula_54 one can determine the multiplicities of the roots of the polynomial. Combining this renormalization with the tangent iteration one can extract directly from the coefficients at the corners of the envelope the roots of the original polynomial.