00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_POLYNOMIAL_SOLVER_H
00026 #define EIGEN_POLYNOMIAL_SOLVER_H
00027
00028 namespace Eigen {
00029
00043 template< typename _Scalar, int _Deg >
00044 class PolynomialSolverBase
00045 {
00046 public:
00047 EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
00048
00049 typedef _Scalar Scalar;
00050 typedef typename NumTraits<Scalar>::Real RealScalar;
00051 typedef std::complex<RealScalar> RootType;
00052 typedef Matrix<RootType,_Deg,1> RootsType;
00053
00054 typedef DenseIndex Index;
00055
00056 protected:
00057 template< typename OtherPolynomial >
00058 inline void setPolynomial( const OtherPolynomial& poly ){
00059 m_roots.resize(poly.size()); }
00060
00061 public:
00062 template< typename OtherPolynomial >
00063 inline PolynomialSolverBase( const OtherPolynomial& poly ){
00064 setPolynomial( poly() ); }
00065
00066 inline PolynomialSolverBase(){}
00067
00068 public:
00070 inline const RootsType& roots() const { return m_roots; }
00071
00072 public:
00083 template<typename Stl_back_insertion_sequence>
00084 inline void realRoots( Stl_back_insertion_sequence& bi_seq,
00085 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00086 {
00087 bi_seq.clear();
00088 for(Index i=0; i<m_roots.size(); ++i )
00089 {
00090 if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold ){
00091 bi_seq.push_back( m_roots[i].real() ); }
00092 }
00093 }
00094
00095 protected:
00096 template<typename squaredNormBinaryPredicate>
00097 inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const
00098 {
00099 Index res=0;
00100 RealScalar norm2 = internal::abs2( m_roots[0] );
00101 for( Index i=1; i<m_roots.size(); ++i )
00102 {
00103 const RealScalar currNorm2 = internal::abs2( m_roots[i] );
00104 if( pred( currNorm2, norm2 ) ){
00105 res=i; norm2=currNorm2; }
00106 }
00107 return m_roots[res];
00108 }
00109
00110 public:
00114 inline const RootType& greatestRoot() const
00115 {
00116 std::greater<Scalar> greater;
00117 return selectComplexRoot_withRespectToNorm( greater );
00118 }
00119
00123 inline const RootType& smallestRoot() const
00124 {
00125 std::less<Scalar> less;
00126 return selectComplexRoot_withRespectToNorm( less );
00127 }
00128
00129 protected:
00130 template<typename squaredRealPartBinaryPredicate>
00131 inline const RealScalar& selectRealRoot_withRespectToAbsRealPart(
00132 squaredRealPartBinaryPredicate& pred,
00133 bool& hasArealRoot,
00134 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00135 {
00136 hasArealRoot = false;
00137 Index res=0;
00138 RealScalar abs2(0);
00139
00140 for( Index i=0; i<m_roots.size(); ++i )
00141 {
00142 if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold )
00143 {
00144 if( !hasArealRoot )
00145 {
00146 hasArealRoot = true;
00147 res = i;
00148 abs2 = m_roots[i].real() * m_roots[i].real();
00149 }
00150 else
00151 {
00152 const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real();
00153 if( pred( currAbs2, abs2 ) )
00154 {
00155 abs2 = currAbs2;
00156 res = i;
00157 }
00158 }
00159 }
00160 else
00161 {
00162 if( internal::abs( m_roots[i].imag() ) < internal::abs( m_roots[res].imag() ) ){
00163 res = i; }
00164 }
00165 }
00166 return internal::real_ref(m_roots[res]);
00167 }
00168
00169
00170 template<typename RealPartBinaryPredicate>
00171 inline const RealScalar& selectRealRoot_withRespectToRealPart(
00172 RealPartBinaryPredicate& pred,
00173 bool& hasArealRoot,
00174 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00175 {
00176 hasArealRoot = false;
00177 Index res=0;
00178 RealScalar val(0);
00179
00180 for( Index i=0; i<m_roots.size(); ++i )
00181 {
00182 if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold )
00183 {
00184 if( !hasArealRoot )
00185 {
00186 hasArealRoot = true;
00187 res = i;
00188 val = m_roots[i].real();
00189 }
00190 else
00191 {
00192 const RealScalar curr = m_roots[i].real();
00193 if( pred( curr, val ) )
00194 {
00195 val = curr;
00196 res = i;
00197 }
00198 }
00199 }
00200 else
00201 {
00202 if( internal::abs( m_roots[i].imag() ) < internal::abs( m_roots[res].imag() ) ){
00203 res = i; }
00204 }
00205 }
00206 return internal::real_ref(m_roots[res]);
00207 }
00208
00209 public:
00224 inline const RealScalar& absGreatestRealRoot(
00225 bool& hasArealRoot,
00226 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00227 {
00228 std::greater<Scalar> greater;
00229 return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
00230 }
00231
00232
00247 inline const RealScalar& absSmallestRealRoot(
00248 bool& hasArealRoot,
00249 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00250 {
00251 std::less<Scalar> less;
00252 return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
00253 }
00254
00255
00270 inline const RealScalar& greatestRealRoot(
00271 bool& hasArealRoot,
00272 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00273 {
00274 std::greater<Scalar> greater;
00275 return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
00276 }
00277
00278
00293 inline const RealScalar& smallestRealRoot(
00294 bool& hasArealRoot,
00295 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00296 {
00297 std::less<Scalar> less;
00298 return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold );
00299 }
00300
00301 protected:
00302 RootsType m_roots;
00303 };
00304
00305 #define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE ) \
00306 typedef typename BASE::Scalar Scalar; \
00307 typedef typename BASE::RealScalar RealScalar; \
00308 typedef typename BASE::RootType RootType; \
00309 typedef typename BASE::RootsType RootsType;
00310
00311
00312
00342 template< typename _Scalar, int _Deg >
00343 class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
00344 {
00345 public:
00346 EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
00347
00348 typedef PolynomialSolverBase<_Scalar,_Deg> PS_Base;
00349 EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
00350
00351 typedef Matrix<Scalar,_Deg,_Deg> CompanionMatrixType;
00352 typedef EigenSolver<CompanionMatrixType> EigenSolverType;
00353
00354 public:
00356 template< typename OtherPolynomial >
00357 void compute( const OtherPolynomial& poly )
00358 {
00359 assert( Scalar(0) != poly[poly.size()-1] );
00360 internal::companion<Scalar,_Deg> companion( poly );
00361 companion.balance();
00362 m_eigenSolver.compute( companion.denseMatrix() );
00363 m_roots = m_eigenSolver.eigenvalues();
00364 }
00365
00366 public:
00367 template< typename OtherPolynomial >
00368 inline PolynomialSolver( const OtherPolynomial& poly ){
00369 compute( poly ); }
00370
00371 inline PolynomialSolver(){}
00372
00373 protected:
00374 using PS_Base::m_roots;
00375 EigenSolverType m_eigenSolver;
00376 };
00377
00378
00379 template< typename _Scalar >
00380 class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1>
00381 {
00382 public:
00383 typedef PolynomialSolverBase<_Scalar,1> PS_Base;
00384 EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
00385
00386 public:
00388 template< typename OtherPolynomial >
00389 void compute( const OtherPolynomial& poly )
00390 {
00391 assert( Scalar(0) != poly[poly.size()-1] );
00392 m_roots[0] = -poly[0]/poly[poly.size()-1];
00393 }
00394
00395 protected:
00396 using PS_Base::m_roots;
00397 };
00398
00399 }
00400
00401 #endif // EIGEN_POLYNOMIAL_SOLVER_H