mirror of
https://gitlab.com/libeigen/eigen.git
synced 2024-12-21 07:19:46 +08:00
263 lines
8.2 KiB
C++
263 lines
8.2 KiB
C++
// This file is part of Eigen, a lightweight C++ template library
|
|
// for linear algebra.
|
|
//
|
|
// Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
|
|
//
|
|
// Eigen is free software; you can redistribute it and/or
|
|
// modify it under the terms of the GNU Lesser General Public
|
|
// License as published by the Free Software Foundation; either
|
|
// version 3 of the License, or (at your option) any later version.
|
|
//
|
|
// Alternatively, you can redistribute it and/or
|
|
// modify it under the terms of the GNU General Public License as
|
|
// published by the Free Software Foundation; either version 2 of
|
|
// the License, or (at your option) any later version.
|
|
//
|
|
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
|
|
// GNU General Public License for more details.
|
|
//
|
|
// You should have received a copy of the GNU Lesser General Public
|
|
// License and a copy of the GNU General Public License along with
|
|
// Eigen. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
#include "main.h"
|
|
#include <unsupported/Eigen/Polynomials>
|
|
#include <iostream>
|
|
#include <algorithm>
|
|
|
|
#ifdef HAS_GSL
|
|
#include "gsl_helper.h"
|
|
#endif
|
|
|
|
using namespace std;
|
|
|
|
namespace Eigen {
|
|
namespace internal {
|
|
template<int Size>
|
|
struct increment_if_fixed_size
|
|
{
|
|
enum {
|
|
ret = (Size == Dynamic) ? Dynamic : Size+1
|
|
};
|
|
};
|
|
}
|
|
}
|
|
|
|
|
|
template<int Deg, typename POLYNOMIAL, typename SOLVER>
|
|
bool aux_evalSolver( const POLYNOMIAL& pols, SOLVER& psolve )
|
|
{
|
|
typedef typename POLYNOMIAL::Index Index;
|
|
typedef typename POLYNOMIAL::Scalar Scalar;
|
|
|
|
typedef typename SOLVER::RootsType RootsType;
|
|
typedef Matrix<Scalar,Deg,1> EvalRootsType;
|
|
|
|
const Index deg = pols.size()-1;
|
|
|
|
psolve.compute( pols );
|
|
const RootsType& roots( psolve.roots() );
|
|
EvalRootsType evr( deg );
|
|
for( int i=0; i<roots.size(); ++i ){
|
|
evr[i] = std::abs( poly_eval( pols, roots[i] ) ); }
|
|
|
|
bool evalToZero = evr.isZero( test_precision<Scalar>() );
|
|
if( !evalToZero )
|
|
{
|
|
cerr << "WRONG root: " << endl;
|
|
cerr << "Polynomial: " << pols.transpose() << endl;
|
|
cerr << "Roots found: " << roots.transpose() << endl;
|
|
cerr << "Abs value of the polynomial at the roots: " << evr.transpose() << endl;
|
|
cerr << endl;
|
|
}
|
|
|
|
#ifdef HAS_GSL
|
|
if (internal::is_same< Scalar, double>::value)
|
|
{
|
|
typedef GslTraits<Scalar> Gsl;
|
|
RootsType gslRoots(deg);
|
|
Gsl::eigen_poly_solve( pols, gslRoots );
|
|
EvalRootsType gslEvr( deg );
|
|
for( int i=0; i<gslRoots.size(); ++i )
|
|
{
|
|
gslEvr[i] = std::abs( poly_eval( pols, gslRoots[i] ) );
|
|
}
|
|
bool gslEvalToZero = gslEvr.isZero( test_precision<Scalar>() );
|
|
if( !evalToZero )
|
|
{
|
|
if( !gslEvalToZero ){
|
|
cerr << "GSL also failed" << endl; }
|
|
else{
|
|
cerr << "GSL did NOT failed" << endl; }
|
|
cerr << "GSL roots found: " << gslRoots.transpose() << endl;
|
|
cerr << "Abs value of the polynomial at the GSL roots: " << gslEvr.transpose() << endl;
|
|
cerr << endl;
|
|
}
|
|
}
|
|
#endif //< HAS_GSL
|
|
|
|
|
|
std::vector<Scalar> rootModuli( roots.size() );
|
|
Map< EvalRootsType > aux( &rootModuli[0], roots.size() );
|
|
aux = roots.array().abs();
|
|
std::sort( rootModuli.begin(), rootModuli.end() );
|
|
bool distinctModuli=true;
|
|
for( size_t i=1; i<rootModuli.size() && distinctModuli; ++i )
|
|
{
|
|
if( internal::isApprox( rootModuli[i], rootModuli[i-1] ) ){
|
|
distinctModuli = false; }
|
|
}
|
|
VERIFY( evalToZero || !distinctModuli );
|
|
|
|
return distinctModuli;
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
template<int Deg, typename POLYNOMIAL>
|
|
void evalSolver( const POLYNOMIAL& pols )
|
|
{
|
|
typedef typename POLYNOMIAL::Scalar Scalar;
|
|
|
|
typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType;
|
|
|
|
PolynomialSolverType psolve;
|
|
aux_evalSolver<Deg, POLYNOMIAL, PolynomialSolverType>( pols, psolve );
|
|
}
|
|
|
|
|
|
|
|
|
|
template< int Deg, typename POLYNOMIAL, typename ROOTS, typename REAL_ROOTS >
|
|
void evalSolverSugarFunction( const POLYNOMIAL& pols, const ROOTS& roots, const REAL_ROOTS& real_roots )
|
|
{
|
|
typedef typename POLYNOMIAL::Scalar Scalar;
|
|
|
|
typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType;
|
|
|
|
PolynomialSolverType psolve;
|
|
if( aux_evalSolver<Deg, POLYNOMIAL, PolynomialSolverType>( pols, psolve ) )
|
|
{
|
|
//It is supposed that
|
|
// 1) the roots found are correct
|
|
// 2) the roots have distinct moduli
|
|
|
|
typedef typename POLYNOMIAL::Scalar Scalar;
|
|
typedef typename REAL_ROOTS::Scalar Real;
|
|
|
|
typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType;
|
|
typedef typename PolynomialSolverType::RootsType RootsType;
|
|
typedef Matrix<Scalar,Deg,1> EvalRootsType;
|
|
|
|
//Test realRoots
|
|
std::vector< Real > calc_realRoots;
|
|
psolve.realRoots( calc_realRoots );
|
|
VERIFY( calc_realRoots.size() == (size_t)real_roots.size() );
|
|
|
|
const Scalar psPrec = internal::sqrt( test_precision<Scalar>() );
|
|
|
|
for( size_t i=0; i<calc_realRoots.size(); ++i )
|
|
{
|
|
bool found = false;
|
|
for( size_t j=0; j<calc_realRoots.size()&& !found; ++j )
|
|
{
|
|
if( internal::isApprox( calc_realRoots[i], real_roots[j] ), psPrec ){
|
|
found = true; }
|
|
}
|
|
VERIFY( found );
|
|
}
|
|
|
|
//Test greatestRoot
|
|
VERIFY( internal::isApprox( roots.array().abs().maxCoeff(),
|
|
internal::abs( psolve.greatestRoot() ), psPrec ) );
|
|
|
|
//Test smallestRoot
|
|
VERIFY( internal::isApprox( roots.array().abs().minCoeff(),
|
|
internal::abs( psolve.smallestRoot() ), psPrec ) );
|
|
|
|
bool hasRealRoot;
|
|
//Test absGreatestRealRoot
|
|
Real r = psolve.absGreatestRealRoot( hasRealRoot );
|
|
VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
|
|
if( hasRealRoot ){
|
|
VERIFY( internal::isApprox( real_roots.array().abs().maxCoeff(), internal::abs(r), psPrec ) ); }
|
|
|
|
//Test absSmallestRealRoot
|
|
r = psolve.absSmallestRealRoot( hasRealRoot );
|
|
VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
|
|
if( hasRealRoot ){
|
|
VERIFY( internal::isApprox( real_roots.array().abs().minCoeff(), internal::abs( r ), psPrec ) ); }
|
|
|
|
//Test greatestRealRoot
|
|
r = psolve.greatestRealRoot( hasRealRoot );
|
|
VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
|
|
if( hasRealRoot ){
|
|
VERIFY( internal::isApprox( real_roots.array().maxCoeff(), r, psPrec ) ); }
|
|
|
|
//Test smallestRealRoot
|
|
r = psolve.smallestRealRoot( hasRealRoot );
|
|
VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
|
|
if( hasRealRoot ){
|
|
VERIFY( internal::isApprox( real_roots.array().minCoeff(), r, psPrec ) ); }
|
|
}
|
|
}
|
|
|
|
|
|
template<typename _Scalar, int _Deg>
|
|
void polynomialsolver(int deg)
|
|
{
|
|
typedef internal::increment_if_fixed_size<_Deg> Dim;
|
|
typedef Matrix<_Scalar,Dim::ret,1> PolynomialType;
|
|
typedef Matrix<_Scalar,_Deg,1> EvalRootsType;
|
|
|
|
cout << "Standard cases" << endl;
|
|
PolynomialType pols = PolynomialType::Random(deg+1);
|
|
evalSolver<_Deg,PolynomialType>( pols );
|
|
|
|
cout << "Hard cases" << endl;
|
|
_Scalar multipleRoot = internal::random<_Scalar>();
|
|
EvalRootsType allRoots = EvalRootsType::Constant(deg,multipleRoot);
|
|
roots_to_monicPolynomial( allRoots, pols );
|
|
evalSolver<_Deg,PolynomialType>( pols );
|
|
|
|
cout << "Test sugar" << endl;
|
|
EvalRootsType realRoots = EvalRootsType::Random(deg);
|
|
roots_to_monicPolynomial( realRoots, pols );
|
|
evalSolverSugarFunction<_Deg>(
|
|
pols,
|
|
realRoots.template cast <
|
|
std::complex<
|
|
typename NumTraits<_Scalar>::Real
|
|
>
|
|
>(),
|
|
realRoots );
|
|
}
|
|
|
|
void test_polynomialsolver()
|
|
{
|
|
for(int i = 0; i < g_repeat; i++)
|
|
{
|
|
CALL_SUBTEST_1( (polynomialsolver<float,1>(1)) );
|
|
CALL_SUBTEST_2( (polynomialsolver<double,2>(2)) );
|
|
CALL_SUBTEST_3( (polynomialsolver<double,3>(3)) );
|
|
CALL_SUBTEST_4( (polynomialsolver<float,4>(4)) );
|
|
CALL_SUBTEST_5( (polynomialsolver<double,5>(5)) );
|
|
CALL_SUBTEST_6( (polynomialsolver<float,6>(6)) );
|
|
CALL_SUBTEST_7( (polynomialsolver<float,7>(7)) );
|
|
CALL_SUBTEST_8( (polynomialsolver<double,8>(8)) );
|
|
|
|
CALL_SUBTEST_9( (polynomialsolver<float,Dynamic>(
|
|
internal::random<int>(9,13)
|
|
)) );
|
|
CALL_SUBTEST_10((polynomialsolver<double,Dynamic>(
|
|
internal::random<int>(9,13)
|
|
)) );
|
|
}
|
|
}
|