blob: cd777a4e2d085e45bda8b775503d85704b8c3268 [file] [log] [blame]
 #include #include #include using namespace Eigen; using namespace std; int main() { typedef Matrix Vector5d; Vector5d roots = Vector5d::Random(); cout << "Roots: " << roots.transpose() << endl; Eigen::Matrix polynomial; roots_to_monicPolynomial( roots, polynomial ); PolynomialSolver psolve( polynomial ); cout << "Complex roots: " << psolve.roots().transpose() << endl; std::vector realRoots; psolve.realRoots( realRoots ); Map mapRR( &realRoots[0] ); cout << "Real roots: " << mapRR.transpose() << endl; cout << endl; cout << "Illustration of the convergence problem with the QR algorithm: " << endl; cout << "---------------------------------------------------------------" << endl; Eigen::Matrix hardCase_polynomial; hardCase_polynomial << -0.957, 0.9219, 0.3516, 0.9453, -0.4023, -0.5508, -0.03125; cout << "Hard case polynomial defined by floats: " << hardCase_polynomial.transpose() << endl; PolynomialSolver psolvef( hardCase_polynomial ); cout << "Complex roots: " << psolvef.roots().transpose() << endl; Eigen::Matrix evals; for( int i=0; i<6; ++i ){ evals[i] = std::abs( poly_eval( hardCase_polynomial, psolvef.roots()[i] ) ); } cout << "Norms of the evaluations of the polynomial at the roots: " << evals.transpose() << endl << endl; cout << "Using double's almost always solves the problem for small degrees: " << endl; cout << "-------------------------------------------------------------------" << endl; PolynomialSolver psolve6d( hardCase_polynomial.cast() ); cout << "Complex roots: " << psolve6d.roots().transpose() << endl; for( int i=0; i<6; ++i ) { std::complex castedRoot( psolve6d.roots()[i].real(), psolve6d.roots()[i].imag() ); evals[i] = std::abs( poly_eval( hardCase_polynomial, castedRoot ) ); } cout << "Norms of the evaluations of the polynomial at the roots: " << evals.transpose() << endl << endl; cout.precision(10); cout << "The last root in float then in double: " << psolvef.roots()[5] << "\t" << psolve6d.roots()[5] << endl; std::complex castedRoot( psolve6d.roots()[5].real(), psolve6d.roots()[5].imag() ); cout << "Norm of the difference: " << std::abs( psolvef.roots()[5] - castedRoot ) << endl; }