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