PolynomialSolver1.cpp
Go to the documentation of this file.
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;
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 }
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
main
int main()
Definition: PolynomialSolver1.cpp:8
Eigen::PolynomialSolverBase::roots
const RootsType & roots() const
Definition: PolynomialSolver.h:55
Eigen::poly_eval
T poly_eval(const Polynomials &poly, const T &x)
Definition: PolynomialUtils.h:46
Eigen::PolynomialSolver
A polynomial solver.
Definition: PolynomialSolver.h:331
Eigen::roots_to_monicPolynomial
void roots_to_monicPolynomial(const RootVector &rv, Polynomial &poly)
Definition: PolynomialUtils.h:127
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
std
Definition: BFloat16.h:88
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
abs
#define abs(x)
Definition: datatypes.h:17
Eigen::PolynomialSolverBase::realRoots
void realRoots(Stl_back_insertion_sequence &bi_seq, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
Definition: PolynomialSolver.h:69
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:34:23