FLINT3.7.0-dev
Algorithm

Certified roots of polynomials

Cost
O(n2)O(n^2) operations at pp bits per iteration for degree nn, with O(n)O(n) iterations typically; validation O(n2)O(n^2) operations; isolation needs a precision of the order of nn times the logarithm of the inverse root separation in the worst case.

Given a polynomial f(x)=∑k=0nakxkf(x) = \sum_{k=0}^{n} a_k x^k of degree nn with complex ball coefficients, acb_poly_find_roots computes balls B1,…,BnB_1, \ldots, B_n and a number ss such that each of B1,…,BsB_1, \ldots, B_s contains exactly one root of every polynomial in the coefficient balls, and each of the remaining balls contains at least one root. If s=ns = n, all roots are isolated. The computation separates the numerical part, which approximates the roots with floating-point arithmetic and no error control, from the certification, which is a short computation in ball arithmetic based on the following classical bound.

Inclusion disks#

Theorem 1 (Inclusion disk). Let ff be a polynomial of degree n≥1n \ge 1, m∈Cm \in \mathbb{C} with f′(m)≠0f'(m) \ne 0, and ρ=n∣f(m)/f′(m)∣\rho = n |f(m)/f'(m)|. Then the closed disk of radius ρ\rho around mm contains a root of ff.

Proof. Let ζ1,…,ζn\zeta_1, \ldots, \zeta_n be the roots of ff. If f(m)=0f(m) = 0 there is nothing to prove. Otherwise f′(m)/f(m)=∑i1/(m−ζi)f'(m)/f(m) = \sum_i 1/(m - \zeta_i). If ∣m−ζi∣>ρ|m - \zeta_i| > \rho for all ii, then

∣f′(m)f(m)∣≤∑i=1n1∣m−ζi∣<nρ=∣f′(m)f(m)∣,\left|\frac{f'(m)}{f(m)}\right| \le \sum_{i=1}^{n} \frac{1}{|m - \zeta_i|} < \frac{n}{\rho} = \left|\frac{f'(m)}{f(m)}\right|,

a contradiction [Kob2010].

Corollary 2 (Isolation). If disks D1,…,DnD_1, \ldots, D_n around points m1,…,mnm_1, \ldots, m_n with the radii of Theorem 1 are pairwise disjoint, then each contains exactly one root of ff, and these are all the roots.

Proof. Each disk contains at least one root, the disks are disjoint, and ff has exactly nn roots counted with multiplicity.

In ball arithmetic, _acb_poly_root_inclusion evaluates f(m)f(m) and f′(m)f'(m) at the exact point mm (the midpoint of an approximation) with the coefficient balls, takes an upper bound u≥∣f(m)∣u \ge |f(m)| and an upper bound v≥1/∣f′(m)∣v \ge 1/|f'(m)| over all polynomials in the balls, and returns the box around the disk of radius nuvn u v. Since uvu v bounds ∣f(m)/f′(m)∣|f(m)/f'(m)| for every polynomial in the coefficient balls, the box contains a root of each of them, and Corollary 2 applies to each of them. If f′(m)f'(m) contains zero, vv is infinite and the box is the whole plane. The evaluation uses Horner's scheme with the error of each partial result bounded by a disk (the larger of the real and imaginary radii would overestimate it by a factor 2\sqrt 2 per step, which matters for high degree).

_acb_poly_validate_roots computes these boxes for all approximations, determines which boxes overlap no other box, moves them to the front and returns their number ss; by Theorem 1 and the argument of Corollary 2, each of these ss boxes contains exactly one root if s=ns = n, and in general each of them contains at least one root and the boxes that overlap no other one are isolated from each other, but two roots may share an overlapping group.

The radius n∣f(m)/f′(m)∣n|f(m)/f'(m)| is about nn times the distance of mm to the nearest root when that root is simple and well separated, so the enclosure is only log⁡2n\log_2 n bits less accurate than the approximation. At a root of multiplicity kk, ∣f/f′∣≈∣m−ζ∣/k|f/f'| \approx |m - \zeta|/k, but mm is only accurate to about p/kp/k bits, so the disks of the kk approximations overlap at every precision: multiple roots can never be isolated, and clusters of roots need a precision that resolves them.

