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 typedef typename POLYNOMIAL::Scalar Scalar; 96 97 typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType; 98 99 PolynomialSolverType psolve; 100 if( aux_evalSolver<Deg, POLYNOMIAL, PolynomialSolverType>( pols, psolve ) ) 101 { 102 //It is supposed that 103 // 1) the roots found are correct 104 // 2) the roots have distinct moduli 105 106 typedef typename POLYNOMIAL::Scalar Scalar; 107 typedef typename REAL_ROOTS::Scalar Real; 108 109 typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType; 110 typedef typename PolynomialSolverType::RootsType RootsType; 111 typedef Matrix<Scalar,Deg,1> EvalRootsType; 112 113 //Test realRoots 114 std::vector< Real > calc_realRoots; 115 psolve.realRoots( calc_realRoots ); 116 VERIFY( calc_realRoots.size() == (size_t)real_roots.size() ); 117 118 const Scalar psPrec = internal::sqrt( test_precision<Scalar>() ); 119 120 for( size_t i=0; i<calc_realRoots.size(); ++i ) 121 { 122 bool found = false; 123 for( size_t j=0; j<calc_realRoots.size()&& !found; ++j ) 124 { 125 if( internal::isApprox( calc_realRoots[i], real_roots[j] ), psPrec ){ 126 found = true; } 127 } 128 VERIFY( found ); 129 } 130 131 //Test greatestRoot 132 VERIFY( internal::isApprox( roots.array().abs().maxCoeff(), 133 internal::abs( psolve.greatestRoot() ), psPrec ) ); 134 135 //Test smallestRoot 136 VERIFY( internal::isApprox( roots.array().abs().minCoeff(), 137 internal::abs( psolve.smallestRoot() ), psPrec ) ); 138 139 bool hasRealRoot; 140 //Test absGreatestRealRoot 141 Real r = psolve.absGreatestRealRoot( hasRealRoot ); 142 VERIFY( hasRealRoot == (real_roots.size() > 0 ) ); 143 if( hasRealRoot ){ 144 VERIFY( internal::isApprox( real_roots.array().abs().maxCoeff(), internal::abs(r), psPrec ) ); } 145 146 //Test absSmallestRealRoot 147 r = psolve.absSmallestRealRoot( hasRealRoot ); 148 VERIFY( hasRealRoot == (real_roots.size() > 0 ) ); 149 if( hasRealRoot ){ 150 VERIFY( internal::isApprox( real_roots.array().abs().minCoeff(), internal::abs( r ), psPrec ) ); } 151 152 //Test greatestRealRoot 153 r = psolve.greatestRealRoot( hasRealRoot ); 154 VERIFY( hasRealRoot == (real_roots.size() > 0 ) ); 155 if( hasRealRoot ){ 156 VERIFY( internal::isApprox( real_roots.array().maxCoeff(), r, psPrec ) ); } 157 158 //Test smallestRealRoot 159 r = psolve.smallestRealRoot( hasRealRoot ); 160 VERIFY( hasRealRoot == (real_roots.size() > 0 ) ); 161 if( hasRealRoot ){ 162 VERIFY( internal::isApprox( real_roots.array().minCoeff(), r, psPrec ) ); } 163 } 164 } 165 166 167 template<typename _Scalar, int _Deg> 168 void polynomialsolver(int deg) 169 { 170 typedef internal::increment_if_fixed_size<_Deg> Dim; 171 typedef Matrix<_Scalar,Dim::ret,1> PolynomialType; 172 typedef Matrix<_Scalar,_Deg,1> EvalRootsType; 173 174 cout << "Standard cases" << endl; 175 PolynomialType pols = PolynomialType::Random(deg+1); 176 evalSolver<_Deg,PolynomialType>( pols ); 177 178 cout << "Hard cases" << endl; 179 _Scalar multipleRoot = internal::random<_Scalar>(); 180 EvalRootsType allRoots = EvalRootsType::Constant(deg,multipleRoot); 181 roots_to_monicPolynomial( allRoots, pols ); 182 evalSolver<_Deg,PolynomialType>( pols ); 183 184 cout << "Test sugar" << endl; 185 EvalRootsType realRoots = EvalRootsType::Random(deg); 186 roots_to_monicPolynomial( realRoots, pols ); 187 evalSolverSugarFunction<_Deg>( 188 pols, 189 realRoots.template cast < 190 std::complex< 191 typename NumTraits<_Scalar>::Real 192 > 193 >(), 194 realRoots ); 195 } 196 197 void test_polynomialsolver() 198 { 199 for(int i = 0; i < g_repeat; i++) 200 { 201 CALL_SUBTEST_1( (polynomialsolver<float,1>(1)) ); 202 CALL_SUBTEST_2( (polynomialsolver<double,2>(2)) ); 203 CALL_SUBTEST_3( (polynomialsolver<double,3>(3)) ); 204 CALL_SUBTEST_4( (polynomialsolver<float,4>(4)) ); 205 CALL_SUBTEST_5( (polynomialsolver<double,5>(5)) ); 206 CALL_SUBTEST_6( (polynomialsolver<float,6>(6)) ); 207 CALL_SUBTEST_7( (polynomialsolver<float,7>(7)) ); 208 CALL_SUBTEST_8( (polynomialsolver<double,8>(8)) ); 209 210 CALL_SUBTEST_9( (polynomialsolver<float,Dynamic>( 211 internal::random<int>(9,13) 212 )) ); 213 CALL_SUBTEST_10((polynomialsolver<double,Dynamic>( 214 internal::random<int>(9,13) 215 )) ); 216 } 217 } 218
