Polynomial Root Finding

Overview

Finding all roots (including complex) of $p(x) = c_0 + c_1 x + \cdots + c_n x^n$ has its own family of dedicated algorithms beyond general 1D root finding. Closed-form formulas, iterative methods, and reduction to eigenvalue problems are used depending on the degree.

Related API: solveLinear, solveQuadratic, solveCubic, solveQuartic, jenkinsTraub, laguerre, durandKernerAberth, bairstow.

Closed-Form Solutions for Degree 1–4

Explicit formulas exist for $n \leq 4$. By Galois theory (Abel-Ruffini), no general closed-form solution exists for $n \geq 5$.

The quadratic case is the classic cancellation example; sangi uses an inverse-Vieta form when $a \approx 0$ or $b^2 \gg 4ac$.

Companion Matrix + QR (Francis Algorithm)

For the monic polynomial $p(x) = x^n + c_{n-1} x^{n-1} + \cdots + c_0$, the companion matrix has $p$ as its characteristic polynomial:

$$C(p) = \begin{pmatrix} 0 & 0 & \cdots & 0 & -c_0 \\ 1 & 0 & \cdots & 0 & -c_1 \\ 0 & 1 & \cdots & 0 & -c_2 \\ \vdots & & \ddots & & \vdots \\ 0 & 0 & \cdots & 1 & -c_{n-1} \end{pmatrix}.$$

$\det(xI - C(p)) = p(x)$, so the eigenvalues of $C(p)$ equal the zeros of $p$. This casts polynomial root finding as an eigenvalue problem.

Francis algorithm (QR iteration)

QR iteration with shift strategies (Wilkinson shift, Francis double shift) achieves quadratic and cubic convergence and computes all eigenvalues in $O(n^3)$. NumPy's np.roots and MATLAB's roots() both use this approach.

Pros: numerically mature, complex roots handled uniformly.

Cons: $O(n^3)$, coefficient scaling needed, accuracy drops near clustered roots.

sangi does not expose companion+QR as a top-level function, but it can be invoked via a linear-algebra backend such as Eigen. The default paths are Jenkins-Traub and Aberth-Ehrlich.

Durand-Kerner (Weierstrass)

Iterate all $n$ root approximations simultaneously. Starting from $z_1^{(0)}, \ldots, z_n^{(0)}$:

$$z_i^{(k+1)} = z_i^{(k)} - \frac{p(z_i^{(k)})}{\prod_{j \neq i} (z_i^{(k)} - z_j^{(k)})}.$$

The denominator is the product of differences from the other roots; at the true roots $z_i = \alpha_i$ the right-hand side vanishes. Discovered by Weierstrass (1891), rediscovered independently by Durand (1960) and Kerner (1966).

Convergence order: quadratic at simple roots.

Standard initial values are equally spaced points on a circle, $z_i^{(0)} = \beta \cdot e^{2\pi i (i-1)/n + \theta}$, where $\beta$ is determined by coefficient magnitudes.

Aberth-Ehrlich

Aberth (1973) and Ehrlich (1967) refined Durand-Kerner by blending in Newton's method:

