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