Skip to content

Polynomial roots - #1452

Merged
NAThompson merged 1 commit into
developfrom
polynomial_roots_imbach_moroz
Oct 12, 2026
Merged

NAThompson merged 1 commit into
developfrom
polynomial_roots_imbach_moroz

Conversation

@NAThompson

@NAThompson NAThompson commented Sep 2, 2026 •

Copy link
Copy Markdown
Collaborator

polynomial_roots: the Imbach-Moroz piecewise root finder (arXiv:2302.06244), whose sector DFTs use tools::fft; polynomial_roots_companion_matrix: the Aurentz-Mach-Vandebril-Watkins companion QR, O(n^2) and normwise backward stable. Both take any random access container of built-in or Boost.Multiprecision reals, with Python bindings.

timing

Max root error / (κ ε), κ the coefficientwise condition number:

Polynomial polynomial_roots polynomial_roots_companion_matrix numpy (LAPACK)
Kac, degree 1000 1.03 483 101
x^1000 - 1 0.71 192 393
Wilkinson 20 0.075 6.7e3 0.99
Chebyshev T_30 0.24 3.4e4 467
∏(x - 10^i), i = -8..8 0.31 4.3e20 166
x^50 - 1e300 0.74 15 1.3e17

🤖 Generated with Claude Code

@NAThompson

Copy link
Copy Markdown
Collaborator Author

https://arxiv.org/abs/2302.06244

@amlel-el-mahrouss amlel-el-mahrouss left a comment •

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have some questions and remarks, so, per void check_residuals in test/polynomial_roots_test.cpp:

  • This could be const auto& for for(auto root : roots).

per void fft in polynomial_roots.hpp:

  • This could benefit of some comments, and more testing as well due to the nested loops.

@jeremy-murphy

Copy link
Copy Markdown
Contributor

Just wondering naming, since there are so many algorithms for finding polynomial roots; you wouldn't call it imbach_moroz_polynomial_roots?

@NAThompson

Copy link
Copy Markdown
Collaborator Author

@jeremy-murphy : I believe this algorithm finishes polynomial rootfinding as it is nlogn and stable.

@jeremy-murphy

Copy link
Copy Markdown
Contributor

@jeremy-murphy : I believe this algorithm finishes polynomial rootfinding as it is nlogn and stable.

Cool, that's pretty exciting! Do you have any results on how accurate it is compared to the existing algorithms we have or against some optimal answer?

@NAThompson

Copy link
Copy Markdown
Collaborator Author

Cool, that's pretty exciting! Do you have any results on how accurate it is compared to the existing algorithms we have or against some optimal answer?

Not yet-I think Aurentz's 2015 work (which is O(n^2)) is probably the reference for stability. I do not believe this algorithm is stable in quite the same sense.

assert(std::count(zero_roots.begin(), zero_roots.end(), std::complex<double>(0, 0)) == 2);

std::mt19937_64 generator(0x5eed);
std::uniform_real_distribution<double> distribution(-2, 2);

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Have you tested how accurate this method is with multiprecision types like boost::multiprecision::cpp_bin_float?

Comment on lines +28 to +29
using std::abs;
using std::log2;

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would adding BOOST_MATH_STD_USING here and below allow multiprecision types?

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JacobHass8 : BOOST_MATH_STD_USING just brings in all math functions. using std::abs allows ADL of abs.

Comment thread test/polynomial_roots_test.cpp Outdated
Comment on lines +91 to +93
try { polynomial_roots(std::vector<double>{0, 0}); }
catch(std::domain_error const&) { threw = true; }
assert(threw);

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should all tests here use boost test functions like BOOST_CHECK_THROW? Additionally, is there any benefit in adding policies to polynomial_roots?

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure why BOOST_CHECK_THROW would matter here, what's your rationale for it?

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JacobHass8 : I eliminated Boost.Test from the library because it takes forever to compile and its error messages are not actionable. We'll use math_unit_test.hpp throughout.

Looks like it needs a CHECK_THROW macro though.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

math_unit_test.hpp does have CHECK_THROW now; the test uses it for the zero polynomial and non-finite coefficients.

