11 #include <unsupported/Eigen/Polynomials> 29 template<
typename PolynomialType>
30 PolynomialType
polyder(
const PolynomialType&
p)
33 PolynomialType
res(p.size());
40 template<
int Deg,
typename POLYNOMIAL,
typename SOLVER>
46 typedef typename SOLVER::RootsType RootsType;
49 const Index deg = pols.size()-1;
52 SOLVER solve_constr (pols);
54 psolve.compute( pols );
55 const RootsType& roots( psolve.roots() );
56 EvalRootsType evr( deg );
57 POLYNOMIAL pols_der =
polyder(pols);
58 EvalRootsType der( deg );
59 for(
int i=0;
i<roots.size(); ++
i ){
67 bool evalToZero = (evr.cwiseQuotient(der)).isZero( test_precision<Scalar>() );
70 cerr <<
"WRONG root: " << endl;
71 cerr <<
"Polynomial: " << pols.transpose() << endl;
72 cerr <<
"Roots found: " << roots.transpose() << endl;
73 cerr <<
"Abs value of the polynomial at the roots: " << evr.transpose() << endl;
77 std::vector<RealScalar> rootModuli( roots.size() );
79 aux = roots.array().abs();
80 std::sort( rootModuli.begin(), rootModuli.end() );
81 bool distinctModuli=
true;
82 for(
size_t i=1;
i<rootModuli.size() && distinctModuli; ++
i )
85 distinctModuli =
false; }
87 VERIFY( evalToZero || !distinctModuli );
89 return distinctModuli;
98 template<
int Deg,
typename POLYNOMIAL>
105 PolynomialSolverType psolve;
106 aux_evalSolver<Deg, POLYNOMIAL, PolynomialSolverType>( pols, psolve );
112 template<
int Deg,
typename POLYNOMIAL,
typename ROOTS,
typename REAL_ROOTS >
121 PolynomialSolverType psolve;
122 if( aux_evalSolver<Deg, POLYNOMIAL, PolynomialSolverType>( pols, psolve ) )
129 std::vector< RealScalar > calc_realRoots;
130 psolve.
realRoots( calc_realRoots, test_precision<RealScalar>());
133 const RealScalar psPrec =
sqrt( test_precision<RealScalar>() );
135 for(
size_t i=0;
i<calc_realRoots.size(); ++
i )
138 for(
size_t j=0;
j<calc_realRoots.size()&& !found; ++
j )
148 abs( psolve.greatestRoot() ), psPrec ) );
152 abs( psolve.smallestRoot() ), psPrec ) );
156 RealScalar r = psolve.absGreatestRealRoot( hasRealRoot );
157 VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
162 r = psolve.absSmallestRealRoot( hasRealRoot );
163 VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
168 r = psolve.greatestRealRoot( hasRealRoot );
169 VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
174 r = psolve.smallestRealRoot( hasRealRoot );
175 VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
182 template<
typename _Scalar,
int _Deg>
186 typedef internal::increment_if_fixed_size<_Deg> Dim;
191 cout <<
"Standard cases" << endl;
192 PolynomialType pols = PolynomialType::Random(deg+1);
193 evalSolver<_Deg,PolynomialType>( pols );
195 cout <<
"Hard cases" << endl;
196 _Scalar multipleRoot = internal::random<_Scalar>();
197 EvalRootsType allRoots = EvalRootsType::Constant(deg,multipleRoot);
199 evalSolver<_Deg,PolynomialType>( pols );
201 cout <<
"Test sugar" << endl;
202 RealRootsType realRoots = RealRootsType::Random(deg);
204 evalSolverSugarFunction<_Deg>(
206 realRoots.template cast <std::complex<RealScalar> >().
eval(),
224 internal::random<int>(9,13)
227 internal::random<int>(9,13)
void evalSolver(const POLYNOMIAL &pols)
#define CALL_SUBTEST_12(FUNC)
#define CALL_SUBTEST_9(FUNC)
#define CALL_SUBTEST_6(FUNC)
#define CALL_SUBTEST_4(FUNC)
T poly_eval(const Polynomials &poly, const T &x)
A matrix or vector expression mapping an existing array of data.
#define CALL_SUBTEST_3(FUNC)
#define CALL_SUBTEST_7(FUNC)
#define CALL_SUBTEST_11(FUNC)
Namespace containing all symbols from the Eigen library.
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T &x, const T &y)
void roots_to_monicPolynomial(const RootVector &rv, Polynomial &poly)
EIGEN_DECLARE_TEST(polynomialsolver)
#define CALL_SUBTEST_10(FUNC)
PolynomialType polyder(const PolynomialType &p)
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
#define VERIFY_IS_EQUAL(a, b)
#define CALL_SUBTEST_1(FUNC)
void evalSolverSugarFunction(const POLYNOMIAL &pols, const ROOTS &roots, const REAL_ROOTS &real_roots)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
#define CALL_SUBTEST_8(FUNC)
bool aux_evalSolver(const POLYNOMIAL &pols, SOLVER &psolve)
NumTraits< Scalar >::Real RealScalar
void realRoots(Stl_back_insertion_sequence &bi_seq, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
#define CALL_SUBTEST_5(FUNC)
#define CALL_SUBTEST_2(FUNC)
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Jet< T, N > sqrt(const Jet< T, N > &f)
void polynomialsolver(int deg)
The matrix class, also used for vectors and row-vectors.
EIGEN_DEVICE_FUNC bool isApprox(const Scalar &x, const Scalar &y, const typename NumTraits< Scalar >::Real &precision=NumTraits< Scalar >::dummy_precision())