tensor_eenv_fcts.hpp

Go to the documentation of this file.

#ifndef realroot_SBDSLV_EENV_C
#define realroot_SBDSLV_EENV_C
#include <realroot/tensor_vctops.hpp>
namespace mmx { 
namespace tensor
{
  inline void eenv::swap (eenv & e) { std::swap (data, e.data); };

  inline int        eenv::nvr () const { return data[1]; };

  inline const int *eenv::szs () const { return data + 2; };
  inline int       *eenv::szs ()       { return cdup ().data + 2; };
  inline int       *eenv::szs_data ()  { return data + 2; };

  inline int    eenv::sz (int v) const { return szs ()[v]; };
  inline int    eenv::sz () const      { return str ()[-1]; };

  inline int       *eenv::vrs ()       { return szs () + nvr (); };
  inline const int *eenv::vrs () const { return szs () + nvr (); };

  inline int       *eenv::str ()       { return vrs () + nvr () + 1; }; 
  inline const int *eenv::str () const { return vrs () + nvr () + 1; };

  inline int eenv::st (int v) const { return str ()[v]; };

  inline int eenv::vr (int v) const { return vrs ()[v]; };

  inline int eenv::msz (int nvr_) { return 1 + 1 + nvr_ + nvr_ + nvr_ + 1 + 2 * nvr_; };
  inline int eenv::msz () const   { return eenv::msz (nvr ()); };

  inline bool eenv::operator< (const eenv & b) const { return sz () < b.sz (); };

  inline int *eenv::alloc (int nvr_) { 
    //    return this->data = new int[msz (nvr_)]; 
    //    data = new int[eenv::msz (nvr_)]; 
    new (data)int[eenv::msz (nvr_)]; 
    return data;
  }

  inline void eenv::construct (int nvr_, int const *szs_, int const *vrs_,
                   int *mem)
  {
#ifdef EENV_DEBUG
    std::cout << "eenv::construct begin\n";
    std::cout <<nvr_<<" "<<mem<<" "<<msz (nvr_)<<std::endl;
#endif
    data = (mem==0)?alloc(nvr_):mem;
    data[0] = 1; 
    data[1] = nvr_;
    if (vrs_)
      {
    std::pair<int,int> *tmp = new std::pair<int,int>[ nvr_ ];
    for ( int i = 0; i < nvr_; i ++ ) { 
      tmp[i].first  = vrs_[i]; 
      tmp[i].second = szs_[i]; 
    }
    std::sort(tmp,tmp+nvr_);
    for ( int i = 0; i < nvr(); i ++ ) { 
      data[2+nvr()+i] = tmp[i].first;    // vrs()[i] = tmp[i].first; 
      data[2+i] = tmp[i].second;   // szs()[i] = tmp[i].second; 
    };
    cstr ();
      }
    else
      if (szs_)
    {
      std::copy (szs_, szs_ + nvr_, data+2); //szs();
      cstr ();
      /*gcc-bug       for (int v = 0; v < nvr (); vrs()[v] = v ++ ); */
      for (int v = 0; v < nvr () ; data[2+nvr()+v] = v, v++ ) ;
    };
  };
  
  inline bool eenv::equal (const eenv & a, const eenv & b)
  {
        if (a.data == b.data)     return true;
    
    for (int i = 1; i < 2 + 2 * a.nvr (); i++)
        if (a.data[i] != b.data[i]) return false;
    return true;
  };
  
  inline void eenv::dealloc () { if (data && (--data[0] == 0) ) delete[] data; };

  inline eenv & eenv::copy (eenv & e, const eenv & b)
  {    
    //    std::cout<<"copy"<<std::endl;
    e.dealloc ();
    e.data = b.data;
    if (e.data)
      e.data[0]++;
    return e;
  };

  inline eenv & eenv::cdup ()
  {
    if (data[0] == 1)
      return *this;
    //    std::cout<<"cdup"<<std::endl;
    data[0]--;
    int *sv = data;
    int _sz = msz();
    data = alloc (nvr ());
    std::copy (sv, sv + _sz, data);
    return *this;
  };

