@@ -31,6 +31,7 @@ build.ninja
 .ninja*
 a.out
 *.json
+*.patch
 
 # Test executables built in place (extensionless), but not the sources or the test_* directories:
 /test/test_*
 
@@ -0,0 +1,87 @@
+[section:polynomial_roots Complex Roots of a Polynomial]
+
+``
+#include <boost/math/tools/polynomial_roots.hpp>
+#include <boost/math/tools/amvw_roots.hpp>
+``
+
+   template <class Real>
+   std::vector<std::complex<Real> >
+   polynomial_roots(std::vector<Real> coefficients,
+                    unsigned precision = std::numeric_limits<Real>::digits);
+
+   template <class Real>
+   std::vector<std::complex<Real> >
+   amvw_roots(std::vector<Real> coefficients);
+
+`polynomial_roots` computes all finite complex roots of a polynomial with real,
+floating-point coefficients. Coefficients are supplied in ascending order, so
+`{a0, a1, ..., an}` represents
+[@https://en.wikipedia.org/wiki/Polynomial a0 + a1 z + ... + an z^n].
+The returned order is unspecified.
+
+The implementation follows the piecewise approximation strategy of Imbach and
+Moroz: a Newton-polygon sweep divides the complex plane into radial rings, the
+rings are divided into angular sectors, and low-degree Taylor polynomials are
+factored locally. Candidate roots are polished and checked against the input
+polynomial. A global Aberth iteration is used as a safety net for unresolved or
+ill-conditioned clusters.
+
+Exact leading zero coefficients are ignored. Exact trailing zero coefficients
+produce the corresponding number of roots at zero. The zero polynomial, a zero
+precision request, and non-finite coefficients cause `std::domain_error` to be
+thrown.
+
+The optional `precision` argument is measured in bits and controls the local
+approximation order. Its default is appropriate for the coefficient type.
+
+[heading The companion QR algorithm]
+
+`amvw_roots` computes the same roots as the eigenvalues of the companion matrix,
+by the core-chasing QR algorithm of Aurentz, Mach, Vandebril and Watkins. It
+stores the companion matrix as O(n) 2x2 unitary rotations, so it takes O(n[super 2])
+time and O(n) memory, against O(n[super 3]) and O(n[super 2]) for a dense
+eigensolver. Before the iteration the variable is scaled so that the constant
+and leading coefficients have equal magnitude.
+
+The two functions make different guarantees. `amvw_roots` is proven normwise
+backward stable: all of its roots are the exact roots of a single polynomial
+whose coefficient vector is within a modest multiple of n[epsilon] of the input,
+relative to the norm of the coefficients. `polynomial_roots` checks a stronger,
+coefficientwise, bound on each root: each one is an exact root of a polynomial
+whose coefficients each differ from the input in only the last few bits.
+It has no a priori proof of convergence, but it returns only roots that
+pass this check. The coefficientwise bound matters when the roots span many
+orders of magnitude, as in [prod](z - 10[super -k]): the small roots are
+determined by the small coefficients, which a normwise perturbation may change
+completely, so `amvw_roots` gets the small roots wrong while `polynomial_roots`
+gets them to nearly full precision. For roots of comparable magnitude the
+two agree.
+
+Exact zero coefficients are handled as for `polynomial_roots`. The zero
+polynomial and non-finite coefficients cause `std::domain_error` to be
+thrown, and `std::runtime_error` is thrown in the unlikely event that the QR
+iteration fails to converge.
+
+[heading Example]
+
+   std::vector<double> p{1, 0, 1}; // 1 + z^2
+   auto roots = boost::math::tools::polynomial_roots(p);
+
+The result contains approximations to `i` and `-i`.
+
+[heading References]
+
+J. Imbach and G. Moroz, "Fast evaluation and root finding for polynomials with
+floating-point coefficients", 2023.
+
+J. L. Aurentz, T. Mach, R. Vandebril and D. S. Watkins, "Fast and backward
+stable computation of roots of polynomials", SIAM J. Matrix Anal. Appl. 36
+(2015), 942-973.
+
+J. L. Aurentz, T. Mach, L. Robol, R. Vandebril and D. S. Watkins, "Fast and
+backward stable computation of roots of polynomials, Part II: backward error
+analysis; companion matrix and companion pencil", SIAM J. Matrix Anal. Appl.
+39 (2018), 1245-1269.
+
+[endsect]
@@ -19,6 +19,7 @@ There are several fully-worked __root_finding_examples, including:
 [include roots.qbk]
 [include cubic_roots.qbk]
 [include quartic_roots.qbk]
+[include polynomial_roots.qbk]
 [include root_finding_examples.qbk]
 [include minima.qbk]
 [include root_comparison.qbk]