CAD.h

Go to the documentation of this file.

/*
 *  GiNaCRA - GiNaC Real Algebra package
 * Copyright (C) 2010-2012  Ulrich Loup, Joachim Redies, Sebastian Junges
 *
 * This file is part of GiNaCRA.
 *
 * GiNaCRA is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * GiNaCRA is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with GiNaCRA.  If not, see <http://www.gnu.org/licenses/>.
 *
 */


#ifndef GINACRA_CAD_H
#define GINACRA_CAD_H

//#define GINACRA_CAD_DEBUG

#include <unordered_map>

#include "tree.h"
#include "settings.h"
#include "Constraint.h"
#include "UnivariatePolynomialSet.h"
#include "RealAlgebraicNumber.h"
#include "RealAlgebraicNumberFactory.h"
#include "RealAlgebraicPoint.h"

namespace GiNaCRA
{
    // TYPES //

    typedef std::unordered_map<RealAlgebraicNumberPtr, RealAlgebraicNumberPtr, RealAlgebraicNumberPtrHasher> SampleSimplification;

    struct SampleList:
        public std::list<RealAlgebraicNumberPtr>
    {
        private:
            list<RealAlgebraicNumberPtr> mQueue;
            pair<list<RealAlgebraicNumberNRPtr>, list<RealAlgebraicNumberIRPtr> > mNRsIRs;
            pair<list<RealAlgebraicNumberPtr>, list<RealAlgebraicNumberPtr> > mNonrootsRoots;

        public:

            pair<iterator, bool> insert( const RealAlgebraicNumberPtr& r )
            {
                RealAlgebraicNumberNRPtr           rNR      = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberNR>( r );
                std::list<RealAlgebraicNumberPtr>::iterator position = this->begin();
                if( !this->empty() )
                {
                    position = std::lower_bound( position, this->end(), r, RealAlgebraicNumberFactory::less );
                    if( position != this->end() && RealAlgebraicNumberFactory::equal( *position, r ))    // already contained in the list
                        return pair<std::list<RealAlgebraicNumberPtr>::iterator, bool>( position, false );    // return iterator to the already contained element
                    // else: append r to the end of the list
                }
                if( rNR != 0 )
                    mNRsIRs.first.push_back( rNR );
                else    // r is numerically represented
                    mNRsIRs.second.push_back( std::tr1::dynamic_pointer_cast<RealAlgebraicNumberIR>( r ));
                if( r->isRoot() )
                    mNonrootsRoots.second.push_back( r );
                else
                    mNonrootsRoots.first.push_back( r );
                mQueue.push_back( r );
                return pair<std::list<RealAlgebraicNumberPtr>::iterator, bool>( std::list<RealAlgebraicNumberPtr>::insert( position, r ), true );    // insert safely and return iterator to the new element
            }

            template<class InputIterator>
            void insert( InputIterator first, InputIterator last )
            {
                for( InputIterator i = first; i != last; ++i )
                    insert( *i );
            }

            void insert( const SampleList& l )
            {
                this->insert( l.begin(), l.end() );
            }

            SampleList::iterator remove( SampleList::iterator position )
            {
                assert( position != this->end() );
                removeFromQueue( *position );
                removeFromNRIR( *position );
                removeFromNonrootRoot( *position );
                return this->erase( position );
            }

            inline const RealAlgebraicNumberPtr next()
            {
                if( this->empty() )
                    assert( false );    // do not call next on empty lists
                return mQueue.front();
            }

            inline const RealAlgebraicNumberPtr nextNR()
            {
                if( this->empty() )
                    assert( false );    // do not call next on empty lists
                if( mNRsIRs.first.empty() )    // only IRs left
                    return mNRsIRs.second.front();
                return mNRsIRs.first.front();    // return NR
            }

            inline const RealAlgebraicNumberPtr nextNonroot()
            {
                if( this->empty() )
                    assert( false );    // do not call next on empty lists
                if( mNonrootsRoots.first.empty() )    // only roots left
                    return mNonrootsRoots.second.front();
                return mNonrootsRoots.first.front();    // return non-root
            }

            inline const RealAlgebraicNumberPtr nextRoot()
            {
                if( this->empty() )
                    assert( false );    // do not call next on empty lists
                if( mNonrootsRoots.second.empty() )    // only non-roots left
                    return mNonrootsRoots.first.front();
                return mNonrootsRoots.second.front();    // return root
            }

            void pop()
            {
#ifdef GINACRA_CAD_DEBUG
            cout << "pop()" << endl;
#endif
                if( this->empty() )
                    return;    // nothing to pop
                RealAlgebraicNumberPtr        r        = this->next();
                list<RealAlgebraicNumberPtr>::iterator position = std::lower_bound( this->begin(), this->end(), r, RealAlgebraicNumberFactory::less );
                if( position != this->end() )    // found in the list
                    this->erase( position );    // remove next()
                else
                    assert( false );    // r should be in this list
                mQueue.pop_front();
                removeFromNRIR(r);
                removeFromNonrootRoot(r);
            }

            void popNR()
            {
#ifdef GINACRA_CAD_DEBUG
            cout << "popNR()" << endl;
#endif
                if( this->empty() )
                    return;    // nothing to pop
                RealAlgebraicNumberPtr        r        = this->nextNR();
                list<RealAlgebraicNumberPtr>::iterator position = std::lower_bound( this->begin(), this->end(), r, RealAlgebraicNumberFactory::less );
                if( position != this->end() )    // found in the list
                    this->erase( position );    // remove nextNR()
                else
                    assert( false );    // r should be in this list
                // remove next also from its bucket
                if( mNRsIRs.first.empty() ) // only IRs left, so pop from them
                    mNRsIRs.second.pop_front();
                else // NRs left, so pop from them
                    mNRsIRs.first.pop_front();
                removeFromNonrootRoot(r);
                removeFromQueue(r);
            }

            void popNonroot()
            {
#ifdef GINACRA_CAD_DEBUG
            cout << "popNonroot()" << endl;
#endif
                if( this->empty() )
                    return;    // nothing to pop
                RealAlgebraicNumberPtr        r        = this->nextNonroot();
                list<RealAlgebraicNumberPtr>::iterator position = std::lower_bound( this->begin(), this->end(), r, RealAlgebraicNumberFactory::less );
                if( position != this->end() )    // found in the list
                    this->erase( position );    // remove nextNonroot()
                else
                    assert( false );    // r should be in this list
                // remove next from its bucket
                if( mNonrootsRoots.first.empty() ) // only roots left
                    mNonrootsRoots.second.pop_front();
                else
                    mNonrootsRoots.first.pop_front();
                removeFromNRIR(r);
                removeFromQueue(r);
            }

            void popRoot()
            {
#ifdef GINACRA_CAD_DEBUG
            cout << "popRoot()" << endl;
#endif
                if( this->empty() )
                    return;    // nothing to pop
                RealAlgebraicNumberPtr        r        = this->nextRoot();
                list<RealAlgebraicNumberPtr>::iterator position = std::lower_bound( this->begin(), this->end(), r, RealAlgebraicNumberFactory::less );
                if( position != this->end() )    // found in the list
                    this->erase( position );    // remove nextNonroot()
                else
                    assert( false );    // r should be in this list
                // remove next from its bucket
                if( mNonrootsRoots.second.empty() ) // only non-roots left
                    mNonrootsRoots.first.pop_front();
                else
                    mNonrootsRoots.second.pop_front();
                removeFromNRIR(r);
                removeFromQueue(r);
            }

            pair< SampleSimplification, bool> simplify()
            {
                pair< SampleSimplification, bool> simplification = pair< SampleSimplification, bool>();
                simplification.second = false;
                for( list<RealAlgebraicNumberIRPtr>::iterator irIter = mNRsIRs.second.begin(); irIter != mNRsIRs.second.end(); )
                {
                    if( !(*irIter)->isNumeric() && (*irIter)->refinementCount() == 0 )    // try at least one refinement
                        (*irIter)->refine();
                    if( (*irIter)->isNumeric() )
                    {
#ifdef GINACRA_CAD_DEBUG
                        cout << "FOUND numeric in ";
                        print( *irIter, cout );
#endif
                        RealAlgebraicNumberPtr r    = std::tr1::dynamic_pointer_cast<RealAlgebraicNumber>( *irIter );
                        RealAlgebraicNumberNRPtr rNR    = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberNR>( *irIter );
                        if( rNR != 0 ) // already simplified!
                            continue;
                        RealAlgebraicNumberNRPtr nr = RealAlgebraicNumberNRPtr( new RealAlgebraicNumberNR( (*irIter)->value(), (*irIter)->isRoot() ));
#ifdef GINACRA_CAD_DEBUG
                        cout << "SIMP: replacing " << **irIter << " by " << *nr << endl;
#endif
                        // store simplification result
                        simplification.first[r] = nr;
                        simplification.second = true;
                        // add to NRs
                        mNRsIRs.first.push_back( nr );
                        // erase in IRs
                        irIter      = mNRsIRs.second.erase( irIter );
                        // replace in basic list
                        list<RealAlgebraicNumberPtr>::iterator position = std::lower_bound( this->begin(), this->end(), r,
                                                                                            RealAlgebraicNumberFactory::less );
                        if( position != this->end() )    // found in the list
                            *position = nr;    // replace ir by nr
                        else
                            assert( false );    // there must be an occurrence in the sample list or there was an error inserting the number
                        // replace in root/non-root lists
                        if( nr->isRoot() )
                        {
                            position = std::find( mNonrootsRoots.second.begin(), mNonrootsRoots.second.end(), r );
                            if( position != mNonrootsRoots.second.end() )    // found in the list
                                *position = nr;    // replace ir by nr
                            else
                                assert( false );    // there must be an occurrence in the sample list or there was an error inserting the number
                        }
                        else
                        {
                            position = std::find( mNonrootsRoots.first.begin(), mNonrootsRoots.first.end(), r );
                            if( position != mQueue.end() )    // found in the list
                                *position = nr;    // replace ir by nr
                            else
                                assert( false );    // there must be an occurrence in the sample list or there was an error inserting the number
                        }
                        // replace in queue
                        position = std::find( mQueue.begin(), mQueue.end(), r );
                        if( position != mQueue.end() )    // found in the list
                            *position = nr;    // replace ir by nr
                        else
                            assert( false );    // there must be an occurrence in the sample list or there was an error inserting the number
                    }
                    else
                        ++irIter;
                }
                return simplification;
            }

            bool contains( RealAlgebraicNumberPtr r ) const
            {
                std::list<RealAlgebraicNumberPtr>::const_iterator position = std::lower_bound( this->begin(), this->end(), r, RealAlgebraicNumberFactory::less );
                return position != this->end();
            }

            bool emptyNR() const
            {
                return mNRsIRs.first.empty();
            }

            bool emptyIR() const
            {
                return mNRsIRs.second.empty();
            }

            bool emptyNonroot() const
            {
                return mNonrootsRoots.first.empty();
            }

            bool emptyRoot() const
            {
                return mNonrootsRoots.second.empty();
            }

    private:

        // AUXILIARY METHODS //

        void removeFromNonrootRoot( RealAlgebraicNumberPtr r )
        {
            if( r->isRoot() )
            {
                list<RealAlgebraicNumberPtr>::iterator pos = std::find( mNonrootsRoots.second.begin(), mNonrootsRoots.second.end(), r );    // find in roots non-root list
                if( pos != mNonrootsRoots.second.end() )
                    mNonrootsRoots.second.erase( pos );
                else
                    assert( false );    // r is marked as root
            }
            else
            {
                list<RealAlgebraicNumberPtr>::iterator pos = std::find( mNonrootsRoots.first.begin(), mNonrootsRoots.first.end(), r );    // find in roots non-root list
                if( pos != mNonrootsRoots.first.end() )
                    mNonrootsRoots.first.erase( pos );
                else
                    assert( false );    // r is marked as root
            }
        }

        void removeFromQueue( RealAlgebraicNumberPtr r )
        {
            // remove from queue
            list<RealAlgebraicNumberPtr>::iterator pos = std::find( mQueue.begin(), mQueue.end(), r );    // find in roots non-root list
            if( pos != mQueue.end() )
                mQueue.erase( pos );
            else
                assert( false );    // r is marked as root
        }

        void removeFromNRIR( RealAlgebraicNumberPtr r )
        {
#ifdef GINACRA_CAD_DEBUG
            cout << "removeFromNRIR: " << r << endl;
            cout << "List: " << endl;
            for( list<RealAlgebraicNumberPtr>::const_iterator i = begin(); i != end(); ++i )
                cout << "  " << *i;
            cout << endl << "NR: " << endl;
            for( list<RealAlgebraicNumberPtr>::const_iterator i = mNRsIRs.first.begin(); i != mNRsIRs.first.end(); ++i )
                cout << "  " << **i;
            cout << endl << "IR: " << endl;
            for( list<RealAlgebraicNumberPtr>::const_iterator i = mNRsIRs.second.begin(); i != mNRsIRs.second.end(); ++i )
                cout << "  " << **i;
#endif
            RealAlgebraicNumberNRPtr rNR = std::tr1::dynamic_pointer_cast<RealAlgebraicNumberNR>( r );
            if( rNR != 0 )
            {
                list<RealAlgebraicNumberNRPtr>::iterator pos = std::find( mNRsIRs.first.begin(), mNRsIRs.first.end(), rNR );
                if( pos != mNRsIRs.first.end() )
                    mNRsIRs.first.erase( pos );
                else
                    assert( false );    // r should be in this list, otherwise it was maybe simplified and moved to the other list
            }
            else
            {
                list<RealAlgebraicNumberIRPtr>::iterator pos = std::find( mNRsIRs.second.begin(), mNRsIRs.second.end(),
                                                                            std::tr1::dynamic_pointer_cast<RealAlgebraicNumberIR>( r ));
                if( pos != mNRsIRs.second.end() )
                    mNRsIRs.second.erase( pos );
                else
                    assert( false );    // r should be in this list
            }
        }

    };

