Polynomial roots via eigenvalues of the companion matrix

doc/internals/polynomial.qbk

@@ -48,7 +48,8 @@
 
       T operator()(T z) const;
 
-
+      // Only available if Eigen library is available:
+      std::vector<std::complex<T>> roots() const;
 
       // modify:
       void set_zero();
@@ -183,11 +184,17 @@ for output
 
 The source code is at [@../../example/polynomial_arithmetic.cpp polynomial_arithmetic.cpp]
 
-[endsect] [/section:polynomials Polynomials]
+Polynomial roots are recovered by the eigenvalues of the companion matrix.
+This requires that the Eigen C++ library to be available; otherwise calling `.roots()` is a compile-time error.
+
+[endsect]
+
+[/section:polynomials Polynomials]
 
 [/
   Copyright 2006 John Maddock and Paul A. Bristow.
   Copyright 2015 Jeremy William Murphy.
+  Copyright 2024 Nick Thompson.
 
   Distributed under the Boost Software License, Version 1.0.
   (See accompanying file LICENSE_1_0.txt or copy at

include/boost/math/tools/polynomial.hpp

@@ -1,5 +1,6 @@
 //  (C) Copyright John Maddock 2006.
 //  (C) Copyright Jeremy William Murphy 2015.
+//  (C) Copyright Nick Thompson 2024.
 
 
 //  Use, modification and distribution are subject to the
@@ -29,6 +30,10 @@
 #include <type_traits>
 #include <iterator>
 
+#if __has_include(<Eigen/Eigenvalues>)
+#include <Eigen/Eigenvalues>
+#endif
+
 namespace boost{ namespace math{ namespace tools{
 
 template <class T>
@@ -575,6 +580,67 @@ class polynomial
       m_data.erase(std::find_if(m_data.rbegin(), m_data.rend(), [](const T& x)->bool { return x != T(0); }).base(), m_data.end());
    }
 
+#if __has_include(<Eigen/Eigenvalues>)
+   /*
+    * Polynomial root recovery by the eigenvalues of the companion matrix.
+    * N.B.: This algorithm is not the state of the art; a faster algorithm is
+    * "Fast and backward stable computation of roots of polynomials" by Aurentz et al.
+    */
+  [[nodiscard]] auto roots() const {
+    // At least as of Eigen 3.4.0, we cannot provide the eigensolver with complex numbers.
+    // Hence this static assert:
+    static_assert(std::is_floating_point<T>::value, "Must be a floating-point type to recover the roots");
+    using std::complex;
+    using Complex = complex<T>;
+    if (m_data.size() == 1) {
+      return std::vector<Complex>();
+    }
+    // There is a temptation to split off the linear and quadratic case, since
+    // they are so easy. Resist the temptation! Your best unit tests will become
+    // tautological.
+    std::size_t n = m_data.size() - 1;
+    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> C(n, n);
+    C << Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>::Zero(n,n);
+    for (std::size_t i = 0; i < n; ++i) {
+      // Remember the class invariant m_data.back() != 0 from the normalize() call?
+      // Reaping blessings right here y'all:
+      C(i, n - 1) = -m_data[i] / m_data.back();
+    }
+    for (std::size_t i = 0; i < n - 1; ++i) {
+      C(i + 1, i) = 1;
+    }
+    Eigen::EigenSolver<decltype(C)> es;
+    es.compute(C, /*computeEigenvectors=*/ false);
+    auto info = es.info();
+    if (info != Eigen::ComputationInfo::Success) {
+      std::ostringstream oss;
+      oss << __FILE__ << ":" << __LINE__ << ":" << __func__ << ": Eigen's eigensolver did not succeed.";
+      switch (info) {
+        case Eigen::ComputationInfo::NumericalIssue:
+          oss << " Problem: numerical issue.";
+          break;
+        case Eigen::ComputationInfo::NoConvergence:
+          oss << " Problem: no convergence.";
+          break;
+        case Eigen::ComputationInfo::InvalidInput:
+          oss << " Problem: Invalid input.";
+          break;
+        default:
+          oss << " Problem: Unknown.";
+      }
+      BOOST_MATH_THROW_EXCEPTION(std::runtime_error(oss.str()));
+    }
+    // Don't want to expose Eigen types to the rest of the world;
+    // Eigen is a detail of this algorithm, so big sad copy:
+    auto eigen_zeros = es.eigenvalues();
+    std::vector<Complex> zeros(eigen_zeros.size());
+    for (std::size_t i = 0; i < zeros.size(); ++i) {
+      zeros[i] = eigen_zeros[i];
+    }
+    return zeros;
+  }
+#endif
+
 private:
     template <class U, class R>
     polynomial& addition(const U& value, R op)

test/polynomial_roots_test.cpp

@@ -0,0 +1,116 @@
+/*
+ * Copyright Nick Thompson, 2024
+ * Use, modification and distribution are subject to the
+ * Boost Software License, Version 1.0. (See accompanying file
+ * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+#include <vector>
+#include <iostream>
+#include <list>
+#include <random>
+#include <cmath>
+#include <complex>
+#include <utility>
+#include <limits>
+#include <algorithm>
+#include <boost/math/tools/polynomial.hpp>
+using boost::math::tools::polynomial;
+#ifdef BOOST_HAS_FLOAT128
+#include <boost/multiprecision/float128.hpp>
+using boost::multiprecision::float128;
+#endif
+#include "math_unit_test.hpp"
+
+#if __has_include(<Eigen/Eigenvalues>)
+
+void test_random_coefficients() {
+    std::random_device rd;
+    uint32_t seed = rd(); 
+    std::mt19937_64 mt(seed);
+    std::uniform_real_distribution<double> unif(-1, 1);
+    std::size_t n = seed % 3 + 3;
+    std::vector<double> coeffs(n, std::numeric_limits<double>::quiet_NaN());
+    for (std::size_t i = 0; i < coeffs.size(); ++i) {
+        coeffs[i] = unif(mt);
+    }
+    coeffs[coeffs.size() -1] = 1.0;
+    auto p = polynomial(std::move(coeffs));
+    auto roots = p.roots();
+    CHECK_EQUAL(roots.size(), p.size() - 1);
+    std::complex<double> root_product = -1;
+    std::complex<double> root_sum = 0.0;
+    for (auto const & root : roots) {
+        root_product *= static_cast<std::complex<double>>(root);
+        root_sum += static_cast<std::complex<double>>(root);
+    }
+    if (p.size() & 1) {
+       root_product *= -1;
+    }
+    CHECK_ULP_CLOSE(root_product.real(), p[0], 1000);
+    CHECK_LE(root_product.imag(), 1e-6);
+
+    CHECK_LE(root_sum.imag(), 1e-7);
+    CHECK_ULP_CLOSE(root_sum.real(), -p[p.size() - 2], 1000);
+}
+
+
+
+void test_wilkinson_polynomial() {
+    // CoefficientList[Expand[(x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x - 7)*(x - 8)*(x - 9)*(x - 10)], x]
+    std::vector<float> coeffs{3628800.0, -10628640.0, 12753576.0, -8409500.0, 3416930.0, -902055.0, 157773.0, -18150.0, 1320.0, -55.0 ,1.0};
+    auto p = polynomial(std::move(coeffs));
+    auto roots = p.roots();
+    CHECK_EQUAL(roots.size(), p.size() - 1);
+    std::complex<double> root_product = -1;
+    std::complex<double> root_sum = 0.0;
+    for (auto const & root : roots) {
+        root_product *= static_cast<std::complex<double>>(root);
+        root_sum += static_cast<std::complex<double>>(root);
+    }
+    if (p.size() & 1) {
+       root_product *= -1;
+    }
+    CHECK_ABSOLUTE_ERROR(root_product.real(), double(p[0]), double(1e-3*p[0]));
+    CHECK_LE(root_product.imag(), 1e-8);
+
+    CHECK_LE(root_sum.imag(), 1e-8);
+    CHECK_ABSOLUTE_ERROR(root_sum.real(), -double(p[p.size() - 2]), 1e-5);
+
+    std::complex<double> c = 0.0;
+    for (std::size_t i = 0; i < roots.size(); ++i) {
+        auto ri = static_cast<std::complex<double>>(roots[i]);
+        for (std::size_t j = i + 1; j < roots.size(); ++j) {
+            c += ri*static_cast<std::complex<double>>(roots[j]);
+        }
+    }
+    CHECK_ULP_CLOSE(p[p.size()-3], static_cast<float>(c.real()), 10);
+    CHECK_ABSOLUTE_ERROR(0.0, c.imag(), 1e-8);
+
+}
+
+template<typename T>
+void test_singular_companion()
+{
+    std::vector<T> coeffs{0.0, 0.0, 1.0}; 
+    auto p = polynomial(std::move(coeffs));
+    auto roots = p.roots();
+    CHECK_EQUAL(roots.size(), p.size() - 1);
+    for (std::size_t i = 0; i < roots.size() - 1; ++i) {
+        CHECK_ABSOLUTE_ERROR(T(0), roots[i].real(), std::numeric_limits<T>::epsilon());
+        CHECK_ABSOLUTE_ERROR(T(0), roots[i].imag(), std::numeric_limits<T>::epsilon());
+    }
+}
+
+
+int main()
+{
+    test_random_coefficients();
+    test_singular_companion<float>();
+    test_singular_companion<double>();
+#if BOOST_HAS_FLOAT128
+    test_singular_companion<float128>();
+#endif
+    test_wilkinson_polynomial();
+    return boost::math::test::report_errors();
+}
+#endif