diff --git a/doc/internals/polynomial.qbk b/doc/internals/polynomial.qbk index b5c6b5fedf..7098b6d535 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,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 diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp index 6f9b9039fd..e249742be8 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,10 @@ #include #include +#if __has_include() +#include +#endif + namespace boost{ namespace math{ namespace tools{ template @@ -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() + /* + * 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::value, "Must be a floating-point type to recover the roots"); + using std::complex; + using Complex = 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..9a0aaca5fe --- /dev/null +++ b/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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +using boost::math::tools::polynomial; +#ifdef BOOST_HAS_FLOAT128 +#include +using boost::multiprecision::float128; +#endif +#include "math_unit_test.hpp" + +#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(); +#if BOOST_HAS_FLOAT128 + test_singular_companion(); +#endif + test_wilkinson_polynomial(); + return boost::math::test::report_errors(); +} +#endif