Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Polynomial roots via eigenvalues of the companion matrix #1131

Open
wants to merge 1 commit into
base: develop
Choose a base branch
from

Conversation

NAThompson
Copy link
Collaborator

No description provided.

@NAThompson NAThompson force-pushed the polynomial_roots_via_companion branch from ebc1f48 to 2344f88 Compare May 4, 2024 17:00
@NAThompson
Copy link
Collaborator Author

@jzmaddock : I did test it!

@mborland , @jzmaddock : Quick question: Given a floating point type T, how to I form the associated complex type in the multiprecision case? To wit, std::complex<T> will fail for float128, cpp_bin_float_100, etc.

@NAThompson NAThompson force-pushed the polynomial_roots_via_companion branch from 2344f88 to 4e6c0ce Compare May 4, 2024 17:05
@jzmaddock
Copy link
Collaborator

Quick question: Given a floating point type T, how to I form the associated complex type in the multiprecision case? To wit, std::complex will fail for float128, cpp_bin_float_100, etc.

Ugh, seems to be a bridge we haven't crossed yet.

However, typeof(polar(T(), T()) would do it, since we never use expression templates for that.

Note that in order to use multiprecision types with Eigen we need to include <boost/multiprecision/eigen.hpp> which specializes a bunch of Eigen traits classes. I guess to avoid nasty surprises this header should include that when available, though it adds an often not needed dependency :(

The static_assert on is_floating_point would have to go as well...

@NAThompson NAThompson force-pushed the polynomial_roots_via_companion branch 2 times, most recently from 3607f34 to 5df7bca Compare May 4, 2024 18:17
@NAThompson
Copy link
Collaborator Author

The static_assert on is_floating_point would have to go as well...

Meh, who cares. Got rid of it.

Note that in order to use multiprecision types with Eigen we need to include <boost/multiprecision/eigen.hpp> which specializes a bunch of Eigen traits classes. I guess to avoid nasty surprises this header should include that when available, though it adds an often not needed dependency :(

I did the following:

#if __has_include(<Eigen/Eigenvalues>)
#include <Eigen/Eigenvalues>
#include <complex> // roots are complex numbers.
   #if __has_include(<boost/multiprecision/eigen.hpp>)
   #include <boost/multiprecision/eigen.hpp>
   #include <boost/math/cstdfloat/cstdfloat_complex_std.hpp>
   #endif
#endif

I don't feel too bad about this . . .

However, typeof(polar(T(), T()) would do it, since we never use expression templates for that.

Suffering just a bit on this one:

In file included from test/polynomial_roots_test.cpp:16:
In file included from ./include/boost/math/tools/polynomial.hpp:38:
./include/boost/math/cstdfloat/cstdfloat_complex_std.hpp:31:11: fatal error: reference to 'complex' is ambiguous
   31 |     class complex<BOOST_CSTDFLOAT_EXTENDED_COMPLEX_FLOAT_TYPE>;
      |           ^
./include/boost/math/cstdfloat/cstdfloat_complex_std.hpp:28:11: note: candidate found by name lookup is 'std::complex'
   28 |     class complex;
      |           ^
/opt/homebrew/opt/llvm/bin/../include/c++/v1/complex:260:28: note: candidate found by name lookup is 'std::__1::complex'
  260 | class _LIBCPP_TEMPLATE_VIS complex {

@jzmaddock
Copy link
Collaborator

#include <boost/math/cstdfloat/cstdfloat_complex_std.hpp>

We really shouldn't need that, and it might get rid of the error?

And of course I should have said decltype since typeof is rather last century (like me!!) :)

@NAThompson
Copy link
Collaborator Author

NAThompson commented May 4, 2024

We really shouldn't need that, and it might get rid of the error?

Got rid of it; now I'm here:

In file included from test/polynomial_roots_test.cpp:12:
/opt/homebrew/opt/llvm/bin/../include/c++/v1/complex:730:62: fatal error: no matching function for call to '__constexpr_fabs'
  730 |   _Tp __logbw  = std::__constexpr_logb(std::__constexpr_fmax(std::__constexpr_fabs(__c), std::__constexpr_fabs(__d)));
      |                                                              ^~~~~~~~~~~~~~~~~~~~~
/opt/homebrew/Cellar/eigen/3.4.0_1/include/eigen3/Eigen/src/Eigenvalues/EigenSolver.h:542:74: note: in instantiation of function template specialization 'std::operator/<boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<100>>>' requested here
  542 |         ComplexScalar cc = ComplexScalar(Scalar(0),-m_matT.coeff(n-1,n)) / ComplexScalar(m_matT.coeff(n-1,n-1)-p,q);

Looks like this is a pure Eigen/Multiprecision interop issue? For example, this:

#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <boost/multiprecision/eigen.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>

using boost::multiprecision::cpp_bin_float_100;
int main() {

    Eigen::Matrix<cpp_bin_float_100, Eigen::Dynamic, Eigen::Dynamic> A = Eigen::Matrix<cpp_bin_float_100, Eigen::Dynamic, Eigen::Dynamic>::Identity(3,3);
    Eigen::EigenSolver<decltype(A)> es;
    es.compute(A, /*computeEigenvectors=*/ false);

    auto eigs = es.eigenvalues();
    for (auto eig : eigs) {
      std::cout << eig << "\n";
    }
}

yields the same error:

In file included from /opt/homebrew/Cellar/eigen/3.4.0_1/include/eigen3/Eigen/Dense:1:
In file included from /opt/homebrew/Cellar/eigen/3.4.0_1/include/eigen3/Eigen/Core:50:
/opt/homebrew/opt/llvm/bin/../include/c++/v1/complex:730:62: fatal error: no matching function for call to '__constexpr_fabs'
  730 |   _Tp __logbw  = std::__constexpr_logb(std::__constexpr_fmax(std::__constexpr_fabs(__c), std::__constexpr_fabs(__d)));
      |                                                              ^~~~~~~~~~~~~~~~~~~~~
/opt/homebrew/Cellar/eigen/3.4.0_1/include/eigen3/Eigen/src/Eigenvalues/EigenSolver.h:542:74: note: in instantiation of function template specialization 'std::operator/<boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<100>>>' requested here
  542 |         ComplexScalar cc = ComplexScalar(Scalar(0),-m_matT.coeff(n-1,n)) / ComplexScalar(m_matT.coeff(n-1,n-1)-p,q);
      |

Naively I feel like I should submit this to Eigen and see how much they want to prioritize multiprecision integration.

@NAThompson
Copy link
Collaborator Author

Ok, upon further investigation is appears that Eigen has the following lines of code everywhere:

typedef std::complex<RealScalar> ComplexScalar;

I think we could probably convince them to do:

using std::complex;
typedef complex<RealScalar> ComplexScalar;

but I doubt we could get them to agree to

using std::complex;
using std::polar;
typedef decltype(polar(RealScalar(), RealScalar()) ComplexScalar;

So we may have to construct the bridge before we can get multiprecision integration.

@NAThompson NAThompson force-pushed the polynomial_roots_via_companion branch from 5df7bca to 04bd67a Compare May 4, 2024 19:28
Copy link

codecov bot commented May 4, 2024

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 93.69%. Comparing base (a53b013) to head (04bd67a).
Report is 471 commits behind head on develop.

Additional details and impacted files

Impacted file tree graph

@@             Coverage Diff             @@
##           develop    #1131      +/-   ##
===========================================
- Coverage    93.71%   93.69%   -0.02%     
===========================================
  Files          771      771              
  Lines        61105    61105              
===========================================
- Hits         57264    57254      -10     
- Misses        3841     3851      +10     
Files with missing lines Coverage Δ
include/boost/math/tools/polynomial.hpp 99.27% <ø> (ø)
test/test_autodiff_5.cpp 100.00% <ø> (ø)

... and 1 file with indirect coverage changes


Continue to review full report in Codecov by Sentry.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a53b013...04bd67a. Read the comment docs.

@jzmaddock
Copy link
Collaborator

Well, I was going to say we should file a bug report again Eigen... and tried but gitlab tied me up in circles and I gave up :(

@NAThompson
Copy link
Collaborator Author

NAThompson commented May 5, 2024

Well, I was going to say we should file a bug report again Eigen... and tried but gitlab tied me up in circles and I gave up :(

Done: https://gitlab.com/libeigen/eigen/-/issues/2818; I do have some worry that they'll set it to WONTFIX until we have some canonical way of providing complex<T> for multiprecision T, but maybe they know some tricks.

@cantonios
Copy link

Well, I was going to say we should file a bug report again Eigen... and tried but gitlab tied me up in circles and I gave up :(

Done: https://gitlab.com/libeigen/eigen/-/issues/2818; I do have some worry that they'll set it to WONTFIX until we have some canonical way of providing complex<T> for multiprecision T,

We won't WONTFIX ;). But you'll need to provide the merge request.

but maybe they know some tricks.

Provided a suggestion on the Eigen bug. Not really a "trick", but we have some traits you can use.

Otherwise, you could explicitly specialize std::complex<your_custom_type> - that should be legal as long as your type satisfies the NumericType requirement.

@jzmaddock
Copy link
Collaborator

Thanks @cantonios for the prompt reply.

Otherwise, you could explicitly specialize std::complex<your_custom_type> - that should be legal as long as your type satisfies the NumericType requirement.

We're looking at that now, as you say specializing complex appears to be legal, however overloading all the non-member functions (in namespace std) is definitely NOT legal and may well have unintended consequences, but I don't see a better way for now.

@cantonios
Copy link

however overloading all the non-member functions (in namespace std) is definitely NOT legal

Why do you need to? You should be able to rig it up so that non-members are found via ADL. e.g. https://godbolt.org/z/fxTb8c8cf

but I don't see a better way for now.

The better way (at least for Eigen) is to add an Eigen::NumTraits<...>::Complex member and use that.

@jzmaddock
Copy link
Collaborator

Why do you need to? You should be able to rig it up so that non-members are found via ADL.

For sure, those overloads can be found via ADL, my expectation is that most users using std::complex, will be calling std::exp. Perhaps I'm wrong?

The better way (at least for Eigen) is to add an Eigen::NumTraits<...>::Complex member and use that.

Yes, that's what we were hoping for ;)

@cantonios
Copy link

For sure, those overloads can be found via ADL, my expectation is that most users using std::complex, will be calling std::exp. Perhaps I'm wrong?

Yes, this is usually true, but arguably a bug on the users' end if they expect it to work with custom types. Besides, you still have this issue if you don't specialize std::complex, opting for your own custom CustomComplex implementation, whatever decltype(polar(...)) is.

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.

3 participants