@NAThompson
NAThompson changed the base branch from develop to fft October 11, 2026 18:20
@NAThompson
NAThompson force-pushed the polynomial_roots_imbach_moroz branch 2 times, most recently from 424a4ca to 553e1b9 Compare October 11, 2026 19:01
@NAThompson
NAThompson force-pushed the polynomial_roots_imbach_moroz branch from 553e1b9 to 6854237 Compare October 11, 2026 19:04
@NAThompson
NAThompson force-pushed the fft branch 3 times, most recently from 7e12e54 to ef22186 Compare October 11, 2026 19:17
@NAThompson
NAThompson force-pushed the polynomial_roots_imbach_moroz branch from 6854237 to e03a9c7 Compare October 11, 2026 19:18
@NAThompson
NAThompson force-pushed the fft branch 2 times, most recently from bcb2e29 to 31fbee6 Compare October 11, 2026 21:24
@NAThompson
NAThompson force-pushed the polynomial_roots_imbach_moroz branch from e03a9c7 to f96a5cc Compare October 11, 2026 21:25
@NAThompson
NAThompson force-pushed the polynomial_roots_imbach_moroz branch from f96a5cc to d93034b Compare October 11, 2026 23:20
@NAThompson
NAThompson force-pushed the polynomial_roots_imbach_moroz branch from d93034b to 29df02a Compare October 11, 2026 23:40
@NAThompson
NAThompson force-pushed the fft branch 3 times, most recently from f374f17 to 79f1c61 Compare October 12, 2026 05:20
@NAThompson
NAThompson force-pushed the polynomial_roots_imbach_moroz branch from 29df02a to 21f5308 Compare October 12, 2026 05:25
@NAThompson
NAThompson force-pushed the polynomial_roots_imbach_moroz branch 3 times, most recently from c3754ee to 4981897 Compare October 12, 2026 12:32
Base automatically changed from fft to develop October 12, 2026 13:21
@NAThompson
NAThompson force-pushed the polynomial_roots_imbach_moroz branch 2 times, most recently from d6da86c to 2b9c7bc Compare October 12, 2026 15:01
Comment thread doc/roots/polynomial_roots.qbk Outdated

``
#include <boost/math/tools/polynomial_roots.hpp>
#include <boost/math/tools/amvw_roots.hpp>

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@claude : I don't think amvw_roots will be scrutable to anyone. We can also do multiple implementations in a single file and dispatch to a detail, rather than have multiple files.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done: both are in polynomial_roots.hpp now (amvw_roots.hpp is gone). The public functions are polynomial_roots(c) and polynomial_roots_companion_matrix(c), each a thin wrapper over polynomial_roots_detail::piecewise / polynomial_roots_detail::companion_qr.

Comment thread doc/roots/polynomial_roots.qbk Outdated
``

template <class Real>
std::vector<std::complex<Real> >

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These should be templated on a RandomAccessContainer, and the return type deduced.

It should work with std::array, std::span, and std::vector.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done: both take RandomAccessContainer const& (anything with value_type, begin and end) and return auto, a std::vector of std::complex<Real> for built-in types or the matching Boost.Multiprecision complex type. The test covers std::vector, std::array, and std::span under C++20.

Comment thread doc/roots/polynomial_roots.qbk Outdated
template <class Real>
std::vector<std::complex<Real> >
polynomial_roots(std::vector<Real> coefficients,
unsigned precision = std::numeric_limits<Real>::digits);

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Precision is not a useful parameter here-how would a user set it?

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed, removed. The local expansion order is now fixed by std::numeric_limits<Real>::digits, which is what the default was.

@NAThompson
NAThompson force-pushed the polynomial_roots_imbach_moroz branch 2 times, most recently from 9f2e824 to ea3689a Compare October 12, 2026 15:37
polynomial_roots: the Imbach-Moroz piecewise root finder (arXiv:2302.06244), whose sector DFTs use tools::fft; polynomial_roots_companion_matrix: the Aurentz-Mach-Vandebril-Watkins companion QR, O(n^2) and normwise backward stable. Both take any random access container of built-in or Boost.Multiprecision reals, with Python bindings.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@NAThompson
NAThompson force-pushed the polynomial_roots_imbach_moroz branch from ea3689a to fadcc22 Compare October 12, 2026 15:40
@NAThompson
NAThompson merged commit 60b4ed1 into develop Oct 12, 2026
47 of 50 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants