system_system.h

Go to the documentation of this file.

/********************************************************************
 *   This file is part of the source code of the realroot library.
 *   Author(s): J.P. Pavone, GALAAD, INRIA
 *   $Id: system.h,v 1.1 2005/07/11 10:03:57 jppavone Exp $
 ********************************************************************/
#ifndef realroot_SOLVE_SBDSLV_SYSTEM_HPP
#define realroot_SOLVE_SBDSLV_SYSTEM_HPP
//--------------------------------------------------------------------
#include <assert.h>
#include <realroot/rounding_mode.hpp>
//#include <realroot/bernstein_eenv.hpp>
#include <realroot/tensor_eenv.hpp>
#include <realroot/tensor_eenv_loops.hpp>
#include <fstream>
#include <realroot/method_base.hpp>
#include <realroot/system_dstack.h>
#include <realroot/system_pstack.h>

//--------------------------------------------------------------------
namespace mmx {
//--------------------------------------------------------------------
namespace realroot
{
template<typename real_t, bool rsafe = false, bool tsafe = false>
struct system : dstack<real_t>, pstack<real_t>
{
    /* types */
    typedef real_t creal_t;
    typedef pstack<creal_t> pstack_t;
    typedef dstack<creal_t> dstack_t;
    typedef typename dstack_t::interval_t interval_t;
    typedef typename pstack_t::vstack_t vstack_t;
    typedef system< creal_t, rsafe, tsafe > self_t ;
    typedef tensor::eenv eenv_t;
    typedef int sz_t;
    
    /* description du systeme */
    sz_t         m_neqs;
    eenv_t     * eenvs;
    
    static const int stRestricted  = 1; /* restriction disponible        */
    static const int stElevated    = 2; /* elevation disponible          */
    static const int stTransformed = 4; /* transformation disponible     */
    static const int stMonomials   = 8; /* version monomiale disponible  */
    /* etats des equations, selection d'equations, selection des variables     */
    sz_t *m_estate,*m_eqssel, *m_vrssel;
    /* taille max, taille totale, addresses des équations       */
    sz_t m_esz,m_ssz,*m_esft;
    /* systeme original, formes monomiales, buffer d'evaluation, restrictions  */
    
    creal_t *m_errors;
    creal_t *m_data,*m_mdata,*m_evchnk,*m_mnrdata,*m_mxrdata;
    creal_t      m_seps;
    binomials< creal_t > m_bzenv; /* cache des coefficients binomiaux (conversions) */
    

    struct rstcmd_t     /* commande de restriction = ordre d'application de de casteljau */
    {
        int m_type;
        static const int RSTCMDRL = 0; /* droite puis gauche */
        static const int RSTCMDLR = 1; /* gauche puis droite */
        creal_t m_aa; creal_t m_ab; /* valeur utilisée à gauche */
        creal_t m_ba; creal_t m_bb; /* valeur utilisée à droite */
        rstcmd_t(){};
        rstcmd_t( int t, const creal_t& a, const creal_t& b ) : m_type(t)
        {
            using namespace numerics::rdw;
            m_ab = a;
            m_ba = b;
            m_aa = upsub((creal_t)1.0,a);
            m_bb = upsub((creal_t)1.0,b);
        };
    };
    
    std::ostream& print( std::ostream& o, const rstcmd_t& rmc )
    {
        o << rmc.m_aa << " " << rmc.m_ab << ", " << rmc.m_ba << ", " << rmc.m_bb;
        return o;

    };
    rstcmd_t * m_rstcmd; /* commandes de restriction pour chaque variable */

