1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2010 Manuel Yguel <[email protected]>
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 template<typename PolynomialType>
polyder(const PolynomialType & p)30 PolynomialType polyder(const PolynomialType& p)
31 {
32 typedef typename PolynomialType::Scalar Scalar;
33 PolynomialType res(p.size());
34 for(Index i=1; i<p.size(); ++i)
35 res[i-1] = p[i]*Scalar(i);
36 res[p.size()-1] = 0.;
37 return res;
38 }
39
40 template<int Deg, typename POLYNOMIAL, typename SOLVER>
aux_evalSolver(const POLYNOMIAL & pols,SOLVER & psolve)41 bool aux_evalSolver( const POLYNOMIAL& pols, SOLVER& psolve )
42 {
43 typedef typename POLYNOMIAL::Scalar Scalar;
44 typedef typename POLYNOMIAL::RealScalar RealScalar;
45
46 typedef typename SOLVER::RootsType RootsType;
47 typedef Matrix<RealScalar,Deg,1> EvalRootsType;
48
49 const Index deg = pols.size()-1;
50
51 // Test template constructor from coefficient vector
52 SOLVER solve_constr (pols);
53
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 ){
60 evr[i] = std::abs( poly_eval( pols, roots[i] ) );
61 der[i] = numext::maxi(RealScalar(1.), std::abs( poly_eval( pols_der, roots[i] ) ));
62 }
63
64 // we need to divide by the magnitude of the derivative because
65 // with a high derivative is very small error in the value of the root
66 // yiels a very large error in the polynomial evaluation.
67 bool evalToZero = (evr.cwiseQuotient(der)).isZero( test_precision<Scalar>() );
68 if( !evalToZero )
69 {
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;
74 cerr << endl;
75 }
76
77 std::vector<RealScalar> rootModuli( roots.size() );
78 Map< EvalRootsType > aux( &rootModuli[0], 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 )
83 {
84 if( internal::isApprox( rootModuli[i], rootModuli[i-1] ) ){
85 distinctModuli = false; }
86 }
87 VERIFY( evalToZero || !distinctModuli );
88
89 return distinctModuli;
90 }
91
92
93
94
95
96
97
98 template<int Deg, typename POLYNOMIAL>
evalSolver(const POLYNOMIAL & pols)99 void evalSolver( const POLYNOMIAL& pols )
100 {
101 typedef typename POLYNOMIAL::Scalar Scalar;
102
103 typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType;
104
105 PolynomialSolverType psolve;
106 aux_evalSolver<Deg, POLYNOMIAL, PolynomialSolverType>( pols, psolve );
107 }
108
109
110
111
112 template< int Deg, typename POLYNOMIAL, typename ROOTS, typename REAL_ROOTS >
evalSolverSugarFunction(const POLYNOMIAL & pols,const ROOTS & roots,const REAL_ROOTS & real_roots)113 void evalSolverSugarFunction( const POLYNOMIAL& pols, const ROOTS& roots, const REAL_ROOTS& real_roots )
114 {
115 using std::sqrt;
116 typedef typename POLYNOMIAL::Scalar Scalar;
117 typedef typename POLYNOMIAL::RealScalar RealScalar;
118
119 typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType;
120
121 PolynomialSolverType psolve;
122 if( aux_evalSolver<Deg, POLYNOMIAL, PolynomialSolverType>( pols, psolve ) )
123 {
124 //It is supposed that
125 // 1) the roots found are correct
126 // 2) the roots have distinct moduli
127
128 //Test realRoots
129 std::vector< RealScalar > calc_realRoots;
130 psolve.realRoots( calc_realRoots, test_precision<RealScalar>());
131 VERIFY_IS_EQUAL( calc_realRoots.size() , (size_t)real_roots.size() );
132
133 const RealScalar psPrec = sqrt( test_precision<RealScalar>() );
134
135 for( size_t i=0; i<calc_realRoots.size(); ++i )
136 {
137 bool found = false;
138 for( size_t j=0; j<calc_realRoots.size()&& !found; ++j )
139 {
140 if( internal::isApprox( calc_realRoots[i], real_roots[j], psPrec ) ){
141 found = true; }
142 }
143 VERIFY( found );
144 }
145
146 //Test greatestRoot
147 VERIFY( internal::isApprox( roots.array().abs().maxCoeff(),
148 abs( psolve.greatestRoot() ), psPrec ) );
149
150 //Test smallestRoot
151 VERIFY( internal::isApprox( roots.array().abs().minCoeff(),
152 abs( psolve.smallestRoot() ), psPrec ) );
153
154 bool hasRealRoot;
155 //Test absGreatestRealRoot
156 RealScalar r = psolve.absGreatestRealRoot( hasRealRoot );
157 VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
158 if( hasRealRoot ){
159 VERIFY( internal::isApprox( real_roots.array().abs().maxCoeff(), abs(r), psPrec ) ); }
160
161 //Test absSmallestRealRoot
162 r = psolve.absSmallestRealRoot( hasRealRoot );
163 VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
164 if( hasRealRoot ){
165 VERIFY( internal::isApprox( real_roots.array().abs().minCoeff(), abs( r ), psPrec ) ); }
166
167 //Test greatestRealRoot
168 r = psolve.greatestRealRoot( hasRealRoot );
169 VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
170 if( hasRealRoot ){
171 VERIFY( internal::isApprox( real_roots.array().maxCoeff(), r, psPrec ) ); }
172
173 //Test smallestRealRoot
174 r = psolve.smallestRealRoot( hasRealRoot );
175 VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
176 if( hasRealRoot ){
177 VERIFY( internal::isApprox( real_roots.array().minCoeff(), r, psPrec ) ); }
178 }
179 }
180
181
182 template<typename _Scalar, int _Deg>
polynomialsolver(int deg)183 void polynomialsolver(int deg)
184 {
185 typedef typename NumTraits<_Scalar>::Real RealScalar;
186 typedef internal::increment_if_fixed_size<_Deg> Dim;
187 typedef Matrix<_Scalar,Dim::ret,1> PolynomialType;
188 typedef Matrix<_Scalar,_Deg,1> EvalRootsType;
189 typedef Matrix<RealScalar,_Deg,1> RealRootsType;
190
191 cout << "Standard cases" << endl;
192 PolynomialType pols = PolynomialType::Random(deg+1);
193 evalSolver<_Deg,PolynomialType>( pols );
194
195 cout << "Hard cases" << endl;
196 _Scalar multipleRoot = internal::random<_Scalar>();
197 EvalRootsType allRoots = EvalRootsType::Constant(deg,multipleRoot);
198 roots_to_monicPolynomial( allRoots, pols );
199 evalSolver<_Deg,PolynomialType>( pols );
200
201 cout << "Test sugar" << endl;
202 RealRootsType realRoots = RealRootsType::Random(deg);
203 roots_to_monicPolynomial( realRoots, pols );
204 evalSolverSugarFunction<_Deg>(
205 pols,
206 realRoots.template cast <std::complex<RealScalar> >().eval(),
207 realRoots );
208 }
209
EIGEN_DECLARE_TEST(polynomialsolver)210 EIGEN_DECLARE_TEST(polynomialsolver)
211 {
212 for(int i = 0; i < g_repeat; i++)
213 {
214 CALL_SUBTEST_1( (polynomialsolver<float,1>(1)) );
215 CALL_SUBTEST_2( (polynomialsolver<double,2>(2)) );
216 CALL_SUBTEST_3( (polynomialsolver<double,3>(3)) );
217 CALL_SUBTEST_4( (polynomialsolver<float,4>(4)) );
218 CALL_SUBTEST_5( (polynomialsolver<double,5>(5)) );
219 CALL_SUBTEST_6( (polynomialsolver<float,6>(6)) );
220 CALL_SUBTEST_7( (polynomialsolver<float,7>(7)) );
221 CALL_SUBTEST_8( (polynomialsolver<double,8>(8)) );
222
223 CALL_SUBTEST_9( (polynomialsolver<float,Dynamic>(
224 internal::random<int>(9,13)
225 )) );
226 CALL_SUBTEST_10((polynomialsolver<double,Dynamic>(
227 internal::random<int>(9,13)
228 )) );
229 CALL_SUBTEST_11((polynomialsolver<float,Dynamic>(1)) );
230 CALL_SUBTEST_12((polynomialsolver<std::complex<double>,Dynamic>(internal::random<int>(2,13))) );
231 }
232 }
233