/Users/mourrain/Devel/mmx/realroot/include/realroot/tensor_monomials_fcts.hpp

Go to the documentation of this file.

#ifndef realroot_MPOLDSE_MONOMIAL_C
#define realroot_MPOLDSE_MONOMIAL_C

#include <realroot/monomial.hpp>
#include <realroot/tensor_monomials.hpp>
#include <realroot/tensor_convert.hpp>
#include <iostream>
namespace mmx { 

#ifndef MMX_WITH_PLUS_SIGN
#define MMX_WITH_PLUS_SIGN
  template<class X> inline bool with_plus_sign(const X& x) {return x>0;}

template<class C> inline bool with_plus_sign(const tensor::monomials<C>& x) {
  return true;
}

#endif  //MMX_WITH_PLUS_SIGN

namespace  tensor 
{
  template < class C > inline
  monomials < C >::monomials () {
    //   int szs = 1; new (this) monomials (1, &szs); data[0] = (C)0; 
    int *szs = new int, *vrs= new int;
    szs[0] = 1;
    vrs[0] = 0;
    env = eenv(1, szs, vrs);
    data = vector_type (env.sz(),(C)0);
    data[0] = (C)0; 
  }
  template < class C > inline
  monomials < C >::monomials (const eenv & e)//:env (e), data (e.sz ()) {}
  {
    env = e;
    data = vector_type (e.sz());
  };
  template < class C > inline
  monomials < C >::monomials (const eenv & e, const C & c)//:env(e),data (e.sz (), c) {};
  {
    env = e;
    data = vector_type (e.sz(),c);
  };
  template < class C > inline
  monomials < C >::monomials (int v, int d)
  {
    int *szs = new int, *vrs= new int;
    szs[0] = d+1;
    vrs[0] = v;
    env = eenv(1, szs, vrs);
    data = vector_type (env.sz(),(C)0);
    data[env.sz()-1]= 1;
  }

  template < class C > inline
  monomials < C >::monomials (const C& c, int d, unsigned v)
  {
    assert(d>=0);
    int *szs = new int, *vrs= new int;
    szs[0] = d+1;
    vrs[0] = v;
    env = eenv(1, szs, vrs);
    data = vector_type (env.sz(),(C)0);
    data[env.sz()-1]= c;
  }

  template < class C > inline
  monomials < C >::monomials (int nvr, const int *szs, const int *vrs) //:env (nvr, szs, vrs),data (env.sz(),(C)0)
  {
    env = eenv(nvr, szs, vrs);
    data = vector_type (env.sz(),(C)0);
  };
  template < class C > inline
  monomials < C >::monomials (const C & x)
  { int szs = 1; new (this) monomials (1, &szs); data[0] = x; }

  template < class C > 
  void realloc (monomials < C > &mpl, const eenv & nenv)
  {
    if (!eenv::equal (mpl.env, nenv ))
      {
        monomials< C > tmp (nenv);
        mpl.swap (tmp);
      };
  };

  template < class C >  inline
  void clear ( monomials < C >& monoms ) { std::fill(monoms.begin(),monoms.end(),C(0)); };
  template < class C >  inline
  void clear ( monomials < C > &mpl, const eenv & nenv)
  {
    if (mpl.env != nenv) {
      monomials < C > tmp (nenv, C (0));
      mpl.swap (tmp);
    } else 
      clear(mpl);
  };

  template < class C > 
  bool varindex (int &lv, const monomials<C>& mpl, int gv) 
  {
    return mpl.env.hasvar (lv, gv);
  }

  template < class C > bool monomials < C >::operator== (const monomials & mpl) const
  {
    if (!eenv::equal (env, mpl.env))
      return false;
    unsigned *oadd = new unsigned[esz ()];
    eenv::oaddress (env, oadd, mpl.env, 0, 0);
    for (int i = 0; i < esz (); i++)
      if (data[i] != mpl[oadd[i]])
        {
          delete[]oadd;
          return false;
        };
    delete[]oadd;
    return true;
  };

  template < class C >
  std::ostream & operator<< (std::ostream & o, const monomials < C > &mpl)
  {
    o << "monomials  object:\n";
    o << "\t(addresses) o: " << &mpl << " e: " << mpl.env.data <<  "d: " << mpl.
      begin () << std::endl;
    o << "\t(eenv)" << mpl.env << std::endl;
    o << "\t(data)";
    array::print (o, mpl.data);
    o << std::endl;
    o << ":object monomials\n";
    return o;
  };
  template < class Result, class Coeff, class Parameters > inline 
  void eval (Result & result, const monomials < Coeff > & monoms , const Parameters & parameters )
  { 
    //    if (monoms.nbvar()==1)
    vct::horner(result, monoms.begin (), monoms.size(), parameters, 1); 
    // levalm( result, monoms.begin (), monoms.env, parameters ); 
  };
  template < class Result, class Coeff, class Parameters > inline 
  void eval (Result & result, const monomials < Coeff > & monoms , const Parameters & parameters , unsigned n)
  { 
      //assert(parameters.size()<n);
    levalm( result, monoms.begin (), monoms.env, parameters ); 
  };
   