    struct CADSettings
    {
        // ATTRIBUTES //

        bool (*mUP_isLess)( const UnivariatePolynomial&, const UnivariatePolynomial& );
        RealAlgebraicNumberSettings::IsolationStrategy mIsolationStrategy;

        // METHODS //

        static const CADSettings getSettings( unsigned setting = DEFAULT_CADSETTING,
                                              RealAlgebraicNumberSettings::IsolationStrategy isolationStrategy = RealAlgebraicNumberSettings::DEFAULT_ISOLATIONSTRATEGY )
        {
            CADSettings cadSettings        = CADSettings();
            cadSettings.mIsolationStrategy = isolationStrategy;
            if( setting & LOWDEG_CADSETTING )
                cadSettings.mUP_isLess = UnivariatePolynomial::univariatePolynomialIsLessLowDeg;
            if( setting & ODDDEG_CADSETTING )
                cadSettings.mUP_isLess = UnivariatePolynomial::univariatePolynomialIsLessOddDeg;
            if( setting & LOWDEG_CADSETTING & ODDDEG_CADSETTING )
                cadSettings.mUP_isLess = UnivariatePolynomial::univariatePolynomialIsLessOddLowDeg;
            if( setting & EVENDEG_CADSETTING )
                cadSettings.mUP_isLess = UnivariatePolynomial::univariatePolynomialIsLessEvenDeg;
            if( setting & LOWDEG_CADSETTING & EVENDEG_CADSETTING )
                cadSettings.mUP_isLess = UnivariatePolynomial::univariatePolynomialIsLessEvenLowDeg;
            if( setting & EAGERLIFTING_CADSETTING )
                cadSettings.mPreferNRSamples = true;
            if( setting & GROEBNER_CADSETTING )
                cadSettings.mSimplifyByGroebner = true;
            if( setting & REALROOTCOUNT_CADSETTING )
                cadSettings.mSimplifyByRootcounting = true;
            if( setting & SQUAREFREEELIMINATION_CADSETTING )
                cadSettings.mSimplifyBySquarefreeing = true;
            return cadSettings;
        }