    /* acces environement global  */
    inline sz_t  nvr() const { return dstack_t::nvr(); };
    inline sz_t  neqs()  const { return m_neqs; };
    /* acces environements locaux */
    inline eenv_t * env   ( sz_t i ) const { return eenvs + i; };
    inline sz_t    nvars  ( sz_t e )         const { return eenvs[e].nvr(); };
    inline sz_t    size   ( sz_t e )         const { return eenvs[e].sz();  };
    inline sz_t    size   ( sz_t e, sz_t v ) const { return eenvs[e].sz(v); };
    inline sz_t    varname( sz_t e, sz_t v ) const { return eenvs[e].vr(v); };
    inline sz_t    stride ( sz_t e, sz_t v ) const { return eenvs[e].st(v); };
    /* precision (securité) */
    inline  const  creal_t & spceps()        const { return m_seps; };
    /* acces aux équations */
    inline         creal_t * data ( sz_t e )       { return    m_data+m_esft[e]; };
    inline  const  creal_t * data ( sz_t e ) const { return    m_data+m_esft[e]; };
    inline         creal_t * rdata( sz_t e )       { return m_mnrdata+m_esft[e]; };
    inline  const  creal_t * rdata( sz_t e ) const { return m_mnrdata+m_esft[e]; };
    inline         creal_t * rmndata(sz_t e)       { return m_mnrdata+m_esft[e]; };
    inline  const  creal_t * rmndata(sz_t e) const { return m_mnrdata+m_esft[e]; };
    inline         creal_t * rmxdata(sz_t e)       { return m_mxrdata+m_esft[e]; };
    inline  const  creal_t * rmxdata(sz_t e) const { return m_mxrdata+m_esft[e]; };

    /* acces aux équations (monomiales) */
    
    void tomonoms( sz_t e )
    {
        std::copy(data(e),data(e)+size(e),m_mdata+m_esft[e]);
        convertb2m(m_mdata+m_esft[e],eenvs[e],m_bzenv);
        m_estate[e] |= stMonomials;
    };
    
    inline creal_t * mdata( sz_t e )
    {
        if ( m_estate[e] & stMonomials ) return m_mdata + m_esft[e]; else tomonoms(e);
        return mdata(e);
    };
    
    void meval( creal_t * dst, sz_t * eqs, sz_t n, creal_t * prms, int nsmp = 1 )
    {
        /* dst+i*n pour i = 0..nsmp pointe sur n valeurs correspondant aux évaluation des équations    */
        /* eqs[j], j = 0..n sous forme monomiales sur les valeurs pointées par prms+i*n, (sens global) */
        creal_t locals[nvr()];
        for ( sz_t i = 0; i < nsmp; i++, dst += n, prms += nvr()  )
            for ( sz_t ie = 0; ie < n; ie ++ )
            {
                sz_t e = eqs[ie]; assert(e<neqs());
                for ( sz_t v = 0; v < nvars(e); locals[v] = prms[varname(e,v)], v++ ) ;
                levalm(dst[e],mdata(e),eenvs[e],locals);
            };
    };
    
    inline void meval( creal_t * dst, sz_t e, creal_t * prms, int nsmp = 1 )
    { sz_t eqs[]={e}; meval(dst,eqs,1,prms,nsmp); };
    
    void set_monomials( creal_t * src, sz_t e, bool scale_ = false )
    {
        m_estate[e] = stMonomials;
        std::copy(src,src+size(e),data(e));
        std::copy(src,src+size(e),m_mdata+m_esft[e]);
        convertm2b(data(e),eenvs[e],m_bzenv);
    };

    void set_bernstein( creal_t * src, sz_t e )
    {
        m_estate[e] = 0;
        std::copy(src,src+size(e),data(e));
    };

    void set_monomials2( const creal_t * src, sz_t e, bool scale_ = false )
    {
        m_estate[e] = stMonomials;
        std::copy(src,src+size(e),data(e));
        std::copy(src,src+size(e),m_mdata+m_esft[e]);
        convertm2b(data(e),eenvs[e],m_bzenv);
    };

    void set_bernstein2( const creal_t * src, sz_t e )
    {
        m_estate[e] = 0;
        std::copy(src,src+size(e),data(e));
    };

