diff options
Diffstat (limited to 'unsupported/Eigen/src/Polynomials')
| -rw-r--r-- | unsupported/Eigen/src/Polynomials/Companion.h | 94 | ||||
| -rw-r--r-- | unsupported/Eigen/src/Polynomials/PolynomialSolver.h | 46 | ||||
| -rw-r--r-- | unsupported/Eigen/src/Polynomials/PolynomialUtils.h | 8 |
3 files changed, 87 insertions, 61 deletions
diff --git a/unsupported/Eigen/src/Polynomials/Companion.h b/unsupported/Eigen/src/Polynomials/Companion.h index b515c2920..59a15b098 100644 --- a/unsupported/Eigen/src/Polynomials/Companion.h +++ b/unsupported/Eigen/src/Polynomials/Companion.h @@ -20,12 +20,6 @@ namespace internal { #ifndef EIGEN_PARSED_BY_DOXYGEN -template <typename T> -T radix(){ return 2; } - -template <typename T> -T radix2(){ return radix<T>()*radix<T>(); } - template<int Size> struct decrement_if_fixed_size { @@ -75,8 +69,7 @@ class companion void setPolynomial( const VectorType& poly ) { const Index deg = poly.size()-1; - m_monic = -1/poly[deg] * poly.head(deg); - //m_bl_diag.setIdentity( deg-1 ); + m_monic = -poly.head(deg)/poly[deg]; m_bl_diag.setOnes(deg-1); } @@ -89,13 +82,13 @@ class companion { const Index deg = m_monic.size(); const Index deg_1 = deg-1; - DenseCompanionMatrixType companion(deg,deg); - companion << + DenseCompanionMatrixType companMat(deg,deg); + companMat << ( LeftBlock(deg,deg_1) << LeftBlockFirstRow::Zero(1,deg_1), BottomLeftBlock::Identity(deg-1,deg-1)*m_bl_diag.asDiagonal() ).finished() , m_monic; - return companion; + return companMat; } @@ -104,20 +97,20 @@ class companion /** 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 + * colB and rowB are respectively 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 ); + bool balanced( RealScalar colNorm, RealScalar rowNorm, + bool& isBalanced, RealScalar& colB, RealScalar& 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 + * colB and rowB are respectively 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 ); + bool balancedR( RealScalar colNorm, RealScalar rowNorm, + bool& isBalanced, RealScalar& colB, RealScalar& rowB ); public: /** @@ -139,10 +132,13 @@ class companion template< typename _Scalar, int _Deg > inline -bool companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm, - bool& isBalanced, Scalar& colB, Scalar& rowB ) +bool companion<_Scalar,_Deg>::balanced( RealScalar colNorm, RealScalar rowNorm, + bool& isBalanced, RealScalar& colB, RealScalar& rowB ) { - if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; } + if( RealScalar(0) == colNorm || RealScalar(0) == rowNorm + || !(numext::isfinite)(colNorm) || !(numext::isfinite)(rowNorm)){ + return true; + } else { //To find the balancing coefficients, if the radix is 2, @@ -150,53 +146,61 @@ bool companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm, // \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 / radix<Scalar>(); - colB = Scalar(1); - const Scalar s = colNorm + rowNorm; - - while (colNorm < rowB) + const RealScalar radix = RealScalar(2); + const RealScalar radix2 = RealScalar(4); + + rowB = rowNorm / radix; + colB = RealScalar(1); + const RealScalar s = colNorm + rowNorm; + + // Find sigma s.t. rowNorm / 2 <= 2^(2*sigma) * colNorm + RealScalar scout = colNorm; + while (scout < rowB) { - colB *= radix<Scalar>(); - colNorm *= radix2<Scalar>(); + colB *= radix; + scout *= radix2; } - - rowB = rowNorm * radix<Scalar>(); - - while (colNorm >= rowB) + + // We now have an upper-bound for sigma, try to lower it. + // Find sigma s.t. 2^(2*sigma) * colNorm / 2 < rowNorm + scout = colNorm * (colB / radix) * colB; // Avoid overflow. + while (scout >= rowNorm) { - colB /= radix<Scalar>(); - colNorm /= radix2<Scalar>(); + colB /= radix; + scout /= radix2; } - //This line is used to avoid insubstantial balancing - if ((rowNorm + colNorm) < Scalar(0.95) * s * colB) + // This line is used to avoid insubstantial balancing. + if ((rowNorm + radix * scout) < RealScalar(0.95) * s * colB) { isBalanced = false; - rowB = Scalar(1) / colB; + rowB = RealScalar(1) / colB; return false; } - else{ - return true; } + else + { + return true; + } } } template< typename _Scalar, int _Deg > inline -bool companion<_Scalar,_Deg>::balancedR( Scalar colNorm, Scalar rowNorm, - bool& isBalanced, Scalar& colB, Scalar& rowB ) +bool companion<_Scalar,_Deg>::balancedR( RealScalar colNorm, RealScalar rowNorm, + bool& isBalanced, RealScalar& colB, RealScalar& rowB ) { - if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; } + if( RealScalar(0) == colNorm || RealScalar(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; + const RealScalar q = colNorm/rowNorm; if( !isApprox( q, _Scalar(1) ) ) { rowB = sqrt( colNorm/rowNorm ); - colB = Scalar(1)/rowB; + colB = RealScalar(1)/rowB; isBalanced = false; return false; @@ -219,8 +223,8 @@ void companion<_Scalar,_Deg>::balance() while( !hasConverged ) { hasConverged = true; - Scalar colNorm,rowNorm; - Scalar colB,rowB; + RealScalar colNorm,rowNorm; + RealScalar colB,rowB; //First row, first column excluding the diagonal //============================================== diff --git a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h index 03198ec8e..5e0ecbb43 100644 --- a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h +++ b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h @@ -99,7 +99,7 @@ class PolynomialSolverBase */ inline const RootType& greatestRoot() const { - std::greater<Scalar> greater; + std::greater<RealScalar> greater; return selectComplexRoot_withRespectToNorm( greater ); } @@ -108,7 +108,7 @@ class PolynomialSolverBase */ inline const RootType& smallestRoot() const { - std::less<Scalar> less; + std::less<RealScalar> less; return selectComplexRoot_withRespectToNorm( less ); } @@ -126,7 +126,7 @@ class PolynomialSolverBase for( Index i=0; i<m_roots.size(); ++i ) { - if( abs( m_roots[i].imag() ) < absImaginaryThreshold ) + if( abs( m_roots[i].imag() ) <= absImaginaryThreshold ) { if( !hasArealRoot ) { @@ -144,10 +144,10 @@ class PolynomialSolverBase } } } - else + else if(!hasArealRoot) { if( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){ - res = i; } + res = i;} } } return numext::real_ref(m_roots[res]); @@ -167,7 +167,7 @@ class PolynomialSolverBase for( Index i=0; i<m_roots.size(); ++i ) { - if( abs( m_roots[i].imag() ) < absImaginaryThreshold ) + if( abs( m_roots[i].imag() ) <= absImaginaryThreshold ) { if( !hasArealRoot ) { @@ -213,7 +213,7 @@ class PolynomialSolverBase bool& hasArealRoot, const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const { - std::greater<Scalar> greater; + std::greater<RealScalar> greater; return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold ); } @@ -236,7 +236,7 @@ class PolynomialSolverBase bool& hasArealRoot, const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const { - std::less<Scalar> less; + std::less<RealScalar> less; return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold ); } @@ -259,7 +259,7 @@ class PolynomialSolverBase bool& hasArealRoot, const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const { - std::greater<Scalar> greater; + std::greater<RealScalar> greater; return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold ); } @@ -282,7 +282,7 @@ class PolynomialSolverBase bool& hasArealRoot, const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const { - std::less<Scalar> less; + std::less<RealScalar> less; return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold ); } @@ -327,7 +327,7 @@ class PolynomialSolverBase * 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 > +template<typename _Scalar, int _Deg> class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> { public: @@ -337,7 +337,10 @@ class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) typedef Matrix<Scalar,_Deg,_Deg> CompanionMatrixType; - typedef EigenSolver<CompanionMatrixType> EigenSolverType; + typedef typename internal::conditional<NumTraits<Scalar>::IsComplex, + ComplexEigenSolver<CompanionMatrixType>, + EigenSolver<CompanionMatrixType> >::type EigenSolverType; + typedef typename internal::conditional<NumTraits<Scalar>::IsComplex, Scalar, std::complex<Scalar> >::type ComplexScalar; public: /** Computes the complex roots of a new polynomial. */ @@ -352,6 +355,25 @@ class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> companion.balance(); m_eigenSolver.compute( companion.denseMatrix() ); m_roots = m_eigenSolver.eigenvalues(); + // cleanup noise in imaginary part of real roots: + // if the imaginary part is rather small compared to the real part + // and that cancelling the imaginary part yield a smaller evaluation, + // then it's safe to keep the real part only. + RealScalar coarse_prec = RealScalar(std::pow(4,poly.size()+1))*NumTraits<RealScalar>::epsilon(); + for(Index i = 0; i<m_roots.size(); ++i) + { + if( internal::isMuchSmallerThan(numext::abs(numext::imag(m_roots[i])), + numext::abs(numext::real(m_roots[i])), + coarse_prec) ) + { + ComplexScalar as_real_root = ComplexScalar(numext::real(m_roots[i])); + if( numext::abs(poly_eval(poly, as_real_root)) + <= numext::abs(poly_eval(poly, m_roots[i]))) + { + m_roots[i] = as_real_root; + } + } + } } else if(poly.size () == 2) { diff --git a/unsupported/Eigen/src/Polynomials/PolynomialUtils.h b/unsupported/Eigen/src/Polynomials/PolynomialUtils.h index 40ba65b7e..394e857ac 100644 --- a/unsupported/Eigen/src/Polynomials/PolynomialUtils.h +++ b/unsupported/Eigen/src/Polynomials/PolynomialUtils.h @@ -20,8 +20,8 @@ namespace Eigen { * 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> + * \note for stability: + * \f$ |x| \le 1 \f$ */ template <typename Polynomials, typename T> inline @@ -67,8 +67,8 @@ T poly_eval( const Polynomials& poly, const T& x ) * 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> + * \pre + * the leading coefficient of the input polynomial poly must be non zero */ template <typename Polynomial> inline |