$$z_i^{(k+1)} = z_i^{(k)} - \frac{N_i}{1 - N_i \cdot \sum_{j \neq i} (z_i^{(k)} - z_j^{(k)})^{-1}}, \quad N_i = \frac{p(z_i^{(k)})}{p'(z_i^{(k)})}.$$

$N_i$ is the standard Newton step; the correction $\sum_{j \neq i} (z_i - z_j)^{-1}$ represents "repulsion" from the other roots.

Convergence order: cubic at simple roots.

Benefits of simultaneous iteration

Deflation-based methods like Jenkins-Traub compute $p(x)/(x - \alpha)$ after each root extraction. In floating-point this perturbs the polynomial slightly, degrading accuracy for later roots. Aberth always evaluates the original $p(x)$, so deflation error never accumulates.

It parallelizes naturally, and modern implementations such as MPSolve use Aberth for large-degree polynomials (degree $\geq 10^4$). sangi's durandKernerAberth implements this method.

Jenkins-Traub

Jenkins and Traub (1970) is a three-stage iterative algorithm, available as RPOLY (real) and CPOLY (complex). Known for excellent numerical stability, distributed as ACM Algorithm 419/493.

Three-stage structure

1. Stage 1 (no-shift): iterate the $H$-polynomial without a shift to collect global information ($M \approx 5$ steps).
2. Stage 2 (fixed-shift): fix a candidate root $s$ and iterate to determine the direction of $H$ ($L \approx 9 + n/2$ steps).
3. Stage 3 (variable-shift): update $s$ each step to lock onto quadratic convergence.

$H$-polynomials update as $H^{(j+1)}(x) = (H^{(j)}(x) - H^{(j)}(s) p(x)/p(s)) / (x - s)$, gradually removing unwanted roots and focusing on the target. After Stage 3 converges to a root, it is deflated out and the algorithm restarts on the remaining polynomial.

Convergence order: about 1.62 (asymptotic across all three stages).

Although the implementation is complex, the method is numerically mature, and sangi's jenkinsTraub is the default path of solvePolynomial for degree ≥ 5.

Laguerre's Method

A polynomial-specific cubically convergent method in the Newton-Halley family:

$$x_{n+1} = x_n - \frac{n \cdot p(x_n)}{G(x_n) \pm \sqrt{(n-1)(n \cdot H(x_n) - G(x_n)^2)}},$$

where $G(x) = p'(x)/p(x)$ and $H(x) = G(x)^2 - p''(x)/p(x)$. The sign of the square root is chosen to maximize the magnitude of the denominator.

Global convergence

Laguerre's method achieves cubic convergence at simple roots and is known for global convergence from almost any initial value (Adams-Smith). Unlike Newton and Halley, which converge only locally, this makes the choice of starting point much less critical when seeking complex roots of a real-coefficient polynomial.

sangi's laguerre is one of the options in solvePolynomial.

Bairstow's Method

Finds complex roots of a real-coefficient polynomial using real arithmetic only. Divide $p$ by a quadratic factor $(x^2 - rx - s)$ and use Newton's iteration on $(r, s)$ to drive the remainder $b_1 x + b_0$ to zero:

$$\begin{pmatrix} \partial b_0/\partial r & \partial b_0/\partial s \\ \partial b_1/\partial r & \partial b_1/\partial s \end{pmatrix} \begin{pmatrix} \Delta r \\ \Delta s \end{pmatrix} = - \begin{pmatrix} b_0 \\ b_1 \end{pmatrix}.$$

The quadratic factor directly yields a complex-conjugate root pair, with no complex arithmetic required. Useful when complex arithmetic is unavailable or undesirable (e.g. embedded targets). On modern general-purpose hardware Jenkins-Traub and Aberth-Ehrlich are usually preferable, but sangi still ships bairstow.

Sturm Sequences and Real-Root Counting

Construct the Sturm sequence $f_0 = f, f_1 = f', f_{i+1} = -(f_{i-1} \bmod f_i)$. The number of distinct real roots of $f$ in $[a, b]$ is exactly

$$\#\{\alpha \in [a, b] : f(\alpha) = 0\} = v(a) - v(b),$$

where $v(x)$ counts sign changes in the Sturm sequence (Sturm's theorem). Multiple roots are counted only once.

Use cases

• Exact count of real roots before searching.
• Isolation of real roots by repeated bisection (one root per subinterval).
• Post-isolation refinement via 1D bracketing (e.g. Brent).

sangi's signVariations(p) returns the Descartes-rule upper bound, and vincentRealRootIsolation uses the Vincent-Akritas-Strzeboński continued-fraction approach (faster than Sturm).

sangi Dispatch Rules

solvePolynomial(p) selects an algorithm by degree:

• Degree 1–4: closed-form formulas (solveLinear/Quadratic/Cubic/Quartic).
• Degree ≥ 5: jenkinsTraub.

Reasons for defaulting to Jenkins-Traub:

• Mature implementation with proven resilience to coefficient scaling and ill-conditioned inputs.
• 50+ years of production use as RPOLY/CPOLY.
• Stable accuracy/speed balance from degree 5 up to a few hundred.

For very large degrees (thousands or more) and parallel environments, Aberth-Ehrlich becomes preferable. For polynomials with many multiple or close roots, Laguerre's global convergence helps. For real-coefficient polynomials where complex arithmetic must be avoided, use Bairstow. Users can call jenkinsTraub, laguerre, durandKernerAberth, or bairstow directly to override the default.

Exact-coefficient (Int / Rational) cases

When roots are algebraically expressible (e.g. rational roots), use the factorization pipeline to extract linear factors before invoking numerical root finding. For floating-point root extraction the path delegates to double or Complex<double> numerical methods.

References

• Jenkins, M. A. and Traub, J. F. (1970). "A three-stage algorithm for real polynomials using quadratic iteration". SIAM Journal on Numerical Analysis, 7(4), 545–566.
• Aberth, O. (1973). "Iteration methods for finding all zeros of a polynomial simultaneously". Mathematics of Computation, 27(122), 339–344.
• Bini, D. A. and Robol, L. (2014). "Solving secular and polynomial equations: A multiprecision algorithm". Journal of Computational and Applied Mathematics, 272, 276–292.
• Adams, D. A. (1967). "A stopping criterion for polynomial root finding". Communications of the ACM, 10(10), 655–658.