    void _rstcmd_() /* mise à jour des commandes de restriction par rapport au domaine courant */
    {
        numerics::rounding<creal_t> rnd( numerics::rnd_dw() );
        using namespace numerics::rdw;
        for ( sz_t v = 0; v < nvr(); v ++ )
        {
            creal_t a = this->current()[v].lower(); // -2
            creal_t b = this->current()[v].upper(); //  0
            assert(a<=b);
            if ( a > (creal_t)1.0 )
                new (m_rstcmd+v) rstcmd_t(rstcmd_t::RSTCMDRL,dwdiv(a,b),b);
            else
                new (m_rstcmd+v) rstcmd_t(rstcmd_t::RSTCMDLR,a,updiv(upsub(b,a),dwsub((creal_t)1.0,a)));
        };
    };
    
    void _rst_( sz_t e, creal_t * base  ) /* restriction de l'equation e dans le buffer de restriction donné par base */
    {
        std::copy(data(e),data(e)+size(e),base+m_esft[e]);

        for ( sz_t v = 0; v < nvars(e); v ++ )  /* pour toute les variables de l'équation e */
        {
            /*   application des commandes de restriction correspondantes (passage local->global) */
            switch( m_rstcmd[ varname(e,v) ].m_type )
            {
            case rstcmd_t::RSTCMDRL: /* restriction à droite puis à gauche */
                brestrictRL(base+m_esft[e],eenvs[e],v,
                            m_rstcmd[varname(e,v)].m_aa,
                            m_rstcmd[varname(e,v)].m_ab,
                            m_rstcmd[varname(e,v)].m_ba,
                            m_rstcmd[varname(e,v)].m_bb );
                break;
            case rstcmd_t::RSTCMDLR: /* restriction à gauche puis à droite */
                brestrictLR(base+m_esft[e],eenvs[e],v,
                            m_rstcmd[varname(e,v)].m_aa,
                            m_rstcmd[varname(e,v)].m_ab,
                            m_rstcmd[varname(e,v)].m_ba,
                            m_rstcmd[varname(e,v)].m_bb );
                break;
            };
        };
    };
    
    
private:
    
    void restrict( sz_t * eqs , sz_t n, creal_t * ) /* restriction des équations données par eqs[i], i = 0..n-1  */
    {
        for ( sz_t e = 0; e < n;  e ++ )
        {
            if ( m_estate[eqs[e]] & stRestricted  ) continue; /* on vérifie que l'équation n'a pas déjà été restreinte */
            _rst_(eqs[e],m_mnrdata);                          /* sinon on calcule sa restriction                       */
            m_estate[eqs[e]] |= stRestricted;                 /* on note que la restriction est maintenant connue      */
        };
    };
    
    inline void restrict( sz_t * eqs, sz_t n, const texp::false_t&  ) /* restriction en arrondi au plus proche */
    {
        _rstcmd_();                                            /* calcul des commandes de restriction      */
        assert(m_mnrdata==m_mxrdata);
        numerics::rounding<creal_t> rnd( numerics::rnd_nr() ); /* passage en mode d'arrondi au plus proche */
        restrict(eqs,n,m_mnrdata);                             /* restriction                              */
    };
    
    void GETERROR( double * mn_, double * mx, int e )
    {
        int sze = size(e);
        creal_t err = 0;
        for ( double * mn = mn_; mn != mn_+sze; mn++, mx++ )
        {
            creal_t d   = *mx-*mn;
            if ( d < 0 ) { *mn=*mx; d = -d; };
            if ( d > err ) err = d;
        };
    };
    
    inline void restrict( sz_t * eqs, sz_t n, const texp::true_t&   ) /* restriction en arrondi inférieur et supérieur */
    {
        _rstcmd_();                                                     /* calcul des commandes de restriction */
        {
            numerics::rounding<creal_t> rnd( numerics::rnd_dw() );        /* passage en mode d'arrondi inférieur */
            restrict(eqs,n,m_mnrdata);                                    /* restriction                         */
        };
        {
            numerics::rounding<creal_t> rnd( numerics::rnd_up() );        /* passage en mode d'arrondi inférieur */
            restrict(eqs,n,m_mxrdata);                                    /* restriction                         */
        };
        for ( int e = 0; e < n; e ++ )
        {
            GETERROR(rmndata(eqs[e]),rmxdata(eqs[e]),eqs[e]);
        };
    };
    
public:
    inline void restrict( sz_t * eqs = 0, sz_t n  = 0 )              /* restrict, dispatch sur les fonctions precedentes */
    {

        restrict( (eqs)?eqs:m_eqssel, (n)?n:neqs(), typename  texp::BoolType< rsafe >::result_t() ); };
    inline void restrict( sz_t e )                                   /* restriction d'une seule  équation */
    { sz_t tmp[] = { e }; restrict( tmp, 1 ); };

