diff --git a/doc/internals/polynomial.qbk b/doc/internals/polynomial.qbk
index b5c6b5fedf..20091cbb42 100644
--- a/doc/internals/polynomial.qbk
+++ b/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,20 @@ for output
The source code is at [@../../example/polynomial_arithmetic.cpp polynomial_arithmetic.cpp]
-[endsect] [/section:polynomials Polynomials]
+[h4 Roots]
+
+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.
+In addition, the polynomial coefficients must be of floating point type, or a static assertion will fail.
+
+[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
diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp
index 6f9b9039fd..c9c48a5d62 100644
--- a/include/boost/math/tools/polynomial.hpp
+++ b/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,11 @@
#include <type_traits>
#include <iterator>
+#if __has_include(<Eigen/Eigenvalues>)
+#include <complex> // roots are complex numbers.
+#include <Eigen/Eigenvalues>
+#endif
+
namespace boost{ namespace math{ namespace tools{
template <class T>
@@ -575,6 +581,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.
+ static_assert(std::is_floating_point<T>::value, "Roots only can be recovered for floating point coefficients.");
+ // We can only support std::complex at this time; refer to the discussion
+ // in pull request #1131:
+ using Complex = std::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)
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index 394f3d36f0..e9895edb2c 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -1068,6 +1068,7 @@ test-suite interpolators :
[ run jso_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ]
[ run random_search_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] ]
[ run cma_es_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../../multiprecision/config//has_eigen : : <build>no ] ]
+ [ run polynomial_roots_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../../multiprecision/config//has_eigen : : <build>no ] ]
[ compile compile_test/random_search_incl_test.cpp : [ requires cxx17_if_constexpr cxx17_std_apply ] ]
[ compile compile_test/differential_evolution_incl_test.cpp : [ requires cxx17_if_constexpr cxx17_std_apply ] ]
[ compile compile_test/jso_incl_test.cpp : [ requires cxx17_if_constexpr cxx17_std_apply ] ]
diff --git a/test/polynomial_roots_test.cpp b/test/polynomial_roots_test.cpp
new file mode 100644
index 0000000000..bd88772a27
--- /dev/null
+++ b/test/polynomial_roots_test.cpp
@@ -0,0 +1,109 @@
+/*
+ * 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>
+#include "math_unit_test.hpp"
+using boost::math::tools::polynomial;
+
+#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>();
+ test_wilkinson_polynomial();
+ return boost::math::test::report_errors();
+}
+#endif
diff --git a/test/test_autodiff_5.cpp b/test/test_autodiff_5.cpp
index 8c875c31f5..16f2c76316 100644
--- a/test/test_autodiff_5.cpp
+++ b/test/test_autodiff_5.cpp
@@ -58,7 +58,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(chebyshev_hpp, T, all_float_types) {
auto x = x_sampler.next();
BOOST_CHECK_CLOSE(
boost::math::chebyshev_t(n, make_fvar<T, m>(x)).derivative(0u),
- boost::math::chebyshev_t(n, x), 40 * test_constants::pct_epsilon());
+ boost::math::chebyshev_t(n, x), 80 * test_constants::pct_epsilon());
// Lower accuracy with clang/apple:
BOOST_CHECK_CLOSE(
boost::math::chebyshev_u(n, make_fvar<T, m>(x)).derivative(0u),