mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-21 07:19:46 +08:00
Code Issues Packages Projects Releases Wiki Activity

Various fixes in polynomial solver and its unit tests:

Browse Source 
- cleanup noise in imaginary part of real roots
 - take into account the magnitude of the derivative to check roots.
 - use <= instead of < at appropriate places
This commit is contained in:
Gael Guennebaud 2018-12-09 22:54:39 +01:00
parent 348bb386d1
commit 450dc97c6b
3 changed files with 46 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
unsupported/Eigen/src/Polynomials/Companion.h
View File

@ -75,8 +75,7 @@ class companion
void setPolynomial( const VectorType& poly )
{
const Index deg = poly.size()-1;
m_monic = Scalar(-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);
}

30
unsupported/Eigen/src/Polynomials/PolynomialSolver.h
View File

@ -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 )
{
@ -340,6 +340,7 @@ class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
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. */
@ -354,6 +355,27 @@ class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
companion.balance();
m_eigenSolver.compute( companion.denseMatrix() );
m_roots = m_eigenSolver.eigenvalues();
MatrixXcd A = companion.denseMatrix();
// 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 = std::pow(4,poly.size()+1)*NumTraits<RealScalar>::epsilon();
std::cout << coarse_prec << "\n";
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)
{

21
unsupported/test/polynomialsolver.cpp
View File

@ -26,6 +26,16 @@ struct increment_if_fixed_size
}
}
template<typename PolynomialType>
PolynomialType polyder(const PolynomialType& p)
{
typedef typename PolynomialType::Scalar Scalar;
PolynomialType res(p.size());
for(Index i=1; i<p.size(); ++i)
res[i-1] = p[i]*Scalar(i);
res[p.size()-1] = 0.;
return res;
}
template<int Deg, typename POLYNOMIAL, typename SOLVER>
bool aux_evalSolver( const POLYNOMIAL& pols, SOLVER& psolve )
@ -44,10 +54,17 @@ bool aux_evalSolver( const POLYNOMIAL& pols, SOLVER& psolve )
psolve.compute( pols );
const RootsType& roots( psolve.roots() );
EvalRootsType evr( deg );
POLYNOMIAL pols_der = polyder(pols);
EvalRootsType der( deg );
for( int i=0; i<roots.size(); ++i ){
evr[i] = std::abs( poly_eval( pols, roots[i] ) ); }
evr[i] = std::abs( poly_eval( pols, roots[i] ) );
der[i] = numext::maxi(RealScalar(1.), std::abs( poly_eval( pols_der, roots[i] ) ));
}
bool evalToZero = evr.isZero( test_precision<Scalar>() );
// we need to divide by the magnitude of the derivative because
// with a high derivative is very small error in the value of the root
// yiels a very large error in the polynomial evaluation.
bool evalToZero = (evr.cwiseQuotient(der)).isZero( test_precision<Scalar>() );
if( !evalToZero )
{
cerr << "WRONG root: " << endl;