  inline void eenv::new_copy (eenv & e, const eenv & b)
  {    
    e.dealloc ();
    unsigned sz = b.msz();
    e.data = new int[sz];
    for(unsigned i=0;i<sz;i++) e.data[i]=b.data[i];
    e.data[0]=1;
  }

  inline void eenv::cstr ()
  {
#ifdef EENV_DEBUG
    std::cout << "eenv::cstr()" << std::endl;
#endif
    int i, s = 2+2*nvr()+1;//* _st(str());
    const int *_sz(szs_data()); //data+2
    for ( i = nvr()-2, data[s+nvr()-1] = 1; i >= -1; data[s+i]=_sz[i+1]*data[s+i+1], i-- ) ;
#ifdef EENV_DEBUG
    std::cout << (*this) << std::endl;
#endif
    //    std::cout << "esz= " << _st[-1] << std::endl;
  };
  
  inline bool eenv::hasvar (int &lv, int gv) const
  {
    int g, d, m;
    const int * _vr = vrs();
    if ((gv < _vr[0]) || (gv > _vr[(d=nvr())-1])) return false;
    if ( gv == _vr[g = 0] ) { lv = 0; return true; };
    if ( gv == _vr[d=nvr()-1]) { lv = nvr()-1; return true;};
    do
      {
    m = (g + d) / 2;
    if ( _vr[m]  == gv ) { lv = m; return true; };
    if ( _vr[m] < gv ) g = m; else d = m;
      }
    while (d - g > 1);
    return false;
  };

  
  inline eenv::~eenv ()  { dealloc (); };
  
  inline eenv::eenv (int nvr_, int const *szs_, int const *vrs_, int *mem)
    :data(new int[msz(nvr_)])
  {
    construct (nvr_, szs_, vrs_, mem);
  };

  inline eenv::eenv (const eenv & e ):data(e.data) { 
    //    data = e.data; 
    if ( data ) data[0]++; 
  };
  inline eenv & eenv::operator= (const eenv & e) { return copy (*this, e); };
  inline eenv eenv::mrgvrs (const eenv & a, const eenv & b)
  {
    
    unsigned c = 0;
    int va, vb;
    int prv = (int) (-1);
    int anvr = a.nvr ();
    int bnvr = b.nvr ();
    const int * avr = a.vrs();
    const int * bvr = b.vrs() ;
    //const eenv::vr_t * agvr = a.gvr ();
    // const eenv::vr_t * bgvr = b.gvr ();

    int tmp[anvr + bnvr];

    c = va = vb = 0;

    do
      {
    if ( avr[va] <= bvr[vb] )
      {
        if ( avr[va] != prv ) tmp[c++] = prv = avr[va];
        va++;
        continue;
      };

    if ( bvr[vb] < avr[va] )
      {
        if ( bvr[vb] != prv ) tmp[c++] = prv = bvr[vb];
        vb++;
        continue;
      };
      }
    while (va < anvr && vb < bnvr);

    if (va < anvr) for (; va < anvr; va++)
      if ( avr[va] != prv ) tmp[c++] = avr[va];
    
    if (vb < bnvr) for (; vb < bnvr; vb++)
      if ( bvr[vb] != prv ) tmp[c++] = bvr[vb];
    
    eenv oenv (c);


    std::copy(tmp,tmp+oenv.nvr(),oenv.vrs());
    return oenv;
  };


  inline eenv eenv::common (const eenv & a, const eenv & b)
  {
      if ( a == b ) return eenv(a);
    int v, lva, lvb;
    eenv oenv (mrgvrs (a, b));

    const int *ovr = oenv.vrs ();
    const int *aszs = a.szs ();
    const int *bszs = b.szs ();
    int *oszs = oenv.szs ();

    for (v = 0; v < oenv.nvr (); v++)
      {
    if (a.hasvar (lva, ovr[v]))
      {
        if (b.hasvar (lvb, ovr[v]))
          {
        oszs[v] = std::max(aszs[lva], bszs[lvb]);
          }
        else
          oszs[v] = aszs[lva];
      }
    else
      {
        if (b.hasvar (lvb, ovr[v]))
          oszs[v] = bszs[lvb];
        else
          std::cout << "erreur dans common\n";
      };
      };

    oenv.cstr ();
    return oenv;
  };
  