    //  private:
    
    eenv_t       m_elenv, *m_elenvs;                                 /* environement(s) d'elevation du systeme   */
    creal_t   ** m_elbins;                                           /* coefficients binomiaux pour l'élévation  */
    unsigned  ** m_eladdr;                                           /* addresses de destination des coefficients*/
    creal_t    * m_mnelldata;                                        /* élévation des restrictions basses        */
    creal_t    * m_mxelldata;                                        /* élévation des restrictions hautes        */
    creal_t    * m_mntrnsfdata;                                      /* transformations basses                   */
    creal_t    * m_mxtrnsfdata;                                      /* transformations hautes                   */
    
    /* precalcul pour l'élévation                                                       */
    /* --------------------------                                                       */
    /* avantages:                                                                       */
    /*              - elevation reduite à un produit scalaire                           */
    /* inconvenients:                                                                   */
    /*              - consomme beaucoup de mémoire                                      */
    /*              - risque d'overflow existe                                          */
    
    void prcinfo( unsigned * address, creal_t * bindata, eenv_t * o, eenv_t * a, eenv_t * b )
    {
        int ia, ib, add, vo, ca, cb;
        sz_t v;
        int vidx = o->mxvr(); /* index de la variable maximum (scotch dans eenv) */
        sz_t adda[vidx+1];
        sz_t addb[vidx+1];
        unsigned c = 0;
        creal_t sumsc = 0.0;

        std::fill( adda, adda + vidx+1, 0 );
        for ( ia = 0; ia < a->sz(); ia ++ )
        {
            std::fill( addb, addb + vidx+1, 0 );
            for ( ib = 0; ib < b->sz(); ib ++ )
            {
                creal_t sc = 1.0;
                /* on risque l'overflow ici */
                for ( v  = 0; v < a->nvr(); sc*=m_bzenv(a->sz(v)-1,adda[a->vr(v)]), v ++ ) ;
                for ( v  = 0; v < b->nvr(); sc*=m_bzenv(b->sz(v)-1,addb[b->vr(v)]), v ++ ) ;
                for ( add = vo = 0; vo < o->nvr(); vo++ )
                {
                    if ( o->sz(vo) > 1 )
                    {
                        sc  /= m_bzenv( o->sz(vo)-1, adda[o->vr(vo)] + addb[o->vr(vo)] );
                        add += (adda[o->vr(vo)]+addb[o->vr(vo)])*o->st(vo);
                    };
                }
                sumsc += sc;
                address[c]   = add;
                bindata[c++] = sc;
                cb = b->nvr()-1;
                addb[b->vr(cb)] ++;
                while ( cb>0 && (addb[b->vr(cb)]==b->sz(cb)) )
                {
                    addb[b->vr(cb)] = 0;
                    cb--;
                    addb[b->vr(cb)]++;
                };
            }
            ca = a->nvr()-1;
            adda[a->vr(ca)]++;
            while ( ca >0 && (adda[a->vr(ca)] == a->szs()[ca]))
            {
                adda[a->vr(ca)] = 0;
                ca--;
                adda[a->vr(ca)]++;
            };
        }
    };
    
