diff options
Diffstat (limited to 'unsupported/test/polynomialsolver.cpp')
| -rw-r--r-- | unsupported/test/polynomialsolver.cpp | 217 |
1 files changed, 217 insertions, 0 deletions
diff --git a/unsupported/test/polynomialsolver.cpp b/unsupported/test/polynomialsolver.cpp new file mode 100644 index 000000000..fefeaff01 --- /dev/null +++ b/unsupported/test/polynomialsolver.cpp @@ -0,0 +1,217 @@ +// 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/. + +#include "main.h" +#include <unsupported/Eigen/Polynomials> +#include <iostream> +#include <algorithm> + +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; + } + + 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) + )) ); + } +} |