cart-elc

Source code for CART-ELC
git clone git://git.laack.co/cart-elc.git
Log | Files | Refs | README | LICENSE

polynomialsolver.cpp (7486B)

      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 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 
    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 }