    /* collecte des informations utiles pour l'élévation (si requise) */
    void _alloc_elldata_()
    {
        numerics::rounding<creal_t> rnd(numerics::rnd_nr());
        assert(m_elenvs==0);
        //int nvrMin = std::min(1000000, nvr());
        int nvrMin = nvr();
        sz_t evrs[nvrMin];

        sz_t eszs[nvrMin];

        m_errors = new creal_t[ neqs() ];

        for ( sz_t v = 0; v < nvr(); evrs[v] = v, eszs[v++] = 0 ) ;
        for ( sz_t e = 0; e < neqs(); e ++ )                          /* recuperation de l'environement élevé */
            for ( sz_t v = 0; v < nvars(e); v ++ )
                if ( eszs[varname(e,v)] < size(e,v) ) eszs[varname(e,v)] = size(e,v);

        m_elenv.~eenv_t();
        new (&m_elenv) eenv_t(nvr(),eszs,evrs);

        m_elenvs = new eenv_t[ neqs() ]; /* allocation des "environements d'élevation". (differences avec m_elenv) */
        assert(m_elenvs);

        unsigned ssz = 0; /* quantité de coefficients necessaires aux calculs d'élevation */
        for ( sz_t e = 0; e < neqs(); e ++ ) /* pour tout nos environements  */
        {
            /* s'il ne contient pas la variable v on initialise sa difference avec l'environement élevé */
            for ( sz_t v = 0; v < nvr(); eszs[v] = m_elenv.szs()[v], v++ ) ;
            /* sinon pour chacune de ses variables on calcule la différence */
            for ( sz_t v = 0; v < nvars(e); v ++ ) eszs[varname(e,v)] = m_elenv.szs()[varname(e,v)]-size(e,v)+1;
            /* on initialise l'environement d'élevation que l'on vient de calculer */
            (m_elenvs+e)->~eenv_t();
            new (m_elenvs+e) eenv_t(nvr(),eszs,evrs);
            ssz += m_elenvs[e].sz() * env(e)->sz(); /* voir def de ssz */
        };

        assert(sizeof(creal_t*)==sizeof(unsigned*)); /* normalement une addresse est une addresse */
        m_elbins = new creal_t* [ 2*neqs() + ssz];
        assert(m_elbins);
        m_eladdr = (unsigned**)(m_elbins + neqs());
        unsigned * achnk = (unsigned*)(m_eladdr + neqs());
        unsigned  ellssz = m_elenv.sz() * neqs();
        creal_t  * cchnk = new creal_t[ ssz + ((rsafe?1:2)+(tsafe?1:2))*ellssz];
        assert(cchnk);
        m_mnelldata = cchnk + ssz;
        m_mxelldata   = rsafe?(m_mnelldata+ellssz):m_mnelldata;
        m_mntrnsfdata = m_mxelldata+ellssz;
        m_mxtrnsfdata = tsafe?(m_mntrnsfdata+ellssz):m_mntrnsfdata;

        for ( sz_t e = 0; e < neqs(); e ++ )
        {
            unsigned isz = m_elenvs[e].sz() * eenvs[e].sz();
            m_elbins[e] = cchnk;
            m_eladdr[e] = achnk;
            cchnk += isz;
            achnk += isz;
            prcinfo( m_eladdr[e], m_elbins[e], &m_elenv, m_elenvs+e, eenvs+e );
        };
    };
    
    void _free_elldata_()
    {
        if ( m_elenvs ) delete[] m_elenvs;
        delete[] m_elbins[0];
        delete[] m_elbins;
    };
    
    inline creal_t * mnelldata  ( sz_t e  )  { return m_mnelldata   + m_elenv.sz() * e; };
    inline creal_t * mxelldata  ( sz_t e  )  { return m_mxelldata   + m_elenv.sz() * e; };
    inline creal_t * elldata( sz_t e )       { return mnelldata(e); };
    inline creal_t * mntrnsfdata( sz_t e )   { return m_mntrnsfdata + m_elenv.sz() * e;};
    inline creal_t * mxtrnsfdata( sz_t e )   { return m_mxtrnsfdata + m_elenv.sz() * e;};
    inline creal_t * trnsfdata  ( sz_t e  )  { return mntrnsfdata(e); };
    