        friend std::ostream& operator <<( std::ostream& os, const CADSettings& settings )
        {
            list<string> settingStrs = list<string>();
            if( settings.mSimplifyByGroebner )
                settingStrs.push_back( "Simplify the input polynomials corresponding to equations by a Groebner basis." );
            if( settings.mSimplifyByRootcounting )
                settingStrs.push_back( "Simplify the base elimination level by real root counting." );
            if( settings.mSimplifyBySquarefreeing )
                settingStrs.push_back( "Simplify all elimination levels by replacing the polynomials by their square-free part." );
            if( settings.mPreferNRSamples )
                settingStrs.push_back( "Prefer numerics to interval representations for sample choice." );
            if( settings.mPreferSamplesByIsRoot && settings.mPreferNonrootSamples )
                settingStrs.push_back( "Prefer non-root to root samples for sample choice." );
            if( settings.mPreferSamplesByIsRoot && !settings.mPreferNonrootSamples )
                settingStrs.push_back( "Prefer root to non-root samples for sample choice." );
            os << "+------------------------------------ CAD Setting -----------------------------------";
            if( settingStrs.empty() )
                os << endl << "| Default ";
            else
            {
                for( list<string>::const_iterator i = settingStrs.begin(); i != settingStrs.end(); ++i )
                    os << endl << "↳ " << *i;
            }
            return os << endl << "+------------------------------------------------------------------------------------";
        }

