1 | #include <unsupported/Eigen/Polynomials>
|
---|
2 | #include <vector>
|
---|
3 | #include <iostream>
|
---|
4 |
|
---|
5 | using namespace Eigen;
|
---|
6 | using namespace std;
|
---|
7 |
|
---|
8 | int main()
|
---|
9 | {
|
---|
10 | typedef Matrix<double,5,1> Vector5d;
|
---|
11 |
|
---|
12 | Vector5d roots = Vector5d::Random();
|
---|
13 | cout << "Roots: " << roots.transpose() << endl;
|
---|
14 | Eigen::Matrix<double,6,1> polynomial;
|
---|
15 | roots_to_monicPolynomial( roots, polynomial );
|
---|
16 |
|
---|
17 | PolynomialSolver<double,5> psolve( polynomial );
|
---|
18 | cout << "Complex roots: " << psolve.roots().transpose() << endl;
|
---|
19 |
|
---|
20 | std::vector<double> realRoots;
|
---|
21 | psolve.realRoots( realRoots );
|
---|
22 | Map<Vector5d> mapRR( &realRoots[0] );
|
---|
23 | cout << "Real roots: " << mapRR.transpose() << endl;
|
---|
24 |
|
---|
25 | cout << endl;
|
---|
26 | cout << "Illustration of the convergence problem with the QR algorithm: " << endl;
|
---|
27 | cout << "---------------------------------------------------------------" << endl;
|
---|
28 | Eigen::Matrix<float,7,1> hardCase_polynomial;
|
---|
29 | hardCase_polynomial <<
|
---|
30 | -0.957, 0.9219, 0.3516, 0.9453, -0.4023, -0.5508, -0.03125;
|
---|
31 | cout << "Hard case polynomial defined by floats: " << hardCase_polynomial.transpose() << endl;
|
---|
32 | PolynomialSolver<float,6> psolvef( hardCase_polynomial );
|
---|
33 | cout << "Complex roots: " << psolvef.roots().transpose() << endl;
|
---|
34 | Eigen::Matrix<float,6,1> evals;
|
---|
35 | for( int i=0; i<6; ++i ){ evals[i] = std::abs( poly_eval( hardCase_polynomial, psolvef.roots()[i] ) ); }
|
---|
36 | cout << "Norms of the evaluations of the polynomial at the roots: " << evals.transpose() << endl << endl;
|
---|
37 |
|
---|
38 | cout << "Using double's almost always solves the problem for small degrees: " << endl;
|
---|
39 | cout << "-------------------------------------------------------------------" << endl;
|
---|
40 | PolynomialSolver<double,6> psolve6d( hardCase_polynomial.cast<double>() );
|
---|
41 | cout << "Complex roots: " << psolve6d.roots().transpose() << endl;
|
---|
42 | for( int i=0; i<6; ++i )
|
---|
43 | {
|
---|
44 | std::complex<float> castedRoot( psolve6d.roots()[i].real(), psolve6d.roots()[i].imag() );
|
---|
45 | evals[i] = std::abs( poly_eval( hardCase_polynomial, castedRoot ) );
|
---|
46 | }
|
---|
47 | cout << "Norms of the evaluations of the polynomial at the roots: " << evals.transpose() << endl << endl;
|
---|
48 |
|
---|
49 | cout.precision(10);
|
---|
50 | cout << "The last root in float then in double: " << psolvef.roots()[5] << "\t" << psolve6d.roots()[5] << endl;
|
---|
51 | std::complex<float> castedRoot( psolve6d.roots()[5].real(), psolve6d.roots()[5].imag() );
|
---|
52 | cout << "Norm of the difference: " << std::abs( psolvef.roots()[5] - castedRoot ) << endl;
|
---|
53 | }
|
---|