    /* boucle d'élevation */
    void elevate_loop( sz_t e, creal_t * dst, creal_t * src, creal_t * bins, unsigned * address )
    {
        int c = 0;
        std::fill( dst, dst + m_elenv.sz(), 0.0 );
        for ( sz_t ia = 0; ia < m_elenvs[e].sz(); ia++ )
            for ( sz_t ib = 0; ib < eenvs[e].sz(); ib ++, c++ )
                dst[address[c]] += bins[c] * src[ib];
        /* pourrait être remplacé par                            */
        /* (1)  dst[obins[ia]+osrc[ib]] += bins[c] * src[ib]     */
        /*      pour des raisons d'économie mémoire              */
        /* (2)  bins[c] par sqbinsa[c]*sqbinsb[c]*src[a]         */
        /*      où sqbinsa[c] = binsa[ia]/sqrt(binsd[ia+ib])     */
        /*      pour éviter les dépassements de capacités        */
        /*      les racines carrées devant être bien-entendues   */
        /* précalculées à la création                            */
        /* la meilleure solution serait d'elever progressivement */
        /* l'équation, mais le code n'est pas encore clair pour  */
        /* l'instant, quand ce sera le cas on sortira l'élevation*/
        /* du systeme comme pour la pile de projection           */
    };

    /* elevation de l'equation e (safe = false) */
    inline void elevate( sz_t e, const texp::false_t& )
    {
        numerics::rounding<creal_t> rnd( numerics::rnd_nr() );
        elevate_loop( e, elldata(e), rdata(e),m_elbins[e],m_eladdr[e] );
    };
    
    /* elevation de l'equation e (safe = true) */
    inline void elevate( sz_t e, const texp::true_t& )
    {
        {
            numerics::rounding<creal_t> rnd( numerics::rnd_dw() );
            elevate_loop( e, mnelldata(e), rmndata(e), m_elbins[e], m_eladdr[e] );
        };
        {
            numerics::rounding<creal_t> rnd( numerics::rnd_up() );
            elevate_loop( e, mxelldata(e), rmxdata(e), m_elbins[e], m_eladdr[e] );
        };
    };
    
    inline void elevate( sz_t * eqs = 0, sz_t n = 0 )
    {
        eqs = eqs?eqs:m_eqssel;
        n   = n?n:neqs();
        if ( m_elenvs == 0 ) _alloc_elldata_();
        for ( sz_t i = 0; i < n; i ++ )
        {
            if ( m_estate[eqs[i]] & stElevated ) continue;
            if ( !(m_estate[eqs[i]] & stRestricted) ) { this->error("elevate: restriction needed"); continue; };
            /* pour le cas ou m_elenv est l'environement d'un équation du système (non-exploité) */
            if ( mnelldata(eqs[i]) != rmndata(eqs[i]) )
                elevate( eqs[i],typename texp::BoolType<rsafe>::result_t() );
            m_estate[eqs[i]] |=  stElevated;
        };
    };
    
    inline void force_elevate( sz_t * eqs = 0, sz_t n  = 0 ) { restrict(eqs,n); elevate(eqs,n); };
    inline void elevate( sz_t e ) { sz_t eqs[] ={ e }; elevate( eqs, 1 ); };
    
