From 9a4a08da46fb88a64550c65c43231df7f5193276 Mon Sep 17 00:00:00 2001
From: Manuel Yguel <manuel.yguel@gmail.com>
Date: Thu, 25 Mar 2010 03:15:05 +0100
Subject: Creation of the Polynomials module with the following features: *
 convenient functions:   - Horner and stabilized Horner evaluation   -
 polynomial coefficients from a set of given roots   - Cauchy bounds * a QR
 based polynomial solver

---
 unsupported/doc/examples/PolynomialSolver1.cpp | 53 ++++++++++++++++++++++++++
 unsupported/doc/examples/PolynomialUtils1.cpp  | 20 ++++++++++
 2 files changed, 73 insertions(+)
 create mode 100644 unsupported/doc/examples/PolynomialSolver1.cpp
 create mode 100644 unsupported/doc/examples/PolynomialUtils1.cpp

(limited to 'unsupported/doc')

diff --git a/unsupported/doc/examples/PolynomialSolver1.cpp b/unsupported/doc/examples/PolynomialSolver1.cpp
new file mode 100644
index 000000000..c875c9361
--- /dev/null
+++ b/unsupported/doc/examples/PolynomialSolver1.cpp
@@ -0,0 +1,53 @@
+#include <unsupported/Eigen/Polynomials>
+#include <vector>
+#include <iostream>
+
+using namespace Eigen;
+using namespace std;
+
+int main()
+{
+  typedef Matrix<double,5,1> Vector5d;
+
+  Vector5d roots = Vector5d::Random();
+  cout << "Roots: " << roots.transpose() << endl;
+  Eigen::Matrix<double,6,1> polynomial;
+  roots_to_monicPolynomial( roots, polynomial );
+
+  PolynomialSolver<double,5> psolve( polynomial );
+  cout << "Complex roots: " << psolve.roots().transpose() << endl;
+
+  std::vector<double> realRoots;
+  psolve.realRoots( realRoots );
+  Map<Vector5d> 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<float,7,1> 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<float,6> psolvef( hardCase_polynomial );
+  cout << "Complex roots: " << psolvef.roots().transpose() << endl;
+  Eigen::Matrix<float,6,1> 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<double,6> psolve6d( hardCase_polynomial.cast<double>() );
+  cout << "Complex roots: " << psolve6d.roots().transpose() << endl;
+  for( int i=0; i<6; ++i )
+  {
+    std::complex<float> 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<float> castedRoot( psolve6d.roots()[5].real(), psolve6d.roots()[5].imag() );
+  cout << "Norm of the difference: " << ei_abs( psolvef.roots()[5] - castedRoot ) << endl;
+}
diff --git a/unsupported/doc/examples/PolynomialUtils1.cpp b/unsupported/doc/examples/PolynomialUtils1.cpp
new file mode 100644
index 000000000..dbfe520b5
--- /dev/null
+++ b/unsupported/doc/examples/PolynomialUtils1.cpp
@@ -0,0 +1,20 @@
+#include <unsupported/Eigen/Polynomials>
+#include <iostream>
+
+using namespace Eigen;
+using namespace std;
+
+int main()
+{
+  Vector4d roots = Vector4d::Random();
+  cout << "Roots: " << roots.transpose() << endl;
+  Eigen::Matrix<double,5,1> polynomial;
+  roots_to_monicPolynomial( roots, polynomial );
+  cout << "Polynomial: ";
+  for( int i=0; i<4; ++i ){ cout << polynomial[i] << ".x^" << i << "+ "; }
+  cout << polynomial[4] << ".x^4" << endl;
+  Vector4d evaluation;
+  for( int i=0; i<4; ++i ){
+    evaluation[i] = poly_eval( polynomial, roots[i] ); }
+  cout << "Evaluation of the polynomial at the roots: " << evaluation.transpose();
+}
-- 
cgit v1.2.3