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 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>
evalSolver(const POLYNOMIAL & pols)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 >
evalSolverSugarFunction(const POLYNOMIAL & pols,const ROOTS & roots,const REAL_ROOTS & 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>
polynomialsolver(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
test_polynomialsolver()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