GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-22 09:39:34 +08:00
Code Issues Packages Projects Releases Wiki Activity

Creation of the Polynomials module with the following features:

Browse Source 
* convenient functions:
  - Horner and stabilized Horner evaluation
  - polynomial coefficients from a set of given roots
  - Cauchy bounds
* a QR based polynomial solver
This commit is contained in:
Manuel Yguel 2010-03-25 03:15:05 +01:00
parent c68098b9be
commit 9a4a08da46
10 changed files with 1058 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
unsupported/Eigen/CMakeLists.txt
View File

@ -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}

132
unsupported/Eigen/Polynomials Normal file
View File

@ -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: */

1
unsupported/Eigen/src/CMakeLists.txt
View File

@ -5,3 +5,4 @@ ADD_SUBDIRECTORY(MoreVectorization)
# ADD_SUBDIRECTORY(FFT)
# ADD_SUBDIRECTORY(Skyline)
ADD_SUBDIRECTORY(MatrixFunctions)
ADD_SUBDIRECTORY(Polynomials)

6
unsupported/Eigen/src/Polynomials/CMakeLists.txt Normal file
View File

@ -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
)

281
unsupported/Eigen/src/Polynomials/Companion.h Normal file
View File

@ -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

395
unsupported/Eigen/src/Polynomials/PolynomialSolver.h Normal file
View File

@ -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

153
unsupported/Eigen/src/Polynomials/PolynomialUtils.h Normal file
View File

@ -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

53
unsupported/doc/examples/PolynomialSolver1.cpp Normal file
View File

@ -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;
}

20
unsupported/doc/examples/PolynomialUtils1.cpp Normal file
View File

@ -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();
}

16
unsupported/test/CMakeLists.txt
View File

@ -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}" )