Zeroes of complex polynomials: Hirano's method

Directly to the input form

• In the following, p denotes the (complex) polynomial
  p(z) = ∑_{i=0 to n} a_i z^n-i .
  
• The minimization of the absolute value of a complex polynomial is a reliable method for finding zeroes and Newton's method combined with stepcontrol, performed in the complex domain, is a means to do that. This method has been enhanced by Hirano by considering the case of multiple zeroes: one considers a choice of up to n correction directions and takes the one with the greatest promise. Given an arbitrary guess z₀ one computes the complete Taylor expansion of p at this point, using the complete Horner's scheme (or Shaw-Traub, should addition be much faster than multiplication) resulting in the n+1 values
  c_k = p^(k)(z₀)/(k!), k=0,...,n
  For c_k ≠ 0 one defines corrections for z₀ by
  0 = μ*c₀ + c_k(dz_k)^k
  where μ is a damping parameter, originally set to one but possibly diminished later. For k=1this is the damped Newton correction. The k-th root is a complex root and we always take the principal value. From these corrections we select (one of) the shortest. Our aim is to diminish |p(z₀+dz_k)| compared to |p(z₀)|: we diminish μ by multiplication with a fixed factor < 1 until
  |p(z₀+dz_k)| <= (1-δ*μ)|p(z₀)|
  δ is a fixed parameter and chosen as 0.1 here. Then z₀+dz_k becomes our next iterate and the process is repeated. One can show that there exists an universal lower bound for μ (depending only on n and the parameters chosen) and hence the method converges always (since the polynomial value tends to zero at least geometrically).
• Having found a zero one can use synthetic division in order to remove this from p and diminish the degree. It is well known that this can introduce severe roundoff effects. Hence each zero found is refined by adding some Newton steps using the original polynomial. Moreover we try to find the zeroes in increasing absolute value, a heuristic first introduced by Wilkinson. Mardens bound
  B =2 max{|a_i/a_n|^1/(n-i)|}, i=0,..,n-1
  for the polynomial with the reflected coefficients and the reciprocal zeroes is used for bounding all zeroes (a_n is the highest coefficient then) and
  1/B = z₀
  is the first iterate.
• In exact arithmetic this is an always convergent method finding all zeroes, for which even complexity bounds are known. In order to cope with roundoff we use an upper bound of
  maxcnt=100*(n+4*n³)
  steps, which is based on these theoretical bounds.
• As a final error estimate we report the magnitude of the last Newton correction, multiplied by 2.
• For details see the paper by Murota in SINUM 19, (1982), 793-799.
  

Input

• The degree n.
• You have the choice: either you specify n zeroes and let the code compute the coefficients and from these back the zeroes or you specify
• n+1 complex coefficients a₀, ... ,a_n (in this order)
• The complex numbers must be written as pairs in parentheses (re,im) denoting re + sqrt(-1)*im.
• In order to avoid unnecessary preparative work we require a_n*a₀ < > 0. a₀ is the coefficient of zⁿ.
• You may require a contour plot of |p(z)|. Then you must decide on a bounding box given by a parameter b, which defines the size of the area to be plotted: |p| will be shown over [re_min-b,re_max+b]x [im_min-b,im_max+b], where re/im_min, re/im_max are the algebraically smallest resp. largest real and imaginary part of a zero.

Output

• You get a list of the zeroes with an error estimation derived from the last Newton correction (using the original polynomial) and, if you specified the zeroes originally, this input.
• If you requested this, a contourplot of |p| around the zeroes. This is adapted internally to get a reasonable representation (for example if all zeroes are real but you did choose a large b.)

Questions ?!

• What can you learn about sensitivity of zeroes with respect to coefficient perturbations, depending on the distribution of the zeroes (linear or geometric, on the unit circle) and what can you observe when there are multiple zeroes, for example (x-1)⁴ and (x-1)⁴+(x²-2x+0.99)*eps?
• Does it indeed ''never fail'' ?
• Did you believe that a ''harmless function'' like a polynomial would look like this?

To the input form