From 04bd67a0f7c240118fb7426fc8b8c785a068576b Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sat, 4 May 2024 09:58:30 -0700 Subject: [PATCH] Polynomial roots via eigenvalues of the companion matrix --- doc/internals/polynomial.qbk | 14 ++- include/boost/math/tools/polynomial.hpp | 67 +++++++++++++++ test/Jamfile.v2 | 1 + test/polynomial_roots_test.cpp | 109 ++++++++++++++++++++++++ test/test_autodiff_5.cpp | 2 +- 5 files changed, 190 insertions(+), 3 deletions(-) create mode 100644 test/polynomial_roots_test.cpp 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> 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 #include +#if __has_include() +#include // roots are complex numbers. +#include +#endif + namespace boost{ namespace math{ namespace tools{ template @@ -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() + /* + * 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::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; + if (m_data.size() == 1) { + return std::vector(); + } + // 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 C(n, n); + C << Eigen::Matrix::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 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 zeros(eigen_zeros.size()); + for (std::size_t i = 0; i < zeros.size(); ++i) { + zeros[i] = eigen_zeros[i]; + } + return zeros; + } +#endif + private: template 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 : : no ] ] + [ run polynomial_roots_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../../multiprecision/config//has_eigen : : 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "math_unit_test.hpp" +using boost::math::tools::polynomial; + +#if __has_include() + +void test_random_coefficients() { + std::random_device rd; + uint32_t seed = rd(); + std::mt19937_64 mt(seed); + std::uniform_real_distribution unif(-1, 1); + std::size_t n = seed % 3 + 3; + std::vector coeffs(n, std::numeric_limits::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 root_product = -1; + std::complex root_sum = 0.0; + for (auto const & root : roots) { + root_product *= static_cast>(root); + root_sum += static_cast>(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 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 root_product = -1; + std::complex root_sum = 0.0; + for (auto const & root : roots) { + root_product *= static_cast>(root); + root_sum += static_cast>(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 c = 0.0; + for (std::size_t i = 0; i < roots.size(); ++i) { + auto ri = static_cast>(roots[i]); + for (std::size_t j = i + 1; j < roots.size(); ++j) { + c += ri*static_cast>(roots[j]); + } + } + CHECK_ULP_CLOSE(p[p.size()-3], static_cast(c.real()), 10); + CHECK_ABSOLUTE_ERROR(0.0, c.imag(), 1e-8); + +} + +template +void test_singular_companion() +{ + std::vector 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::epsilon()); + CHECK_ABSOLUTE_ERROR(T(0), roots[i].imag(), std::numeric_limits::epsilon()); + } +} + + +int main() +{ + test_random_coefficients(); + test_singular_companion(); + test_singular_companion(); + 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(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(x)).derivative(0u),