companion-matrix root-finding

How does a numerical library find ALL the roots of a polynomial at once — say all five roots of a degree-5 polynomial, including complex ones? Not by iterating Newton five times and hoping. The standard trick is surprising and elegant: turn the polynomial into a matrix whose EIGENVALUES are exactly the polynomial's roots, then use a mature, reliable eigenvalue algorithm to find them all in one go.

Every monic polynomial p(x) = x^n + c_{n-1} x^{n-1} + ... + c_1 x + c_0 has a COMPANION MATRIX — an n-by-n matrix specially built from its coefficients (with 1's on the subdiagonal and the negated coefficients in the last column or top row). Its defining property is that the characteristic polynomial of this matrix is exactly p(x). Therefore the eigenvalues of the companion matrix ARE the roots of p, real and complex alike. To find the roots you hand the companion matrix to the QR algorithm (the workhorse eigenvalue solver in LAPACK), which delivers all n roots robustly. This is precisely how widely used 'roots' functions in numerical software work. An alternative all-roots iteration, the Durand-Kerner (Weierstrass) method, refines all root estimates simultaneously using a clever fixed-point update rather than a matrix.

Two honest points round this out. First, it is delightfully reliable because it offloads the hard work onto the extremely well-tested QR eigenvalue machinery rather than juggling many separate Newton iterations with their convergence headaches. Second, there is a genuine conditioning warning: the roots of a polynomial can be exquisitely sensitive to tiny changes in its coefficients. Wilkinson's famous example — a degree-20 polynomial with roots at 1, 2, ..., 20 — has roots that move enormously when a coefficient is perturbed in the 10th decimal place, so passing through the coefficient form can be ill-conditioned. When you can, it is safer to work with a polynomial's roots or factored form directly than to reconstruct it from coefficients.