1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2010 Manuel Yguel <manuel.yguel (at) gmail.com> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #include "main.h" 11 #include <unsupported/Eigen/Polynomials> 12 #include <iostream> 13 #include <algorithm> 14 15 using namespace std; 16 17 namespace Eigen { 18 namespace internal { 19 template<int Size> 20 struct increment_if_fixed_size 21 { 22 enum { 23 ret = (Size == Dynamic) ? Dynamic : Size+1 24 }; 25 }; 26 } 27 } 28 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 psolve.compute( pols ); 42 const RootsType& roots( psolve.roots() ); 43 EvalRootsType evr( deg ); 44 for( int i=0; i<roots.size(); ++i ){ 45 evr[i] = std::abs( poly_eval( pols, roots[i] ) ); } 46 47 bool evalToZero = evr.isZero( test_precision<Scalar>() ); 48 if( !evalToZero ) 49 { 50 cerr << "WRONG root: " << endl; 51 cerr << "Polynomial: " << pols.transpose() << endl; 52 cerr << "Roots found: " << roots.transpose() << endl; 53 cerr << "Abs value of the polynomial at the roots: " << evr.transpose() << endl; 54 cerr << endl; 55 } 56 57 std::vector<Scalar> rootModuli( roots.size() ); 58 Map< EvalRootsType > aux( &rootModuli[0], roots.size() ); 59 aux = roots.array().abs(); 60 std::sort( rootModuli.begin(), rootModuli.end() ); 61 bool distinctModuli=true; 62 for( size_t i=1; i<rootModuli.size() && distinctModuli; ++i ) 63 { 64 if( internal::isApprox( rootModuli[i], rootModuli[i-1] ) ){ 65 distinctModuli = false; } 66 } 67 VERIFY( evalToZero || !distinctModuli ); 68 69 return distinctModuli; 70 } 71 72 73 74 75 76 77 78 template<int Deg, typename POLYNOMIAL> 79 void evalSolver( const POLYNOMIAL& pols ) 80 { 81 typedef typename POLYNOMIAL::Scalar Scalar; 82 83 typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType; 84 85 PolynomialSolverType psolve; 86 aux_evalSolver<Deg, POLYNOMIAL, PolynomialSolverType>( pols, psolve ); 87 } 88 89 90 91 92 template< int Deg, typename POLYNOMIAL, typename ROOTS, typename REAL_ROOTS > 93 void evalSolverSugarFunction( const POLYNOMIAL& pols, const ROOTS& roots, const REAL_ROOTS& real_roots ) 94 { 95 using std::sqrt; 96 typedef typename POLYNOMIAL::Scalar Scalar; 97 98 typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType; 99 100 PolynomialSolverType psolve; 101 if( aux_evalSolver<Deg, POLYNOMIAL, PolynomialSolverType>( pols, psolve ) ) 102 { 103 //It is supposed that 104 // 1) the roots found are correct 105 // 2) the roots have distinct moduli 106 107 typedef typename POLYNOMIAL::Scalar Scalar; 108 typedef typename REAL_ROOTS::Scalar Real; 109 typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType; 110 111 //Test realRoots 112 std::vector< Real > calc_realRoots; 113 psolve.realRoots( calc_realRoots ); 114 VERIFY( calc_realRoots.size() == (size_t)real_roots.size() ); 115 116 const Scalar psPrec = sqrt( test_precision<Scalar>() ); 117 118 for( size_t i=0; i<calc_realRoots.size(); ++i ) 119 { 120 bool found = false; 121 for( size_t j=0; j<calc_realRoots.size()&& !found; ++j ) 122 { 123 if( internal::isApprox( calc_realRoots[i], real_roots[j] ), psPrec ){ 124 found = true; } 125 } 126 VERIFY( found ); 127 } 128 129 //Test greatestRoot 130 VERIFY( internal::isApprox( roots.array().abs().maxCoeff(), 131 abs( psolve.greatestRoot() ), psPrec ) ); 132 133 //Test smallestRoot 134 VERIFY( internal::isApprox( roots.array().abs().minCoeff(), 135 abs( psolve.smallestRoot() ), psPrec ) ); 136 137 bool hasRealRoot; 138 //Test absGreatestRealRoot 139 Real r = psolve.absGreatestRealRoot( hasRealRoot ); 140 VERIFY( hasRealRoot == (real_roots.size() > 0 ) ); 141 if( hasRealRoot ){ 142 VERIFY( internal::isApprox( real_roots.array().abs().maxCoeff(), abs(r), psPrec ) ); } 143 144 //Test absSmallestRealRoot 145 r = psolve.absSmallestRealRoot( hasRealRoot ); 146 VERIFY( hasRealRoot == (real_roots.size() > 0 ) ); 147 if( hasRealRoot ){ 148 VERIFY( internal::isApprox( real_roots.array().abs().minCoeff(), abs( r ), psPrec ) ); } 149 150 //Test greatestRealRoot 151 r = psolve.greatestRealRoot( hasRealRoot ); 152 VERIFY( hasRealRoot == (real_roots.size() > 0 ) ); 153 if( hasRealRoot ){ 154 VERIFY( internal::isApprox( real_roots.array().maxCoeff(), r, psPrec ) ); } 155 156 //Test smallestRealRoot 157 r = psolve.smallestRealRoot( hasRealRoot ); 158 VERIFY( hasRealRoot == (real_roots.size() > 0 ) ); 159 if( hasRealRoot ){ 160 VERIFY( internal::isApprox( real_roots.array().minCoeff(), r, psPrec ) ); } 161 } 162 } 163 164 165 template<typename _Scalar, int _Deg> 166 void polynomialsolver(int deg) 167 { 168 typedef internal::increment_if_fixed_size<_Deg> Dim; 169 typedef Matrix<_Scalar,Dim::ret,1> PolynomialType; 170 typedef Matrix<_Scalar,_Deg,1> EvalRootsType; 171 172 cout << "Standard cases" << endl; 173 PolynomialType pols = PolynomialType::Random(deg+1); 174 evalSolver<_Deg,PolynomialType>( pols ); 175 176 cout << "Hard cases" << endl; 177 _Scalar multipleRoot = internal::random<_Scalar>(); 178 EvalRootsType allRoots = EvalRootsType::Constant(deg,multipleRoot); 179 roots_to_monicPolynomial( allRoots, pols ); 180 evalSolver<_Deg,PolynomialType>( pols ); 181 182 cout << "Test sugar" << endl; 183 EvalRootsType realRoots = EvalRootsType::Random(deg); 184 roots_to_monicPolynomial( realRoots, pols ); 185 evalSolverSugarFunction<_Deg>( 186 pols, 187 realRoots.template cast < 188 std::complex< 189 typename NumTraits<_Scalar>::Real 190 > 191 >(), 192 realRoots ); 193 } 194 195 void test_polynomialsolver() 196 { 197 for(int i = 0; i < g_repeat; i++) 198 { 199 CALL_SUBTEST_1( (polynomialsolver<float,1>(1)) ); 200 CALL_SUBTEST_2( (polynomialsolver<double,2>(2)) ); 201 CALL_SUBTEST_3( (polynomialsolver<double,3>(3)) ); 202 CALL_SUBTEST_4( (polynomialsolver<float,4>(4)) ); 203 CALL_SUBTEST_5( (polynomialsolver<double,5>(5)) ); 204 CALL_SUBTEST_6( (polynomialsolver<float,6>(6)) ); 205 CALL_SUBTEST_7( (polynomialsolver<float,7>(7)) ); 206 CALL_SUBTEST_8( (polynomialsolver<double,8>(8)) ); 207 208 CALL_SUBTEST_9( (polynomialsolver<float,Dynamic>( 209 internal::random<int>(9,13) 210 )) ); 211 CALL_SUBTEST_10((polynomialsolver<double,Dynamic>( 212 internal::random<int>(9,13) 213 )) ); 214 } 215 } 216
