| 29 | |
| 30 | template<int Deg, typename POLYNOMIAL, typename SOLVER> |
| 31 | bool aux_evalSolver( const POLYNOMIAL& pols, SOLVER& psolve ) |
| 32 | { |
| 33 | typedef typename POLYNOMIAL::Index Index; |
| 34 | typedef typename POLYNOMIAL::Scalar Scalar; |
| 35 | |
| 36 | typedef typename SOLVER::RootsType RootsType; |
| 37 | typedef Matrix<Scalar,Deg,1> EvalRootsType; |
| 38 | |
| 39 | const Index deg = pols.size()-1; |
| 40 | |
| 41 | // Test template constructor from coefficient vector |
| 42 | SOLVER solve_constr (pols); |
| 43 | |
| 44 | psolve.compute( pols ); |
| 45 | const RootsType& roots( psolve.roots() ); |
| 46 | EvalRootsType evr( deg ); |
| 47 | for( int i=0; i<roots.size(); ++i ){ |
| 48 | evr[i] = std::abs( poly_eval( pols, roots[i] ) ); } |
| 49 | |
| 50 | bool evalToZero = evr.isZero( test_precision<Scalar>() ); |
| 51 | if( !evalToZero ) |
| 52 | { |
| 53 | cerr << "WRONG root: " << endl; |
| 54 | cerr << "Polynomial: " << pols.transpose() << endl; |
| 55 | cerr << "Roots found: " << roots.transpose() << endl; |
| 56 | cerr << "Abs value of the polynomial at the roots: " << evr.transpose() << endl; |
| 57 | cerr << endl; |
| 58 | } |
| 59 | |
| 60 | std::vector<Scalar> rootModuli( roots.size() ); |
| 61 | Map< EvalRootsType > aux( &rootModuli[0], roots.size() ); |
| 62 | aux = roots.array().abs(); |
| 63 | std::sort( rootModuli.begin(), rootModuli.end() ); |
| 64 | bool distinctModuli=true; |
| 65 | for( size_t i=1; i<rootModuli.size() && distinctModuli; ++i ) |
| 66 | { |
| 67 | if( internal::isApprox( rootModuli[i], rootModuli[i-1] ) ){ |
| 68 | distinctModuli = false; } |
| 69 | } |
| 70 | VERIFY( evalToZero || !distinctModuli ); |
| 71 | |
| 72 | return distinctModuli; |
| 73 | } |
| 74 | |
| 75 | |
| 76 | |