  template < class Result, class Coeff, class Parameters > inline 
  void heval (Result & result, const monomials < Coeff > & monoms , const Parameters & parameters)
  { 
    hevalm( result, monoms.begin (), monoms.env, parameters ); 
  };
  
  template < class Result, class Coeff, class Parameters > inline 
  void subs0 (Result *result, const monomials < Coeff > & monoms , const Parameters & parameters )
  { 
    subs0m( result, monoms.begin (), monoms.env, parameters ); 
  };

  template < class C > 
  void extend (monomials < C > &mpl, const eenv & nenv)
  {
    if ( mpl.env != nenv )
      {
        unsigned *oadd = new unsigned[mpl.esz ()];
        eenv::oaddress (nenv, oadd, mpl.env, 0, 0);
        monomials < C > tmp (nenv);
        vct::icopy (tmp.begin (), oadd, mpl.esz (), mpl.begin ());
        delete[]oadd;
        mpl.swap (tmp);
      };
  };

  // addition optimisable ...(mais pas maintenant)
  template < class X, class Y >
  void waddm (monomials < X > &mpla, const monomials < Y > &mplb)
  {
    assert (eenv::subset (mplb.env, mpla.env));
    unsigned *oadd = new unsigned[mplb.esz ()];
    eenv::oaddress (mpla.env, oadd, mplb.env, 0, 0);
    vct::ipadd (mpla.begin (), oadd, mplb.esz (), mplb.begin ());
    delete[]oadd;
  };

  template < class X, class C > inline void 
  add (monomials < X > &r, const monomials < X > &a, const C& c ) {
    add(r=a,c);
  }

  template < class C > inline void 
  add (monomials < C > &a, const C & x) {
    a[0] += x;
  }

  template < class C > void
  add (monomials < C > &r, const C & c, const monomials < C > &a) {
    add(r,a,c);
  }

  template < class C > inline
  void add (monomials < C > &r, const monomials < C > &a)
  {
    extend (r, eenv::common (r.env, a.env));
    waddm (r, a);
  }

  template < class C > inline
  void add (monomials < C > &r, const monomials < C > &a, const monomials < C > &b)
  {
    if ( &r == &a  ) {  add(r,b);  return;  };
    if ( &r  == &b ) {  add(r,a); return; };
    realloc(r,eenv::common(a.env,b.env));
    unsigned *oadd = new unsigned[std::max(a.esz (),b.esz())];
    eenv::oaddress (r.env, oadd,a.env, 0, 0);
    vct::icopy(r.begin(),oadd,a.esz(),a.begin());
    eenv::oaddress(r.env,oadd,b.env,0,0);
    vct::ipadd (r.begin (), oadd,b.esz (), b.begin ());
  };
  
  template < class C > inline
  void sub (monomials < C > &r, const monomials < C > &a, const monomials < C > &b)
  {
    if ( &r == &a  ) {  sub(r,b);  return;  };
    if ( &r == &b  ) {  monomials<C> tmp = b; sub(r,a,tmp); return; };
    
    realloc(r,eenv::common(a.env,b.env));
    unsigned *oadd = new unsigned[std::max(a.esz (),b.esz())];
    eenv::oaddress (r.env, oadd,a.env, 0, 0);
    vct::icopy(r.begin(),oadd,a.esz(),a.begin());
    eenv::oaddress(r.env,oadd,b.env,0,0);
    vct::ipsub (r.begin (), oadd,b.esz (), b.begin ());
  };

  template < class C > inline void
  sub (monomials < C > &mpl, const C & c, const monomials < C > &a) {
    add(mpl,a,-c);
  }
  template < class X, class Y >
  void wsubm (monomials < X > &mpla, const monomials < Y > &mplb)
  {
    assert (eenv::subset (mplb.env, mpla.env));
    unsigned *oadd = new unsigned[mplb.esz ()];
    eenv::oaddress (mpla.env, oadd, mplb.env, 0, 0);
    vct::ipsub (mpla.begin (), oadd, mplb.esz (), mplb.begin ());
    delete[]oadd;
  };