        bool simplifyByGroebner() const
        {
            return mSimplifyByGroebner;
        }

        void setSimplifyByGroebner( bool b )
        {
            mSimplifyByGroebner = b;
        }

        bool simplifyByRootcounting() const
        {
            return mSimplifyByRootcounting;
        }

        void setSimplifyByRootcounting( bool b )
        {
            mSimplifyByRootcounting = b;
        }

        bool simplifyBySquarefreeing() const
        {
            return mSimplifyBySquarefreeing;
        }

        void setSimplifyBySquarefreeing( bool b )
        {
            mSimplifyBySquarefreeing = b;
        }

        bool preferNRSamples() const
        {
            return mPreferNRSamples;
        }

        void setPreferNRSamples( bool b )
        {
            mPreferNRSamples = b;
            // dependency: PreferNRSamples excludes the use of PreferSamplesByIsRoot and vice versa
            mPreferSamplesByIsRoot = b ? false : mPreferSamplesByIsRoot;
        }

        bool preferSamplesByIsRoot() const
        {
            return mPreferSamplesByIsRoot;
        }

        void setPreferSamplesByIsRoot( bool b )
        {
            mPreferSamplesByIsRoot = b;
            // dependency: PreferSamplesByIsRoot excludes the use of PreferNRSamples and vice versa
            mPreferNRSamples = b ? false : mPreferNRSamples;
        }

