PolynomialSolver.h 13.7 KB
Newer Older
xuebingbing's avatar
xuebingbing committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_POLYNOMIAL_SOLVER_H
#define EIGEN_POLYNOMIAL_SOLVER_H

namespace Eigen { 

/** \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;

    typedef DenseIndex Index;

  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
    {
      using std::abs;
      bi_seq.clear();
      for(Index i=0; i<m_roots.size(); ++i )
      {
        if( 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
    {
      Index res=0;
      RealScalar norm2 = numext::abs2( m_roots[0] );
      for( Index i=1; i<m_roots.size(); ++i )
      {
        const RealScalar currNorm2 = numext::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
    {
      using std::abs;
      hasArealRoot = false;
      Index res=0;
      RealScalar abs2(0);

      for( Index i=0; i<m_roots.size(); ++i )
      {
        if( 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( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){
            res = i; }
        }
      }
      return numext::real_ref(m_roots[res]);
    }


    template<typename RealPartBinaryPredicate>
    inline const RealScalar& selectRealRoot_withRespectToRealPart(
        RealPartBinaryPredicate& pred,
        bool& hasArealRoot,
        const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
    {
      using std::abs;
      hasArealRoot = false;
      Index res=0;
      RealScalar val(0);

      for( Index i=0; i<m_roots.size(); ++i )
      {
        if( 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( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){
            res = i; }
        }
      }
      return numext::real_ref(m_roots[res]);
    }

  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 )
    {
      eigen_assert( Scalar(0) != poly[poly.size()-1] );
      internal::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 )
    {
      eigen_assert( Scalar(0) != poly[poly.size()-1] );
      m_roots[0] = -poly[0]/poly[poly.size()-1];
    }

  protected:
    using                   PS_Base::m_roots;
};

} // end namespace Eigen

#endif // EIGEN_POLYNOMIAL_SOLVER_H