  template < class X, class Y > inline
  void sub (monomials < X > &r, const monomials < Y > &a)
  {
    extend (r, eenv::common (r.env, a.env)), wsubm (r, a);
  }


  template < class X, class C > inline void 
  sub (monomials < X > &r, const monomials < X > &a, const C& c )
  {
    r=a;sub(r,c);
  }
  template < class C , class X> inline 
  void sub (monomials < C > &a, const X & x)
  {
    a[0] -= (C)x;
  }

  template < class C > void mul (monomials < C > &r, const C & c)
  {
    vct::scmul (r.begin (), c, r.esz (), 1);
  };

  template < class C > void div (monomials < C > &r, const C & c)
  {
    vct::scdiv (r.begin (), c, r.esz (), 1);
  };

  template < class C > void mul (monomials < C > &r, const monomials < C > &a)
  {
    //    monomials<C> s(r);
    conv(r,a);
  };
  
  
  template < class X, class Y, class Z > inline
  void mul (monomials < X > &r, const monomials < Y > &a, const monomials < Z > &b)
  {
    conv(r,a,b);
  };

  template < class C > inline
  void mul (monomials < C > &a, const monomials < C > &b, const C & c)
  {
    mul (a = b, c);
  };

  template < class C > inline void
  mul (monomials < C > &r, const C &c, const monomials < C > &a) {
    mul(r=a,c);
  };

  template < class C > inline void
  div (monomials < C > &r, const monomials < C > &a, const C &c) {
    div(r=a,c);
  }

  template<class OSTREAM, class C> OSTREAM& 
  print( OSTREAM& os, const monomials<C>& mpl, const variables& Var = monom<C>::var) { 
    typedef monom<C>             monom_t;
    typedef typename monom<C>::container_t exponent_t;
    exponent_t tmp (mpl.env.mxvr()+1);
    const int *szs = mpl.szs ();
    const int *vrs = mpl.vrs ();
    bool notfirst = false;
    for (int i = 0; i < mpl.size(); i++)
      if (mpl[i] != 0)
        {
          int a, v;
          for (a = i, v = mpl.nvr () - 1; a; a /= szs[v], v--)
            tmp[vrs[v]] = a % szs[v];
          std::ostringstream sm;
          print(sm, monom_t (mpl[i], tmp),Var);
          //bool is_pol = is_polynomial(sm);
         if ( sm.str()[0]!='-' && sm.str()[0]!='+' && notfirst) { 
            os <<'+'; 
          }
          os<<sm.str();
          notfirst=true;
        };
    if(!notfirst) os << "0";
    return os;
  }

  template<class SYNTAX, class C> void
  print_flatten( SYNTAX& out, const monomials<C>& mpl, const variables& Var = monom<C>::var) { 
    typedef monom<C>             monom_t;
    typedef typename monom<C>::container_t exponent_t;
    exponent_t tmp (mpl.env.mxvr()+1);
    const int *szs = mpl.szs ();
    const int *vrs = mpl.vrs ();
    for (int i = 0; i < mpl.size(); i++)
      if (mpl[i] != 0)
        {
          int a, v;
          for (a = i, v = mpl.nvr () - 1; a; a /= szs[v], v--)
            tmp[vrs[v]] = a % szs[v];

          SYNTAX m = flatten(mpl[i]);
          for(int i=0;i<mpl.env.mxvr()+1;i++)
            m = m * pow(SYNTAX(Var[i].data()),flatten(tmp[i]));
          out+=m;             
        }
  }


  template<class Coeff> inline int
  degree(const monomials<Coeff>& p) {
    int d=0;
    for(int i=0;i<p.nbvar();i++) d+= *(p.env.szs()+i)-1;
    return d;
  }

  template<class Coeff> inline int
  degree(const monomials<Coeff>& p, int v) {
    return p.env.sz(v)-1;
  }


  template<class Coeff> inline int
  leading_variable(const monomials<Coeff>& p) {
    return 0;
  }

  template<class Coeff> inline void 
  diff( monomials<Coeff>& res, const monomials<Coeff>& src, int v ) {
    clear( res, eenv::diff(src.env,v));
    //    mdiff( res.begin(), res.env, src.begin(), src.env, v );
    mdiff( res.begin(), res.env, src.begin(), src.env,  v );
  }

  template < class C > 
  template < class X, class O>
  monomials < C >::monomials (const sparse::monomial_seq< X, O > &imp)
  {
    env= eenv(imp);
    data= vector_type(env.sz (), (C)0);
    //    const int *vr = this->vrs ();
    //    const int *st = this->str ();
    //    const int  nvr = this->nvr();
    //    C *data = this->begin ();
    mpolfill(this->begin(),imp,env);
  };
  