  inline eenv eenv::mul (const eenv & a, const eenv & b)
  {

    int v, lva, lvb;
    eenv oenv (mrgvrs (a, b));
    //  std::cout << " mrgvrs = " << oenv << std::endl;
    int *oszs = oenv.szs ();
    const int *ovrs = oenv.vrs ();
    const int *aszs = a.szs ();
    const int *bszs = b.szs ();

    for (v = 0; v < oenv.nvr (); v++)
      {
    oszs[v] = 0;
    if (a.hasvar (lva, ovrs[v]))
      oszs[v] += aszs[lva] - 1;
    if (b.hasvar (lvb, ovrs[v]))
      oszs[v] += bszs[lvb] - 1;
    oszs[v]++;
      };
    oenv.cstr ();
    return oenv;
  };

  inline eenv eenv::diff (const eenv & b, int lv)
  {
    //    const int * bszs = b.szs();
    eenv oenv; eenv::new_copy(oenv,b);
    int *oszs = oenv.szs();
    oszs[lv] = std::max (oszs[lv] - 1, 1);
    oenv.cstr ();
    return oenv;
  };

  inline bool eenv::vmap (int *vmap, const eenv & o, const eenv & i)
  {
    for (int v = 0; v < i.nvr (); v++)
      {
    int lv;
    if (!o.hasvar (lv, i.vr (v)))
      {
        return false;
      };
    vmap[v] = lv;
      };
    return true;
  };


inline  bool eenv::oaddress (const eenv & oenv, unsigned *osupp,
               const eenv & ienv, unsigned *isupp, unsigned nsp)
  {
    //  std::cout << oenv << std::endl;
    //  std::cout << ienv << std::endl;
    int v;
    unsigned c;
    int vmap_[ienv.nvr ()];
    unsigned addi, addo;
    if (!vmap (vmap_, oenv, ienv))
      return false;
    if (isupp)
      for (c = 0; c < nsp; osupp[c] = addo, c++)
    for (addi = isupp[c], addo = 0, v = ienv.nvr () - 1; addi;
         addi /= ienv.sz (v), v--)
      addo += (addi % ienv.sz (v)) * oenv.st (vmap_[v]);
    else
      {
    const int *istr = ienv.str ();
    const int *ostr = oenv.str ();
    osupp[0] = 0;
    for (v = 0; v < ienv.nvr (); v++)
      for (int j = 0; j < istr[-1]; j += istr[v - 1])
        for (int i = j + istr[v]; i < j + istr[v - 1]; i += istr[v])
          osupp[i] = osupp[i - istr[v]] + ostr[vmap_[v]];
      };
    return true;
  };
  
  inline
  std::ostream & operator<< (std::ostream & out, const eenv & env)
  {
    out << "*eenv*\n";
    out << "\tpos: " << env.data    << "\n";
    out << "\tref: " << env.data[0] << "\n";
    out << "\tnvr: " << env.nvr()   << "\n";
    out << "\tesz: " << env.sz ()   << "\n";
    out << "\tszs: ";
    vct::print (env.szs (), env.nvr (), 1, out);
    out << "\n";
    out << "\tvrs: ";
    vct::print (env.vrs (), env.nvr (), 1, out);
    out << "\n";
    out << "\tstr: ";
    vct::print (env.str (), env.nvr (), 1, out);
    return out;
  };
  
