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>());
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)