The Weierstrass–Durand–Kerner iteration#

Given approximations z1,…,znz_1, \ldots, z_n of the roots of ff, the iteration

zk←zk−f(zk)an∏j≠k(zk−zj)z_k \leftarrow z_k - \frac{f(z_k)}{a_n \prod_{j \ne k} (z_k - z_j)}

of Weierstrass [Wei1891], Durand [Dur1960] and Kerner [Ker1966] is Newton's method for the system of equations that expresses the coefficients of an∏k(x−zk)a_n \prod_k (x - z_k) in terms of the zkz_k (Vieta's formulas). It converges quadratically to simple roots from good starting values, and in practice from almost all starting values on a circle around the roots. The Aberth iteration [Abe1973], which converges cubically, is implemented in _gr_poly_refine_roots_aberth but not used by acb_poly_find_roots. Each iteration costs O(n2)O(n^2) operations.

The starting values are chosen by the Newton polygon technique of Bini [Bin1996]. Let k1<k2<⋯<kqk_1 < k_2 < \cdots < k_q be the indices of the vertices of the upper convex hull of the points (k,log⁡2∣ak∣)(k, \log_2 |a_k|). For each edge, mi=ki−ki−1m_i = k_i - k_{i-1} points are placed on the circle of radius ∣aki−1/aki∣1/mi|a_{k_{i-1}}/a_{k_i}|^{1/m_i}, which is where the polynomial has mim_i roots of approximately that modulus when the coefficients dominate as the polygon indicates. The angles are spread with an irrational offset so that the configuration has no symmetry that the iteration would preserve.

Implementation#

_acb_poly_find_roots proceeds as follows for degree nn:

  1. If the leading coefficient contains zero, the roots are unbounded: all outputs are set to [±∞]+[±∞]i[\pm\infty] + [\pm\infty]i and 0 is returned. Degree 1 is solved directly.
  2. If p≤53p \le 53, the polynomial is first rescaled so that its coefficients fit in double precision: with the line y=αk+βy = \alpha k + \beta that best fits the points (k,log⁡2∣ak∣)(k, \log_2 |a_k|) in the minimax sense (computed from the upper convex hull and the endpoints, a rotating calipers argument), it substitutes x→2−αxx \to 2^{-\alpha} x and divides by 2β2^{\beta}. cd_poly_find_roots then runs the iteration in double precision, with steps limited in size (between 2−102^{-10} and 2−32^{-3}, adapted to progress) and roots of modulus above 1 updated through the reversed polynomial to avoid overflow. If the last relative correction is below 1/81/8, the result goes to validation.
  3. Otherwise the iteration runs at pp bits with floating-point complex numbers (nfloat_complex if available, acf otherwise) through the generic _gr_poly_refine_roots_wdk, for at most 2n+p2n + \sqrt{p} iterations by default. When the largest correction, relative to the largest root, falls below 2−p/22^{-p/2} (respectively 2−p/32^{-p/3}, 2−p/42^{-p/4}), the iteration is stopped two (three, four) iterations later, since quadratic convergence then reaches full precision.
  4. The approximations are validated by _acb_poly_validate_roots.

The working precision must exceed the number of bits lost in evaluating ff near the roots, which for a root ζ\zeta is about log⁡2\log_2 of the condition number ∑k∣ak∣∣ζ∣k/∣f′(ζ)∣\sum_k |a_k| |\zeta|^k / |f'(\zeta)|. Users therefore call the function in a loop with doubling precision, reusing the previous approximations as starting values. For polynomials with exact coefficients, the radius of the inclusion disks shrinks proportionally to 2−p2^{-p} once the precision suffices, and every squarefree polynomial is eventually isolated.

Real roots#

For a real polynomial, the roots in the upper and the lower half-plane come in conjugate pairs, and a root whose isolating box meets the real axis might still be nonreal. acb_poly_validate_real_roots proves that it is real: let rr be the number of boxes that meet the real axis. The test requires r≡n(mod2)r \equiv n \pmod 2 and checks, at the midpoints between the consecutive candidate intervals (sorted by midpoint), that ff has signs that alternate starting from the sign at −∞-\infty and ending with the sign at +∞+\infty. Then ff changes sign within each candidate interval, which by the intermediate value theorem contains a real root; since each box contains exactly one root, that root is real, and the other boxes contain the nonreal roots.