        bool preferNonrootSamples() const
        {
            return mPreferNonrootSamples;
        }

        void setPreferNonrootSamples( bool b )
        {
            setPreferSamplesByIsRoot( true );
            mPreferNonrootSamples = b;
        }

        protected:

            // ATTRIBUTES //

            bool mPreferNRSamples;
            bool mPreferSamplesByIsRoot;
            bool mPreferNonrootSamples;
            bool mSimplifyByGroebner;
            bool mSimplifyByRootcounting;
            bool mSimplifyBySquarefreeing;

        private:

            CADSettings():
                mUP_isLess( UnivariatePolynomial::univariatePolynomialIsLess ),
                mIsolationStrategy( RealAlgebraicNumberSettings::DEFAULT_ISOLATIONSTRATEGY ),
                mPreferNRSamples( false ),
                mPreferSamplesByIsRoot( false ),
                mPreferNonrootSamples( false ),
                mSimplifyByGroebner( false ),
                mSimplifyByRootcounting( false ),
                mSimplifyBySquarefreeing( false )
            {}

    };

    struct isLessInLiftingPositions
    {
        vector<UnivariatePolynomial> mEliminationSet;
        bool (*mIsLess)( const UnivariatePolynomial&, const UnivariatePolynomial& );

        isLessInLiftingPositions( const vector<UnivariatePolynomial>& eliminationSet, bool (*isLess)( const UnivariatePolynomial&, const UnivariatePolynomial& ) ) :
            mEliminationSet(eliminationSet),
            mIsLess( isLess )
        {}

