polynomialsolver.cpp
Go to the documentation of this file.
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>
21 {
22  enum {
23  ret = (Size == Dynamic) ? Dynamic : Size+1
24  };
25 };
26 }
27 }
28 
29 template<typename PolynomialType>
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>
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>
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 >
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>
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 
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 }
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
VERIFY_IS_EQUAL
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:386
Eigen::poly_eval
T poly_eval(const Polynomials &poly, const T &x)
Definition: PolynomialUtils.h:46
ret
DenseIndex ret
Definition: level1_cplx_impl.h:44
Eigen::PolynomialSolver
A polynomial solver.
Definition: PolynomialSolver.h:331
aux_evalSolver
bool aux_evalSolver(const POLYNOMIAL &pols, SOLVER &psolve)
Definition: polynomialsolver.cpp:41
Eigen::internal::isApprox
EIGEN_DEVICE_FUNC bool isApprox(const Scalar &x, const Scalar &y, const typename NumTraits< Scalar >::Real &precision=NumTraits< Scalar >::dummy_precision())
Definition: Eigen/src/Core/MathFunctions.h:1947
res
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Definition: PartialRedux_count.cpp:3
CALL_SUBTEST_11
#define CALL_SUBTEST_11(FUNC)
Definition: split_test_helper.h:64
CALL_SUBTEST_9
#define CALL_SUBTEST_9(FUNC)
Definition: split_test_helper.h:52
CALL_SUBTEST_4
#define CALL_SUBTEST_4(FUNC)
Definition: split_test_helper.h:22
CALL_SUBTEST_3
#define CALL_SUBTEST_3(FUNC)
Definition: split_test_helper.h:16
CALL_SUBTEST_1
#define CALL_SUBTEST_1(FUNC)
Definition: split_test_helper.h:4
evalSolverSugarFunction
void evalSolverSugarFunction(const POLYNOMIAL &pols, const ROOTS &roots, const REAL_ROOTS &real_roots)
Definition: polynomialsolver.cpp:113
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
polyder
PolynomialType polyder(const PolynomialType &p)
Definition: polynomialsolver.cpp:30
Eigen::Dynamic
const int Dynamic
Definition: Constants.h:22
CALL_SUBTEST_10
#define CALL_SUBTEST_10(FUNC)
Definition: split_test_helper.h:58
Eigen::roots_to_monicPolynomial
void roots_to_monicPolynomial(const RootVector &rv, Polynomial &poly)
Definition: PolynomialUtils.h:127
CALL_SUBTEST_5
#define CALL_SUBTEST_5(FUNC)
Definition: split_test_helper.h:28
Eigen::g_repeat
static int g_repeat
Definition: main.h:169
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(polynomialsolver)
Definition: polynomialsolver.cpp:210
CALL_SUBTEST_6
#define CALL_SUBTEST_6(FUNC)
Definition: split_test_helper.h:34
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition: split_test_helper.h:10
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
evalSolver
void evalSolver(const POLYNOMIAL &pols)
Definition: polynomialsolver.cpp:99
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
main.h
std
Definition: BFloat16.h:88
p
float * p
Definition: Tutorial_Map_using.cpp:9
CALL_SUBTEST_12
#define CALL_SUBTEST_12(FUNC)
Definition: split_test_helper.h:70
Eigen::internal::increment_if_fixed_size
Definition: polynomialsolver.cpp:20
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
abs
#define abs(x)
Definition: datatypes.h:17
internal
Definition: BandTriangularSolver.h:13
Eigen::numext::maxi
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T &x, const T &y)
Definition: Eigen/src/Core/MathFunctions.h:1093
CALL_SUBTEST_7
#define CALL_SUBTEST_7(FUNC)
Definition: split_test_helper.h:40
CALL_SUBTEST_8
#define CALL_SUBTEST_8(FUNC)
Definition: split_test_helper.h:46
polynomialsolver
void polynomialsolver(int deg)
Definition: polynomialsolver.cpp:183
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:232
Eigen::PolynomialSolverBase::realRoots
void realRoots(Stl_back_insertion_sequence &bi_seq, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
Definition: PolynomialSolver.h:69
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Eigen::internal::cast
EIGEN_DEVICE_FUNC NewType cast(const OldType &x)
Definition: Eigen/src/Core/MathFunctions.h:460
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
VERIFY
#define VERIFY(a)
Definition: main.h:380
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:34:23