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

1 files changed, 132 insertions, 0 deletions

diff --git a/unsupported/Eigen/Polynomials b/unsupported/Eigen/Polynomials
new file mode 100644
index 000000000..9e1f6b759
--- /dev/null
+++ b/unsupported/Eigen/Polynomials
@@ -0,0 +1,132 @@
+#ifndef EIGEN_POLYNOMIALS_MODULE_H
+#define EIGEN_POLYNOMIALS_MODULE_H
+
+#include <Eigen/Core>
+
+#include <Eigen/src/Core/util/DisableMSVCWarnings.h>
+
+#include <Eigen/QR>
+
+// Note that EIGEN_HIDE_HEAVY_CODE has to be defined per module
+#if (defined EIGEN_EXTERN_INSTANTIATIONS) && (EIGEN_EXTERN_INSTANTIATIONS>=2)
+  #ifndef EIGEN_HIDE_HEAVY_CODE
+  #define EIGEN_HIDE_HEAVY_CODE
+  #endif
+#elif defined EIGEN_HIDE_HEAVY_CODE
+  #undef EIGEN_HIDE_HEAVY_CODE
+#endif
+
+namespace Eigen {
+
+/** \ingroup Unsupported_modules
+  * \defgroup Polynomials_Module Polynomials module
+  *
+  * \nonstableyet
+  *
+  * \brief This module provides a QR based polynomial solver.
+	*
+  * To use this module, add
+  * \code
+  * #include <unsupported/Eigen/Polynomials>
+  * \endcode
+	* at the start of your source file.
+  */
+
+#include "src/Polynomials/PolynomialUtils.h"
+#include "src/Polynomials/Companion.h"
+#include "src/Polynomials/PolynomialSolver.h"
+
+/**
+	\page polynomials Polynomials defines functions for dealing with polynomials
+	and a QR based polynomial solver.
+	\ingroup Polynomials_Module
+
+	The remainder of the page documents first the functions for evaluating, computing
+	polynomials, computing estimates about polynomials and next the QR based polynomial
+	solver.
+
+	\section polynomialUtils convenient functions to deal with polynomials
+	\subsection roots_to_monicPolynomial
+	The function
+	\code
+	void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly )
+	\endcode
+	computes the coefficients \f$ a_i \f$ of
+
+	\f$ p(x) = a_0 + a_{1}x + ... + a_{n-1}x^{n-1} + x^n \f$
+
+	where \f$ p \f$ is known through its roots i.e. \f$ p(x) = (x-r_1)(x-r_2)...(x-r_n) \f$.
+
+	\subsection poly_eval
+	The function
+	\code
+	T poly_eval( const Polynomials& poly, const T& x )
+	\endcode
+	evaluates a polynomial at a given point using stabilized H&ouml;rner method.
+
+	The following code computes the coefficients in the monomial basis of the monic polynomial given through its roots then evaluate it at those roots.
+
+	\include PolynomialUtils1.cpp
+  Output: \verbinclude PolynomialUtils1.out
+
+	\subsection Cauchy bounds
+	The function
+	\code
+	Real cauchy_max_bound( const Polynomial& poly )
+	\endcode
+	provides a maximum bound (the Cauchy one: \f$C(p)\f$) for the absolute value of a root of the given polynomial i.e.
+	\f$ \forall r_i \f$ root of \f$ p(x) = \sum_{k=0}^d a_k x^k \f$,
+	\f$ |r_i| \le C(p) = \sum_{k=0}^{d} \left | \frac{a_k}{a_d} \right | \f$
+	The leading coefficient \f$ p \f$: should be non zero \f$a_d \neq 0\f$.
+
+
+	The function
+	\code
+	Real cauchy_min_bound( const Polynomial& poly )
+	\endcode
+	provides a minimum bound (the Cauchy one: \f$c(p)\f$) for the absolute value of a non zero root of the given polynomial i.e.
+	\f$ \forall r_i \neq 0 \f$ root of \f$ p(x) = \sum_{k=0}^d a_k x^k \f$,
+	\f$ |r_i| \ge c(p) = \left( \sum_{k=0}^{d} \left | \frac{a_k}{a_0} \right | \right)^{-1} \f$
+
+
+
+
+	\section QR polynomial solver class
+	Computes the complex roots of a polynomial by computing the eigenvalues of the associated companion matrix with the QR algorithm.
+	
+	The roots of \f$ p(x) = a_0 + a_1 x + a_2 x^2 + a_{3} x^3 + x^4 \f$ are the eigenvalues of
+	\f$
+	\left [
+	\begin{array}{cccc}
+	0 & 0 &  0 & a_0 \\
+	1 & 0 &  0 & a_1 \\
+	0 & 1 &  0 & a_2 \\
+	0 & 0 &  1 & a_3
+	\end{array} \right ]
+	\f$
+
+	However, the QR algorithm is not guaranteed to converge when there are several eigenvalues with same modulus.
+
+	Therefore the current polynomial solver is guaranteed to provide a correct result only when the complex roots \f$r_1,r_2,...,r_d\f$ have distinct moduli i.e.
+	
+	\f$ \forall i,j \in [1;d],~ \| r_i \| \neq \| r_j \| \f$.
+
+	With 32bit (float) floating types this problem shows up frequently.
+  However, almost always, correct accuracy is reached even in these cases for 64bit
+  (double) floating types and small polynomial degree (<20).
+
+	\include PolynomialSolver1.cpp
+	In the example a polynomial with almost conjugate roots is provided to the solver.
+	Those roots have almost same module therefore the QR algorithm failed to converge: the accuracy
+	of the last root is bad.
+
+	This problem is less visible with double.
+  Output: \verbinclude PolynomialSolver1.out
+*/
+
+} // namespace Eigen
+
+#include <Eigen/src/Core/util/EnableMSVCWarnings.h>
+
+#endif // EIGEN_POLYNOMIALS_MODULE_H
+/* vim: set filetype=cpp et sw=2 ts=2 ai: */