GiNaCRA  0.6.4
RealAlgebraicNumberFactory.cpp
Go to the documentation of this file.
00001 /*
00002  * GiNaCRA - GiNaC Real Algebra package
00003  * Copyright (C) 2010-2012  Ulrich Loup, Joachim Redies, Sebastian Junges
00004  *
00005  * This file is part of GiNaCRA.
00006  *
00007  * GiNaCRA is free software: you can redistribute it and/or modify
00008  * it under the terms of the GNU General Public License as published by
00009  * the Free Software Foundation, either version 3 of the License, or
00010  * (at your option) any later version.
00011  *
00012  * GiNaCRA is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU General Public License
00018  * along with GiNaCRA.  If not, see <http://www.gnu.org/licenses/>.
00019  *
00020  */
00021 
00022 
00032 #include <assert.h>
00033 
00034 #include "RealAlgebraicNumberFactory.h"
00035 #include "utilities.h"
00036 #include "operators.h"
00037 
00038 namespace GiNaCRA
00039 {
00040     using std::cout;
00041     using std::endl;
00042 
00043     struct polynomial_has_nonzero_sign:
00044         public unary_function<RationalUnivariatePolynomial, bool>
00045     {
00046         RationalUnivariatePolynomial p;
00047 
00048         polynomial_has_nonzero_sign( const RationalUnivariatePolynomial& p ):
00049             p( p )
00050         {}
00051 
00056         bool operator ()( const RealAlgebraicNumberPtr& a ) const
00057         {
00058             return a->sgn( p ) != GiNaC::ZERO_SIGN;
00059         }
00060     };
00061 
00063     // Common //
00065 
00066     bool RealAlgebraicNumberFactory::equal( const RealAlgebraicNumberPtr& a, const RealAlgebraicNumberPtr& b )
00067     {
00068         RealAlgebraicNumberIRPtr aIR = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberIR>( a );
00069         RealAlgebraicNumberNRPtr aNR = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberNR>( a );
00070         RealAlgebraicNumberIRPtr bIR = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberIR>( b );
00071         RealAlgebraicNumberNRPtr bNR = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberNR>( b );
00072 
00079         if( aIR != 0 )
00080         {
00081             if( bIR != 0 )    // compare two interval representations
00082                 return *aIR == *bIR;
00083             else    // compare one numeric with an interval representation
00084             {
00085                 if( !aIR->refineAvoiding( *bNR ))    // try to refine the isolating interval of irA to avoid nrB
00086                     return false;    // i.e. nrB is not in the isolating interval or not root of the polynomial of irA
00087             }
00088         }
00089         else if( aNR != 0 )
00090         {
00091             if( bNR != 0 )
00092                 return static_cast<numeric>(*aNR) == static_cast<numeric>(*bNR);
00093             else    // compare one numerical with an interval representation
00094             {
00095                 if( !bIR->refineAvoiding( *aNR ))    // try to refine the isolating interval of irB to avoid nrA
00096                     return false;    // i.e. nrA is not in the isolating interval or not root of the polynomial of irB
00097             }
00098         }
00099         return true;    // nrA must be the exact numeric representation of irB OR nrB must be the exact numeric representation of irA
00100     }
00101 
00102     bool RealAlgebraicNumberFactory::less( const RealAlgebraicNumberPtr& a, const RealAlgebraicNumberPtr& b )
00103     {
00104         if( RealAlgebraicNumberFactory::equal( a, b ))
00105             return false;
00106 
00111         RealAlgebraicNumberIRPtr aIR = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberIR>( a );
00112         RealAlgebraicNumberNRPtr aNR = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberNR>( a );
00113         RealAlgebraicNumberIRPtr bIR = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberIR>( b );
00114         RealAlgebraicNumberNRPtr bNR = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberNR>( b );
00115 
00116         if( aIR != 0 )
00117         {
00118             if( bIR != 0 )
00119                 return aIR->isLessWhileUnequal( *bIR );
00120             else
00121                 return aIR->interval().right() <= b->value();
00122         }
00123         else if( aNR != 0 )
00124         {
00125             if( bNR != 0 )
00126                 return a->value() <= b->value();
00127             else
00128                 return a->value() <= bIR->interval().left();
00129         }
00130         return true;
00131     }
00132 
00133     bool RealAlgebraicNumberFactory::isRealAlgebraicNumberNR( const RealAlgebraicNumberPtr& A )
00134     {
00135         RealAlgebraicNumberNRPtr nrA = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberNR>( A );
00136         if( nrA != 0 )
00137             return true;
00138         else
00139             return false;
00140     }
00141 
00142     bool RealAlgebraicNumberFactory::isRealAlgebraicNumberIR( const RealAlgebraicNumberPtr& A )
00143     {
00144         RealAlgebraicNumberIRPtr irA = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberIR>( A );
00145         if( irA != 0 )
00146             return true;
00147         else
00148             return false;
00149     }
00150 
00152     // Real Roots //
00154 
00155     list<RealAlgebraicNumberPtr> RealAlgebraicNumberFactory::realRoots( const RationalUnivariatePolynomial& p,
00156                                                                         RealAlgebraicNumberSettings::IsolationStrategy pivoting )
00157     {
00158         /*
00159          Annotations to the algorithm:
00160          (1) Use Sturm's theorem:
00161          # real roots of p in ]l, r[ = signVariations<RationalUnivariatePolynomial>(standardSturmSequence<RationalUnivariatePolynomial>(p, p'), l) - signVariations<RationalUnivariatePolynomial>(standardSturmSequence<RationalUnivariatePolynomial>(p, p'), r)
00162          (2) All real roots of p are in
00163          (a)    ]-1-p.maximumNorm(), 1+p.maximumNorm()[,
00164          (b)    ]-p.oneNorm(), p.oneNorm()[,
00165          (c)    ]-p.twoNorm(), p.twoNorm()[, or
00166          (d)    ]-p.cauchyBound(), p.cauchyBound()[ (see Lemmas 10.2 and 10.3 in ISBN 0-387-94090-1).
00167          (3) Perform a binary (or more ~> parallel??) search on initial interval by divide & conquer.
00168          */
00169 
00170         list<RealAlgebraicNumberPtr> roots = list<RealAlgebraicNumberPtr>();    // list of p's roots
00171         if( p.isConstant() )
00172             return roots;
00173         list<RationalUnivariatePolynomial> seq = RationalUnivariatePolynomial::standardSturmSequence( p, p.diff() );
00174         // determine two initial intervals as minimal representatives of the above mentioned bounds, excluding 0 (yields normalized intervals in the first place)
00175         numeric l    = -1 - p.maximumNorm();
00176         numeric r    = 1 + p.maximumNorm();
00177         numeric norm = p.cauchyBound();    // twoNorm and oneNorm are upper bounds for the Cauchy bound; thus, we neglect them
00178         if( r > norm )
00179             r = norm;
00180         norm = -norm;
00181         if( l < norm )
00182             l = norm;
00183         // PREPROCESSING:
00184         // check whether 0 is a root and remove the respective monomial
00185         bool zeroRoot = p.hasZeroRoot();
00186         RationalUnivariatePolynomial q = zeroRoot ? RationalUnivariatePolynomial( p.nonzeropart() ) : p;
00187         if( zeroRoot )    // 0 is a root (which is added in the end)
00188             seq = RationalUnivariatePolynomial::standardSturmSequence( q, q.diff() );    // reduce Sturm sequence
00189         // MAIN-SEARCH:
00190         // recursive divide & conquer search of non-zero roots
00191         const unsigned varMinLeft = RationalUnivariatePolynomial::signVariations( seq, l );    // for root order computations
00192         searchRealRoots( varMinLeft, q, seq, OpenInterval( l, 0 ), &roots, 0, pivoting );
00193         searchRealRoots( varMinLeft, q, seq, OpenInterval( 0, r ), &roots, 0, pivoting );
00194         if( zeroRoot )
00195             roots.push_back( RealAlgebraicNumberNRPtr( new RealAlgebraicNumberNR( 0, true )));    // mark as root
00196         return roots;
00197     }
00198 
00199     list<RealAlgebraicNumberPtr> RealAlgebraicNumberFactory::realRootsEval( const UnivariatePolynomial& p,
00200                                                                             const evalmap& m,
00201                                                                             RealAlgebraicNumberSettings::IsolationStrategy pivoting )
00202             throw ( invalid_argument )
00203     {
00204         list<RealAlgebraicNumberPtr> roots = list<RealAlgebraicNumberPtr>();
00205         if( p.isConstant() )
00206             return roots;
00207 
00208         /* Evaluating
00209          */
00210         symbol y = p.variable();    // variable to be solved for
00211         if( m.find( y ) != m.end() )
00212             throw invalid_argument( "The main variable of the polynomial may not occur in the evaluation map." );
00213         evalmap::const_iterator i = m.begin();
00214         ex currentResultant = GiNaC::resultant( i->second->polynomial(), p, i->first );
00215         //        ex currentResultant = GiNaC::resultant( p, i->second->polynomial(), i->first );
00216         //        ex currentResultant = UnivariatePolynomial(p, i->first).resultant(i->second->polynomial());
00217         map<symbol, OpenInterval, GiNaC::ex_is_less> varToInterval;
00218         varToInterval[i->first] = i->second->interval();
00219         // compute the result polynomial and the interval map over all variables
00220         ++i;
00221         for( ; i != m.end(); ++i )
00222         {
00223             currentResultant = GiNaC::resultant( i->second->polynomial(), currentResultant, i->first );
00224             //            currentResultant = UnivariatePolynomial(currentResultant, i->first).resultant(i->second->polynomial());
00225             varToInterval[i->first] = i->second->interval();
00226         }
00227         RationalUnivariatePolynomial res = RationalUnivariatePolynomial( currentResultant, y );
00228         list<RationalUnivariatePolynomial> seq = RationalUnivariatePolynomial::standardSturmSequence( res, res.diff() );
00229 
00230         /* Root-finding PREPROCESSING:
00231          */
00232         // check whether 0 is a root and remove the respective monomial
00233         bool zeroRoot = res.hasZeroRoot();
00234         RationalUnivariatePolynomial q = zeroRoot ? res.nonzeropart() : res;
00235         if( zeroRoot )    // 0 is a root (which is added in the end)
00236             seq = RationalUnivariatePolynomial::standardSturmSequence( q, q.diff() );    // reduce Sturm sequence
00237         // compute the Cauchy bound of p
00238         OpenInterval cauchyBoundInterval = OpenInterval();
00239         OpenInterval lcfInterval         = OpenInterval::evaluate( p.lcoeff(), varToInterval ).abs();    // we have to perform the conversion of coefficients because it is not clear whether we have a numeric or a RealAlgebraicNumberIR
00240         if( !lcfInterval.isZero() )
00241         {
00242             for( int d = p.ldegree(); d <= p.degree(); ++d )
00243                 cauchyBoundInterval = cauchyBoundInterval.add( OpenInterval::evaluate( p.coeff( d ), varToInterval ).abs().div( lcfInterval ));
00244         }
00245         numeric l = -cauchyBoundInterval.right();
00246         numeric r = cauchyBoundInterval.right();
00247         // possibly optimize bounds by maximum norm of resultant
00248         numeric norm = 1 + res.maximumNorm();
00249         if( r > norm )
00250             r = norm;
00251         norm = -norm;
00252         if( l < norm )
00253             l = norm;
00254         //        cout << "realRootsEval: search root of " << res << " in " << l << " and "  << r << endl;
00255         // Root-finding MAIN-SEARCH:
00256         // recursive divide & conquer search of non-zero roots
00257         const unsigned varMinLeft = RationalUnivariatePolynomial::signVariations( seq, l );    // for root order computations
00258         searchRealRoots( varMinLeft, q, seq, OpenInterval( l, 0 ), &roots, 0, pivoting );
00259         searchRealRoots( varMinLeft, q, seq, OpenInterval( 0, r ), &roots, 0, pivoting );
00260         if( zeroRoot )
00261             roots.push_back( RealAlgebraicNumberNRPtr( new RealAlgebraicNumberNR( 0, true )));    // mark as root
00262         return roots;
00263     }
00264 
00265     list<RealAlgebraicNumberPtr> RealAlgebraicNumberFactory::realRootsEval( const UnivariatePolynomial& p,
00266                                                                             const vector<RealAlgebraicNumberIRPtr>& a,
00267                                                                             const vector<symbol>& v,
00268                                                                             RealAlgebraicNumberSettings::IsolationStrategy pivoting )
00269             throw ( invalid_argument )
00270     {
00271         if( a.size() != v.size() )
00272             throw invalid_argument( "The number of specified variables does not match the number of specified numbers." );
00273         evalmap m = evalmap();
00274         for( unsigned i = 0; i != v.size(); ++i )
00275             m[v.at( i )] = a.at( i );
00276         return RealAlgebraicNumberFactory::realRootsEval( p, m );
00277     }
00278 
00279     list<RealAlgebraicNumberPtr> RealAlgebraicNumberFactory::commonRealRoots( const list<RationalUnivariatePolynomial>& l )
00280     {
00281         // heuristics: determine polynomial p with lowest degree (chance of less roots to test)
00282         list<RationalUnivariatePolynomial>::const_iterator q = l.begin();
00283         RationalUnivariatePolynomial p = *q;
00284         ++q;
00285         for( ; q != l.end(); ++q )
00286             if( q->degree() < p.degree() )
00287                 p = *q;
00288         // determine real roots of p and remove any root which is no root of another polynomial in l
00289         list<RealAlgebraicNumberPtr> commonRoots = realRoots( p );
00290         for( q = l.begin(); q != l.end(); ++q )
00291             commonRoots.remove_if( polynomial_has_nonzero_sign( *q ));
00292         return commonRoots;
00293     }
00294 
00296     // Auxiliary Functions //
00298 
00299     void RealAlgebraicNumberFactory::searchRealRoots( const unsigned varMinLeft,
00300                                                       const RationalUnivariatePolynomial& p,
00301                                                       const list<RationalUnivariatePolynomial>& seq,
00302                                                       const OpenInterval& i,
00303                                                       list<RealAlgebraicNumberPtr>* roots,
00304                                                       unsigned offset,
00305                                                       RealAlgebraicNumberSettings::IsolationStrategy pivoting )
00306     {
00307         //    cout << "Search roots of " << p << " in " << i << endl;
00308         // common block
00309         unsigned varRight  = RationalUnivariatePolynomial::signVariations( seq, i.right() );
00310         int      rootCount = RationalUnivariatePolynomial::signVariations( seq, i.left() ) - varRight;
00311         //    cout << "Found " << rootCount << " root(s)!" << endl;
00312 
00317         rootCount -= offset;
00318         if( rootCount <= 0 )    // Tests showed, that in some cases the offset is not needed what results in a negative root count when subtracting the offset anyway.
00319             return;
00320         bool middleIsRoot = false;
00321         // midpoint should be checked for both cases, there is a root inside i (as a heuristics for a numeric representation) or i has to be split by the middle
00322         numeric pivot = i.midpoint();
00323         if( p.sgn( pivot ) == GiNaC::ZERO_SIGN )
00324         {
00325             if( pivoting != RealAlgebraicNumberSettings::SIMPLE_ISOLATIONSTRATEGY )
00326                 roots->push_back( RealAlgebraicNumberNRPtr( new RealAlgebraicNumberNR( pivot, true )));    // mark as root
00327             middleIsRoot = true;
00328         }
00329 
00330         switch( pivoting )
00331         {
00332             case RealAlgebraicNumberSettings::SIMPLE_ISOLATIONSTRATEGY:
00333             {
00334                 if( rootCount == 1 )
00335                 {    // no dissection needed
00336                     roots->push_back( RealAlgebraicNumberIRPtr( new RealAlgebraicNumberIR( p, i, seq, false )));    // prohibit interval normalization
00337                     return;
00338                 }
00339                 if( middleIsRoot )    // in this case, pivot is a root itself what requires a correction in real root counting
00340                     ++offset;
00341                 // split interval into two parts by the pivot element
00342                 unsigned allRootCount = roots->size();
00343                 numeric middleBoundLeft  = i.left();
00344                 numeric middleBoundRight = i.right();
00345                 searchRealRoots( varMinLeft, p, seq, OpenInterval( i.left(), pivot ), roots, offset, pivoting );    // search left
00346                 if( middleIsRoot && allRootCount < roots->size() )
00347                 {    // found roots at the left
00348                     allRootCount                      = roots->size();
00349                     RealAlgebraicNumberIRPtr lastRoot = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberIR>( roots->back() );
00350                     assert( lastRoot != 0 );
00351                     lastRoot->refineAvoiding( pivot );
00352                     middleBoundLeft = lastRoot->interval().right();
00353                 }
00354                 searchRealRoots( varMinLeft, p, seq, OpenInterval( pivot, i.right() ), roots, offset, pivoting );    // search right
00355                 if( middleIsRoot && allRootCount < roots->size() )
00356                 {    // found roots at the right
00357                     RealAlgebraicNumberIRPtr lastRoot = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberIR>( roots->back() );
00358                     assert( lastRoot != 0 );
00359                     lastRoot->refineAvoiding( pivot );
00360                     middleBoundRight = lastRoot->interval().left();
00361                     // add middle
00362                     roots->push_back( RealAlgebraicNumberIRPtr( new RealAlgebraicNumberIR( p, OpenInterval( middleBoundLeft, middleBoundRight ), seq,
00363                                                                                            false )));    // prohibit interval normalization
00364                 }
00365                 return;
00366             }
00367             case RealAlgebraicNumberSettings::GENERIC_ISOLATIONSTRATEGY:
00368                 if( rootCount == 1 )
00369                 {    // no dissection needed
00370                     if( middleIsRoot )
00371                         return;
00372                     roots->push_back( RealAlgebraicNumberIRPtr( new RealAlgebraicNumberIR( p, i, seq, false )));    // prohibit interval normalization
00373                     return;
00374                 }
00375                 if( middleIsRoot )    // in this case, pivot is a root itself what requires a correction in real root counting
00376                     ++offset;
00377                 // split interval into two parts by the pivot element
00378                 searchRealRoots( varMinLeft, p, seq, OpenInterval( i.left(), pivot ), roots, offset, pivoting );    // search left
00379                 searchRealRoots( varMinLeft, p, seq, OpenInterval( pivot, i.right() ), roots, offset, pivoting );    // search right
00380                 return;
00381             case RealAlgebraicNumberSettings::BINARYSAMPLE_ISOLATIONSTRATEGY:
00382                 if( rootCount == 1 )
00383                 {
00384                     if( middleIsRoot )
00385                         return;
00386                     pivot = i.sample();    // try sample as root
00387                     if( p.sgn( pivot ) == GiNaC::ZERO_SIGN )
00388                         roots->push_back( RealAlgebraicNumberNRPtr( new RealAlgebraicNumberNR( pivot, true )));    // mark as root
00389                     else
00390                         roots->push_back( RealAlgebraicNumberIRPtr( new RealAlgebraicNumberIR( p, i, seq, false )));    // prohibit interval normalization
00391                     return;
00392                 }
00393                 pivot = i.sample();    // try sample as separating element
00394                 if( p.sgn( pivot ) == GiNaC::ZERO_SIGN )    // first check it for being a root
00395                 {
00396                     roots->push_back( RealAlgebraicNumberNRPtr( new RealAlgebraicNumberNR( pivot, true )));    // mark as root
00397                     ++offset;    // because pivot is a root itself and will serve as separating element, a correction in real root counting is required
00398                 }
00399                 // split interval into two parts by the pivot element
00400                 searchRealRoots( varMinLeft, p, seq, OpenInterval( i.left(), pivot ), roots, offset, pivoting );    // search left
00401                 searchRealRoots( varMinLeft, p, seq, OpenInterval( pivot, i.right() ), roots, offset, pivoting );    // search right
00402                 return;
00403             case RealAlgebraicNumberSettings::TERNARYSAMPLE_ISOLATIONSTRATEGY:
00404             case RealAlgebraicNumberSettings::TERNARYNEWTON_ISOLATIONSTRATEGY:
00405                 if( rootCount == 1 )
00406                 {    // perform bisection
00407                     if( middleIsRoot )
00408                         return;
00409                     if( pivoting == RealAlgebraicNumberSettings::TERNARYSAMPLE_ISOLATIONSTRATEGY )
00410                         pivot = i.sample();    // try sample as root
00411                     else    // case RealAlgebraicNumberSettings::TERNARYNEWTON_ISOLATIONSTRATEGY
00412                         pivot = p.approximateRealRoot( pivot, i, 1 );    // try Newton's solution as root
00413                     if( p.sgn( pivot ) == GiNaC::ZERO_SIGN )
00414                         roots->push_back( RealAlgebraicNumberNRPtr( new RealAlgebraicNumberNR( pivot, true )));    // mark as root
00415                     else
00416                         roots->push_back( RealAlgebraicNumberIRPtr( new RealAlgebraicNumberIR( p, i, seq, false )));    // prohibit interval normalization
00417                     return;
00418                 }
00419                 numeric pivot2 = pivot;    // init: change nothing
00420                 if( pivoting == RealAlgebraicNumberSettings::TERNARYSAMPLE_ISOLATIONSTRATEGY )
00421                     pivot2 = i.sample();    // try sample as second separating element
00422                 else    // case RealAlgebraicNumberSettings::TERNARYNEWTON_ISOLATIONSTRATEGY
00423                     pivot2 = p.approximateRealRoot( pivot, i, 1 );    // try Newton's solution as second separating element
00424                 if( pivot == pivot2 )
00425                 {
00426                     // split interval into two parts by the pivot element
00427                     searchRealRoots( varMinLeft, p, seq, OpenInterval( i.left(), pivot ), roots, offset, pivoting );    // search left
00428                     searchRealRoots( varMinLeft, p, seq, OpenInterval( pivot, i.right() ), roots, offset, pivoting );    // search right
00429                 }
00430                 else
00431                 {
00432                     if( p.sgn( pivot2 ) == GiNaC::ZERO_SIGN )    // first check it for being a root
00433                     {
00434                         roots->push_back( RealAlgebraicNumberNRPtr( new RealAlgebraicNumberNR( pivot2, true )));    // mark number as root
00435                         ++offset;    // because pivot2 is a root itself and will serve as separating element, a correction in real root counting is required
00436                     }
00437                     if( middleIsRoot )    // in this case, pivot is also a root requiring a correction in real root counting (maximum offset of 2)
00438                         ++offset;
00439                     numeric pivotMin = std::min( pivot, pivot2 );
00440                     numeric pivotMax = std::max( pivot, pivot2 );
00441                     // split interval into three parts by the two pivot elements
00442                     searchRealRoots( varMinLeft, p, seq, OpenInterval( i.left(), pivotMin ), roots, offset, pivoting );
00443                     searchRealRoots( varMinLeft, p, seq, OpenInterval( pivotMin, pivotMax ), roots, offset, pivoting );
00444                     searchRealRoots( varMinLeft, p, seq, OpenInterval( pivotMax, i.right() ), roots, offset, pivoting );
00445                 }
00446                 return;
00447         }
00448     }
00449 
00450     const RealAlgebraicNumberPtr RealAlgebraicNumberFactory::evaluateIR( const UnivariatePolynomial& p,
00451                                                                          const vector<RealAlgebraicNumberIRPtr>& a,
00452                                                                          const vector<symbol>& v )
00453             throw ( invalid_argument )
00454     {
00455         evalmap m = evalmap();
00456         for( unsigned i = 0; i < v.size(); ++i )
00457             m[v.at( i )] = a.at( i );
00458         return RealAlgebraicNumberFactory::evaluateIR( p, m );
00459     }
00460 
00461     const RealAlgebraicNumberPtr RealAlgebraicNumberFactory::evaluateIR( const UnivariatePolynomial& p, const evalmap m ) throw ( invalid_argument )
00462     {
00463         //        cout << "call evalIR( " << p << "( " << p.variable() <<  " ) , [";
00464         //        for( evalmap::const_iterator iter = m.begin(); iter != m.end(); ++iter )
00465         //            cout << "  " << iter->first << " -> " << *iter->second;
00466         //        cout << "] )" << endl;
00467         //        cout << "p rational: " << GiNaC::is_rational_polynomial(p, p.variable()) <<  endl;
00468         //        cout << "evalmap size: " << m.size() <<  endl;
00469         //
00470         //        if( m.size() == 1 && m.begin()->second->sgn( RationalUnivariatePolynomial( p )) == GiNaC::ZERO_SIGN )
00471         //            return RealAlgebraicNumberIRPtr( RealAlgebraicNumberIR::zero( p.variable() ) );
00472         evalmap::const_iterator i = m.begin();
00473         symbol y                              = symbol();    // local auxiliary variable
00474         UnivariatePolynomial currentResultant = UnivariatePolynomial( y - p, i->first ).resultant( i->second->polynomial().inVariable( i->first ));
00475         map<symbol, OpenInterval, GiNaC::ex_is_less> mInterval;
00476         mInterval[i->first] = i->second->interval();
00477         // compute the result polynomial and the initial result interval
00478         ++i;
00479         for( ; i != m.end(); ++i )
00480         {
00481             currentResultant    = UnivariatePolynomial( currentResultant, i->first ).resultant( i->second->polynomial().inVariable( i->first ));
00482             mInterval[i->first] = i->second->interval();
00483         }
00484         //        cout << "current resultant: " << currentResultant << endl;
00485         RationalUnivariatePolynomial r = RationalUnivariatePolynomial( currentResultant, y );    // r in y??
00486         //        cout << "current resultant poly: " << r << endl;
00487         list<RationalUnivariatePolynomial> seq = RationalUnivariatePolynomial::standardSturmSequence( r, r.diff() );
00488         OpenInterval interval = OpenInterval::evaluate( p, mInterval );
00489         // refine the result interval until it isolates exactly one real root of the result polynomial
00490         //        cout << "p = " << r << endl;
00491         while( RationalUnivariatePolynomial::countRealRoots( seq, interval ) != 1 )
00492         {
00493             //                        cout << "i = " << interval << endl;
00494             for( i = m.begin(); i != m.end(); ++i )
00495             {
00496                 i->second->refine();
00497                 mInterval[i->first] = i->second->interval();
00498             }
00499             interval = OpenInterval::evaluate( p, mInterval );
00500         }
00501         //        cout << "evalIR Result: " << std::tr1::shared_ptr<RealAlgebraicNumber>( new RealAlgebraicNumberIR( r, interval )) << endl;
00502         return std::tr1::shared_ptr<RealAlgebraicNumber>( new RealAlgebraicNumberIR( r, interval ));
00503     }
00504 }