  template < class C, class O>
  void convert (monomials < C > &mpl, const sparse::monomial_seq< C, O> &imp)
  {
    mpl. ~ monomials < C > ();
    new (&mpl) monomials<C>(imp);
  };

  template < class C, class O> void 
  convert (sparse::monomial_seq< C, O> &pol, const monomials < C > &mpl) {
    typedef typename sparse::monomial_seq< C, O>::monom_t monom_t;
    typedef typename monom_t::container_t container_t;
    add(pol,C(0));
    container_t tmp (mpl.env.mxvr()+1);
    const int *szs = mpl.szs ();
    const int *vrs = mpl.vrs ();
    for (int i = 0; i < mpl.size(); i++)
      if (mpl[i] != 0)
        {
          int a, v;
          for (a = i, v = mpl.nvr () - 1; a; a /= szs[v], v--)
          tmp[vrs[v]] = a % szs[v];
          add(pol, monom_t (mpl[i], tmp));
        };
  }

// inline eenv face( const eenv& e, int lv ) 
//  {
//    const int * eszs = e.szs();
//    const int * evrs = e.vrs();
//    const int   envr = e.nvr();
//    int szs[ envr ];
//    int vrs[ envr ];
//    int c = 0;
//    for ( int v = 0; v < envr; v ++ ) if ( v != lv ) { szs[c] = eszs[v]; vrs[c++] = evrs[v]; };//evrs[v]; };
//    return eenv(c,szs,vrs);
//  };

inline eenv face_env( const eenv& e, int lv ) {
  const int * eszs = e.szs();
  const int * evrs = e.vrs();
  const int   envr = e.nvr();
  int szs[ envr ];
  int vrs[ envr ];
  int c = 0;
  for ( int v = 0; v < envr; v ++ ) if ( v != lv ) { szs[c] = eszs[v]; vrs[c++] = evrs[v]; };//evrs[v]; };
  return eenv(c,szs,vrs);
}

template<class C> inline
void lface ( monomials<C>& f, const monomials<C>& src, int v )
{
  lface( f.begin(), src.env, src.begin(), v );
};

template<class C> inline
void rface ( monomials<C>& f, const monomials<C>& src, int v )
{
  rface( f.begin(), src.env, src.begin(), v );
};


template<class C> inline
void face ( monomials<C>& f, const monomials<C>& src, int v, int n )
{
  assert(n==0||n==1);
  eenv fenv(face_env(src.env,src.env.localv(v)));
  if ( f.env != fenv ) realloc( f, fenv );
  if ( n ) rface( f, src, src.env.localv(v)  );
  else     lface( f, src, src.env.localv(v)  );
};
 

  template<class C> inline
  void islice( monomials<C>& f, const monomials<C>& src, int v, int i )
  {
    slice( f.begin(), src.env, src.begin(), v );
  }

  template<class C> inline
  void slice( monomials<C>& f, const monomials<C>& src, int v, int n )
  {
    assert(n==0||n==1);
    eenv fenv(slice(src.env,src.env.localv(v)));
    if ( f.env != fenv ) realloc( f, fenv );
    islice( f, src, src.env.localv(v),n);
  }

    template<class C> inline
    const C monomials<C>::entry(std::vector<int> deg)
    {
        int * vvrs = vrs();
        int * sstr = str();
        int * sszs = szs() ;

        int i,pos=0;
        for (i = 0; i < nvr() ; ++i)
        {
            std::cout<<"entry:" <<deg[i] << std::endl;
            if (deg[i] >= sszs[ vvrs[i] ] || deg[i]<0 ) return C(0); 
            pos += deg[i]*sstr[ vvrs[i] ];
        }
        return data[pos];
    }

   template<class C> //inline
   void shift( monomials<C> & f, const C & t, const int & v )
   {
       std::vector<int> ind = std::vector<int>( f.nvr() );
       int s,k,j, i,pos =0;

        int * vars = f.vrs();
        int * sstr = f.str();
        int * sszs = f.szs() ;
       
       for (;;) 
       {
           for ( s = sszs[v], j = 0; j <= s-2; j++ )
               for( k = s-2; k >= j; k-- ) 
                   f.data[sstr[v]*k+ pos ] += 
                   t*f.data[ sstr[v]*(1+k)+ pos ];

           /* next row */
           for (i = 0; i < f.nvr() ; ++i) 
           {
               if ( vars[i] != v )
               {
                   ind[i] += 1;
                   pos += sstr[i];
               } else continue;

               if (ind[i] < sszs[i])
                   break;
               ind[i] = 0;
               pos -= sszs[i]*sstr[i];
           }
           if (i == f.nvr() )
               break;
       }
   }