Root bounds#

arb_poly_root_bound_fujiwara and acb_poly_root_bound_fujiwara bound the absolute values of all roots by Fujiwara's bound [Fuj1916]

∣ζ∣≤2max⁡{∣an−1an∣,∣an−2an∣1/2,…,∣a1an∣1/(n−1),∣a02an∣1/n},|\zeta| \le 2 \max\left\{\left|\frac{a_{n-1}}{a_n}\right|, \left|\frac{a_{n-2}}{a_n}\right|^{1/2}, \ldots, \left|\frac{a_1}{a_n}\right|^{1/(n-1)}, \left|\frac{a_0}{2 a_n}\right|^{1/n}\right\},

computed with upper bounds for the numerators and a lower bound for ∣an∣|a_n|, so that it holds for every polynomial in the balls. The bound is within a factor 2n2n of the largest root modulus.

Integer polynomials#

For a squarefree integer polynomial, arb_fmpz_poly_complex_roots always succeeds, because it controls the precision itself:

  1. It writes f(x)=xeg(xd)f(x) = x^e g(x^d) with the largest dd, so that the roots of ff are 0 (if e=1e = 1) and the dd-th roots of the roots of gg.
  2. It isolates the real roots of gg exactly by Descartes' rule of signs (Counting and isolating real roots over Z) and refines them by interval Newton iteration (arb_fmpz_poly_refine_root_arb, see Isolating and refining real roots).
  3. If the real roots are numerous or clustered, it divides them out in ball arithmetic, increasing the precision until the quotient (whose roots are the nonreal roots) is accurate to the working precision, so that the Weierstrass–Durand–Kerner iteration does not have to resolve the clusters.
  4. It calls acb_poly_find_roots at 53, 128, 256, ... bits until all roots are isolated and accurate to the target, proves which roots are real with acb_poly_validate_real_roots when necessary, and takes dd-th roots.

The output lists the real roots in increasing order, then the nonreal roots in conjugate pairs; the second member of each pair is the conjugate of the first, so that the pairs are exactly symmetric.

Use in FLINT#

Referenced from#

References#

  • [Kob2010] A. Kobel, "Certified Complex Numerical Root Finding", Seminar on Computational Geometry and Geometric Computing (2010), http://www.mpi-inf.mpg.de/departments/d1/teaching/ss10/Seminar_CGGC/Slides/02_Kobel_NRS.pdf
  • [Wei1891] K. Weierstrass. Neuer Beweis des Satzes, dass jede ganze rationale Function einer Veränderlichen dargestellt werden kann als ein Product aus linearen Functionen derselben Veränderlichen. Sitzungsberichte der Königlich Preussischen Akademie der Wissenschaften zu Berlin, 1891, 1085–1101.
  • [Dur1960] E. Durand. Solutions numériques des équations algébriques, tome I. Masson, Paris, 1960.
  • [Ker1966] I. O. Kerner. Ein Gesamtschrittverfahren zur Berechnung der Nullstellen von Polynomen. Numerische Mathematik 8 (1966), 290–294. https://doi.org/10.1007/BF02162564
  • [Abe1973] O. Aberth. Iteration methods for finding all zeros of a polynomial simultaneously. Mathematics of Computation 27 (1973), 339–344. https://doi.org/10.1090/S0025-5718-1973-0329236-7
  • [Bin1996] D. A. Bini, "Numerical computation of polynomial zeros by means of Aberth's method". Numerical algorithms 13 (1996): 179-200. https://doi.org/10.1007/BF02207694
  • [Fuj1916] M. Fujiwara. Über die obere Schranke des absoluten Betrages der Wurzeln einer algebraischen Gleichung. Tôhoku Mathematical Journal 10 (1916), 167–171.
Cite this page

The FLINT team, “Certified roots of polynomials,” FLINT: Fast Library for Number Theory, documentation, version 3.7.0-dev, 2026. https://flintlib.org/doc/algorithms/ball-arithmetic/polynomial-roots.html (generated 2026-11-02).