        bool operator ()( unsigned i, unsigned j ) const
        {
            return mIsLess( mEliminationSet[i], mEliminationSet[j] );
        }
    };

    class CAD
    {
        public:

            // Constructors and Destructors //

            /*
             * Standard constructor doing nothing.
             */
            CAD();

            CAD( const UnivariatePolynomialSet& s, const vector<symbol>& v, CADSettings setting = CADSettings::getSettings() );

            /*
             * Copy constructor.
             */
            CAD( const CAD& cad ):
                mVariables( cad.mVariables ),
                mSampleTree( cad.mSampleTree ),
                mSampleListIncrements( cad.mSampleListIncrements ),
                mEliminationSets( cad.mEliminationSets ),
                mLiftingPositions( cad.mLiftingPositions ),
                mIsComplete( cad.mIsComplete ),
                mSetting( cad.mSetting )
            {}

            // Selectors //

            /*
             * @return list of main variables of the polynomials of this cad
             */
            const vector<symbol> variables()
            {
                return mVariables;
            }

            /*
             * Returns the input set of polynomials underlying this cad as a list.
             * @return the set of polynomials underlying this cad as a list
             */
            const vector<UnivariatePolynomial> polynomials()
            {
                return mEliminationSets[0];
            }

            /*
             * Sets with successively eliminated variables due to a CAD projection. Index i corresponds to the set where i variables were eliminated.
             * <br/ >For i>0, the i-th set was obtained by eliminating the variable i-1.
             * @return all eliminations of the polynomials and the polynomials themselves (at position 0) computed so far
             */
            const vector<vector<UnivariatePolynomial> > eliminationSets()
            {
                return mEliminationSets;
            }

            /*
             * @return true if the cad is computed completely, false if there are still samples to compute
             */
            bool isComplete() const
            {
                return mIsComplete;
            }

            // Operations on CAD object //

            void complete();

            const vector<RealAlgebraicPoint> samples();

            bool check( const vector<Constraint>& constraints, RealAlgebraicPoint& r );

            void printSampleTree( std::ostream& os = std::cout );

            void addPolynomial( const UnivariatePolynomial& p, const vector<GiNaC::symbol>& v, unsigned setting = DEFAULT_CADSETTING )
            {
                list<UnivariatePolynomial> l = list<UnivariatePolynomial>( 1, p );
                addPolynomials( l.begin(), l.end(), v, setting );
            }