    template<class C>
    void set_variable(monomials<C> & f, C t, int v )
    {
       std::vector<int> ind = std::vector<int>( f.nvr() );
       int k, i, pos =0;

        int * vars = f.vrs();
        int * sstr = f.str();
        int * sszs = f.szs() ;
        C tmp;
       
       for (;;) 
       {
           tmp=0;
           for( k = sszs[v]-1; k !=-1; --k )
           { 
               tmp *= t;
               tmp += f.data[sstr[v]*k + pos];
               f.data[sstr[v]*k + pos] = 0 ;
           }
           f.data[pos]=tmp;
           
           /* next row */
           for (i = 0; i < f.nvr() ; ++i) 
           {
               if ( vars[i] != v )
               {
                   ind[i] += 1;
                   pos += sstr[i];
               } else continue;

               if (ind[i] < sszs[i])
                   break;
               ind[i] = 0;
               pos -= sszs[i]*sstr[i];
           }
           if (i == f.nvr() )
               break;
       }
   }


    template<class C> //inline
    void reciprocal( monomials<C> & f, const int & v)
   {
       std::vector<int> ind = std::vector<int>( f.nvr() );
       int s,i,pos =0;

        int * vars = f.vrs();
        int * sstr = f.str();
        int * sszs = f.szs();
        int l= sszs[v]-1;
       
       for (;;) 
       {
           for ( s = 0 ; s <= l/2; s++ )
               std::swap(f.data[sstr[v]*s+pos], f.data[sstr[v]*(l-s)+pos] );

           // next row
           for (i = 0; i < f.nvr() ; ++i) 
           {
               if ( vars[i] != v )
               {
                   ind[i] += 1;
                   pos += sstr[i];
               } else continue;

               if (ind[i] < sszs[i])
                   break;
               ind[i] = 0;
               pos -= sszs[i]*sstr[i];
           }
           if (i == f.nvr() )
               break;
       }
   }

    template<class C> //inline
    void contraction( monomials<C> & f, const C & t, const int & v )
    {
        std::vector<int> ind = std::vector<int>( f.nvr() );
        int s,i,pos =0;
        
        int * vars = f.vrs();
        int * sstr = f.str();
        int * sszs = f.szs();
        
       for (;;) 
       {
           for ( s = 0 ; s < sszs[v]; ++s )
               f.data[sstr[v]*s+pos] *= pow(t,s) ;

           // next row
           for (i = 0; i < f.nvr() ; ++i) 
           {
               if ( vars[i] != v )
               {
                   ind[i] += 1;
                   pos += sstr[i];
               } else continue;

               if (ind[i] < sszs[i])
                   break;
               ind[i] = 0;
               pos -= sszs[i]*sstr[i];
           }
           if (i == f.nvr() )
               break;
       }
   };

    template<class C>
     void mins(  monomials<C> & mm,  monomials<C> & f, int v )
     {
        int p,k,i,j;
        int * sszs = f.szs() ;
        int * sstr = f.str() ;

        for ( p = 0, k = 0; k < sszs[v]; k++, p += sstr[v] )
        {
            mm.data[k] = f.data[p];
          for ( i = 0; i < f.size(); i += sstr[v-1] )
            for ( j = i; j < i + sstr[v]; j ++ )
              if ( f.data[j+p] < mm[k] ) mm[k] = f.data[j+p];
        }
     }

    template<class C>
    void maxs(  monomials<C> & mm,  monomials<C> & f, int v )
    {

       int p,k,i,j;
       int * sszs = f.szs() ;
       int * sstr = f.str() ;

       for ( p = 0, k = 0; k < sszs[v]; k++, p += sstr[v] )
        {
            mm.data[k] = f.data[p];
          for ( i = 0; i < f.size(); i += sstr[v-1] )
            for ( j = i; j < i + sstr[v]; j ++ )
              if ( f.data[j+p] > mm[k] ) mm[k] = f.data[j+p];
        }
    }
    
    //rename variable v to n
    template<class C>
    void rename_var(monomials<C> & f, const int & v, const int & n )
    {
        int * vr = f.vrs();
        
        for (int i=0; i< f.nvr(); ++i)
            if (vr[i]==v) { 
                vr[i]=n;
                break; }
    }

};//namespace tensor
} //namespace mmx
#endif  /* realroot_MPOLDSE_MONOMIAL_C */