aboutsummaryrefslogtreecommitdiffhomepage
path: root/unsupported
diff options
context:
space:
mode:
authorGravatar Manuel Yguel <manuel.yguel@gmail.com>2010-03-25 03:15:05 +0100
committerGravatar Manuel Yguel <manuel.yguel@gmail.com>2010-03-25 03:15:05 +0100
commit9a4a08da46fb88a64550c65c43231df7f5193276 (patch)
treea27cf42f2fbf8faaef1aaac373d7ec6e039560ec /unsupported
parentc68098b9be036b94bc603bf31646af464810ccc2 (diff)
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
Diffstat (limited to 'unsupported')
-rw-r--r--unsupported/Eigen/CMakeLists.txt2
-rw-r--r--unsupported/Eigen/Polynomials132
-rw-r--r--unsupported/Eigen/src/CMakeLists.txt1
-rw-r--r--unsupported/Eigen/src/Polynomials/CMakeLists.txt6
-rw-r--r--unsupported/Eigen/src/Polynomials/Companion.h281
-rw-r--r--unsupported/Eigen/src/Polynomials/PolynomialSolver.h395
-rw-r--r--unsupported/Eigen/src/Polynomials/PolynomialUtils.h153
-rw-r--r--unsupported/doc/examples/PolynomialSolver1.cpp53
-rw-r--r--unsupported/doc/examples/PolynomialUtils1.cpp20
-rw-r--r--unsupported/test/CMakeLists.txt16
10 files changed, 1058 insertions, 1 deletions
diff --git a/unsupported/Eigen/CMakeLists.txt b/unsupported/Eigen/CMakeLists.txt
index f9d79c51a..87cc4be1e 100644
--- a/unsupported/Eigen/CMakeLists.txt
+++ b/unsupported/Eigen/CMakeLists.txt
@@ -1,4 +1,4 @@
-set(Eigen_HEADERS AdolcForward BVH IterativeSolvers MatrixFunctions MoreVectorization AutoDiff AlignedVector3)
+set(Eigen_HEADERS AdolcForward BVH IterativeSolvers MatrixFunctions MoreVectorization AutoDiff AlignedVector3 Polynomials)
install(FILES
${Eigen_HEADERS}
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: */
diff --git a/unsupported/Eigen/src/CMakeLists.txt b/unsupported/Eigen/src/CMakeLists.txt
index a0870d3af..035a84ca8 100644
--- a/unsupported/Eigen/src/CMakeLists.txt
+++ b/unsupported/Eigen/src/CMakeLists.txt
@@ -5,3 +5,4 @@ ADD_SUBDIRECTORY(MoreVectorization)
# ADD_SUBDIRECTORY(FFT)
# ADD_SUBDIRECTORY(Skyline)
ADD_SUBDIRECTORY(MatrixFunctions)
+ADD_SUBDIRECTORY(Polynomials)
diff --git a/unsupported/Eigen/src/Polynomials/CMakeLists.txt b/unsupported/Eigen/src/Polynomials/CMakeLists.txt
new file mode 100644
index 000000000..51f13f3cb
--- /dev/null
+++ b/unsupported/Eigen/src/Polynomials/CMakeLists.txt
@@ -0,0 +1,6 @@
+FILE(GLOB Eigen_Polynomials_SRCS "*.h")
+
+INSTALL(FILES
+ ${Eigen_Polynomials_SRCS}
+ DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/Polynomials COMPONENT Devel
+ )
diff --git a/unsupported/Eigen/src/Polynomials/Companion.h b/unsupported/Eigen/src/Polynomials/Companion.h
new file mode 100644
index 000000000..7c9e9c518
--- /dev/null
+++ b/unsupported/Eigen/src/Polynomials/Companion.h
@@ -0,0 +1,281 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_COMPANION_H
+#define EIGEN_COMPANION_H
+
+// This file requires the user to include
+// * Eigen/Core
+// * Eigen/src/PolynomialSolver.h
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+
+template <typename T>
+T ei_radix(){ return 2; }
+
+template <typename T>
+T ei_radix2(){ return ei_radix<T>()*ei_radix<T>(); }
+
+template<int Size>
+struct ei_decrement_if_fixed_size
+{
+ enum {
+ ret = (Size == Dynamic) ? Dynamic : Size-1 };
+};
+
+#endif
+
+template< typename _Scalar, int _Deg >
+class ei_companion
+{
+ public:
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
+
+ enum {
+ Deg = _Deg,
+ Deg_1=ei_decrement_if_fixed_size<Deg>::ret
+ };
+
+ typedef _Scalar Scalar;
+ typedef typename NumTraits<Scalar>::Real RealScalar;
+ typedef Matrix<Scalar, Deg, 1> RightColumn;
+ //typedef DiagonalMatrix< Scalar, Deg_1, Deg_1 > BottomLeftDiagonal;
+ typedef Matrix<Scalar, Deg_1, 1> BottomLeftDiagonal;
+
+ typedef Matrix<Scalar, Deg, Deg> DenseCompanionMatrixType;
+ typedef Matrix< Scalar, _Deg, Deg_1 > LeftBlock;
+ typedef Matrix< Scalar, Deg_1, Deg_1 > BottomLeftBlock;
+ typedef Matrix< Scalar, 1, Deg_1 > LeftBlockFirstRow;
+
+ public:
+ EIGEN_STRONG_INLINE const _Scalar operator()( int row, int col ) const
+ {
+ if( m_bl_diag.rows() > col )
+ {
+ if( 0 < row ){ return m_bl_diag[col]; }
+ else{ return 0; }
+ }
+ else{ return m_monic[row]; }
+ }
+
+ public:
+ template<typename VectorType>
+ void setPolynomial( const VectorType& poly )
+ {
+ const int deg = poly.size()-1;
+ m_monic = -1/poly[deg] * poly.head(deg);
+ //m_bl_diag.setIdentity( deg-1 );
+ m_bl_diag.setOnes(deg-1);
+ }
+
+ template<typename VectorType>
+ ei_companion( const VectorType& poly ){
+ setPolynomial( poly ); }
+
+ public:
+ DenseCompanionMatrixType denseMatrix() const
+ {
+ const int deg = m_monic.size();
+ const int deg_1 = deg-1;
+ DenseCompanionMatrixType companion(deg,deg);
+ companion <<
+ ( LeftBlock(deg,deg_1)
+ << LeftBlockFirstRow::Zero(1,deg_1),
+ BottomLeftBlock::Identity(deg-1,deg-1)*m_bl_diag.asDiagonal() ).finished()
+ , m_monic;
+ return companion;
+ }
+
+
+
+ protected:
+ /** Helper function for the balancing algorithm.
+ * \returns true if the row and the column, having colNorm and rowNorm
+ * as norms, are balanced, false otherwise.
+ * colB and rowB are repectively the multipliers for
+ * the column and the row in order to balance them.
+ * */
+ bool balanced( Scalar colNorm, Scalar rowNorm,
+ bool& isBalanced, Scalar& colB, Scalar& rowB );
+
+ /** Helper function for the balancing algorithm.
+ * \returns true if the row and the column, having colNorm and rowNorm
+ * as norms, are balanced, false otherwise.
+ * colB and rowB are repectively the multipliers for
+ * the column and the row in order to balance them.
+ * */
+ bool balancedR( Scalar colNorm, Scalar rowNorm,
+ bool& isBalanced, Scalar& colB, Scalar& rowB );
+
+ public:
+ /**
+ * Balancing algorithm from B. N. PARLETT and C. REINSCH (1969)
+ * "Balancing a matrix for calculation of eigenvalues and eigenvectors"
+ * adapted to the case of companion matrices.
+ * A matrix with non zero row and non zero column is balanced
+ * for a certain norm if the i-th row and the i-th column
+ * have same norm for all i.
+ */
+ void balance();
+
+ protected:
+ RightColumn m_monic;
+ BottomLeftDiagonal m_bl_diag;
+};
+
+
+
+template< typename _Scalar, int _Deg >
+inline
+bool ei_companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm,
+ bool& isBalanced, Scalar& colB, Scalar& rowB )
+{
+ if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; }
+ else
+ {
+ //To find the balancing coefficients, if the radix is 2,
+ //one finds \f$ \sigma \f$ such that
+ // \f$ 2^{2\sigma-1} < rowNorm / colNorm \le 2^{2\sigma+1} \f$
+ // then the balancing coefficient for the row is \f$ 1/2^{\sigma} \f$
+ // and the balancing coefficient for the column is \f$ 2^{\sigma} \f$
+ rowB = rowNorm / ei_radix<Scalar>();
+ colB = Scalar(1);
+ const Scalar s = colNorm + rowNorm;
+
+ while (colNorm < rowB)
+ {
+ colB *= ei_radix<Scalar>();
+ colNorm *= ei_radix2<Scalar>();
+ }
+
+ rowB = rowNorm * ei_radix<Scalar>();
+
+ while (colNorm >= rowB)
+ {
+ colB /= ei_radix<Scalar>();
+ colNorm /= ei_radix2<Scalar>();
+ }
+
+ //This line is used to avoid insubstantial balancing
+ if ((rowNorm + colNorm) < Scalar(0.95) * s * colB)
+ {
+ isBalanced = false;
+ rowB = Scalar(1) / colB;
+ return false;
+ }
+ else{
+ return true; }
+ }
+}
+
+template< typename _Scalar, int _Deg >
+inline
+bool ei_companion<_Scalar,_Deg>::balancedR( Scalar colNorm, Scalar rowNorm,
+ bool& isBalanced, Scalar& colB, Scalar& rowB )
+{
+ if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; }
+ else
+ {
+ /**
+ * Set the norm of the column and the row to the geometric mean
+ * of the row and column norm
+ */
+ const _Scalar q = colNorm/rowNorm;
+ if( !ei_isApprox( q, _Scalar(1) ) )
+ {
+ rowB = ei_sqrt( colNorm/rowNorm );
+ colB = Scalar(1)/rowB;
+
+ isBalanced = false;
+ return false;
+ }
+ else{
+ return true; }
+ }
+}
+
+
+template< typename _Scalar, int _Deg >
+void ei_companion<_Scalar,_Deg>::balance()
+{
+ EIGEN_STATIC_ASSERT( 1 < Deg, YOU_MADE_A_PROGRAMMING_MISTAKE );
+ const int deg = m_monic.size();
+ const int deg_1 = deg-1;
+
+ bool hasConverged=false;
+ while( !hasConverged )
+ {
+ hasConverged = true;
+ Scalar colNorm,rowNorm;
+ Scalar colB,rowB;
+
+ //First row, first column excluding the diagonal
+ //==============================================
+ colNorm = ei_abs(m_bl_diag[0]);
+ rowNorm = ei_abs(m_monic[0]);
+
+ //Compute balancing of the row and the column
+ if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) )
+ {
+ m_bl_diag[0] *= colB;
+ m_monic[0] *= rowB;
+ }
+
+ //Middle rows and columns excluding the diagonal
+ //==============================================
+ for( int i=1; i<deg_1; ++i )
+ {
+ // column norm, excluding the diagonal
+ colNorm = ei_abs(m_bl_diag[i]);
+
+ // row norm, excluding the diagonal
+ rowNorm = ei_abs(m_bl_diag[i-1]) + ei_abs(m_monic[i]);
+
+ //Compute balancing of the row and the column
+ if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) )
+ {
+ m_bl_diag[i] *= colB;
+ m_bl_diag[i-1] *= rowB;
+ m_monic[i] *= rowB;
+ }
+ }
+
+ //Last row, last column excluding the diagonal
+ //============================================
+ const int ebl = m_bl_diag.size()-1;
+ VectorBlock<RightColumn,Deg_1> headMonic( m_monic, 0, deg_1 );
+ colNorm = headMonic.array().abs().sum();
+ rowNorm = ei_abs( m_bl_diag[ebl] );
+
+ //Compute balancing of the row and the column
+ if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) )
+ {
+ headMonic *= colB;
+ m_bl_diag[ebl] *= rowB;
+ }
+ }
+}
+
+
+#endif // EIGEN_COMPANION_H
diff --git a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h
new file mode 100644
index 000000000..49a5ffffa
--- /dev/null
+++ b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h
@@ -0,0 +1,395 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_POLYNOMIAL_SOLVER_H
+#define EIGEN_POLYNOMIAL_SOLVER_H
+
+/** \ingroup Polynomials_Module
+ * \class PolynomialSolverBase.
+ *
+ * \brief Defined to be inherited by polynomial solvers: it provides
+ * convenient methods such as
+ * - real roots,
+ * - greatest, smallest complex roots,
+ * - real roots with greatest, smallest absolute real value,
+ * - greatest, smallest real roots.
+ *
+ * It stores the set of roots as a vector of complexes.
+ *
+ */
+template< typename _Scalar, int _Deg >
+class PolynomialSolverBase
+{
+ public:
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
+
+ typedef _Scalar Scalar;
+ typedef typename NumTraits<Scalar>::Real RealScalar;
+ typedef std::complex<RealScalar> RootType;
+ typedef Matrix<RootType,_Deg,1> RootsType;
+
+ protected:
+ template< typename OtherPolynomial >
+ inline void setPolynomial( const OtherPolynomial& poly ){
+ m_roots.resize(poly.size()); }
+
+ public:
+ template< typename OtherPolynomial >
+ inline PolynomialSolverBase( const OtherPolynomial& poly ){
+ setPolynomial( poly() ); }
+
+ inline PolynomialSolverBase(){}
+
+ public:
+ /** \returns the complex roots of the polynomial */
+ inline const RootsType& roots() const { return m_roots; }
+
+ public:
+ /** Clear and fills the back insertion sequence with the real roots of the polynomial
+ * i.e. the real part of the complex roots that have an imaginary part which
+ * absolute value is smaller than absImaginaryThreshold.
+ * absImaginaryThreshold takes the dummy_precision associated
+ * with the _Scalar template parameter of the PolynomialSolver class as the default value.
+ *
+ * \param[out] bi_seq : the back insertion sequence (stl concept)
+ * \param[in] absImaginaryThreshold : the maximum bound of the imaginary part of a complex
+ * number that is considered as real.
+ * */
+ template<typename Stl_back_insertion_sequence>
+ inline void realRoots( Stl_back_insertion_sequence& bi_seq,
+ const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
+ {
+ bi_seq.clear();
+ for( int i=0; i<m_roots.size(); ++i )
+ {
+ if( ei_abs( m_roots[i].imag() ) < absImaginaryThreshold ){
+ bi_seq.push_back( m_roots[i].real() ); }
+ }
+ }
+
+ protected:
+ template<typename squaredNormBinaryPredicate>
+ inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const
+ {
+ int res=0;
+ RealScalar norm2 = ei_abs2( m_roots[0] );
+ for( int i=1; i<m_roots.size(); ++i )
+ {
+ const RealScalar currNorm2 = ei_abs2( m_roots[i] );
+ if( pred( currNorm2, norm2 ) ){
+ res=i; norm2=currNorm2; }
+ }
+ return m_roots[res];
+ }
+
+ public:
+ /**
+ * \returns the complex root with greatest norm.
+ */
+ inline const RootType& greatestRoot() const
+ {
+ std::greater<Scalar> greater;
+ return selectComplexRoot_withRespectToNorm( greater );
+ }
+
+ /**
+ * \returns the complex root with smallest norm.
+ */
+ inline const RootType& smallestRoot() const
+ {
+ std::less<Scalar> less;
+ return selectComplexRoot_withRespectToNorm( less );
+ }
+
+ protected:
+ template<typename squaredRealPartBinaryPredicate>
+ inline const RealScalar& selectRealRoot_withRespectToAbsRealPart(
+ squaredRealPartBinaryPredicate& pred,
+ bool& hasArealRoot,
+ const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
+ {
+ hasArealRoot = false;
+ int res=0;
+ RealScalar abs2;
+
+ for( int i=0; i<m_roots.size(); ++i )
+ {
+ if( ei_abs( m_roots[i].imag() ) < absImaginaryThreshold )
+ {
+ if( !hasArealRoot )
+ {
+ hasArealRoot = true;
+ res = i;
+ abs2 = m_roots[i].real() * m_roots[i].real();
+ }
+ else
+ {
+ const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real();
+ if( pred( currAbs2, abs2 ) )
+ {
+ abs2 = currAbs2;
+ res = i;
+ }
+ }
+ }
+ else
+ {
+ if( ei_abs( m_roots[i].imag() ) < ei_abs( m_roots[res].imag() ) ){
+ res = i; }
+ }
+ }
+ return m_roots[res].real();
+ }
+
+
+ template<typename RealPartBinaryPredicate>
+ inline const RealScalar& selectRealRoot_withRespectToRealPart(
+ RealPartBinaryPredicate& pred,
+ bool& hasArealRoot,
+ const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
+ {
+ hasArealRoot = false;
+ int res=0;
+ RealScalar val;
+
+ for( int i=0; i<m_roots.size(); ++i )
+ {
+ if( ei_abs( m_roots[i].imag() ) < absImaginaryThreshold )
+ {
+ if( !hasArealRoot )
+ {
+ hasArealRoot = true;
+ res = i;
+ val = m_roots[i].real();
+ }
+ else
+ {
+ const RealScalar curr = m_roots[i].real();
+ if( pred( curr, val ) )
+ {
+ val = curr;
+ res = i;
+ }
+ }
+ }
+ else
+ {
+ if( ei_abs( m_roots[i].imag() ) < ei_abs( m_roots[res].imag() ) ){
+ res = i; }
+ }
+ }
+ return m_roots[res].real();
+ }
+
+ public:
+ /**
+ * \returns a real root with greatest absolute magnitude.
+ * A real root is defined as the real part of a complex root with absolute imaginary
+ * part smallest than absImaginaryThreshold.
+ * absImaginaryThreshold takes the dummy_precision associated
+ * with the _Scalar template parameter of the PolynomialSolver class as the default value.
+ * If no real root is found the boolean hasArealRoot is set to false and the real part of
+ * the root with smallest absolute imaginary part is returned instead.
+ *
+ * \param[out] hasArealRoot : boolean true if a real root is found according to the
+ * absImaginaryThreshold criterion, false otherwise.
+ * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide
+ * whether or not a root is real.
+ */
+ inline const RealScalar& absGreatestRealRoot(
+ bool& hasArealRoot,
+ const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
+ {
+ std::greater<Scalar> greater;
+ return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
+ }
+
+
+ /**
+ * \returns a real root with smallest absolute magnitude.
+ * A real root is defined as the real part of a complex root with absolute imaginary
+ * part smallest than absImaginaryThreshold.
+ * absImaginaryThreshold takes the dummy_precision associated
+ * with the _Scalar template parameter of the PolynomialSolver class as the default value.
+ * If no real root is found the boolean hasArealRoot is set to false and the real part of
+ * the root with smallest absolute imaginary part is returned instead.
+ *
+ * \param[out] hasArealRoot : boolean true if a real root is found according to the
+ * absImaginaryThreshold criterion, false otherwise.
+ * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide
+ * whether or not a root is real.
+ */
+ inline const RealScalar& absSmallestRealRoot(
+ bool& hasArealRoot,
+ const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
+ {
+ std::less<Scalar> less;
+ return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
+ }
+
+
+ /**
+ * \returns the real root with greatest value.
+ * A real root is defined as the real part of a complex root with absolute imaginary
+ * part smallest than absImaginaryThreshold.
+ * absImaginaryThreshold takes the dummy_precision associated
+ * with the _Scalar template parameter of the PolynomialSolver class as the default value.
+ * If no real root is found the boolean hasArealRoot is set to false and the real part of
+ * the root with smallest absolute imaginary part is returned instead.
+ *
+ * \param[out] hasArealRoot : boolean true if a real root is found according to the
+ * absImaginaryThreshold criterion, false otherwise.
+ * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide
+ * whether or not a root is real.
+ */
+ inline const RealScalar& greatestRealRoot(
+ bool& hasArealRoot,
+ const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
+ {
+ std::greater<Scalar> greater;
+ return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
+ }
+
+
+ /**
+ * \returns the real root with smallest value.
+ * A real root is defined as the real part of a complex root with absolute imaginary
+ * part smallest than absImaginaryThreshold.
+ * absImaginaryThreshold takes the dummy_precision associated
+ * with the _Scalar template parameter of the PolynomialSolver class as the default value.
+ * If no real root is found the boolean hasArealRoot is set to false and the real part of
+ * the root with smallest absolute imaginary part is returned instead.
+ *
+ * \param[out] hasArealRoot : boolean true if a real root is found according to the
+ * absImaginaryThreshold criterion, false otherwise.
+ * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide
+ * whether or not a root is real.
+ */
+ inline const RealScalar& smallestRealRoot(
+ bool& hasArealRoot,
+ const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
+ {
+ std::less<Scalar> less;
+ return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold );
+ }
+
+ protected:
+ RootsType m_roots;
+};
+
+#define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE ) \
+ typedef typename BASE::Scalar Scalar; \
+ typedef typename BASE::RealScalar RealScalar; \
+ typedef typename BASE::RootType RootType; \
+ typedef typename BASE::RootsType RootsType;
+
+
+
+/** \ingroup Polynomials_Module
+ *
+ * \class PolynomialSolver
+ *
+ * \brief A polynomial solver
+ *
+ * Computes the complex roots of a real polynomial.
+ *
+ * \param _Scalar the scalar type, i.e., the type of the polynomial coefficients
+ * \param _Deg the degree of the polynomial, can be a compile time value or Dynamic.
+ * Notice that the number of polynomial coefficients is _Deg+1.
+ *
+ * This class implements a polynomial solver and provides convenient methods such as
+ * - real roots,
+ * - greatest, smallest complex roots,
+ * - real roots with greatest, smallest absolute real value.
+ * - greatest, smallest real roots.
+ *
+ * WARNING: this polynomial solver is experimental, part of the unsuported Eigen modules.
+ *
+ *
+ * Currently a QR algorithm is used to compute the eigenvalues of the companion matrix of
+ * the polynomial to compute its roots.
+ * This supposes that the complex moduli of the roots are all distinct: e.g. there should
+ * be no multiple roots or conjugate roots for instance.
+ * 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).
+ */
+template< typename _Scalar, int _Deg >
+class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
+{
+ public:
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
+
+ typedef PolynomialSolverBase<_Scalar,_Deg> PS_Base;
+ EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
+
+ typedef Matrix<Scalar,_Deg,_Deg> CompanionMatrixType;
+ typedef EigenSolver<CompanionMatrixType> EigenSolverType;
+
+ public:
+ /** Computes the complex roots of a new polynomial. */
+ template< typename OtherPolynomial >
+ void compute( const OtherPolynomial& poly )
+ {
+ assert( Scalar(0) != poly[poly.size()-1] );
+ ei_companion<Scalar,_Deg> companion( poly );
+ companion.balance();
+ m_eigenSolver.compute( companion.denseMatrix() );
+ m_roots = m_eigenSolver.eigenvalues();
+ }
+
+ public:
+ template< typename OtherPolynomial >
+ inline PolynomialSolver( const OtherPolynomial& poly ){
+ compute( poly ); }
+
+ inline PolynomialSolver(){}
+
+ protected:
+ using PS_Base::m_roots;
+ EigenSolverType m_eigenSolver;
+};
+
+
+template< typename _Scalar >
+class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1>
+{
+ public:
+ typedef PolynomialSolverBase<_Scalar,1> PS_Base;
+ EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
+
+ public:
+ /** Computes the complex roots of a new polynomial. */
+ template< typename OtherPolynomial >
+ void compute( const OtherPolynomial& poly )
+ {
+ assert( Scalar(0) != poly[poly.size()-1] );
+ m_roots[0] = -poly[0]/poly[poly.size()-1];
+ }
+
+ protected:
+ using PS_Base::m_roots;
+};
+
+#endif // EIGEN_POLYNOMIAL_SOLVER_H
diff --git a/unsupported/Eigen/src/Polynomials/PolynomialUtils.h b/unsupported/Eigen/src/Polynomials/PolynomialUtils.h
new file mode 100644
index 000000000..d78821f90
--- /dev/null
+++ b/unsupported/Eigen/src/Polynomials/PolynomialUtils.h
@@ -0,0 +1,153 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_POLYNOMIAL_UTILS_H
+#define EIGEN_POLYNOMIAL_UTILS_H
+
+/** \ingroup Polynomials_Module
+ * \returns the evaluation of the polynomial at x using Horner algorithm.
+ *
+ * \param[in] poly : the vector of coefficients of the polynomial ordered
+ * by degrees i.e. poly[i] is the coefficient of degree i of the polynomial
+ * e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$.
+ * \param[in] x : the value to evaluate the polynomial at.
+ *
+ * <i><b>Note for stability:</b></i>
+ * <dd> \f$ |x| \le 1 \f$ </dd>
+ */
+template <typename Polynomials, typename T>
+inline
+T poly_eval_horner( const Polynomials& poly, const T& x )
+{
+ T val=poly[poly.size()-1];
+ for( int i=poly.size()-2; i>=0; --i ){
+ val = val*x + poly[i]; }
+ return val;
+}
+
+/** \ingroup Polynomials_Module
+ * \returns the evaluation of the polynomial at x using stabilized Horner algorithm.
+ *
+ * \param[in] poly : the vector of coefficients of the polynomial ordered
+ * by degrees i.e. poly[i] is the coefficient of degree i of the polynomial
+ * e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$.
+ * \param[in] x : the value to evaluate the polynomial at.
+ */
+template <typename Polynomials, typename T>
+inline
+T poly_eval( const Polynomials& poly, const T& x )
+{
+ typedef typename NumTraits<T>::Real Real;
+
+ if( ei_abs2( x ) <= Real(1) ){
+ return poly_eval_horner( poly, x ); }
+ else
+ {
+ T val=poly[0];
+ T inv_x = T(1)/x;
+ for( int i=1; i<poly.size(); ++i ){
+ val = val*inv_x + poly[i]; }
+
+ return std::pow(x,(T)(poly.size()-1)) * val;
+ }
+}
+
+/** \ingroup Polynomials_Module
+ * \returns a maximum bound for the absolute value of any root of the polynomial.
+ *
+ * \param[in] poly : the vector of coefficients of the polynomial ordered
+ * by degrees i.e. poly[i] is the coefficient of degree i of the polynomial
+ * e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$.
+ *
+ * <i><b>Precondition:</b></i>
+ * <dd> the leading coefficient of the input polynomial poly must be non zero </dd>
+ */
+template <typename Polynomial>
+inline
+typename NumTraits<typename Polynomial::Scalar>::Real cauchy_max_bound( const Polynomial& poly )
+{
+ typedef typename Polynomial::Scalar Scalar;
+ typedef typename NumTraits<Scalar>::Real Real;
+
+ assert( Scalar(0) != poly[poly.size()-1] );
+ const Scalar inv_leading_coeff = Scalar(1)/poly[poly.size()-1];
+ Real cb(0);
+
+ for( int i=0; i<poly.size()-1; ++i ){
+ cb += ei_abs(poly[i]*inv_leading_coeff); }
+ return cb + Real(1);
+}
+
+/** \ingroup Polynomials_Module
+ * \returns a minimum bound for the absolute value of any non zero root of the polynomial.
+ * \param[in] poly : the vector of coefficients of the polynomial ordered
+ * by degrees i.e. poly[i] is the coefficient of degree i of the polynomial
+ * e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$.
+ */
+template <typename Polynomial>
+inline
+typename NumTraits<typename Polynomial::Scalar>::Real cauchy_min_bound( const Polynomial& poly )
+{
+ typedef typename Polynomial::Scalar Scalar;
+ typedef typename NumTraits<Scalar>::Real Real;
+
+ int i=0;
+ while( i<poly.size()-1 && Scalar(0) == poly(i) ){ ++i; }
+ if( poly.size()-1 == i ){
+ return Real(1); }
+
+ const Scalar inv_min_coeff = Scalar(1)/poly[i];
+ Real cb(1);
+ for( int j=i+1; j<poly.size(); ++j ){
+ cb += ei_abs(poly[j]*inv_min_coeff); }
+ return Real(1)/cb;
+}
+
+/** \ingroup Polynomials_Module
+ * Given the roots of a polynomial compute the coefficients in the
+ * monomial basis of the monic polynomial with same roots and minimal degree.
+ * If RootVector is a vector of complexes, Polynomial should also be a vector
+ * of complexes.
+ * \param[in] rv : a vector containing the roots of a polynomial.
+ * \param[out] poly : the vector of coefficients of the polynomial ordered
+ * by degrees i.e. poly[i] is the coefficient of degree i of the polynomial
+ * e.g. \f$ 3 + x^2 \f$ is stored as a vector \f$ [ 3, 0, 1 ] \f$.
+ */
+template <typename RootVector, typename Polynomial>
+void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly )
+{
+
+ typedef typename Polynomial::Scalar Scalar;
+
+ poly.setZero( rv.size()+1 );
+ poly[0] = -rv[0]; poly[1] = Scalar(1);
+ for( int i=1; i<(int)rv.size(); ++i )
+ {
+ for( int j=i+1; j>0; --j ){ poly[j] = poly[j-1] - rv[i]*poly[j]; }
+ poly[0] = -rv[i]*poly[0];
+ }
+}
+
+
+#endif // EIGEN_POLYNOMIAL_UTILS_H
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();
+}
diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt
index 9618a5a56..f01d99c36 100644
--- a/unsupported/test/CMakeLists.txt
+++ b/unsupported/test/CMakeLists.txt
@@ -24,3 +24,19 @@ if(FFTW_FOUND)
ei_add_test(FFTW "-DEIGEN_FFTW_DEFAULT " "-lfftw3 -lfftw3f -lfftw3l" )
endif(FFTW_FOUND)
+find_package(GSL)
+if(GSL_FOUND AND GSL_VERSION_MINOR LESS 9)
+ set(GSL_FOUND "")
+endif(GSL_FOUND AND GSL_VERSION_MINOR LESS 9)
+if(GSL_FOUND)
+ add_definitions("-DHAS_GSL" ${GSL_DEFINITIONS})
+ include_directories(${GSL_INCLUDE_DIR})
+ ei_add_property(EIGEN_TESTED_BACKENDS "GSL, ")
+else(GSL_FOUND)
+ ei_add_property(EIGEN_MISSING_BACKENDS "GSL, ")
+ set(GSL_LIBRARIES " ")
+endif(GSL_FOUND)
+
+ei_add_test(polynomialutils)
+ei_add_test(polynomialsolver " " "${GSL_LIBRARIES}" )
+