            template<class InputIterator>
            void addPolynomials( InputIterator first, InputIterator last, const vector<GiNaC::symbol>& v, unsigned setting = DEFAULT_CADSETTING )
            {
                this->alterSetting( setting );
                // settings (need to be set before every other process because of VariableListPool)
                if( mSetting.simplifyByGroebner() )
                {
                    MultivariatePolynomialSettings::InitializeGiNaCRAMultivariateMR();
                    for( vector<symbol>::const_iterator i = mVariables.begin(); i != mVariables.end(); ++i )
                        VariableListPool::addVariable( *i );
                }

                /* Algorithm overview:
                 * (1) Determine the variables differing from mVariables and add them to the front of the existing variables.
                 * (2) Add as many levels to the front of eliminationSets as new variables were determined.
                 * (3) Compute the elimination for the new elimination sets with the new polynomials, add the prevailing ones, and perform pairwise elimination with the prevailing ones.
                 * (4) Re-compute elimination for the newly arriving polynomials in the old elimination sets.
                 */
                // (1)
                vector<GiNaC::symbol> newVariables     = vector<GiNaC::symbol>();
                unsigned     newVariableCount = 0;
                for( vector<GiNaC::symbol>::const_iterator i = v.begin(); i != v.end(); ++i )
                {
                    if( std::find( mVariables.begin(), mVariables.end(), *i ) == mVariables.end() )
                    {    // found a new variable
                        newVariables.push_back( *i );
                        ++newVariableCount;
                    }
                }
                for( vector<GiNaC::symbol>::const_iterator i = mVariables.begin(); i != mVariables.end(); ++i )
                    newVariables.push_back( *i );
                if( newVariableCount != 0 )
                {    // extend the other data structures depending on level
                    vector<SampleList> newSampleListIncrements = vector<SampleList>( newVariableCount, SampleList() );
                    for( vector<SampleList>::const_iterator i = mSampleListIncrements.begin(); i != mSampleListIncrements.end(); ++i )
                        newSampleListIncrements.push_back( *i );
                    mSampleListIncrements = newSampleListIncrements;
                }
                // (2)
                vector<vector<UnivariatePolynomial> > newEliminationSets =
                    vector<vector<UnivariatePolynomial> >( mEliminationSets.size() + newVariableCount, vector<UnivariatePolynomial>() );
                for( unsigned i = newVariableCount; i < newEliminationSets.size(); ++i )    // copy previous elimination sets
                    newEliminationSets[i] = mEliminationSets[i - newVariableCount];
                UnivariatePolynomialSet currentEliminationSet = UnivariatePolynomialSet( newEliminationSets[0].begin(), newEliminationSets[0].end() );
                UnivariatePolynomialSet s                     = UnivariatePolynomialSet( first, last );    // collect the new polynomials
                for( UnivariatePolynomialSet::const_iterator i = s.begin(); i != s.end(); ++i )
                {    // add new polynomials to level 0, unifying their variables
                    UnivariatePolynomial pNewVar( *i, newVariables.front() );
                    newEliminationSets[0].push_back( pNewVar );
                    currentEliminationSet.insert( pNewVar );
                }
                // (3) and (4) [can be made more efficient by making use of the atomic elimination operators]
                // Caution: The elimination is recomputed completely in case we have multivariate polynomials.
                // this loop does nothing if we have univariate polynomials
                for( unsigned i = 1; i != newVariables.size(); ++i )
                {    // perform elimination of level i-1
                    // position i of mEliminationSets corresponds to variable i-1 (current main variable) eliminated in position i-1 of mEliminationSets
                    currentEliminationSet = eliminationSet( currentEliminationSet, newVariables[i] );
                    for( vector<UnivariatePolynomial>::const_iterator j = newEliminationSets[i].begin(); j != newEliminationSets[i].end();
                            ++j )    // insert possibly existing polynomials of the current level
                        currentEliminationSet.insert( *j );
                    vector<UnivariatePolynomial> currentEliminationList( currentEliminationSet.begin(), currentEliminationSet.end() );
                    // *** CADSettings: simplifyBySquarefreeing
                    if( mSetting.simplifyBySquarefreeing() )
                    {
                        for( unsigned k = 0; k < currentEliminationList.size(); ++k )
                            currentEliminationList[k] = currentEliminationList[k].sepapart();
                    }
                    std::sort( currentEliminationList.begin(), currentEliminationList.end(), mSetting.mUP_isLess );
                    if( mSetting.simplifyBySquarefreeing() )
                        std::unique( currentEliminationList.begin(), currentEliminationList.end() );
                    newEliminationSets[i] = vector<UnivariatePolynomial>( currentEliminationList.begin(), currentEliminationList.end() );
                    // *** /CADSettings: simplifyBySquarefreeing
                }
                // apply heuristics
                // *** CADSettings: simplifyByRootcounting
                if( mSetting.simplifyByRootcounting() )
                {
                    vector<UnivariatePolynomial> baseLevel = mEliminationSets.back();
                    for( unsigned i = 0; i < baseLevel.size(); )
                    {
                        if( baseLevel[i].degree() % 2 == 0 )
                        {    // there could only be complex roots if the degree is even
                            RationalUnivariatePolynomial p( baseLevel[i] );
                            if( p.countRealRoots() == 0 )
                            {
                                baseLevel.erase( baseLevel.begin() + i );
                                continue;    // do not increase i
                            }
                        }
                        ++i;
                    }
                    mEliminationSets.back() = baseLevel;
                }
                // *** CADSettings: /simplifyByRootcounting
                mVariables       = newVariables;
                mEliminationSets = newEliminationSets;
                mIsComplete      = false;
                // reset all lifting positions
                mLiftingPositions = vector< list<unsigned> >( mVariables.size(), list<unsigned>() );
                for( unsigned i = 0; i < mVariables.size(); ++i )
                {
                    for( unsigned j = 0; j < mEliminationSets[i].size(); ++j )
                        mLiftingPositions[i].push_back(j);
        //            l.sort( isLessInLiftingPositions( mEliminationSets[i], mSetting.mUP_isLess ) );
                }
            }

