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@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>
aux_evalSolver(const POLYNOMIAL & pols,SOLVER & psolve)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
77
78
79
80
81 template<int Deg, typename POLYNOMIAL>
evalSolver(const POLYNOMIAL & pols)82 void evalSolver( const POLYNOMIAL& pols )
83 {
84 typedef typename POLYNOMIAL::Scalar Scalar;
85
86 typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType;
87
88 PolynomialSolverType psolve;
89 aux_evalSolver<Deg, POLYNOMIAL, PolynomialSolverType>( pols, psolve );
90 }
91
92
93
94
95 template< int Deg, typename POLYNOMIAL, typename ROOTS, typename REAL_ROOTS >
evalSolverSugarFunction(const POLYNOMIAL & pols,const ROOTS & roots,const REAL_ROOTS & real_roots)96 void evalSolverSugarFunction( const POLYNOMIAL& pols, const ROOTS& roots, const REAL_ROOTS& real_roots )
97 {
98 using std::sqrt;
99 typedef typename POLYNOMIAL::Scalar Scalar;
100
101 typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType;
102
103 PolynomialSolverType psolve;
104 if( aux_evalSolver<Deg, POLYNOMIAL, PolynomialSolverType>( pols, psolve ) )
105 {
106 //It is supposed that
107 // 1) the roots found are correct
108 // 2) the roots have distinct moduli
109
110 typedef typename POLYNOMIAL::Scalar Scalar;
111 typedef typename REAL_ROOTS::Scalar Real;
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 = 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 abs( psolve.greatestRoot() ), psPrec ) );
134
135 //Test smallestRoot
136 VERIFY( internal::isApprox( roots.array().abs().minCoeff(),
137 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(), 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(), 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>
polynomialsolver(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
test_polynomialsolver()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 CALL_SUBTEST_11((polynomialsolver<float,Dynamic>(1)) );
217 }
218 }
219