  inline 
  bool eenv::subset (const eenv & a, const eenv & b)
  {
    const int *avrs = a.vrs ();
    const int *aszs = a.szs ();
    const int *bszs = b.szs ();
    int lv;
    for (int i = 0; i < a.nvr (); i++)
      if ((!b.hasvar (lv, avrs[i])) || aszs[i] > bszs[lv])
    return false;
    return true;
  };

  inline bool eenv::operator== (const eenv & e) const
  {
    return eenv::equal (*this, e);
  };

  inline bool eenv::operator!= (const eenv & e) const
  {
    return !eenv::equal (*this, e);
  };

  inline eenv eenv::elevation (const eenv & o, const eenv & i)
  {
    int nszs[o.nvr ()];
    const int * ovrs = o.vrs();
    const int * oszs = o.szs ();
    const int * iszs = i.szs ();
    int lv;
    for (int v = 0; v < o.nvr (); v++)
      if (i.hasvar (lv, ovrs[v]))
    nszs[v] = oszs[v] - iszs[lv] + 1;
      else
    nszs[v] = oszs[v];
    eenv tmp(o.nvr (), nszs, o.vrs ());
    //    std::cout << " tmp = " << tmp << std::endl;
    return tmp;
  };


// mpoldst Connection 

  struct vd
  {
    int n, d;
      vd ():n (-1), d (0)
    {
    };
  };

  inline std::ostream & operator<< (std::ostream & o, const vd & vd_)
  {
    o << "vd(" << vd_.n << ", " << vd_.d << ")";
    return o;
  };

  template < class MONOM > void mescan (int &last, vd * vdeg, const MONOM & m)
  {
    for (unsigned i = 0; i < m.size (); i++)
      {
    if (m[i] != 0)
      {
        if (vdeg[i].n != -1)
          vdeg[i].d = std::max ((int) m[i], vdeg[i].d);
        else
          {
        if (last != -1)
          {
            vdeg[i].n = vdeg[last].n;
            vdeg[i].d = m[i];
            vdeg[last].n = i;
            last = i;
          }
        else
          {
            last = i;
            vdeg[i].n = i;
            vdeg[i].d = m[i];
          };
          };
      };
      };
  }
  
  
  struct oulala { 
    int initial; int name; 
    bool operator<( const oulala& b ) const { return name < b.name; }; 
  };

  template < class C, class O >
  inline eenv::eenv (const sparse::monomial_seq< C, O > &mpol)
  {
    int nvr = nbvar (mpol);
    vd *lenv = new vd[nvr];
    int last = -1;
    typename sparse::monomial_seq < C, O >::const_iterator it;
    for (it = mpol.begin (); it != mpol.end (); mescan (last, lenv, *it++)) ;
    if ( last == -1 ) 
      { 
    nvr = 1; int szs = 1; 
    new (this) eenv(1,&szs,0);
    return; 
      };
    int c = last;
    nvr = 0;
    do
      {
    c = lenv[c].n;
    nvr++;
      }
    while (c != last);
    int vrs[nvr];
    int szs[nvr];
    c = last;

    int k = 0;
    do
      {
    vrs[k] = c;
    szs[k] = lenv[c].d + 1;
    c = lenv[c].n;
    lenv[vrs[k]].n = -1;
    lenv[vrs[k]].d = 0;
    k++;
      }
    while (c != last);
    new (this) eenv (nvr, szs, vrs );
  };

  inline bool eenv::equiv (const eenv & a, const eenv & b)
  {
    return equal(a,b);
  };


  inline eenv eenv::rewrite( const eenv& o, const eenv& i )
  {
    int vmap_[ o.nvr() ];
    if ( ! vmap( vmap_, o, i ) ) { std::cerr << "vmap error(rewrite)\n"; exit(1); };
    eenv tmp(o);
    for ( int v = 0; v < o.nvr(); v ++ ) tmp.szs()[v] = 1;
    for ( int v = 0; v < i.nvr(); v ++ ) tmp.szs()[vmap_[v]] = i.szs()[v];
    tmp.cstr();
    return tmp;
  };
  
}

} //namespace mmx
#endif