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