1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2010 Manuel Yguel <manuel.yguel (at) gmail.com> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_POLYNOMIAL_SOLVER_H 11 #define EIGEN_POLYNOMIAL_SOLVER_H 12 13 namespace Eigen { 14 15 /** \ingroup Polynomials_Module 16 * \class PolynomialSolverBase. 17 * 18 * \brief Defined to be inherited by polynomial solvers: it provides 19 * convenient methods such as 20 * - real roots, 21 * - greatest, smallest complex roots, 22 * - real roots with greatest, smallest absolute real value, 23 * - greatest, smallest real roots. 24 * 25 * It stores the set of roots as a vector of complexes. 26 * 27 */ 28 template< typename _Scalar, int _Deg > 29 class PolynomialSolverBase 30 { 31 public: 32 EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg) 33 34 typedef _Scalar Scalar; 35 typedef typename NumTraits<Scalar>::Real RealScalar; 36 typedef std::complex<RealScalar> RootType; 37 typedef Matrix<RootType,_Deg,1> RootsType; 38 39 typedef DenseIndex Index; 40 41 protected: 42 template< typename OtherPolynomial > 43 inline void setPolynomial( const OtherPolynomial& poly ){ 44 m_roots.resize(poly.size()); } 45 46 public: 47 template< typename OtherPolynomial > 48 inline PolynomialSolverBase( const OtherPolynomial& poly ){ 49 setPolynomial( poly() ); } 50 51 inline PolynomialSolverBase(){} 52 53 public: 54 /** \returns the complex roots of the polynomial */ 55 inline const RootsType& roots() const { return m_roots; } 56 57 public: 58 /** Clear and fills the back insertion sequence with the real roots of the polynomial 59 * i.e. the real part of the complex roots that have an imaginary part which 60 * absolute value is smaller than absImaginaryThreshold. 61 * absImaginaryThreshold takes the dummy_precision associated 62 * with the _Scalar template parameter of the PolynomialSolver class as the default value. 63 * 64 * \param[out] bi_seq : the back insertion sequence (stl concept) 65 * \param[in] absImaginaryThreshold : the maximum bound of the imaginary part of a complex 66 * number that is considered as real. 67 * */ 68 template<typename Stl_back_insertion_sequence> 69 inline void realRoots( Stl_back_insertion_sequence& bi_seq, 70 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 71 { 72 using std::abs; 73 bi_seq.clear(); 74 for(Index i=0; i<m_roots.size(); ++i ) 75 { 76 if( abs( m_roots[i].imag() ) < absImaginaryThreshold ){ 77 bi_seq.push_back( m_roots[i].real() ); } 78 } 79 } 80 81 protected: 82 template<typename squaredNormBinaryPredicate> 83 inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const 84 { 85 Index res=0; 86 RealScalar norm2 = numext::abs2( m_roots[0] ); 87 for( Index i=1; i<m_roots.size(); ++i ) 88 { 89 const RealScalar currNorm2 = numext::abs2( m_roots[i] ); 90 if( pred( currNorm2, norm2 ) ){ 91 res=i; norm2=currNorm2; } 92 } 93 return m_roots[res]; 94 } 95 96 public: 97 /** 98 * \returns the complex root with greatest norm. 99 */ 100 inline const RootType& greatestRoot() const 101 { 102 std::greater<Scalar> greater; 103 return selectComplexRoot_withRespectToNorm( greater ); 104 } 105 106 /** 107 * \returns the complex root with smallest norm. 108 */ 109 inline const RootType& smallestRoot() const 110 { 111 std::less<Scalar> less; 112 return selectComplexRoot_withRespectToNorm( less ); 113 } 114 115 protected: 116 template<typename squaredRealPartBinaryPredicate> 117 inline const RealScalar& selectRealRoot_withRespectToAbsRealPart( 118 squaredRealPartBinaryPredicate& pred, 119 bool& hasArealRoot, 120 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 121 { 122 using std::abs; 123 hasArealRoot = false; 124 Index res=0; 125 RealScalar abs2(0); 126 127 for( Index i=0; i<m_roots.size(); ++i ) 128 { 129 if( abs( m_roots[i].imag() ) < absImaginaryThreshold ) 130 { 131 if( !hasArealRoot ) 132 { 133 hasArealRoot = true; 134 res = i; 135 abs2 = m_roots[i].real() * m_roots[i].real(); 136 } 137 else 138 { 139 const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real(); 140 if( pred( currAbs2, abs2 ) ) 141 { 142 abs2 = currAbs2; 143 res = i; 144 } 145 } 146 } 147 else 148 { 149 if( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){ 150 res = i; } 151 } 152 } 153 return numext::real_ref(m_roots[res]); 154 } 155 156 157 template<typename RealPartBinaryPredicate> 158 inline const RealScalar& selectRealRoot_withRespectToRealPart( 159 RealPartBinaryPredicate& pred, 160 bool& hasArealRoot, 161 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 162 { 163 using std::abs; 164 hasArealRoot = false; 165 Index res=0; 166 RealScalar val(0); 167 168 for( Index i=0; i<m_roots.size(); ++i ) 169 { 170 if( abs( m_roots[i].imag() ) < absImaginaryThreshold ) 171 { 172 if( !hasArealRoot ) 173 { 174 hasArealRoot = true; 175 res = i; 176 val = m_roots[i].real(); 177 } 178 else 179 { 180 const RealScalar curr = m_roots[i].real(); 181 if( pred( curr, val ) ) 182 { 183 val = curr; 184 res = i; 185 } 186 } 187 } 188 else 189 { 190 if( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){ 191 res = i; } 192 } 193 } 194 return numext::real_ref(m_roots[res]); 195 } 196 197 public: 198 /** 199 * \returns a real root with greatest absolute magnitude. 200 * A real root is defined as the real part of a complex root with absolute imaginary 201 * part smallest than absImaginaryThreshold. 202 * absImaginaryThreshold takes the dummy_precision associated 203 * with the _Scalar template parameter of the PolynomialSolver class as the default value. 204 * If no real root is found the boolean hasArealRoot is set to false and the real part of 205 * the root with smallest absolute imaginary part is returned instead. 206 * 207 * \param[out] hasArealRoot : boolean true if a real root is found according to the 208 * absImaginaryThreshold criterion, false otherwise. 209 * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide 210 * whether or not a root is real. 211 */ 212 inline const RealScalar& absGreatestRealRoot( 213 bool& hasArealRoot, 214 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 215 { 216 std::greater<Scalar> greater; 217 return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold ); 218 } 219 220 221 /** 222 * \returns a real root with smallest absolute magnitude. 223 * A real root is defined as the real part of a complex root with absolute imaginary 224 * part smallest than absImaginaryThreshold. 225 * absImaginaryThreshold takes the dummy_precision associated 226 * with the _Scalar template parameter of the PolynomialSolver class as the default value. 227 * If no real root is found the boolean hasArealRoot is set to false and the real part of 228 * the root with smallest absolute imaginary part is returned instead. 229 * 230 * \param[out] hasArealRoot : boolean true if a real root is found according to the 231 * absImaginaryThreshold criterion, false otherwise. 232 * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide 233 * whether or not a root is real. 234 */ 235 inline const RealScalar& absSmallestRealRoot( 236 bool& hasArealRoot, 237 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 238 { 239 std::less<Scalar> less; 240 return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold ); 241 } 242 243 244 /** 245 * \returns the real root with greatest value. 246 * A real root is defined as the real part of a complex root with absolute imaginary 247 * part smallest than absImaginaryThreshold. 248 * absImaginaryThreshold takes the dummy_precision associated 249 * with the _Scalar template parameter of the PolynomialSolver class as the default value. 250 * If no real root is found the boolean hasArealRoot is set to false and the real part of 251 * the root with smallest absolute imaginary part is returned instead. 252 * 253 * \param[out] hasArealRoot : boolean true if a real root is found according to the 254 * absImaginaryThreshold criterion, false otherwise. 255 * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide 256 * whether or not a root is real. 257 */ 258 inline const RealScalar& greatestRealRoot( 259 bool& hasArealRoot, 260 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 261 { 262 std::greater<Scalar> greater; 263 return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold ); 264 } 265 266 267 /** 268 * \returns the real root with smallest value. 269 * A real root is defined as the real part of a complex root with absolute imaginary 270 * part smallest than absImaginaryThreshold. 271 * absImaginaryThreshold takes the dummy_precision associated 272 * with the _Scalar template parameter of the PolynomialSolver class as the default value. 273 * If no real root is found the boolean hasArealRoot is set to false and the real part of 274 * the root with smallest absolute imaginary part is returned instead. 275 * 276 * \param[out] hasArealRoot : boolean true if a real root is found according to the 277 * absImaginaryThreshold criterion, false otherwise. 278 * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide 279 * whether or not a root is real. 280 */ 281 inline const RealScalar& smallestRealRoot( 282 bool& hasArealRoot, 283 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 284 { 285 std::less<Scalar> less; 286 return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold ); 287 } 288 289 protected: 290 RootsType m_roots; 291 }; 292 293 #define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE ) \ 294 typedef typename BASE::Scalar Scalar; \ 295 typedef typename BASE::RealScalar RealScalar; \ 296 typedef typename BASE::RootType RootType; \ 297 typedef typename BASE::RootsType RootsType; 298 299 300 301 /** \ingroup Polynomials_Module 302 * 303 * \class PolynomialSolver 304 * 305 * \brief A polynomial solver 306 * 307 * Computes the complex roots of a real polynomial. 308 * 309 * \param _Scalar the scalar type, i.e., the type of the polynomial coefficients 310 * \param _Deg the degree of the polynomial, can be a compile time value or Dynamic. 311 * Notice that the number of polynomial coefficients is _Deg+1. 312 * 313 * This class implements a polynomial solver and provides convenient methods such as 314 * - real roots, 315 * - greatest, smallest complex roots, 316 * - real roots with greatest, smallest absolute real value. 317 * - greatest, smallest real roots. 318 * 319 * WARNING: this polynomial solver is experimental, part of the unsuported Eigen modules. 320 * 321 * 322 * Currently a QR algorithm is used to compute the eigenvalues of the companion matrix of 323 * the polynomial to compute its roots. 324 * This supposes that the complex moduli of the roots are all distinct: e.g. there should 325 * be no multiple roots or conjugate roots for instance. 326 * With 32bit (float) floating types this problem shows up frequently. 327 * However, almost always, correct accuracy is reached even in these cases for 64bit 328 * (double) floating types and small polynomial degree (<20). 329 */ 330 template< typename _Scalar, int _Deg > 331 class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> 332 { 333 public: 334 EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg) 335 336 typedef PolynomialSolverBase<_Scalar,_Deg> PS_Base; 337 EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) 338 339 typedef Matrix<Scalar,_Deg,_Deg> CompanionMatrixType; 340 typedef EigenSolver<CompanionMatrixType> EigenSolverType; 341 342 public: 343 /** Computes the complex roots of a new polynomial. */ 344 template< typename OtherPolynomial > 345 void compute( const OtherPolynomial& poly ) 346 { 347 eigen_assert( Scalar(0) != poly[poly.size()-1] ); 348 internal::companion<Scalar,_Deg> companion( poly ); 349 companion.balance(); 350 m_eigenSolver.compute( companion.denseMatrix() ); 351 m_roots = m_eigenSolver.eigenvalues(); 352 } 353 354 public: 355 template< typename OtherPolynomial > 356 inline PolynomialSolver( const OtherPolynomial& poly ){ 357 compute( poly ); } 358 359 inline PolynomialSolver(){} 360 361 protected: 362 using PS_Base::m_roots; 363 EigenSolverType m_eigenSolver; 364 }; 365 366 367 template< typename _Scalar > 368 class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1> 369 { 370 public: 371 typedef PolynomialSolverBase<_Scalar,1> PS_Base; 372 EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) 373 374 public: 375 /** Computes the complex roots of a new polynomial. */ 376 template< typename OtherPolynomial > 377 void compute( const OtherPolynomial& poly ) 378 { 379 eigen_assert( Scalar(0) != poly[poly.size()-1] ); 380 m_roots[0] = -poly[0]/poly[poly.size()-1]; 381 } 382 383 protected: 384 using PS_Base::m_roots; 385 }; 386 387 } // end namespace Eigen 388 389 #endif // EIGEN_POLYNOMIAL_SOLVER_H 390