            // PUBLIC STATIC METHODS //

            static const UnivariatePolynomialSet eliminationSet( const UnivariatePolynomialSet& P,
                                                                 const symbol& nextVariable )
                    throw ( invalid_argument );

            static const SampleList samples( const list<RealAlgebraicNumberPtr>& roots, SampleList& currentSamples ) throw ( invalid_argument );

            static const SampleList samples( const list<RationalUnivariatePolynomial>& polynomials,
                                             SampleList& currentSamples )
                    throw ( invalid_argument );

            // ATOMIC METHODS //

            static void elimination( const UnivariatePolynomial& p, const symbol& variable, UnivariatePolynomialSet& eliminated ) throw ( invalid_argument );

            static void elimination( const UnivariatePolynomial& p,
                                                                 const UnivariatePolynomial& q,
                                                                 const symbol& variable,
                                                                 UnivariatePolynomialSet& eliminated )
                    throw ( invalid_argument );

            static const SampleList samples( const RationalUnivariatePolynomial& p,
                                             SampleList& currentSamples,
                                             CADSettings settings = CADSettings::getSettings() )
                    throw ( invalid_argument );

            static const SampleList samples( const UnivariatePolynomial& p,
                                             const list<RealAlgebraicNumberPtr>& sample,
                                             const list<symbol>& variables,
                                             SampleList& currentSamples,
                                             CADSettings settings = CADSettings::getSettings() )
                    throw ( invalid_argument );

        private:

            // Heuristics //

            void alterSetting( unsigned setting = DEFAULT_CADSETTING,
                               RealAlgebraicNumberSettings::IsolationStrategy isolationStrategy = RealAlgebraicNumberSettings::DEFAULT_ISOLATIONSTRATEGY )
            {
                mSetting = CADSettings::getSettings( setting, isolationStrategy );
            }

            // ATTRIBUTES //

            vector<symbol> mVariables;
            tree<RealAlgebraicNumberPtr> mSampleTree;
            vector<SampleList> mSampleListIncrements;
            vector<vector<UnivariatePolynomial> > mEliminationSets;
            vector<list<unsigned> > mLiftingPositions;
            bool mIsComplete;
            CADSettings mSetting;

            // Auxiliary methods //

            inline const RealAlgebraicPoint constructSampleAt( tree<RealAlgebraicNumberPtr>::iterator node,
                                                               const tree<RealAlgebraicNumberPtr>::iterator& root ) const;

            inline const bool liftCheck( tree<RealAlgebraicNumberPtr>::iterator node,
                                         const list<RealAlgebraicNumberPtr>& sample,
                                         unsigned level,
                                         const list<symbol>& variables,
                                         const vector<Constraint>& constraints,
                                         RealAlgebraicPoint& r )
                    throw ( invalid_argument );

            inline const bool satisfys( const RealAlgebraicPoint& r, const vector<Constraint>& constraints ) const;

            // STATIC AUXILIARY METHODS //

            static const UnivariatePolynomialSet truncation( const UnivariatePolynomial& P );

            static const UnivariatePolynomialSet truncation( const UnivariatePolynomialSet& P );
    };
}    // namespace GiNaC
#endif