Accelerated Approximation of the Complex Roots of a Univariate Polynomial (Extended Abstract)

Highly efficient and even nearly optimal algorithms have been developed for the classical problem of univariate polynomial root-finding (see, e.g., \cite{P95}, \cite{P02}, \cite{MNP13}, and the bibliography therein), but this is still an area of active research. By combining some powerful techniques developed in this area we devise new nearly optimal algorithms, whose substantial merit is their simplicity, important for the implementation.

💡 Research Summary

The paper addresses the classic problem of finding all complex roots of a univariate polynomial, focusing on achieving near‑optimal computational complexity while keeping the algorithmic structure simple enough for practical implementation. The authors build on a series of powerful techniques that have appeared in the literature—Newton’s method, isolation ratios, power‑sum computations, and fast Fourier transform (FFT) based multipoint evaluation—to construct a three‑stage procedure that can isolate, refine, and certify each root with provably low cost.

The first concept introduced is the isolation ratio of a disc D(c,r) that contains a single root of the polynomial p(x). If the disc is (1+η)‑isolated (η>0), the authors show that by increasing the isolation ratio to at least 5d² (where d is the degree) one can guarantee quadratic convergence of Newton’s iteration started at the disc’s centre. This is formalised in Theorem 1, which essentially states that a disc whose radius is reduced to Δ = O(r η/d²) becomes 5d²‑isolated, and Newton’s method converges from the centre in a single step.

To achieve such a reduction without expensive root‑finding, the paper exploits the fact that the power sum s₁ = Σ z_j of the roots inside the disc equals the root itself when there is only one root. The authors approximate s₁ by evaluating the rational function p(ω)/p′(ω) at the q‑th roots of unity ω_j = exp(2πi j/q) and forming the weighted average

 s₁* = (1/q) ∑_{j=0}^{q−1} ω_j · p(ω_j)/p′(ω_j).

Choosing q = Θ(log d) guarantees that the error |s₁*−s₁| is smaller than the target Δ. The crucial observation is that evaluating p and p′ at all ω_j can be reduced to three FFTs of size q: two FFTs compute the values of p and p′, and a third multiplies the resulting vectors elementwise. This costs O(q log q) arithmetic operations, i.e. O(log d · log log d). The authors also discuss how to construct a reduced‑degree polynomial p_q(x) (degree ≤ q−1) from the original coefficients using only d additions, which further streamlines the evaluation.

Having obtained a refined centre ĉ = s₁*, the algorithm shrinks the disc to radius Δ and thereby obtains a 5d²‑isolated sub‑disc. Newton’s iteration x_{k+1}=x_k−p(x_k)/p′(x_k) is then applied starting from ĉ. Because the disc is now sufficiently isolated, the iteration converges quadratically, and the number of Newton steps needed to achieve a prescribed absolute error ε is O(log log (1/ε)). Consequently, for a single root the total arithmetic cost is

 O(log d · log log d) + O(log log (1/ε))

operations, which is essentially optimal up to polylogarithmic factors.

The paper extends the analysis to all d roots. Assuming we are given d discs, each (1+η)‑isolated and containing a distinct simple root, the same procedure can be run in parallel on each disc. By replacing the uniform roots of unity with equally spaced points on each disc’s boundary and employing the Moenck–Borodin multipoint evaluation algorithm, the authors achieve a total cost of

 O(d · log² d · log log (1/ε))

arithmetic operations for approximating all roots within a fixed absolute error. This matches the best known bounds for the problem while avoiding the heavy machinery of earlier near‑optimal algorithms.

A substantial portion of the paper is devoted to Boolean (bit) complexity analysis. Assuming the input coefficients are known exactly, the authors introduce a working precision λ that depends on the desired output precision ℓ, the coefficient magnitude τ, and the degree d. They show that the dominant Boolean cost comes from the three FFTs and the subsequent Newton steps, yielding an overall bound of

 eOB(d²·τ + d·ℓ)

where eOB hides polylogarithmic factors. For refining a root to L bits of precision, the cost becomes eOB(d²·τ + d·L). When all roots are refined simultaneously, the authors employ fast multipoint evaluation (again via Moenck–Borodin) to keep the Boolean cost within the same asymptotic range.

The authors also discuss extensions. By combining their method with the initial isolation algorithm of Mehlhorn et al. (