    void _safe_transform_( creal_t * mat, sz_t * eqs, sz_t n,
                           creal_t * trsf_up, creal_t * trsf_dw,
                           creal_t * ell_up,  creal_t * ell_dw )
    {
        numerics::rounding<creal_t> rnd( numerics::rnd_dw() );
        creal_t *dst_up,*dst_dw,*src_up,*src_dw,cm;
        for ( sz_t l = 0; l < n;  l ++ )
        {
            m_estate[eqs[l]] |= stTransformed;
            dst_up = trsf_up + eqs[l] * m_elenv.sz();
            dst_dw = trsf_dw + eqs[l] * m_elenv.sz();
            std::fill( dst_up, dst_up + m_elenv.sz(), (creal_t)0.0 );
            std::fill( dst_dw, dst_dw + m_elenv.sz(), (creal_t)0.0 );
            for ( sz_t c = 0; c < n; c ++ )
            {
                cm = mat[l*n+c];
                src_up = ell_up + eqs[c] * m_elenv.sz();
                src_dw = ell_dw + eqs[c] * m_elenv.sz();
                if ( cm > 0 )
                {
                    for ( sz_t k = 0; k < m_elenv.sz(); dst_dw[k] += cm*src_dw[k], k ++ ) 
                        ;
                    {
                        numerics::rounding<creal_t> rnd(numerics::rnd_up());
                        for ( sz_t k = 0; k < m_elenv.sz(); k ++ )
                            dst_up[k] += cm*src_up[k];
                    };
                };
                if ( cm < 0 )
                {
                    for ( sz_t k = 0; k < m_elenv.sz(); dst_dw[k] += cm*src_up[k], k ++ ) 
                        ;
                    {
                        numerics::rounding<creal_t> rnd(numerics::rnd_up());
                        for ( sz_t k = 0; k < m_elenv.sz(); k ++ )
                            dst_up[k] += cm*src_dw[k];
                    };
                };
            };
        };
    };
    
    void _transform_( creal_t * mat, sz_t * eqs, sz_t n, creal_t * basetrsf, creal_t * baseell  )
    {
        creal_t *dst,*src,cm;
        for ( sz_t l = 0; l < n;  l ++ )
        {
            m_estate[eqs[l]] |= stTransformed;
            dst = basetrsf + eqs[l] * m_elenv.sz();
            std::fill( dst, dst + m_elenv.sz(), (creal_t)0.0 );
            for ( sz_t c = 0; c < n; c ++ )
            {
                cm = mat[l*n+c];
                src = baseell + eqs[c] * m_elenv.sz();
                for ( sz_t k = 0; k < m_elenv.sz(); dst[k] += cm*src[k], k ++ ) ;
            };
        };
    };
    
    inline void transform( creal_t * mat, sz_t * eqs, sz_t n  )
    {
        elevate( eqs, n );
        if ( tsafe )
        {
            { numerics::rounding< creal_t > rnd( numerics::rnd_dw() );
                _transform_( mat, eqs, n, m_mntrnsfdata, m_mnelldata ); };
            {numerics::rounding< creal_t > rnd( numerics::rnd_up() );
                _transform_( mat, eqs, n, m_mxtrnsfdata, m_mxelldata );};
            //_safe_transform_(mat,eqs,n,m_mxtrnsfdata,m_mntrnsfdata,m_mxelldata, m_mnelldata );
        }
        else
        {
            numerics::rounding< creal_t > rnd( numerics::rnd_nr() );
            _transform_( mat, eqs, n, m_mntrnsfdata, m_mnelldata );
        };
    };
    
    bool project( sz_t * eqs = 0 , sz_t ne = 0 )
    {
        ne = ne?ne:neqs();
        eqs=eqs?eqs:m_eqssel;
        unsigned nobj = 0;
        creal_t * pmn, * pmx;
        for ( sz_t i = 0; i < ne; i ++ )
        {
            sz_t e = eqs[i];
            if ( m_estate[e] & stTransformed )
            {
                //        std::cout << "project " << i << " transformed \n";
                nobj ++;
                for ( sz_t v = 0; v < nvr(); v ++ )
                {
                    pstack_t::push(v,pmn,pmx,m_elenv.sz(v));
                    mins(pmn,mntrnsfdata(e), m_elenv, v );
                    maxs(pmx,mxtrnsfdata(e), m_elenv, v );
                };
                continue;
            };

            if ( m_estate[e] & stRestricted )
            {
                nobj ++;
                for ( sz_t lv = 0; lv < env(e)->nvr(); lv ++ )
                {
                    pstack_t::push(varname(e,lv),pmn,pmx,size(e,lv));
                    mins(pmn,rmndata(e),eenvs[e],lv);
                    maxs(pmx,rmxdata(e),eenvs[e],lv);
                };
                continue;
            };
        };

        if ( nobj == 0 ) this->error("no restricted equations to project!\n");
        return true;
    };
    
