synaps/mpol/numexp.h

/***********************************************************************
History: Exponents for monomials are represented by element of E.
    Writen by Ph. Trebuchet (10/12/99).
$Id: numexp.h,v 1.1 2005/07/11 11:17:56 mourrain Exp $
**********************************************************************/
// A verifier le D (decallage).
#ifndef SYNAPS_MPOL_NUMEXP_H
#define SYNAPS_MPOL_NUMEXP_H


#include <iostream>
#include <synaps/init.h>
#include <synaps/mpol/Monom.h>

__BEGIN_NAMESPACE_SYNAPS


// Numeric exponent. To be used with caution for the moment.
static int mask=0;
template<char X, int D=2, class T=int> 
struct numexp {
  typedef int size_type;
  typedef T         exponent_t;
  typedef int         degree_t;
  typedef numexp<X,D,T> self_t;
  // static T mask;
  //Constructors
  numexp():ind(0) {}
  numexp(int t):ind(0) {}
  numexp(const self_t &A): ind(A.ind) {}
  numexp(int s, int *t)
    {
      for(int i=s-1;i>=0;i--)
        {
          ind+=t[i];
          ind=(ind<<D);
        }
    }

  //destructors // rien a faire ...
  int operator != (const numexp<X,D,T>&m)
  {
    return ind!=m.ind;
  }
  int operator != (const numexp<X,D,T>&m) const 
  {
    return ind!=m.ind;
  }
  exponent_t operator[](int i) const 
  {
    //  static T mask=0;
    if (mask==0)
      {
        for(int j=0;j<D;j++)
          mask=(mask<<1)+1;
      }
    return ((mask)&(ind>>(D*i)));
  }

  void reserve(int d) {};
  void resize(int d) {};
  void setExponent(int i,int d) 
    {
      T tmp1=(ind>>(D*i+D))<<(D*i+D);
      T tmp2=0;
       if(mask==0)
         {
           for(int j=0;j<D;j++)
             mask=(mask<<1)+1;
         }
      for(int j=0;j<i;j++) tmp2=(tmp2<<D)+mask;
      tmp2&=ind;
      ind=tmp1+tmp2+(((T)d)<<(D*i));
    }
  
  
  //operators
  self_t & operator =(const numexp<X,D,T> &A)
  {
      ind=A.ind; return *this;
  }

  self_t & operator += (const self_t & A)
  {
    ind+=A.ind; return *this;
  }
  self_t & operator *=(const self_t & A)
    {
      ind+=A.ind; return *this;
    }
  self_t & operator /=(const numexp<X,D,T> &A)
    {
      ind-=A.ind; return *this;
    }

  bool operator == (const numexp<X,D,T> &A) const
    {
      return (ind==A.ind);
    }
  
  //data
  T ind;
  friend void add(numexp<X,D,T> &n,const numexp<X,D,T> &m ){
    n.ind+=m.ind;
  };
};

template<char X,int D, class T> 
int Degree(const numexp<X,D,T> & m)
{
  int d=0;
  for(int i=0;i<lvar(m)+1;i++)  d+= m[i];
  return d;
}

template<char X,int D,class T> 
int GDegree(const numexp<X,D,T> &A,int i)
{
  return (3&(A.ind>>(D*i)));
}

template<char X,int D,class T> 
int SDegree(numexp<X,D,T> &A,int i,int d)
{
  int tmp1=(A.ind>>(D*i+D))<<(D*i+D);
  int tmp2=1;
  for(int j=0;j<i;j++) tmp2=(tmp2<<1)+1;
  tmp2&=A.ind;
  A.ind=tmp1+tmp2+(d<<(D*i));
  return A.ind;
}

template<char X,int D,class T> 
int lvar(const numexp<X,D,T> &A)
{
  int i=0;
  for(;(A.ind>>(D*i))>0;i++);
  return i-1;
}

template<char X,int D,class T> 
void erase(numexp<X,D,T> &A) {}

template<char X,int D,class T> 
void add(numexp<X,D,T> &r,const numexp<X,D,T> &A,const numexp<X,D,T> &B)
{
  r.ind=A.ind+B.ind;
}

template<char X,int D,class T> 
std::ostream & operator<<(std::ostream &os, const numexp<X,D,T> & rep)
{
  //  T mask=0;//3 tjrs pareil
  T tmp;
  int flag=0;
  if (mask==0){
    for(int j=0;j<D;j++) mask=(mask<<1)+1;
  }
  for(int i=0;i<lvar(rep)+1;i++)
     {
       if((rep.ind>>(D*i))&mask)
        {
          if (flag) 
            {
              os<<"*";
            }
          if (!flag) flag =1;
          os<<X<<i;
          if ((tmp=(rep.ind>>(D*i))&mask)>1)
            os<<'^'<<tmp;
        }
     }
  return os;
}

//----------------------------------------------------------------------//
template <class C,char X,int D, typename T> inline 
void assign(Monom<C, 
            numexp<X,D,T> > & a, 
            const OP<'*',Monom<C,numexp<X,D,T> >,Monom<C,numexp<X,D,T> > > & M)
{
  a.coeff() = (M.op1.coeff()*M.op2.coeff());
  a.rep().ind =  M.op1.rep().ind;
  a.rep()*=M.op2.rep();
}
//----------------------------------------------------------------------//
template <class R,char X,int D ,typename T> inline 
void assign(Monom<R,numexp<X,D,T> > & res, 
            const OP<'/',Monom<R,numexp<X,D,T> >,Monom<R,numexp<X,D,T> > > & M)
{
  if(M.op1.nvars()< M.op2.nvars())
    res= Monom<R,numexp<X,D,T> >(0);
  else {
    res= Monom<R,numexp<X,D,T> >(M.op1);
    res.coeff() /=M.op2.coeff();
    for (int i = 0; i <= M.op1.nvars() && res.coeff() !=0 ; ++i)
      if(M.op1[i] <M.op2[i])
        {
          res= Monom<R,numexp<X,D,T> >(0);
        }
      else
        res.rep().setExponent(i,M.op1[i] - M.op2[i]);
  }
}
//begining of partial specialisation of Monom


//template< char X, int D, typename T> inline
//void add(numexp<X,D,T> &n,numexp<X,D,T> &m)
//{
  //  n+=m;
//}
//template <class C, char X, int D, typename T> inline
//Monom<C,numexp<X,D,T> > & Monom<C,numexp<X,D,T> >::operator *= 
//           (const Monom<C,numexp<X,D,T> > & m) {
//      coeff() *= m.coeff();
//      if (coeff()!=0) 
//        rep()+=m.rep();
//      else
//        rep()=0;
//      return *this;
//}

__END_NAMESPACE_SYNAPS

#endif // SYNAPS_MPOL_NUMEXP_H