GiNaCRA  0.6.4
CAD.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 "settings.h"
00033 #include "utilities.h"
00034 #include "operators.h"
00035 #include "CAD.h"
00036 #include "Groebner.h"
00037 
00038 using std::sort;
00039 using GiNaC::signedSubresultantsCoefficients;
00040 using GiNaC::ZERO_SIGN;
00041 
00042 namespace GiNaCRA
00043 {
00045     // Constructors and Destructors //
00047 
00048     CAD::CAD():
00049         mVariables(),
00050         mSampleTree(),
00051         mSampleListIncrements(),
00052         mEliminationSets(),
00053         mLiftingPositions(),
00054         mIsComplete( false ),
00055         mSetting( CADSettings::getSettings() )
00056     {
00057         // initialize root
00058         mSampleTree.insert( mSampleTree.begin(), RealAlgebraicNumberPtr() );    // empty root node!!
00059     }
00060 
00061     CAD::CAD( const UnivariatePolynomialSet& s, const vector<symbol>& v, CADSettings setting ):
00062         mVariables( v ),
00063         mSampleTree( tree<RealAlgebraicNumberPtr>() ),
00064         mSampleListIncrements( v.size(), SampleList() ),
00065         mEliminationSets( v.size() ),
00066         mLiftingPositions( v.size(), list<unsigned>() ),
00067         mIsComplete( false ),
00068         mSetting( setting )
00069     {
00070         // settings (need to be set before every other process because of VariableListPool)
00071         if( mSetting.simplifyByGroebner() )
00072         {
00073             MultivariatePolynomialSettings::InitializeGiNaCRAMultivariateMR();
00074             for( vector<symbol>::const_iterator i = mVariables.begin(); i != mVariables.end(); ++i )
00075                 VariableListPool::addVariable( *i );
00076         }
00077 
00081         unsigned dim = mVariables.size();
00082         // initialize the elimination sets
00083         // normalization: the first main variable shall be the first variable in mVariables (for elimination order)
00084         UnivariatePolynomialSet currentEliminationSet = UnivariatePolynomialSet();
00085         mEliminationSets[0]                           = vector<UnivariatePolynomial>();    // position 0 corresponds to 0 variables eliminated
00086         for( UnivariatePolynomialSet::const_iterator i = s.cbegin(); i != s.cend(); ++i )
00087         {    // add new polynomials to level 0, unifying their variables
00088             UnivariatePolynomial pNewVar( *i, mVariables.front() );
00089             mEliminationSets[0].push_back( pNewVar );
00090             currentEliminationSet.insert( pNewVar );
00091         }
00092         sort( mEliminationSets[0].begin(), mEliminationSets[0].end(), mSetting.mUP_isLess );
00093         // this loop does nothing if we have univariate polynomials
00094         for( unsigned i = 1; i != dim; ++i )
00095         {    // perform elimination of level i-1
00096             // position i of mEliminationSets corresponds to variable i-1 (current main variable) eliminated in position i-1 of mEliminationSets
00097             currentEliminationSet = eliminationSet( currentEliminationSet, mVariables[i] );
00098             vector<UnivariatePolynomial> currentEliminationList( currentEliminationSet.begin(), currentEliminationSet.end() );
00099             // *** CADSettings: simplifyBySquarefreeing
00100             if( mSetting.simplifyBySquarefreeing() )
00101             {
00102                 for( unsigned i = 0; i < currentEliminationList.size(); ++i )
00103                     currentEliminationList[i] = currentEliminationList[i].sepapart();
00104             }
00105             std::sort( currentEliminationList.begin(), currentEliminationList.end(), mSetting.mUP_isLess );
00106             if( mSetting.simplifyBySquarefreeing() )
00107                 std::unique( currentEliminationList.begin(), currentEliminationList.end() );
00108             // *** /CADSettings: simplifyBySquarefreeing
00109             mEliminationSets[i] = vector<UnivariatePolynomial>( currentEliminationList.begin(), currentEliminationList.end() );
00110         }
00111         if( mEliminationSets.back().empty() )
00112             mIsComplete = true;
00113         // apply heuristics
00114         // *** CADSettings: simplifyByRootcounting
00115         if( mSetting.simplifyByRootcounting() )
00116         {
00117             vector<UnivariatePolynomial> baseLevel = mEliminationSets.back();
00118             for( unsigned i = 0; i < baseLevel.size(); )
00119             {
00120                 if( baseLevel[i].degree() % 2 == 0 )
00121                 {    // there could only be complex roots if the degree is even
00122                     RationalUnivariatePolynomial p( baseLevel[i] );
00123                     if( p.countRealRoots() == 0 )
00124                     {
00125                         baseLevel.erase( baseLevel.begin() + i );
00126                         continue;    // do not increase i
00127                     }
00128                 }
00129                 ++i;
00130             }
00131             mEliminationSets.back() = baseLevel;
00132         }
00133         // *** CADSettings: /simplifyByRootcounting
00134         // initialize root
00135         mSampleTree.insert( mSampleTree.begin(), RealAlgebraicNumberPtr() );    // empty root node!!
00136         // initialize lifting positions
00137         for( unsigned i = 0; i < dim; ++i )
00138         {
00139             list<unsigned> l = list<unsigned>();
00140             for( unsigned j = 0; j < mEliminationSets[i].size(); ++j )
00141                 l.push_back( j );
00142             //            l.sort( isLessInLiftingPositions( mEliminationSets[i], mSetting.mUP_isLess ) );
00143             mLiftingPositions[i] = l;
00144         }
00145     }
00146 
00148     // Operations on CAD object //
00150 
00151     void CAD::complete()
00152     {
00153         RealAlgebraicPoint r = RealAlgebraicPoint();
00154         check( vector<Constraint>( 1, Constraint( Polynomial( 1 ), ZERO_SIGN, mVariables )), r );
00155     }
00156 
00157     inline const vector<RealAlgebraicPoint> CAD::samples()
00158     {
00159         unsigned                               dim  = mVariables.size();
00160         tree<RealAlgebraicNumberPtr>::iterator root = mSampleTree.begin( mSampleTree.head );
00161         vector<RealAlgebraicPoint>             s    = vector<RealAlgebraicPoint>();
00162         for( tree<RealAlgebraicNumberPtr>::leaf_iterator leaf = mSampleTree.begin_leaf(); leaf != mSampleTree.end_leaf(); ++leaf )
00163         {
00164             RealAlgebraicPoint sample = constructSampleAt( leaf, root );    // for each leaf construct the path by iterating back to the root
00165             if( sample.dim() == dim )    // discard points which are ill-formed (possible by intermediate nodes which did not yield valid child nodes)
00166                 s.push_back( sample );
00167         }
00168         return s;
00169     }
00170 
00171     bool CAD::check( const vector<Constraint>& constraints, RealAlgebraicPoint& r )
00172     {
00173         tree<RealAlgebraicNumberPtr>::iterator root       = mSampleTree.begin( mSampleTree.head );
00174         int                                    dim        = mVariables.size();
00175         vector<RealAlgebraicPoint>             sampleList = samples();    // extract all valid samples
00176         // traverse the current sample tree for satisfying samples
00177         for( vector<RealAlgebraicPoint>::const_iterator a = sampleList.begin(); a != sampleList.end(); ++a )
00178         {    // check each sample point
00179             if( satisfys( *a, constraints ))
00180             {
00181                 r = *a;
00182                 return true;
00183             }
00184         }
00185         if( mIsComplete )    // there are no more samples producible
00186             return false;
00187         // collect all constraints which do have a polynomial in common with the polynomials of this CAD
00188 //        for( vector<UnivariatePolynomial>::const_iterator pol = mEliminationSets[0].begin(); pol != mEliminationSets[0].end(); ++pol )
00189 //
00190         // if still not complete, construct new samples starting at the base level mVariables.size( )-1
00191         return liftCheck( root, list<RealAlgebraicNumberPtr>(), dim, list<symbol>(), constraints, r );
00192     }
00193 
00194     void CAD::printSampleTree( std::ostream& os )
00195     {
00196         for( tree<RealAlgebraicNumberPtr>::iterator i = mSampleTree.begin(); i != mSampleTree.end(); ++i )
00197         {
00198             for( int d = 0; d != mSampleTree.depth( i ); ++d )
00199                 os << "  [";
00200             print( *i, os );
00201             os << endl;
00202         }
00203     }
00204 
00206     // PUBLIC STATIC METHODS //
00208 
00209     const UnivariatePolynomialSet CAD::eliminationSet( const UnivariatePolynomialSet& polynomials,
00210                                                        const symbol& nextVariable )
00211             throw ( invalid_argument )
00212     {
00213         if( polynomials.empty() )
00214             return polynomials;
00215         UnivariatePolynomialSet eliminatedPolynomials = UnivariatePolynomialSet();
00216         symbol x                                      = polynomials.variable();
00217         assert( (ex)nextVariable != x );    // Next variable must not equal the main variable of the set!
00218         // !PAIRED:
00219         for( UnivariatePolynomialSet::const_iterator it1 = polynomials.begin(); it1 != polynomials.end(); ++it1 )
00220             elimination( *it1, nextVariable, eliminatedPolynomials );
00221         // PAIRED:
00222         for( UnivariatePolynomialSet::const_iterator pol_it1 = polynomials.begin(); pol_it1 != polynomials.end(); ++pol_it1 )
00223         {
00224             UnivariatePolynomialSet::const_iterator it2 = pol_it1;
00225             ++it2;    // guarantees pol_it1 != it2
00226             for( ; it2 != polynomials.end(); ++it2 )
00227                 elimination( *pol_it1, *it2, nextVariable, eliminatedPolynomials );
00228         }
00229         eliminatedPolynomials = eliminatedPolynomials.makePrimitive();
00230         eliminatedPolynomials.removeNumbers();
00231         return eliminatedPolynomials;
00232     }
00233 
00234     const SampleList CAD::samples( const RationalUnivariatePolynomial& p,
00235                                    SampleList& currentSamples,
00236                                    CADSettings settings )
00237             throw ( invalid_argument )
00238     {
00239         return CAD::samples( RealAlgebraicNumberFactory::realRoots( p, settings.mIsolationStrategy ), currentSamples );
00240     }
00241 
00242     const SampleList CAD::samples( const UnivariatePolynomial& p,
00243                                    const list<RealAlgebraicNumberPtr>& sample,
00244                                    const list<symbol>& variables,
00245                                    SampleList& currentSamples,
00246                                    CADSettings settings )
00247             throw ( invalid_argument )
00248     {
00249         assert( variables.size() == sample.size() );
00250 #ifdef GINACRA_CAD_DEBUG
00251         cout << "samples of " << p << " in ";
00252         for( auto i = variables.begin(); i!= variables.end(); ++i )
00253             cout << " " << *i;
00254         cout << " eval. at ";
00255         for( auto i = sample.begin(); i!= sample.end(); ++i )
00256             cout << " " << *i;
00257         cout << endl;
00258 #endif
00259         vector<RealAlgebraicNumberIRPtr> numbersIR = vector<RealAlgebraicNumberIRPtr>( variables.size() );    // shall contain only interval-represented components after the preprocessing
00260         vector<symbol> variablesIR = vector<symbol>( variables.size() );    // shall contain the variable indices corresponding to the components of rInterval
00261         int j = 0;
00262         ex pEx = p;
00263         list<RealAlgebraicNumberPtr>::const_iterator sampleValue = sample.begin();
00264         list<symbol>::const_iterator                 variable    = variables.begin();
00265         // Preprocessing: substitute all NumericRepresentation occurrences of r directly
00266         while( variable != variables.end() )
00267         {
00268             RealAlgebraicNumberNRPtr rNumeric = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberNR>( *sampleValue );
00269             if( rNumeric == 0 )
00270             {    // store interval representation for later substitution
00271                 numbersIR[j]     = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberIR>( *sampleValue );    // cast safe here
00272                 variablesIR[j++] = *variable;
00273             }
00274             else
00275                 pEx = pEx.subs( *variable == static_cast<numeric>(*rNumeric) );
00276             ++sampleValue;
00277             ++variable;
00278         }
00279         try
00280         {    // case: all components of r were rational numeric
00281 #ifdef GINACRA_CAD_DEBUG
00282             cout << "Use roots of " << pEx << endl;
00283 #endif
00284             return samples( RationalUnivariatePolynomial( pEx, p.variable() ), currentSamples );
00285         }
00286         catch( ... )
00287         {    // case: there are interval-represented components
00288             numbersIR.resize( j );
00289             variablesIR.resize( j );
00290             // go through the remaining interval components recursively
00291             --j;
00292 #ifdef GINACRA_CAD_DEBUG
00293             cout << "Roots of " << pEx << " in ";
00294             for( auto k = numbersIR.begin(); k != numbersIR.end(); ++k )
00295                 cout << " " << **k;
00296             cout << ": " << endl;
00297 #endif
00298             list<RealAlgebraicNumberPtr> roots = RealAlgebraicNumberFactory::realRootsEval( UnivariatePolynomial( pEx, p.variable() ), numbersIR, variablesIR,
00299                                                                        settings.mIsolationStrategy );
00300 #ifdef GINACRA_CAD_DEBUG
00301             for( auto k = roots.begin(); k != roots.end(); ++k )
00302                 cout << " " << *k;
00303             cout << "." << endl;
00304 #endif
00305             return samples( roots, currentSamples );
00306         }
00307     }
00308 
00309     // ATOMIC METHODS //
00311 
00312     void CAD::elimination( const UnivariatePolynomial& p, const symbol& variable, UnivariatePolynomialSet& eliminated ) throw ( invalid_argument )
00313     {
00314         UnivariatePolynomialSet truncations = CAD::truncation( p );
00315         for( UnivariatePolynomialSet::const_iterator it2 = truncations.begin(); it2 != truncations.end(); ++it2 )
00316         {
00317             eliminated.insert( UnivariatePolynomial( it2->lcoeff(), variable ));
00318             UnivariatePolynomial it2Diff = it2->diff();
00319             if( !it2Diff.isZero() )
00320             {
00321                 vector<ex> subresultants = UnivariatePolynomial::principalSubresultantCoefficients( *it2, it2Diff );
00322                 for( int i = 0; i <= it2->degree() - 2 || i == 0; ++i )
00323                     for( vector<ex>::const_iterator i = subresultants.begin(); i != subresultants.end(); ++i )
00324                         eliminated.insert( UnivariatePolynomial( *i, variable ));
00325             }
00326         }
00327     }
00328 
00329     void CAD::elimination( const UnivariatePolynomial& p,
00330                            const UnivariatePolynomial& q,
00331                            const symbol& variable,
00332                            UnivariatePolynomialSet& eliminated )
00333                            throw ( invalid_argument )
00334     {
00335         UnivariatePolynomialSet truncations = CAD::truncation( p );
00336         for( UnivariatePolynomialSet::const_iterator it1 = truncations.begin(); it1 != truncations.end(); ++it1 )
00337         {
00338             vector<ex> subresultants = UnivariatePolynomial::principalSubresultantCoefficients( *it1, q );
00339             for( vector<ex>::const_iterator i = subresultants.begin(); i != subresultants.end(); ++i )
00340                 eliminated.insert( UnivariatePolynomial( *i, variable ));
00341         }
00342     }
00343 
00344     const SampleList CAD::samples( const list<RationalUnivariatePolynomial>& polynomials, SampleList& currentSamples ) throw ( invalid_argument )
00345     {
00346         list<RealAlgebraicNumberPtr> roots = list<RealAlgebraicNumberPtr>();
00347         for( list<RationalUnivariatePolynomial>::const_iterator it = polynomials.begin(); it != polynomials.end(); ++it )
00348         {
00349             list<RealAlgebraicNumberPtr> tmp = RealAlgebraicNumberFactory::realRoots( *it );
00350             roots.insert( roots.end(), tmp.begin(), tmp.end() );
00351         }
00352         return CAD::samples( roots, currentSamples );
00353     }
00354 
00356     // Auxiliary methods //
00358 
00359     inline const RealAlgebraicPoint CAD::constructSampleAt( tree<RealAlgebraicNumberPtr>::iterator node,
00360                                                             const tree<RealAlgebraicNumberPtr>::iterator& root ) const
00361     {
00362         list<RealAlgebraicNumberPtr> v = list<RealAlgebraicNumberPtr>();
00363         while( node != root )
00364         {    // proceed from the deepest node up to the root while the children of root represent the last component of the sample point and the deepest node the first
00365             v.push_back( *node );
00366             node = mSampleTree.parent( node );
00367         }
00368         return RealAlgebraicPoint( v );
00369     }
00370 
00371     inline const bool CAD::liftCheck( tree<RealAlgebraicNumberPtr>::iterator node,
00372                                       const list<RealAlgebraicNumberPtr>& sample,
00373                                       unsigned level,
00374                                       const list<symbol>& variables,
00375                                       const vector<Constraint>& constraints,
00376                                       RealAlgebraicPoint& r )
00377             throw ( invalid_argument )
00378     {
00379         // base level: zero variables left to substitute => evaluate the constraint
00380         if( level == 0 )
00381         {
00382             RealAlgebraicPoint t = RealAlgebraicPoint( sample );
00383             if( satisfys( t, constraints ))
00384             {
00385                 r = t;
00386                 return true;
00387             }
00388             return false;
00389         }
00390         // level > 0: lifting
00391         --level;    // previous variable will be substituted next
00392         list<RealAlgebraicNumberPtr> extSample    = list<RealAlgebraicNumberPtr>( sample );    // initial sample point
00393         list<symbol>                 newVariables = list<symbol>( variables );    // initial variable list
00394         newVariables.push_front( mVariables[level] );    // the first variable is always the last one lifted
00395 
00396         /*
00397          * Main loop: performs all operations possible in one level > 0, in particular, 2 phases.
00398          * Phase 1: Choose a lifting position and construct the corresponding samples.
00399          * Phase 2: Choose a sample and trigger liftCheck for the next level with the chosen sample.
00400          */
00401         bool computeMoreSamples = true;    
00402         SampleList currentSamples = SampleList();  
00403         currentSamples.insert(mSampleTree.begin( node ), mSampleTree.end( node ));
00404         while( true )
00405         {
00406             /* Phase 1
00407              * Lifting position choice and corresponding sample construction.
00408              */
00409             unsigned liftingPosition;
00410             // loop if no samples are present at all or heuristics according to the respective setting demand to continue with a new lifting position, construct new samples
00411             while( computeMoreSamples || mSampleListIncrements[level].empty()
00412                     || (mSetting.preferNRSamples() && mSampleListIncrements[level].emptyNR())
00413                     || (mSetting.preferSamplesByIsRoot() && mSetting.preferNonrootSamples() && mSampleListIncrements[level].emptyNonroot())
00414                     || (mSetting.preferSamplesByIsRoot() &&!mSetting.preferNonrootSamples() && mSampleListIncrements[level].emptyRoot()))
00415             {
00416                 if( computeMoreSamples )    // disable blind sample construction
00417                     computeMoreSamples = false;
00418                 if( !mLiftingPositions[level].empty() )    // choose next lifting position, if possible
00419                     liftingPosition = mLiftingPositions[level].front();
00420                 else
00421                     break;
00422 #ifdef GINACRA_CAD_DEBUG
00423                 cout << "New samples for " << mEliminationSets[level][liftingPosition] << " at " << sample << endl << " given {";
00424                 for( auto k = currentSamples.begin(); k != currentSamples.end(); ++k )
00425                     cout << " " << *k << ( (*k)->isRoot() ? "r" : "" ) << ( (*k)->isNumeric() ? "n" : "" );
00426                 cout << " }:" << endl;
00427 #endif
00428                 SampleList sampls = samples( mEliminationSets[level][liftingPosition], sample, variables, currentSamples );
00429 #ifdef GINACRA_CAD_DEBUG
00430                 for( auto k = sampls.begin(); k != sampls.end(); ++k )
00431                     cout << " " << *k << ( (*k)->isRoot() ? "r" : "" ) << ( (*k)->isNumeric() ? "n" : "" ) << endl;
00432                 cout << "Current samples are " << endl;
00433                 for( auto k = currentSamples.begin(); k != currentSamples.end(); ++k )
00434                     cout << " " << *k << ( (*k)->isRoot() ? "r" : "" ) << ( (*k)->isNumeric() ? "n" : "" ) << endl;
00435 #endif
00436                 mSampleListIncrements[level].insert( sampls );
00437                 mLiftingPositions[level].pop_front();    // discard lifting position just used for sample construction
00438                 if( mSetting.preferSamplesByIsRoot() || mSetting.preferNRSamples() )
00439                 {
00440                     pair<SampleSimplification, bool> simplification = mSampleListIncrements[level].simplify();
00441                     if( simplification.second )
00442                     { // simplification took place => replace all
00443                         for( SampleSimplification::const_iterator i = simplification.first.begin(); i != simplification.first.end(); ++i )
00444                         {
00445                             SampleList::iterator simplEntry = std::lower_bound( currentSamples.begin( ), currentSamples.end( ), i->first, RealAlgebraicNumberFactory::less );
00446                             if( simplEntry != currentSamples.end( ) )
00447                                 *simplEntry = i->second;
00448                             // it could happen that the node is not yet in the list of children
00449                         }
00450                     }
00451                 }
00452             }
00453 #ifdef GINACRA_CAD_DEBUG
00454             if( !mLiftingPositions.empty() )
00455                 cout << "Current lifting position in level " << level << ": " << mEliminationSets[level][liftingPosition] << endl;
00456             else
00457                 cout << "Reached final lifting position in level " << level << "." << endl;
00458             cout << "Current sample point in level " << level << ": " << sample << endl;
00459             //            cout << "Current sample list increment in level " << level << " has NRs: " << !mSampleListIncrements[level].emptyNR() << endl;
00460             //            cout << "Current sample list increment in level " << level << " has IRs: " << !mSampleListIncrements[level].emptyIR() << endl;
00461             if( mSampleListIncrements[level].empty() )
00462                 cout << "Current sample list increment in level " << level << " empty.";
00463             else
00464             {
00465                 cout << "Current next() of sample list increment in level " << level << ": ";
00466                 print( mSampleListIncrements[level].next() );
00467             }
00468             cout << endl;
00469             cout << "Sample list increment for lifting position " << level << ": " << endl;
00470             for( SampleList::const_iterator s = mSampleListIncrements[level].begin(); s != mSampleListIncrements[level].end(); ++s )
00471                 print( *s, cout << " " );
00472             cout << endl;
00473 #endif
00474 
00475             /* Phase 2
00476              * Lifting of the current level.
00477              */
00478             while( !mSampleListIncrements[level].empty() )    // iterate through all samples found by the next() method
00479             {
00480                 tree<RealAlgebraicNumberPtr>::iterator newNode;
00481 
00482                 /*
00483                  * Sample choice
00484                  */
00485                 RealAlgebraicNumberPtr newSample;
00486                 if( mSetting.preferNRSamples() )
00487                 {
00488                     if( mSampleListIncrements[level].emptyNR() &&!mLiftingPositions[level].empty() )
00489                     {
00490                         computeMoreSamples = true;
00491                         break;    // construct new samples until there are lifting positions
00492                     }
00493                     // otherwise take also IR samples (as implemented in nextNR
00494                     newSample = mSampleListIncrements[level].nextNR();
00495                 }
00496                 else if( mSetting.preferSamplesByIsRoot() )
00497                 {
00498                     if( mSetting.preferNonrootSamples() )
00499                     {
00500                         if( mSampleListIncrements[level].emptyNonroot() &&!mLiftingPositions[level].empty() )
00501                         {
00502                             computeMoreSamples = true;
00503                             break;    // construct new samples until there are lifting positions
00504                         }
00505                         // otherwise take also IR samples (as implemented in nextNonroot
00506                         newSample = mSampleListIncrements[level].nextNonroot();
00507                     }
00508                     else
00509                     {
00510                         if( mSampleListIncrements[level].emptyRoot() &&!mLiftingPositions[level].empty() )
00511                         {
00512                             computeMoreSamples = true;
00513                             break;    // construct new samples until there are lifting positions
00514                         }
00515                         // otherwise take also IR samples (as implemented in nextNonroot
00516                         newSample = mSampleListIncrements[level].nextRoot();
00517                     }
00518                 }
00519                 else    // use FCFS
00520                     newSample = mSampleListIncrements[level].next();
00521 
00522                 /*
00523                  * Sample storage
00524                  */
00525 
00526                 newNode = std::lower_bound( mSampleTree.begin( node ), mSampleTree.end( node ), newSample, RealAlgebraicNumberFactory::less );
00527                 if( newNode == mSampleTree.end( node ) ) // the new sample is either contained in the children nor any child is greater than it
00528                     newNode = mSampleTree.append_child( node, newSample );
00529                 else if( RealAlgebraicNumberFactory::equal( *newNode, newSample ) ) // && newNode != mSampleTree.end( node )
00530                     newNode = mSampleTree.replace( newNode, newSample );    // update existing data in tree (could be helpful in case an interval representation was refined or converted to a numeric)
00531                 else // newNode is a child being greater than newSample
00532                     newNode = mSampleTree.insert( newNode, newSample );
00533 
00534                 extSample.push_front( newSample );    // insert at the first position in order to meet the correct variable order
00535 
00536                 // Lifting
00537 
00538                 if( liftCheck( newNode, extSample, level, newVariables, constraints, r ))
00539                     return true;    // current lifting position remains in mLiftingPositions.front()
00540 
00541                 /*
00542                  * Sample pop if lifting unsuccessful
00543                  */
00544                 if( mSetting.preferNRSamples() )
00545                     mSampleListIncrements[level].popNR();    // remove sample from increment list because it was completely lifted (pop() uses heuristics already used in next())
00546                 else if( mSetting.preferSamplesByIsRoot() )
00547                 {
00548                     if( mSetting.preferNonrootSamples() )
00549                         mSampleListIncrements[level].popNonroot();
00550                     else
00551                         mSampleListIncrements[level].popRoot();
00552                 }
00553                 else
00554                     mSampleListIncrements[level].pop();
00555 
00556                 extSample.pop_front();    // clean sample point component again
00557             }
00558             if( mLiftingPositions[level].empty() )
00559                 break;    // all samples tried and no lifting possible any more
00560         }
00561         // restore lifting positions
00562         list<unsigned> l = list<unsigned>();
00563         for( unsigned j = 0; j < mEliminationSets[level].size(); ++j )
00564             l.push_back( j );
00565         mLiftingPositions[level] = l;
00566 
00567         if( sample.empty() )    // no satisfying samples found in initial level => complete CAD was computed
00568             mIsComplete = true;    // store completeness result
00569         return false;
00570     }
00571 
00572     inline const bool CAD::satisfys( const RealAlgebraicPoint& r, const vector<Constraint>& constraints ) const
00573     {
00574 #ifdef GINACRA_CAD_DEBUG
00575         cout << "Checking whether " << r << " satisfies ";
00576         for( vector<Constraint>::const_iterator i = constraints.begin(); i != constraints.end(); ++i )
00577             cout << *i << " ";
00578         cout << endl;
00579 #endif
00580         for( vector<Constraint>::const_iterator c = constraints.begin(); c != constraints.end(); ++c )
00581         {
00582             if( !c->satisfiedBy( r ))
00583                 return false;
00584         }
00585         return true;
00586     }
00588     // STATIC AUXILIARY METHODS //
00590 
00591     const SampleList CAD::samples( const list<RealAlgebraicNumberPtr>& roots, SampleList& currentSamples ) throw ( invalid_argument )
00592     {
00593         SampleList newSampleSet = SampleList();    // new samples to return - roots and maybe intermediate points
00594         if( roots.empty() )
00595         { // if there is no real root any value results in a positive sign, choose ZERO by default
00596             if( currentSamples.empty() )
00597             {
00598                 RealAlgebraicNumberPtr r = RealAlgebraicNumberPtr( new RealAlgebraicNumberNR( GiNaC::ZERO, false ) );
00599                 newSampleSet.insert( r );
00600                 currentSamples.insert( r );
00601             }
00602 //            else: regardless of where the lifting position is evaluated, it has the same sign, so any point in currentSamples would be appropriate
00603             return newSampleSet;
00604         }
00605         for( list<RealAlgebraicNumberPtr>::const_iterator i = roots.begin(); i != roots.end(); ++i )
00606         {
00607             pair<SampleList::iterator, bool> insertValue = currentSamples.insert( *i );    // now the intervals of possible interval representations are disjoint (by equality checks)
00608             if( !insertValue.second )
00609             {    // value already in the list
00610                 if( !(*insertValue.first)->isRoot() )
00611                 {    // the new root is already contained, but only as sample value => take the root and start sample construction from scratch
00612                     currentSamples.remove( insertValue.first );
00613                     assert( (*i)->isRoot() );
00614                     insertValue = currentSamples.insert( *i );
00615                 }
00616                 else if( !(*insertValue.first)->isNumeric() && (*i)->isNumeric() )
00617                 { // there is already an interval-represented root with the same value present and it can be replaced by a numeric
00618                     currentSamples.remove( insertValue.first );
00619                     insertValue = currentSamples.insert( RealAlgebraicNumberPtr( new RealAlgebraicNumberNR( (*i)->value(), true ) ) );
00620                 }
00621                 else
00622                     continue;
00623             }
00624             SampleList currentSamplesIncrement = SampleList();    // local set storing the elements which shall be added to currentSampleSet and newSampleSet in the end
00625 
00631             //
00632             // next: right neighbor
00633             SampleList::iterator neighbor = SampleList::iterator( insertValue.first );
00634             ++neighbor;    // -> next (safe here, but need to check for end() later)
00635             if( neighbor == currentSamples.end() )    // rightmost position
00636             {    // insert one rightmost sample (by adding 1 or taking the rightmost interval bound)
00637                 RealAlgebraicNumberNRPtr insertValueNR = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberNR>( *insertValue.first );
00638                 if( insertValueNR != 0 )
00639                 {
00640                     currentSamplesIncrement.insert( RealAlgebraicNumberPtr( new RealAlgebraicNumberNR( static_cast<numeric>(*insertValueNR)
00641                                                                                                                            + numeric( 1 ), false )));    // add as a no root
00642                 }
00643                 else    // interval representation
00644                     currentSamplesIncrement.insert( RealAlgebraicNumberPtr(
00645                         new RealAlgebraicNumberNR(
00646                             std::tr1::dynamic_pointer_cast<RealAlgebraicNumberIR>( *insertValue.first )->interval().right(), false )));    // add as a no root
00647             }
00648             else if( (*neighbor)->isRoot() )
00649             {    // sample between neighbor and insertValue.first needed and will be added to newSampleSet
00650                 RealAlgebraicNumberNRPtr insertValueNR = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberNR>( *insertValue.first );
00651                 RealAlgebraicNumberNRPtr neighborNR    = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberNR>( *neighbor );
00652                 if( insertValueNR != 0 )
00653                 {
00654                     if( neighborNR != 0 )
00655                         currentSamplesIncrement.insert( RealAlgebraicNumberPtr(
00656                             new RealAlgebraicNumberNR( OpenInterval( *insertValueNR, *neighborNR ).sampleFast(), false )));    // add as no root
00657                     else
00658                         currentSamplesIncrement.insert( RealAlgebraicNumberPtr(
00659                             new RealAlgebraicNumberNR(
00660                                 std::tr1::dynamic_pointer_cast<RealAlgebraicNumberIR>( *neighbor )->interval().left(), false )));    // add as no root
00661                 }
00662                 else    // interval representation, take right bound of insertValue.first which must be strictly between insertValue.first and neighbor
00663                     currentSamplesIncrement.insert( RealAlgebraicNumberPtr(
00664                         new RealAlgebraicNumberNR(
00665                             std::tr1::dynamic_pointer_cast<RealAlgebraicNumberIR>( *insertValue.first )->interval().right(), false )));    // add as no root
00666             }
00667             //
00668             // previous: left neighbor
00669             neighbor = SampleList::iterator( insertValue.first );
00670             if( neighbor == currentSamples.begin() )    // leftmost position
00671             {    // insert one leftmost sample (by subtracting 1 or taking the leftmost interval bound)
00672                 RealAlgebraicNumberNRPtr insertValueNR = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberNR>( *insertValue.first );
00673                 if( insertValueNR != 0 )
00674                     currentSamplesIncrement.insert( RealAlgebraicNumberPtr( new RealAlgebraicNumberNR( static_cast<numeric>(*insertValueNR)
00675                                                                                                                            - numeric( 1 ), false )));    // add as no root
00676                 else    // interval representation
00677                     currentSamplesIncrement.insert( RealAlgebraicNumberPtr(
00678                         new RealAlgebraicNumberNR(
00679                             std::tr1::dynamic_pointer_cast<RealAlgebraicNumberIR>( *insertValue.first )->interval().left(), false )));    // add as no root
00680             }
00681             else
00682             {
00683                 --neighbor;    // now neighbor is the left bound (can be safely determined now)
00684                 if( (*neighbor)->isRoot() )
00685                 {    // sample between neighbor and insertValue.first needed and will be added to newSampleSet
00686                     RealAlgebraicNumberNRPtr insertValueNR = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberNR>( *insertValue.first );
00687                     RealAlgebraicNumberNRPtr neighborNR    = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberNR>( *neighbor );
00688                     if( insertValueNR != 0 )
00689                     {
00690                         if( neighborNR != 0 )
00691                             currentSamplesIncrement.insert( RealAlgebraicNumberPtr(
00692                                 new RealAlgebraicNumberNR( OpenInterval( *neighborNR, *insertValueNR ).sampleFast(), false )));    // add as no root
00693                         else
00694                             currentSamplesIncrement.insert( RealAlgebraicNumberPtr(
00695                                 new RealAlgebraicNumberNR(
00696                                     std::tr1::dynamic_pointer_cast<RealAlgebraicNumberIR>( *neighbor )->interval().right(), false )));    // add as no root
00697                     }
00698                     else    // interval representation, take left bound of insertValue.first which must be strictly between insertValue.first and neighbor
00699                         currentSamplesIncrement.insert( RealAlgebraicNumberPtr(
00700                             new RealAlgebraicNumberNR(
00701                                 std::tr1::dynamic_pointer_cast<RealAlgebraicNumberIR>( *insertValue.first )->interval().left(), false )));    // add as no root
00702                 }
00703             }
00704             newSampleSet.insert( *insertValue.first );    // add the root to new samples (with root switch on)
00705             newSampleSet.insert( currentSamplesIncrement.begin(), currentSamplesIncrement.end() );
00706             currentSamples.insert( currentSamplesIncrement.begin(), currentSamplesIncrement.end() );
00707         }
00708         return newSampleSet;
00709     }
00710 
00711     const UnivariatePolynomialSet CAD::truncation( const UnivariatePolynomial& p )
00712     {
00713         UnivariatePolynomialSet ret;
00714         ret.insert( p );
00715         UnivariatePolynomial truncation( p );
00716         while( !truncation.isConstant() )
00717         {
00718             truncation = truncation.trunc();
00719             ret.insert( truncation );
00720         }
00721         return ret;
00722     }
00723 
00724     const UnivariatePolynomialSet CAD::truncation( const UnivariatePolynomialSet& P )
00725     {
00726         UnivariatePolynomialSet ret;
00727         for( UnivariatePolynomialSet::const_iterator it1 = P.begin(); it1 != P.end(); ++it1 )
00728         {
00729             UnivariatePolynomialSet M = truncation( *it1 );
00730             for( UnivariatePolynomialSet::const_iterator it2 = M.begin(); it2 != M.end(); ++it2 )
00731                 ret.insert( *it2 );
00732         }
00733         return ret;
00734     }
00735 
00736 }