    /* reinitialisation de l'état des équations */
    inline void dreset() {       dstack_t::reset(); };
    bool thickness()
    {
        return pstack_t::thickness( m_seps);
    };
    inline void reset()
    {
        dstack_t::reset();
        pstack_t::reset();
        /* la seule information indépendante du domaine est la forme monomiale */
        for ( sz_t e = 0; e < neqs(); m_estate[e] &= stMonomials, e ++ ) ;
    };
    
    inline void init( creal_t * inits = 0 )
    {
        dstack_t::init( inits ); /* initialisation du domaine de recherche */
        reset();
    };
    
    
    system() { eenvs = 0; };

    system( sz_t neqs_, sz_t nvars_, sz_t const * nvrs, sz_t const ** vars, sz_t const ** szs,
            const creal_t& prmeps = (creal_t)1e-6,
            const creal_t& spceps = (creal_t)1e-12, const creal_t& isoleps = (creal_t)1e-6, bool use_prcnd = true ) :
        dstack_t(nvars_,prmeps,isoleps),
        pstack_t(nvars_)
    {
        assert(use_prcnd || !use_prcnd );
        m_seps = spceps;
        sz_t i;
        m_neqs = neqs_;
        m_rstcmd = new rstcmd_t[ nvr() ];
        eenvs  = new eenv_t[ neqs() ];
        for (i=0;i<neqs();i++) eenvs[i] = eenv_t( nvrs[i], szs[i], vars[i] );
        m_ssz = 0;
        m_esft  = new sz_t [ neqs()  /* offsets   */+
                neqs()  /* états     */+
                neqs()  /* séléction */+
                nvr() /* séléction */ ];
        assert(m_esft);
        m_estate = m_esft   + neqs(), m_eqssel = m_estate + neqs(), m_vrssel = m_eqssel + neqs();
        for ( i = 0; i < nvr(); m_vrssel[i] = i, i ++ ) ;
        for ( i = m_esz = 0; i < neqs(); m_ssz += size(i), i ++ )
        {
            m_estate[i] = 0; m_eqssel[i] = i; m_esft[i]  = m_ssz;
            m_esz = std::max(m_esz,size(i));
        };

        m_elenvs  = 0;
        m_data    = new creal_t[ ((rsafe)?4:3)*m_ssz + m_esz ];
        assert( m_data );
        m_mnrdata = m_data + m_ssz;
        m_mxrdata = rsafe?(m_mnrdata+m_ssz):m_mnrdata;
        m_mdata   = m_mxrdata  + m_ssz;
        m_evchnk  = m_mdata    + m_ssz;
        _alloc_elldata_();
    };
    
    ~system()
    {
        if ( m_errors ) delete[] m_errors;
        if ( m_rstcmd ) delete[] m_rstcmd;
        if ( eenvs  ) delete[] eenvs;
        if ( m_esft   ) delete[] m_esft;
        if ( m_data   ) delete[] m_data;
        _free_elldata_();
    };

    
    // gloutonneries
    struct cmp_sze
    {
        self_t * m_slv;
        bool operator()( int a, int b ) const
        { return m_slv->env(a)->sz() < m_slv->env(b)->sz(); };
        cmp_sze( self_t * slv ) : m_slv(slv){};
    };
    struct cmp_szv
    {
        self_t * m_slv;
        bool operator()( int a, int b ) const
        { return m_slv->nvars(a) < m_slv->nvars(b); };
    };
    /* tri sur la taille              */
    void szssort()  { std::sort( m_eqssel, m_eqssel + neqs(), cmp_sze(this) ); };
    /* tri sur le nombre de variables */
    void szvsort()  { std::sort( m_eqssel, m_eqssel + neqs(), cmp_szv(this) ); };


};

};



//--------------------------------------------------------------------
} //namespace mmx
/********************************************************************/
#endif //