/* -*- mode:C++  -*-  */
/*
 *  Copyright (C) 2000,2014 B. Parisse, Institut Fourier, 38402 St Martin d'Heres
 *
 *  This program 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.
 *
 *  This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
 */
#ifndef _GIAC_POLY_H_
#define _GIAC_POLY_H_
#include "first.h"
//#include <fstream>
#include "vector.h"
#include <string>
#include "fraction.h"
#include "index.h"
#include "monomial.h"
#include "threaded.h"
#include <algorithm>

#ifndef NO_NAMESPACE_GIAC
namespace giac {
#endif // ndef NO_NAMESPACE_GIAC

  // the tensor class to represent polynomial
  template <class T> class tensor{
  public:
    // members
    int dim; // number of indices, this is the size of monomial.index
    std::vector< monomial<T> > coord; // sorted list of monomials
    // T zero;
    // functional object sorting function for monomial ordering
    bool (* is_strictly_greater)( const index_m &, const index_m &);
    std::pointer_to_binary_function < const monomial<T> &, const monomial<T> &, bool> m_is_strictly_greater ;
    // constructors
    tensor(const tensor<T> & t) : dim(t.dim), coord(t.coord), is_strictly_greater(t.is_strictly_greater), m_is_strictly_greater(t.m_is_strictly_greater) { }
    tensor(const tensor<T> & t, const std::vector< monomial<T> > & v) : dim(t.dim), coord(v), is_strictly_greater(t.is_strictly_greater), m_is_strictly_greater(t.m_is_strictly_greater) { }
    // warning: this constructor prohibits construction of tensor from a value
    // of type T if this value is an int, except by using tensor<T>(T(int))
    tensor() : dim(0), is_strictly_greater(i_lex_is_strictly_greater), m_is_strictly_greater(std::ptr_fun<const monomial<T> &, const monomial<T> &, bool>(m_lex_is_strictly_greater<T>)) { }
    explicit tensor(int d) : dim(d), is_strictly_greater(i_lex_is_strictly_greater), m_is_strictly_greater(std::ptr_fun<const monomial<T> &, const monomial<T> &, bool>(m_lex_is_strictly_greater<T>)) { }
    explicit tensor(int d,const tensor<T> & t) : dim(d),is_strictly_greater(t.is_strictly_greater), m_is_strictly_greater(t.m_is_strictly_greater)  { }
    tensor(const monomial<T> & v) : dim(int(v.index.size())), is_strictly_greater(i_lex_is_strictly_greater), m_is_strictly_greater(std::ptr_fun<const monomial<T> &, const monomial<T> &, bool>(m_lex_is_strictly_greater<T>)) { 
      coord.push_back(v);
    }
    tensor(const T & v, int d) : dim(d), is_strictly_greater(i_lex_is_strictly_greater), m_is_strictly_greater(std::ptr_fun<const monomial<T> &, const monomial<T> &, bool>(m_lex_is_strictly_greater<T>)) {
      if (!is_zero(v))
	coord.push_back(monomial<T>(v,0,d));
    }
    tensor(int d,const std::vector< monomial<T> > & c) : dim(d), coord(c), is_strictly_greater(i_lex_is_strictly_greater),m_is_strictly_greater(std::ptr_fun<const monomial<T> &, const monomial<T> &, bool>(m_lex_is_strictly_greater<T>)) { }
    ~tensor() { coord.clear(); }
    // member functions
    // ordering monomials in the tensor
    void tsort(){ 
#if 1 // def NSPIRE
      sort_helper<T> M(m_is_strictly_greater);
      sort(coord.begin(),coord.end(),M); 
#else
      sort(coord.begin(),coord.end(),m_is_strictly_greater); 
#endif
    }
    int lexsorted_degree() const{ 
      if (!dim)
	return 0;
      if (coord.empty())
	return 0;
      else
	return coord.front().index.front(); 
    }
    int degree(int n) const ;
    int valuation(int n) const ;
    index_t degree() const ;
    int total_degree() const ;
    int partial_degree(int nvars) const ; // total degree wrt to vars 0..nvars-1
    void reverse() ; // reverse variable ordering
    void append(const tensor<T> &);
    void Tcoeffs(std::vector< tensor<T> > & v) const;
    std::vector< tensor<T> > Tcoeffs() const;
    tensor<T> coeff(int deg) const;
    tensor<T> multiplydegrees(int d) const ;
    tensor<T> dividedegrees(int d) const ;
    tensor<T> dividealldegrees(int d) const ;
    tensor<T> homogeneize() const;
    void reorder(const std::vector<int> & permutation) ;
    // shift and multiply, shift and divide, shift only
    tensor<T> shift (const index_m & ishift,const T & fois) const ;
    tensor<T> shift (const T & fois,const index_m & ishift) const ;
    tensor<T> shift (const index_m & ishift) const ;
    // tensor<T> operator + (const T & other) const ;
    void TAdd(const tensor<T> &other,tensor<T> & result) const;
    // tensor<T> operator - (const tensor<T> & other) const ;
    void TSub(const tensor<T> &other,tensor<T> & result) const;
    // tensor<T> operator * (const tensor<T> & other) const ;
    // tensor<T> operator * (const T & fact) const ;
    tensor<T> & operator *= (const T & fact) ;
    // tensor<T> operator - () const ;
    // tensor<T> operator / (const tensor<T> & other) const ;
    // tensor<T> operator / (const T & fact) const ;
    tensor<T> & operator /= (const T & fact) ;
    // tensor<T> operator % (const tensor<T> & other) const ;
    bool TDivRem (const tensor<T> & other, tensor<T> & quo, tensor<T> & rem, bool allowrational = true ) const ; // this=quo*other+rem
    bool TDivRemHash(const tensor<T> & b,tensor<T> & quo,tensor<T> & r,bool allowrational=false,int exactquo=0,double qmax=0.0) const ; // same as TDivRem but allowrationnal=false *and* poly with 1 main variable
    bool TDivRem1(const tensor<T> & b,tensor<T> & quo,tensor<T> & r,bool allowrational=false,int exactquo=0) const ; 
    bool Texactquotient (const tensor<T> & other, tensor<T> & quo,bool allowrational=false ) const ; // this=quo*other+rem, rem must be 0
    bool TPseudoDivRem (const tensor<T> & other, tensor<T> & quo, tensor<T> & rem, tensor<T> & a) const ; // a*this=quo*other+rem
    // bool TDivRem (const T & x0, tensor<T> & quo, tensor<T> & rem) const ;
    tensor<T> operator () (const T & x0 ) const;
    // T operator () (const std::vector<T> & x0) const;
    T constant_term() const ;
    T norm() const;
    tensor<T> derivative() const ;
    tensor<T> integrate() const ;
    // insertion of a monomial and reordering
    void insert_monomial(const monomial<T> & c);
    // position corresponding to an index, return -1 if not found
    int position(const index_m & v) const;
    index_m vector_int(int position) const;
    T & operator () ( const index_m & v) ;
    const T & operator () ( const index_m & v) const ;
    tensor<T> untrunc1 (int j=0) const {
      std::vector< monomial<T> > v;
      Untrunc1(this->coord,j,v);
      return tensor<T>(dim+1,v);
    }
    void untruncn (int j=0) {
      Untruncn(this->coord,j);
      ++dim;
    }
    tensor<T> untrunc (int j,int dim) const {
      std::vector< monomial<T> > v;
      Untrunc(this->coord,j,dim,v);
      return tensor<T>(dim,v);
    }
    tensor<T> trunc1 () const {
      assert(dim);
      std::vector< monomial<T> > v;
      Trunc1(this->coord,v);
      return tensor<T>(dim-1,v);
    }
    int gcddeg(int k) const;
    index_t gcddeg() const;
    // printing
    std::string print() const {
      if (coord.empty())
        return "";
      std::string s;
      typename std::vector< monomial<T> > :: const_iterator it=coord.begin(),itend=coord.end();
      for (;;){
        s += it->print();
        ++it;
        if (it==itend)
	  return s;
        s +=  '+' ;
      }
    };
    const char * dbgprint() const {
#if 0
      static std::string s;
      s=print();
      return s.c_str();
#else
      static std::string * sptr=0;
      if (!sptr) sptr=new std::string;
      *sptr=print();
      return sptr->c_str();
#endif
    }
    void high_order_degree_truncate(int n){
      // suppress terms of order >= n
      typename std::vector< monomial<T> >::iterator it=coord.begin(),itend=coord.end();
      for (;it!=itend;++it){
	if (it->index.front()<n)
	  break;
      }
      if (it!=coord.begin() && it!=itend)
	coord.erase(coord.begin(),it);
    }
    tensor<T> total_degree_truncate(int n) const {
      tensor<T> res(dim);
      // suppress terms of total degree > n
      typename std::vector< monomial<T> >::const_iterator it=coord.begin(),itend=coord.end();
      for (;it!=itend;++it){
	if (sum_degree(it->index)<=n)
	  res.coord.push_back(*it);
      }
      return res;
    }
  };
  template <class T> class Tref_tensor{
  public:
    ref_count_t ref_count;
    tensor<T> t;
    Tref_tensor<T>(const tensor<T> & P): ref_count(1),t(P) {}
    Tref_tensor<T>(int dim): ref_count(1),t(dim) {}
  };

  // convert p to monomial represented by unsigned integers
  // using the rule [a1 a2 .. an] [deg1 ... degn ] -> 
  // (... (a1*deg2+a2)*deg3 +...)*degn+an
  template<class T,class U>
  void convert(const tensor<T> & p,const index_t & deg,std::vector< T_unsigned<T,U> >  & v){
    typename std::vector< monomial<T> >::const_iterator it=p.coord.begin(),itend=p.coord.end(),itstop;
    v.clear();
    v.reserve(itend-it);
    T_unsigned<T,U> gu;
    U u;
    int nterms;
    index_t::const_iterator itit,itcur,idxcur,idxend,ditbeg=deg.begin(),ditend=deg.end(),dit;
    for (;it!=itend;++it){
      u=0;
      itcur=itit=it->index.begin();
      for (dit=ditbeg;dit!=ditend;++itit,++dit)
	u=u*unsigned(*dit)+unsigned(*itit);
      gu.u=u;
      gu.g=it->value;
      v.push_back(gu);
      // dense poly check 
      --itit;
      nterms=*itit;
      if (nterms<2 || nterms>=itend-it)
	continue;
      itstop=it+nterms;
      if (itstop->index.back())
	continue;
      for (idxcur=itstop->index.begin(),idxend=itstop->index.end()-1;idxcur!=idxend;++idxcur,++itcur){
	if (*idxcur!=*itcur)
	  break;
      }
      if (idxcur!=idxend)
	continue;
      // this part is dense
      for (;it!=itstop;){
	++it;
	gu.g=it->value;
	--gu.u;
	v.push_back(gu);
      }
    }
  }

  template<class T,class U>
  void convert(const std::vector< T_unsigned<T,U> > & v,const index_t & deg,tensor<T> & p){
    typename std::vector< T_unsigned<T,U> >::const_iterator it=v.begin(),itend=v.end();
    index_t::const_reverse_iterator ditbeg=deg.rbegin(),ditend=deg.rend(),dit;
    p.dim=ditend-ditbeg;
    p.coord.clear();
    p.coord.reserve(itend-it);
    U u;
    index_t i(p.dim);
    int k;
    for (;it!=itend;++it){
      u=it->u;
      for (k=p.dim-1,dit=ditbeg;dit!=ditend;++dit,--k){
	i[k]=u % unsigned(*dit);
	u = u/unsigned(*dit);
      }
      p.coord.push_back(monomial<T>(it->g,i));
    }
  }

  template <class T>
  bool operator == (const tensor<T> & p,const tensor<T> & q){
    return p.dim==q.dim && p.coord.size()==q.coord.size() && p.coord==q.coord;
  }

  template <class T> bool tensor_is_strictly_greater(const tensor<T> & p,const tensor<T> & q){
    if (q.coord.empty())
      return true;
    if (p.coord.empty())
      return false;
    return p.m_is_strictly_greater(p.coord.front(),q.coord.front());
  }

  template <class T>
  void tensor<T>::insert_monomial(const monomial<T> & c){
    coord.push_back(c);
    this->tsort(); // sort(coord.begin(),coord.end(),m_is_strictly_greater);
  }

  template <class T>
  int tensor<T>::degree(int n) const {
    typename std::vector< monomial<T> >::const_iterator it=this->coord.begin();
    typename std::vector< monomial<T> >::const_iterator it_end=this->coord.end();
    int res=0;
    for (;it!=it_end;++it){
      int temp=(it->index)[n];
      if (res<temp)
	res=temp;
    }
    return res;
  }

  template <class T>
  int tensor<T>::total_degree() const {
    typename std::vector< monomial<T> >::const_iterator it=this->coord.begin();
    typename std::vector< monomial<T> >::const_iterator it_end=this->coord.end();
    int res=0;
    for (;it!=it_end;++it){
      int temp=sum_degree(it->index);
      if (res<temp)
	res=temp;
    }
    return res;
  }

  template <class T>
  tensor<T> tensor<T>::homogeneize() const {
    int d=total_degree();
    tensor<T> res(*this);
    typename std::vector< monomial<T> >::iterator it=res.coord.begin(),itend=res.coord.end();
    for (;it!=itend;++it){
      index_t i=it->index.iref();
      i.push_back(d-sum_degree(i));
      it->index=i;
    }
    return res;
  }

  template <class T>
  int tensor<T>::partial_degree(int vars) const {
    typename std::vector< monomial<T> >::const_iterator it=this->coord.begin();
    typename std::vector< monomial<T> >::const_iterator it_end=this->coord.end();
    int res=0;
    for (;it!=it_end;++it){
      int temp=0;
      index_t::const_iterator jt=it->index.begin(),jtend=jt+vars;
      for (;jt!=jtend;++jt){
	temp += *jt;
      }
      if (res<temp)
	res=temp;
    }
    return res;
  }

  template <class T>
  int tensor<T>::valuation(int n) const {
    typename std::vector< monomial<T> >::const_iterator it=this->coord.begin();
    typename std::vector< monomial<T> >::const_iterator it_end=this->coord.end();
    if (it==it_end)
      return 0;
    int res=(it->index)[n];
    for (;it!=it_end;++it){
      int temp=(it->index)[n];
      if (res>temp)
	res=temp;
    }
    return res;
  }

  template <class T>
  index_t tensor<T>::degree() const {
    typename std::vector< monomial<T> >::const_iterator it=this->coord.begin(),it2;
    typename std::vector< monomial<T> >::const_iterator it_end=this->coord.end();
    index_t res(dim);
    if (!dim) return res;
    index_t::iterator itresbeg=res.begin(),itresend=res.end(),itres;
    index_t::const_iterator ittemp,ittemp2,ittempend;
    if (// false &&
	is_strictly_greater==i_lex_is_strictly_greater){
      for (;it!=it_end;++it){
	ittemp=it->index.begin();
	for (itres=itresbeg;itres!=itresend;++itres,++ittemp){
	  if (*itres<*ittemp)
	    *itres=*ittemp;
	}
	// check if the polynomial is dense -> skip ?
	--ittemp; // point to xn power
	if (*ittemp<3 || *ittemp>=it_end-it)
	  continue;
	it2=it+(*ittemp); // if dense, point to the last monomial with same x1..xn-1 
	if (it2->index.back()) // last power in xn must be 0
	  continue;
	ittemp=it->index.begin();
	ittemp2=it2->index.begin();
	ittempend=ittemp+dim-1; // check all other powers
	for (;ittemp!=ittempend;++ittemp2,++ittemp){
	  if (*ittemp!=*ittemp2)
	    break;
	}
	if (ittemp!=ittempend)
	  continue;
	it=it2;
      }
    }
    else {
      for (;it!=it_end;++it){
	ittemp=it->index.begin();
	for (itres=itresbeg;itres!=itresend;++itres,++ittemp){
	  if (*itres<*ittemp)
	    *itres=*ittemp;
	}
      }
    }
    return res;
  }

  template <class T>
  void tensor<T>::append(const tensor<T> & p){
    if (p.coord.empty())
      return;
    if (this->coord.empty()){
      this->dim=p.dim;
      this->coord=p.coord;
      return;
    }
    if (is_strictly_greater(this->coord.back().index , p.coord.front().index)){
      this->coord.reserve(this->coord.size()+p.coord.size());
      typename std::vector< monomial<T> >::const_iterator it=p.coord.begin();
      typename std::vector< monomial<T> >::const_iterator it_end=p.coord.end();
      for (;it!=it_end;++it)
	this->coord.push_back(*it);
    }
    else
      TAdd(p,*this); // *this=*this+p;
  }

  template <class T>
  tensor<T> tensor<T>::multiplydegrees(int d) const {
    tensor<T> res(dim);
    typename std::vector< monomial<T> >::const_iterator it=coord.begin(),it_end=coord.end();
    for (;it!=it_end;++it){
      index_t i=it->index.iref();
      i.front() *= d;
      res.coord.push_back(monomial<T>(it->value,i));
    }
    return res;
  }

  template <class T>
  tensor<T> tensor<T>::dividedegrees(int d) const {
    tensor<T> res(dim);
    typename std::vector< monomial<T> >::const_iterator it=coord.begin(),it_end=coord.end();
    for (;it!=it_end;++it){
      index_t i=it->index.iref();
      i.front() /= d;
      res.coord.push_back(monomial<T>(it->value,i));
    }
    return res;
  }

  template <class T>
  tensor<T> tensor<T>::dividealldegrees(int d) const {
    tensor<T> res(dim);
    typename std::vector< monomial<T> >::const_iterator it=coord.begin(),it_end=coord.end();
    for (;it!=it_end;++it){
      index_t i=it->index.iref();
      i = i/d;
      res.coord.push_back(monomial<T>(it->value,i));
    }
    return res;
  }

  template <class T> 
  int tensor<T>::position(const index_m & v) const {
    int smax=int(coord.size())-1;
    int smin=0;
    int s;
    for (;smin<smax;){
      s=(smax+smin)/2; // smin <= s < smax
      index_m vs=coord[s].index;
      if (v==vs)
	break;
      if (is_strictly_greater(v,vs)) // if v > v[s] must start above smin+1
	smax=s-1; // same
      else 
	smin=s+1; // keeps smin <=smax
    }
    s=(smax+smin)/2; // if loop breaked return the correct s, else smin=smax
    index_m vs=coord[s].index;
    if (v==vs)
      return(s);
    else
      return(-1);
  }

  template <class T>  
  index_m tensor<T>::vector_int(int position) const {
    return(coord[position].index);
  }

  template <class T>
  T & tensor<T>::operator () ( const index_m & v) {
    int p=position(v);
    if (p!=-1)
      return coord[p].value;
    coord.push_back(T(0),v);
    p=position(v);
    return coord[p].value;
  }

  template<class T>  
  const T & tensor<T>::operator () ( const index_m & v) const{
    static T const myzero(0);
    int p=position(v);
    if (p==-1) {
      return myzero;
    }
    else
      return (coord[p].value);
  }

  template <class T>
  std::ostream & operator << (std::ostream & os, const tensor<T> & t)
  {
    return os << t.print();
  }


  template <class T>
  void lexsort(std::vector < monomial<T> > & v){
#if 1 // def NSPIRE
    sort_helper<T> M(std::ptr_fun<const monomial<T> &, const monomial<T> &, bool>(m_lex_is_strictly_greater<T>));
    sort(v.begin(),v.end(),M);
#else
    sort(v.begin(),v.end(),std::ptr_fun<const monomial<T> &, const monomial<T> &, bool>(m_lex_is_strictly_greater<T>));
#endif
  }



  template <class T>
  int tensor<T>::gcddeg(int k) const {
    assert(k<dim);
    int res=0;
    typename std::vector< monomial<T> >::const_iterator it=coord.begin(),itend=coord.end();
    for (;it!=itend;++it){
      res=mygcd((it->index)[k],res);
      if (res==1)
	break;
    }
    return res;
  }

  template <class T>
  index_t tensor<T>::gcddeg() const {
    typename std::vector< monomial<T> >::const_iterator it=coord.begin(),itend=coord.end();
    assert(itend!=it);
    index_t res(it->index.iref());
    index_t zero(res.size());
    for (;it!=itend;++it){
      res=index_gcd(it->index.iref(),res);
      if (res==zero)
	break;
    }
    return res;
  }

  // univariate Horner evaluation with remainder, rem is the evaluation result
  // as a n-1-dimensional tensor and quo the quotient as a n-dim tensor
  /*
    template <class T>
    bool tensor<T>::TDivRem (const T & x0, tensor<T> & quo, tensor<T> & rem) const
    {
    if (coord.empty()){
    rem=*this;
    quo=*this;
    return true;
    }
    std::vector< monomial<T> > horner_coord(coord);
    tensor<T> add_rem(dim-1), add_quo(dim-1);
    lexsort(horner_coord);
    tensor<T> quotient(*this,new_seq);
    rem.coord.clear();
    quo.coord.clear();
    typename std::vector< monomial<T> >::const_iterator it = horner_coord.begin();
    index_m pui=it->index;
    for (;it!=horner_coord.end();++it){
    if (pui.front()==it->index.front()) { 
    // same external power, add
    add_rem.coord.push_back(it->trunc1());
    add_quo.coord.push_back(*it);
    }
    else {      // different power do an Horner *
    rem=(rem+add_rem)*pow(x0,pui.front()-it->index.front());
    add_rem.coord.clear();
    add_rem.coord.push_back(it->trunc1());
    for (;pui.front()> it->index.front();){
    pui[0]--;
    add_quo.divided_by_x();
    add_quo=add_quo*x0;
    quo=quo+add_quo;
    }
    add_quo.coord.clear();
    pui=it->index;
    }
    }
    rem=(rem+add_rem)*pow(x,pui.front());
    for (;pui.front();){
    pui[0]--;
    add_quo.divided_by_x();
    add_quo=add_quo*x0;
    quo=quo+add_quo;
    }
    // sort remainder and quotient
    rem.tsort();
    quo.tsort();
    return true;
    }
  */

  // univariate Horner evaluation with remainder, rem is the evaluation result
  // as a n-1-dimensional tensor, the quotient is not computed
  template <class T>
  tensor<T> tensor<T>::operator () (const T & x0 ) const
  {
    if (coord.empty())
      return *this;
    if (!dim){
      return *this;
    }
    std::vector< monomial<T> > horner_coord(coord);
    tensor<T> rem(dim-1),add_rem(dim-1) ;
    lexsort(horner_coord);
    typename std::vector< monomial<T> >::const_iterator it = horner_coord.begin();
    index_m pui=(*it).index;
    for (;it!=horner_coord.end();++it){
      if (pui.front()==it->index.front()) { 
	// same external power, add
	add_rem.coord.push_back(it->trunc1());
      }
      else {      // different power do an Horner *
#if !defined NSPIRE && !defined(FXCG)// GIAC_VECTOR
	rem.TAdd(add_rem,rem); rem *= pow(x0,pui.front()-it->index.front()); 
#else
	rem =(add_rem+rem)*pow(x0,pui.front()-it->index.front());
#endif
	add_rem.coord.clear();
	add_rem.coord.push_back(it->trunc1());
	pui=it->index;
      }
    }
    rem.TAdd(add_rem,rem); // rem=(add_rem+rem);
    if (pui.front()){
#if !defined (NSPIRE) && !defined(FXCG) // GIAC_VECTOR
      rem *= pow(x0,pui.front());
#else
      rem = rem*pow(x0,pui.front());
#endif
    }
    rem.tsort();
    return rem;
  }


  template <class T>
  void tensor<T>::TAdd (const tensor<T> & other,tensor<T> & result) const {  
    typename std::vector< monomial<T> >::const_iterator a=coord.begin();
    typename std::vector< monomial<T> >::const_iterator a_end=coord.end();
    if (a == a_end) {
      result=other;
      return ;
    }
    typename std::vector< monomial<T> >::const_iterator b=other.coord.begin();
    typename std::vector< monomial<T> >::const_iterator b_end=other.coord.end();
    if (b==b_end){
      result= *this;
      return ;
    }
    Add(a,a_end,b,b_end,result.coord,is_strictly_greater);
  }

  /*
    template <class T>
    tensor<T> tensor<T>::operator + (const tensor<T> & other) const {  
    // Tensor addition
    typename std::vector< monomial<T> >::const_iterator a=coord.begin();
    typename std::vector< monomial<T> >::const_iterator a_end=coord.end();
    if (a == a_end) {
    return other;
    }
    typename std::vector< monomial<T> >::const_iterator b=other.coord.begin();
    typename std::vector< monomial<T> >::const_iterator b_end=other.coord.end();
    if (b==b_end){
    return *this;
    }
    std::vector< monomial<T> > new_coord;
    Add(a,a_end,b,b_end,new_coord,is_strictly_greater);
    return tensor<T>(*this,new_coord);
    }

  */

  template <class T>
  void tensor<T>::TSub (const tensor<T> & other,tensor<T> & result) const {
    typename std::vector< monomial<T> >::const_iterator a=coord.begin();
    typename std::vector< monomial<T> >::const_iterator a_end=coord.end();
    typename std::vector< monomial<T> >::const_iterator b=other.coord.begin();
    typename std::vector< monomial<T> >::const_iterator b_end=other.coord.end();
    if (b == b_end) {
      result= *this;
      return;
    }
    Sub(a,a_end,b,b_end,result.coord,is_strictly_greater);
  }

  /*
    template <class T>
    tensor<T> tensor<T>::operator - (const tensor<T> & other) const {  
    // Tensor addition
    typename std::vector< monomial<T> >::const_iterator a=coord.begin();
    typename std::vector< monomial<T> >::const_iterator a_end=coord.end();
    typename std::vector< monomial<T> >::const_iterator b=other.coord.begin();
    typename std::vector< monomial<T> >::const_iterator b_end=other.coord.end();
    if (b == b_end) {
    return *this;
    }
    std::vector< monomial<T> > new_coord;
    Sub(a,a_end,b,b_end,new_coord,is_strictly_greater);
    return tensor<T>(*this,new_coord);
    }

    template <class T>
    tensor<T> tensor<T>::operator - () const {  
    // Tensor addition
    std::vector< monomial<T> > new_coord;
    typename std::vector< monomial<T> >::const_iterator a = coord.begin();
    typename std::vector< monomial<T> >::const_iterator a_end = coord.end();
    new_coord.reserve(((int) a_end - (int) a )/(sizeof(monomial<T>)));  
    for (;a!=a_end;++a){
    new_coord.push_back(monomial<T>(-(*a).value,(*a).index));
    }
    return tensor<T>(*this,new_coord);
    }

    template <class T>
    tensor<T> tensor<T>::operator * (const tensor<T> & other) const {  
    // Multiplication
    typename std::vector< monomial<T> >::const_iterator ita = coord.begin();
    typename std::vector< monomial<T> >::const_iterator ita_end = coord.end();
    typename std::vector< monomial<T> >::const_iterator itb = other.coord.begin();
    typename std::vector< monomial<T> >::const_iterator itb_end = other.coord.end();
    //  COUT << coord.size() << " " << (int) ita_end - (int) ita << " " << sizeof(monomial<T>) << '\n' ;
    // first some trivial cases
    if (ita==ita_end)
    return(*this);
    if (itb==itb_end)
    return(other);
    if (is_one(*this))
    return other;
    if (is_one(other))
    return *this;
    // Now look if length a=1 or length b=1, happens frequently
    // think of x^3*y^2*z translated to internal form
    int c1=coord.size();
    if (c1==1)
    return(other.shift(coord.front().index,coord.front().value));
    int c2=other.coord.size();
    if (c2==1)
    return(this->shift(other.coord.front().index,other.coord.front().value));
    std::vector< monomial<T> > new_coord;
    new_coord.reserve(c1+c2); // assumes dense poly (would be c1+c2-1)
    Mul(ita,ita_end,itb,itb_end,new_coord,m_is_strictly_greater);
    return tensor<T>(*this,new_coord);
    }

  */

  template <class T>
  tensor<T> Tpow(const tensor<T> & x,int n){
#ifndef NO_STDEXCEPT
    if (n<0)
      setsizeerr("poly.h/Tpow n<0");
#endif
    if (!n)
      return tensor<T>(T(1),x.dim);
    if (n==1)
      return x;
    if (n==2)
      return x*x;
    if (x.coord.size()==1)
      return tensor<T>(monomial<T>(pow(x.coord.front().value,n),x.coord.front().index*n));
    tensor<T> res(x);
    for (int j=1;j<n;j++)
      res=res*x;
    return res;
    /* 
       Note: contrary to univariate polynomials or integers
       the "fast" powering algorithm is *slower* than the above
       loop for multivariate polynomials (with the current implementation of Mul)
       Indeed a dense poly of deg. aa and d variables may have 
       binomial(aa+d,d) monomials
       hence the last multiplication in the "fast" powering algorithm
       is O(binomial(n/2*deg+d,d)^2)=O(n^(2d))
       as the n multiplications of the above loop are
       O(n*binomial(n*deg+d,d)) = O(n^(d+1))
       tensor<T> temp=Tpow(x,n/2);
       if (n%2)
       return tensor<T>(temp*temp*x);
       else
       return tensor<T>(temp*temp);
    */
  }

  /*
    template <class T>
    tensor<T> tensor<T>::operator / (const tensor<T> & other) const {  
    tensor<T> rem(*this),quo(*this);
    assert( (*this).TDivRem(other,quo,rem) );
    return(quo);
    }

  template <class T>
  tensor<T> tensor<T>::operator % (const tensor<T> & other) const {  
    tensor<T> rem(*this),quo(*this);
    assert( (*this).TDivRem(other,quo,rem) );
    return(rem);
  }
  
  */

  /*
  template <class T>
  tensor<T> tensor<T>::operator * (const T & fact ) const {  
    // Tensor constant multiplication
    if (is_one(fact))
      return *this;
    std::vector< monomial<T> > new_coord;
    if (is_zero(fact))
      return tensor<T>(*this,new_coord);
    typename std::vector< monomial<T> >::const_iterator a = coord.begin();
    typename std::vector< monomial<T> >::const_iterator a_end = coord.end();
    Mul(a,a_end,fact,new_coord);
    return tensor<T>(*this,new_coord);
  }

  inline template <class T>
  tensor<T> operator * (const T & fact ,const tensor<T> & p){
    return p*fact;
  }
  */

  template <class T>
  tensor<T> & tensor<T>::operator *= (const T & fact ) {  
    // Tensor constant multiplication
    if (is_one(fact))
      return *this;
    if (is_zero(fact)){
      coord.clear();
      return *this;
    }
    typename std::vector< monomial<T> >::const_iterator a = coord.begin();
    typename std::vector< monomial<T> >::const_iterator a_end = coord.end();
    Mul<T>(a,a_end,fact,coord);
    return *this;
  }

  /*
  template <class T>
  tensor<T> tensor<T>::operator / (const T & fact ) const {  
    if (is_one(fact))
      return *this;
    std::vector< monomial<T> > new_coord;
    typename std::vector< monomial<T> >::const_iterator a = coord.begin();
    typename std::vector< monomial<T> >::const_iterator a_end = coord.end();
    Div(a,a_end,fact,new_coord);
    return tensor<T>(*this,new_coord);
  }
  */

  template <class T>
  tensor<T> & tensor<T>::operator /= (const T & fact ) {  
    if (is_one(fact))
      return *this;
    typename std::vector< monomial<T> >::const_iterator a = coord.begin();
    typename std::vector< monomial<T> >::const_iterator a_end = coord.end();
    Div<T>(a,a_end,fact,coord);
    return *this;
  }

  template<class T>
  tensor<T> Tnextcoeff(typename std::vector< monomial<T> >::const_iterator & it,const typename std::vector< monomial<T> >::const_iterator & itend){
    if (it==itend)
      return tensor<T>(0);
    int n=it->index.front();
    int d=int(it->index.size());
    tensor<T> res(d-1);
    for (;(it!=itend) && (it->index.front()==n);++it)
      res.coord.push_back(it->trunc1());
    return res;
  }

  template<class T>
  tensor<T> Tlastcoeff(const typename std::vector< monomial<T> >::const_iterator & itbeg,const typename std::vector< monomial<T> >::const_iterator & itend){
    assert(itbeg!=itend);
    typename std::vector< monomial<T> >::const_iterator it=itend;
    --it;
    int n=it->index.front();
    int d=int(it->index.size());
    tensor<T> res(d-1);
    for (;;){
      if (it==itbeg)
	break;
      --it;
      if (it->index.front()!=n){
	++it;
	break;
      }
    }
    for (;it!=itend;++it)
      res.coord.push_back(it->trunc1());
    return res;
  }

  template <class T>
  bool tensor<T>::TDivRem1(const tensor<T> & b,tensor<T> & quo,tensor<T> & r,bool allowrational,int exactquo) const {
    const tensor<T> & a=*this;
    quo.coord.clear();
    quo.dim=a.dim;
    r.dim=a.dim;
    r.coord.clear();
    if ( b.dim<=1 || b.coord.size()==1 || a.coord.empty() ){
      return a.TDivRem(b,quo,r,allowrational) && (exactquo?r.coord.empty():true) ;
    }
    /* alternative code
       std::vector< tensor<T> > R=Tcoeffs();
       std::vector< tensor<T> > B=b.Tcoeffs();
       tensor<T> Q;
       int rs=R.size(),bs=B.size(),qs=rs-bs;
       if (rs<bs){
       r=a;
       return true;
       }
       for (int i=0;i<=qs;++i){
       if (!R[i].Texactquotient(B[0],Q,allowrational))
       return false;
       for (int j=i+1;j<i+bs;++j){
       if (!B[j-i].coord.empty())
       R[j]=R[j]-Q*B[j-i];
       }
       Q=Q.untrunc1(qs-i);
       quo.append(Q);
       }
       // convert R[qs+1..rs] to r
       for (int i=qs+1;i<rs;++i){
       R[i]=R[i].untrunc1(rs-i-1);
       r.append(R[i]);
       }
       return true;
    */
    // old code
    typename std::vector< monomial<T> >::const_iterator it=b.coord.begin(),itend=b.coord.end();
    r=a;
    if (b.coord.empty()){
      if (exactquo)
	return r.coord.empty();
      return true;
    }
    int bdeg=it->index.front(),rdeg=r.lexsorted_degree(),qdeg=rdeg-bdeg; // int bval=(itend-1)->index.front();
    tensor<T> b0(Tnextcoeff<T>(it,b.coord.end()));
    tensor<T> q(b0.dim);
    if (it==b.coord.end()){ // bn==b0, b=b0*main var to a power
      it=r.coord.begin();
      itend=r.coord.end();
      while ( it!=itend){
	rdeg=it->index.front();
	if (rdeg<bdeg)
	  break;
	tensor<T> a0(Tnextcoeff<T>(it,itend)),tmp(a0.dim);
	if (!a0.Texactquotient(b0,q,allowrational))
	  return false;
	q=q.untrunc1(rdeg-bdeg);
	quo.append(q);
      }
      r.coord=std::vector< monomial<T> >(it,itend);
      if (exactquo)
	return it==itend;
      return true;
    }
    tensor<T> q1,b1(b0.dim); // subprincipal coeff
    if (it->index.front()==bdeg-1)
      b1=Tnextcoeff<T>(it,b.coord.end()); 
    // here it might be improved by checking that lastcoeff divides
#if 1
    if (exactquo){
      tensor<T> bn=Tlastcoeff<T>(b.coord.begin(),b.coord.end());
      tensor<T> an(Tlastcoeff<T>(a.coord.begin(),a.coord.end())),qn;
      if (!an.Texactquotient(bn,qn,allowrational))
	return false;
    }
#endif
    for ( ;(rdeg=r.lexsorted_degree()) >=bdeg;){
      it=r.coord.begin();
      itend=r.coord.end();
      // compute 2 terms of the quotient in one iteration,
      // this will save one long subtraction
      tensor<T> a0(Tnextcoeff<T>(it,itend)),a1(a0.dim);
      if (exactquo && it==itend)
	return false;
      if (!a0.Texactquotient(b0,q,allowrational))
	return false;
      qdeg=rdeg-bdeg;
      if (qdeg){
	if (it!=itend && it->index.front()==rdeg-1)
	  a1=Tnextcoeff<T>(it,itend);
#if defined GIAC_VECTOR || defined NSPIRE // fix for * on arm compiler
	tensor<T> tmp(q.dim);
	tmp.coord=q.coord*b1.coord;
	a1.TSub(tmp,a1);
#else
	a1.TSub(q*b1,a1); // a1=a1-q*b1; 
#endif
	q=q.untrunc1(qdeg);
	if (!a1.Texactquotient(b0,q1,allowrational))
	  return false;
	q.TAdd(q1.untrunc1(qdeg-1),q); // q=q+q1.untrunc1(qdeg-1);
      }
      else
	q=q.untrunc1(qdeg);	
      /*
      if (exactquo){
	tensor<T> an(Tlastcoeff<T>(it,itend));
	int rval=(itend-1)->index.front();
	if (rval-bval<rdeg-bdeg){
	  if (rval<bval+count) // we must gain one valuation per loop
	    return false;
	  if (!an.Texactquotient(bn,qn,allowrational))
	    return false;
	  q=q+qn.untrunc1(rval-bval);
	  ++count;
	}
      }
      */
      quo.TAdd(q,quo); // quo=quo+q;
#if defined GIAC_VECTOR || defined NSPIRE
      tensor<T> tmp(q.dim);
      tmp.coord=q.coord*b.coord;
      r.TSub(tmp,r); // r=r-q*b;       
#else
      r.TSub(q*b,r); // r=r-q*b; 
#endif
      if (r.coord.empty())
	return true;
    }
    if (exactquo) return r.coord.empty();
    return true;
  }


  // hashing seems slower than DivRem1 if computation should be done with int
  template <class T>
  bool tensor<T>::TDivRemHash(const tensor<T> & b,tensor<T> & quo,tensor<T> & r,bool allowrational,int exactquo,double qmax) const {
    const tensor<T> & a=*this;
    quo.coord.clear();
    quo.dim=a.dim;
    r.dim=a.dim;
    r.coord.clear();
    int bs=b.coord.size();
    if ( b.dim<=1 || bs==1 || a.coord.empty() ){
      return a.TDivRem(b,quo,r,allowrational) && (exactquo?r.coord.empty():true) ;
    }
    int bdeg=b.coord.front().index.front(),rdeg=lexsorted_degree(),ddeg=rdeg-bdeg;
    if (ddeg>2 && bs>10){
      index_t d1=degree(),d2=b.degree(),d3=b.coord.front().index.iref(),d(dim);
      // i-th degrees of th / other in quotient and remainder
      // are <= i-th degree of th + ddeg*(i-th degree of other - i-th degree of lcoeff of other) 
      double ans=1;
      for (int i=0;i<dim;++i){
	d[i]=d1[i]+(ddeg+1)*(d2[i]-d3[i])+1;
	int j=1;
	// round to newt power of 2
	for (;;j++){
	  if (!(d[i] >>= 1))
	    break;
	}
	d[i] = 1 << j;
	ans = ans*unsigned(d[i]);
	if (ans/RAND_MAX>RAND_MAX)
	  break;
      }
      if (ans<RAND_MAX){
	std::vector< T_unsigned<T,unsigned> > p1,p2,quot,remain;
	std::vector<unsigned> vars(dim);
	vars[dim-1]=1;
	for (int i=dim-2;i>=0;--i){
	  vars[i]=d[i+1]*vars[i+1];
	}
	convert(*this,d,p1);
	convert(b,d,p2);
	if (hashdivrem<T,unsigned>(p1,p2,quot,remain,vars,/* reduce */0,qmax,allowrational)){
	  convert(quot,d,quo);
	  convert(remain,d,r);
	  return true;
	}
	else
	  ans=1e200; // don't make another unsuccessful division!
      }
      if (ans/RAND_MAX<RAND_MAX){
	std::vector< T_unsigned<T,ulonglong> > p1,p2,quot,remain;
	std::vector<ulonglong> vars(dim);
	vars[dim-1]=1;
	for (int i=dim-2;i>=0;--i){
	  vars[i]=d[i+1]*vars[i+1];
	}
	convert(*this,d,p1);
	convert(b,d,p2);
	if (hashdivrem<T,ulonglong>(p1,p2,quot,remain,vars,/* reduce */0,qmax,allowrational)){
	  convert(quot,d,quo);
	  convert(remain,d,r);
	  return true;
	}
      }
    }
    return TDivRem1(b,quo,r,allowrational,exactquo);
  }

  template<class T> tensor<T> tensor<T>::coeff(int deg) const {
    tensor<T> res(dim-1);
    typename std::vector< monomial<T> >::const_iterator it=coord.begin(),itend=coord.end();
    for (;it!=itend;++it){
      int d=it->index.front();
      if (d>deg)
	continue;
      if (d==deg)
	return Tnextcoeff<T>(it,itend);
      else
	return res;
    }
    return res;
  }

  template<class T>
  void tensor<T>::Tcoeffs(std::vector< tensor<T> > & v) const{
    int current_deg=lexsorted_degree();
    std::vector< tensor<T> > w(current_deg+1,dim-1);
    typename std::vector< monomial<T> >::const_iterator it=coord.begin(),itend=coord.end();
    for (;it!=itend;++it){
      w[current_deg-it->index.front()].coord.push_back(it->trunc1());
    }
    w.swap(v);
  }

  template<class T>
  std::vector< tensor<T> > tensor<T>::Tcoeffs() const{
    int current_deg=lexsorted_degree();
    std::vector< tensor<T> > v;
    v.reserve(current_deg+1);
    typename std::vector< monomial<T> >::const_iterator it=coord.begin(),itend=coord.end();
    for (;it!=itend;--current_deg){
      if (it->index.front()==current_deg){
	v.push_back(Tnextcoeff<T>(it,itend));
      }
      else 
	v.push_back(tensor<T>(dim-1));
    }
    for (;current_deg>=0;--current_deg) v.push_back(tensor<T>(dim-1));
    return v;
  }

  template <class T>
  bool tensor<T>::Texactquotient(const tensor<T> & b,tensor<T> & quo,bool allowrational) const {
    if (coord.empty()){
      quo.dim=dim; quo.coord.clear();
      return true;
    }
    if (*this==b){
      quo=tensor<T>(T(1),dim);
      return true;
    }
    if (dim>1 && !allowrational && lexsorted_degree()==b.lexsorted_degree()){
      if (!Tfirstcoeff(*this).trunc1().Texactquotient(Tfirstcoeff(b).trunc1(),quo,allowrational))
	return false;
      quo=quo.untrunc1();
      if (is_one(quo))
	return false; // already tested above
      return *this==quo*b;
    }
    tensor<T> r(b.dim);
    return this->TDivRem1(b,quo,r,allowrational,1);
  }

  template <class T>
  bool tensor<T>::TDivRem (const tensor<T> & other, tensor<T> & quo, tensor<T> & rem, bool allowrational ) const {  
    int asize=int((*this).coord.size());
    if (!asize){
      quo=*this;
      rem=*this; 
      return true;
    }
    int bsize=int(other.coord.size());
    if (!bsize){
      quo.dim=dim; quo.coord.clear();
      rem=*this;
      return true;
    }
    index_m a_max = (*this).coord.front().index;
    index_m b_max = other.coord.front().index;
    quo.coord.clear();
    quo.dim=this->dim;
    rem.dim=this->dim;
    if (bsize==1){
      rem.coord.clear();
      T b=other.coord.front().value;
      if (b_max==b_max*0){
	if (is_one(b))
	  quo = *this ;
	else {
	  typename std::vector< monomial<T> >::const_iterator itend=coord.end();
	  for (typename std::vector< monomial<T> >::const_iterator it=coord.begin();it!=itend;++it){
	    T q=rdiv(it->value,b);
	    if (!allowrational && has_denominator(q))
	      return false;
	    quo.coord.push_back(monomial<T>(q,it->index)); 
	  }
	}
	return true;
      }
      typename std::vector< monomial<T> >::const_iterator itend=coord.end();
      typename std::vector< monomial<T> >::const_iterator it=coord.begin();
      for (;it!=itend;++it){
	if (!(it->index>=b_max))
	  break;
	T q=rdiv(it->value,b);
	if (!allowrational && has_denominator(q))
	  return false;
	quo.coord.push_back(monomial<T>(q,it->index-b_max)); 
      }
      rem.coord=std::vector< monomial<T> >(it,itend);
      return true;
    }
    rem=*this;
    if ( ! (a_max>=b_max) ){
      // test that the first power of a_max is < to that of b_max
      return (a_max.front()<b_max.front());
    }
    T b=other.coord.front().value;
    while (a_max >= b_max){
      // errors should be trapped here and false returned if error occured
      T q=rdiv(rem.coord.front().value,b);
      if (!allowrational){
	if ( has_denominator(q) || 
	     (!is_zero(q*b - rem.coord.front().value)) )
	  return false;
      }
      // end error trapping
      quo.coord.push_back(monomial<T>(q,a_max-b_max));
      const tensor<T> & temp=other.shift(a_max-b_max,q);
      rem.TSub(temp,rem); // rem = rem-temp;
      if (rem.coord.size())
	a_max=rem.coord.front().index;
      else
	break;
    }
    return(true);    
  }

  template <class T>
  bool tensor<T>::TPseudoDivRem (const tensor<T> & other, tensor<T> & quo, tensor<T> & rem, tensor<T> & a) const {
    int m=this->lexsorted_degree();
    int n=other.lexsorted_degree();
    a.coord.clear();
    a.coord.push_back(monomial<T>(T(1),a.dim));
    rem=*this;
    quo.coord.clear();
    if (m<n)
      return true;
    // a=Tpow(Tfirstcoeff(other),m-n+1);
    // return (a*(*this)).TDivRem1(other,quo,rem);
    index_m ishift(dim);
    tensor<T> b0(Tfirstcoeff(other));
    for (int i=m;i>=n;--i){
#if defined NSPIRE || defined(FXCG)
      a=a*b0;
      quo=quo*b0;
#else
#ifdef GIAC_VECTOR
      a.coord *= b0.coord; // a.coord = a.coord*b0.coord; // a=a*b0;
      quo.coord *= b0.coord; // quo.coord = quo.coord*b0.coord; // quo=quo*b0;
#else
      a *= b0; // a=a*b0
      quo *= b0; // quo=quo*b0
#endif // GIAC_VECTOR
#endif // NSPIRE
      typename std::vector< monomial<T> >::const_iterator it=rem.coord.begin(),itend=rem.coord.end();
      if (it==itend || it->index.front()!=i){
#if defined NSPIRE || defined(FXCG)
	rem=rem*b0;
#else
#ifdef GIAC_VECTOR
	rem.coord *= b0.coord; // rem.coord = rem.coord*b0.coord; // rem=rem*b0;
#else
	rem *= b0; // rem=rem*b0;
#endif
#endif // NSPIRE
	continue;
      }
      ishift.front()=i-n;
      const tensor<T> & rem0 = Tfirstcoeff(rem).shift(ishift);
      quo.append(rem0);
#if defined GIAC_VECTOR || defined NSPIRE
      rem.coord = rem.coord*b0.coord; 
      tensor<T> tmp(rem0.dim);
      tmp.coord=rem0.coord*other.coord;
      rem.TSub(tmp,rem); // rem=rem*b0-rem0*other;
#else
      rem=rem*b0-rem0*other;
#endif
    }
    return true;
  }

  template <class T>
  T tensor<T>::constant_term () const {  
    if (!((*this).coord.size()))
      return T(0);
#if 1
    if (sum_degree((*this).coord.back().index)==0)
      return (*this).coord.back().value;
    return T(0);
#else
    index_m i=(*this).coord.front().index*0;
    return (*this)(i);
#endif
  }

  template <class T>
  tensor<T> tensor<T>::shift (const index_m & i,const T & fois) const {
    tensor<T> res(dim,*this);
    res.coord.reserve(coord.size());
    Shift(this->coord,i,fois,res.coord);
    return res;
  }

  template <class T>
  tensor<T> tensor<T>::shift (const T & fois,const index_m & i) const {
    tensor<T> res(dim,*this);
    res.coord.reserve(coord.size());
    Shift(this->coord,fois,i,res.coord);
    return res;
  }

  template <class T>
  tensor<T> tensor<T>::shift (const index_m & i) const {
    tensor<T> res(dim);
    res.coord.reserve(coord.size());
    Shift(this->coord,i,res.coord);
    return res;
  }

  template <class T>
  void tensor<T>::reverse() {
    typename std::vector< monomial<T> >::const_iterator itend=coord.end();
    for (typename std::vector< monomial<T> >::iterator it=coord.begin();it!=itend;++it)
      (*it).reverse();
    this->tsort();
  }

  template <class T>
  void tensor<T>::reorder(const std::vector<int> & permutation) {
    typename std::vector< monomial<T> >::const_iterator itend=coord.end();
    for (typename std::vector< monomial<T> >::iterator it=coord.begin();it!=itend;++it)
      it->reorder(permutation);
    this->tsort();
  }

  template <class T>
  T tensor<T>::norm() const{
    T res( 0);
    typename std::vector< monomial<T> >::const_iterator itend=coord.end();
    for (typename std::vector< monomial<T> >::const_iterator it=coord.begin();it!=itend;++it)
      res=max(res,it->norm());
    return res;
  }

  template <class T>
  T Tcontent(const tensor<T> & p){
    return Content(p.coord);
  }

  template<class T>
  T Tppz(tensor<T> & p,bool divide=true){
    T n=Tcontent(p);
    if (divide)
      p /= n;
    return T(n);
  }

  template<class T>
  bool Tis_one(const tensor<T> &p){
    if (p.coord.size()!=1)
      return false;
    if (!is_one(p.coord.front().value))
      return false;
    const index_m & i = p.coord.front().index;
    index_t::const_iterator it=i.begin(),itend=i.end();
    for (;it!=itend ;++it){
      if ((*it)!=0)
	return false;
    }
    return true;
  }

  template <class T>
  bool Tis_constant(const tensor<T> & p){
    if (p.coord.size()!=1)
      return false;
    const index_m & i = p.coord.front().index;
    index_t::const_iterator it=i.begin(),itend=i.end();
    for (;it!=itend ;++it){
      if ((*it)!=0)
	return false;
    }
    return true;  
  }

  template<class T>
  bool Tlistmax(const tensor<T> &p,T & n){
    n=T(  1);
    for (typename std::vector< monomial<T> >::const_iterator it=p.coord.begin();it!=p.coord.end() ;++it){
      if (!(it->value.is_cinteger()))
	return false;
      n=max(n,linfnorm(it->value));
    }
    return true;
  }

  template<class T>
  void Tapply(const tensor<T> & p,T (*f)(const T &),tensor<T> & pgcd){
    typename std::vector< monomial<T> >::const_iterator it=p.coord.begin(),itend=p.coord.end();
    Apply(it,itend,f,pgcd.coord);
  }

  template<class T>
  tensor<T> Tapply(const tensor<T> &p,T (*f)(const T &)){
    tensor<T> res(p.dim);
    Tapply(p,f,res);
    return res;
  }

  template<class T>
  void Tlgcd(const tensor<T> & p,tensor<T> & pgcd){
    if (!p.dim){
      pgcd=p;
      return ;
    } 
    if (Tis_one(pgcd))
      return;
    pgcd=pgcd.trunc1();
    typename std::vector< monomial<T> >::const_iterator it=p.coord.begin();
    typename std::vector< monomial<T> >::const_iterator itend=p.coord.end();
    // typename std::vector< monomial<T> >::const_iterator itbegin=it;
    for (;it!=itend;){
      if (Tis_one(pgcd))
	break;
      pgcd=gcd(pgcd,Tnextcoeff<T>(it,itend));
    }
    if (pgcd.coord.empty()){
      index_m i;
      for (int j=0;j<p.dim;j++)
	i.push_back(0);
      pgcd.coord.push_back(monomial<T>(T(1),i));
      ++pgcd.dim;
    }
    else
      pgcd=pgcd.untrunc1();
  }

  template<class T>
  tensor<T> Tlgcd(const tensor<T> & p){
    if (p.dim==1){
      const T & c=Tcontent(p);
      return tensor<T>(c,1);
    }
    tensor<T> pgcd(p.dim); // pgcd=0
    Tlgcd(p,pgcd);
    return pgcd;
  }

  template<class T>
  void Tcommonlgcd(const tensor<T> & p,const tensor<T> & q,tensor<T> & pgcd){
    if (!p.dim){
      pgcd=p;
      return ;
    } 
    pgcd=pgcd.trunc1();
    typename std::vector< monomial<T> >::const_iterator it=p.coord.begin();
    typename std::vector< monomial<T> >::const_iterator itend=p.coord.end();
    typename std::vector< monomial<T> >::const_iterator jt=q.coord.begin();
    typename std::vector< monomial<T> >::const_iterator jtend=q.coord.end();
    for (;(it!=itend) && (jt!=jtend);){
      if (Tis_one(pgcd))
	break;
      pgcd=gcd(pgcd,Tnextcoeff<T>(it,itend));
      pgcd=gcd(pgcd,Tnextcoeff<T>(jt,jtend));
    }
    for (;(it!=itend);){
      if (Tis_one(pgcd))
	break;
      pgcd=gcd(pgcd,Tnextcoeff<T>(it,itend));
    }
    for (;(jt!=jtend);){
      if (Tis_one(pgcd))
	break;
      pgcd=gcd(pgcd,Tnextcoeff<T>(jt,jtend));
    }
    if (pgcd.coord.empty()){
      index_m i;
      for (int j=0;j<p.dim;j++)
	i.push_back(0);
      pgcd.coord.push_back(monomial<T>(T(1),i));
    }
    else
      pgcd=pgcd.untrunc1();
  }

  // collect all terms with same 1st Tpower as the 1st power of p
  template<class T>
  tensor<T> Tcarcomp(const tensor<T> & p){
    typename std::vector< monomial<T> >::const_iterator it=p.coord.begin();
    int n=it->index.front();
    tensor<T> res(p.dim);
    for (;n==it->index.front();++it)
      res.coord.push_back(monomial<T>(it->value,it->index));
    return res;
  }

  template<class T>
  tensor<T> Tfirstcoeff(const tensor<T> & p){
    typename std::vector< monomial<T> >::const_iterator it=p.coord.begin();
    typename std::vector< monomial<T> >::const_iterator itend=p.coord.end();
    if (it==itend)
      return p;
    int n=it->index.front();
    tensor<T> res(p.dim);
    for (;(it!=itend) && (n==it->index.front());++it)
      res.coord.push_back(monomial<T>(it->value,it->index.set_first_zero()));
    return res;
  }

  template<class T>
  void Tswap(tensor<T> * & a, tensor<T> * & b){
    tensor<T> * c =  a;
    a = b ;
    b = c ;
  }

  // polynomial subresultant: compute content and primitive part of the GCD
  // with respect to the other vars
  // gcddeg!=0 if we know the probable degree of the gcd
  template<class T>
  void Tcontentgcd(const tensor<T> &p, const tensor<T> & q, tensor<T> & cont,tensor<T> & prim,int gcddeg){
    if (p.coord.empty()){
      cont=Tlgcd(q);
      prim=q/Tlgcd(q);
      return ;
    }
    if (q.coord.empty()){
      cont=Tlgcd(p);
      prim=p/Tlgcd(p);
      return;
    }
    assert(p.dim==q.dim);
    // set auxiliary polynomials g and h to 1
    tensor<T> g(T(1),p.dim);
    tensor<T> h(g);
    // dp and dq are the "content" of p and q w.r.t. other variables
    tensor<T> dp(Tlgcd(p));
    tensor<T> dq(Tlgcd(q));
    cont=gcd(dp.trunc1(),dq.trunc1()).untrunc1();
    if (!p.dim){
      prim=tensor<T>(T(1),0);
      return ;
    }
    // COUT << "Cont" << cont << '\n'; 
    tensor<T> a(p.dim),b(p.dim),quo(p.dim),r(p.dim),tmp(p.dim);
    // a and b are the primitive part of p and q
    p.TDivRem1(dp,a,r,true);
    q.TDivRem1(dq,b,r,true);
    while (
	   !a.coord.empty()
	   && !ctrl_c && !interrupted
	   ){
      int n=b.lexsorted_degree();
      int m=a.lexsorted_degree();
      if (!n) {// if b is constant (then b!=0), gcd=original Tlgcd
	prim=tensor<T>(T(1),p.dim);
	return ;
      }
      if (n==gcddeg){
	b.TDivRem1(Tlgcd(b),prim,r,true);
	if (a.Texactquotient(prim,r,true))
	  return;
      }
      int ddeg=m-n;
      if (ddeg<0){
#if defined RTOS_THREADX || defined BESTA_OS || defined USTL
	tensor<T> t(a); a=b; b=t;
#else
	swap(a,b); // exchange a<->b may occur only at the beginning
#endif
      }
      else {
	tensor<T> b0(Tfirstcoeff(b));
	a.TPseudoDivRem(b,quo,r,tmp);
	// (a*Tpow(b0,ddeg+1)).TDivRem1(b,quo,r); // division works always
	if (r.coord.empty())
	  break;
	// remainder is non 0, loop continue: a <- b
	a=b;
	tensor<T> temp(Tpow(h,ddeg));
	// now divides r by g*h^(m-n), result is the new b
	r.TDivRem1(g*temp,b,quo); // quo is the remainder here, not used
	// new g=b0 and new h=b0^(m-n)*h/temp
	if (ddeg==1) // the normal case, remainder deg. decreases by 1 each time
	  h=b0;
	else // not sure if it's better to keep temp or divide by h^(m-n+1)
	  (Tpow(b0,ddeg)*h).TDivRem1(temp,h,quo);
	g=b0;
      }
    }
    // COUT << "Prim" << b << '\n';
    b.TDivRem1(Tlgcd(b),prim,r,true);
  }

  template<class T>
  void Tsturm_seq(const giac::tensor<T> &p, giac::tensor<T> & cont,std::vector< tensor<T> > & sturm_seq){
    sturm_seq=std::vector< tensor<T> > (1,p);
    if (p.coord.empty()){
      cont=p;
      return ;
    }
    // set auxiliary polynomials g and h to 1
    tensor<T> g(T(1),p.dim);
    tensor<T> h(g);
    // dp and dq are the "content" of p and q w.r.t. other variables
    cont=Tlgcd(p);
    if (!p.dim)
      return ;
    // COUT << "Cont" << cont << '\n'; 
    tensor<T> a(p.dim),b(p.dim),quo(p.dim),r(p.dim),tmp(p.dim);
    tensor<T> b0(g);
    std::vector< tensor<T> > sign_error(2,g);
    // a is the primitive part of p and q
    p.TDivRem1(cont,a,r);
    b=p.derivative();
    sturm_seq.push_back(b);
    for (int loop_counter=0;!a.coord.empty();++loop_counter){
      int m=a.lexsorted_degree();
      int n=b.lexsorted_degree();
      int ddeg=m-n; // should be 1 generically
      if (!n) {// if b is constant (then b!=0), gcd=original Tlgcd
	return ;
      }
      b0=Tfirstcoeff(b); 
      a.TPseudoDivRem(b,quo,r,tmp);
      // (a*Tpow(b0,ddeg+1)).TDivRem1(b,quo,r); // division works always
      if (r.coord.empty())
	break;
      // remainder is non 0, loop continue: a <- b
      a=b;
      tensor<T> temp(Tpow(h,ddeg));
      // now divides r by g*h^(m-n) and change sign, result is the new b
      r.TDivRem1(-g*temp,b,quo); // quo is the remainder here, not used
      // Now save b in the Sturm sequence *with sign adjustement*
      // Since we have -g*h^ddeg*r=b0^(ddeg+1)*a % previous b
      // instead of -r=a % previous b for correct sign
      // the sign error on b does not change the sign error on r
      // if ddeg is odd, the sign error on r is sign error on a*sign(g*h)
      // if ddeg is even, the signe error on r is sign error on a*sign(b0*g)
      // Note that if ddeg is odd and the previous ddeg was generic = 1
      // we do not change the sign error at all
      // Note that the sign errors propagate modulo 2 on the indices
      // of the Sturm sequence hence the vector of sign errors of length 2
      if ( ( (ddeg %2)!=0 ) && (h!=g) )
	sign_error[loop_counter %2] = sign_error[loop_counter %2] * h * g ;
      if ( ( (ddeg %2)==0) && (b0!=g) )
	sign_error[loop_counter %2] = sign_error[loop_counter % 2] * b0 * g;
      sturm_seq.push_back(b*sign_error[loop_counter %2]);
      // new g=b0 and new h=b0^(m-n)*h/temp
      if (ddeg==1) // the normal case, remainder deg. decreases by 1 each time
	h=b0;
      else // not sure if it's better to keep temp or divide by h^(m-n+1)
	(Tpow(b0,ddeg)*h).TDivRem1(temp,h,quo);
      g=b0;
    } // end while loop
  }

  // polynomial sub-resultant: return the gcd only
  template<class T>
  tensor<T> Tgcdpsr(const tensor<T> & p, const tensor<T> & q,int gcddeg=0){
    tensor<T> prim(p.dim),cont(p.dim);
    Tcontentgcd(p,q,prim,cont,gcddeg);
    // COUT << "Prim" << prim << "Cont" << cont << '\n';
    return prim*cont;
  }

  template<class T>
  tensor<T> Tresultant(const tensor<T> &p, const tensor<T> & q){
    assert(p.dim==q.dim);
    if (p.coord.empty())
      return p;
    if (q.coord.empty())
      return q;
    int m=p.lexsorted_degree();
    int n=q.lexsorted_degree();
    int sign=1;
    tensor<T> ptmp(p),qtmp(q);
    if (n > m) {
      if ( (n*m) % 2)
	sign=-1;
      int tmpint=n; n=m; m=tmpint;
#if defined RTOS_THREADX || defined BESTA_OS || defined USTL
      ptmp=q; qtmp=p; // swap(ptmp,qtmp);
#else
      swap(ptmp,qtmp);
#endif
    }
    // degree(qtmp)=n <= degree(ptmp)=m
    if (!n) // q is cst
      return pow(qtmp,m)*T(sign);
    int ddeg;
    T cp(ppz(ptmp)),cq(ppz(qtmp));
    T res(pow(cq,m)*pow(cp,n));
    tensor<T> g(T(1),p.dim), h(T(1),p.dim);
    while (n){
      if (m*n %2)
	sign=-sign;
      ddeg=m-n;
      if (debug_infolevel)
	CERR << CLOCK() << "Tresultant n,m,ddeg: " << n << " ," << m << " ," << ddeg << '\n';
      tensor<T> tmp1(Tfirstcoeff(qtmp)),tmp2(p.dim),tmp3(pow(h,ddeg)),rem(p.dim),a(p.dim);
      ptmp.TPseudoDivRem(qtmp,tmp2,rem,a); // (ptmp*pow(tmp1,ddeg+1)).TDivRem1(qtmp,tmp2,rem,false);
      rem.high_order_degree_truncate(n);
      ptmp=qtmp;
      m=n;
      qtmp=(rem/g)/tmp3; // qtmp=rem/(g*tmp3);
      n=qtmp.lexsorted_degree();
      if (ddeg==1)
	h=tmp1;
      else
	h=(h*pow(tmp1,ddeg))/tmp3;
      g=tmp1;
    }
    // q is cst or zero
    if (qtmp.coord.empty())
      return tensor<T>(T(0),p.dim);
    m=ptmp.lexsorted_degree();
    return (pow(qtmp,m)/pow(h,m-1))*(res*T(sign));
  }

  // Bézout identity
  // given p and q, find u and v s.t. u*p+v*q=d where d=gcd(p,q) using PSR algo
  // Iterative algorithm to find u and d, then q=(d-u*p)/v
  template<class T>
  void Tegcdpsr(const tensor<T> &p1, const tensor<T> & p2, tensor<T> & u,tensor<T> & v,tensor<T> & d){
    assert(p1.dim==p2.dim);
    // set auxiliary polynomials g and h to 1
    tensor<T> g(T(1),p1.dim);
    tensor<T> h(g);
    tensor<T> a(p1.dim),b(p1.dim),q(p1.dim),r(p1.dim);
    const tensor<T> & cp1=Tlgcd(p1);
    const tensor<T> & cp2=Tlgcd(p2);
    bool Tswapped=false;
    if (p1.lexsorted_degree()<p2.lexsorted_degree())
      Tswapped=true;
    // initializes a and b to p1, p2
    const tensor<T> & pp1=Tis_one(cp1)?p1:p1/cp1;
    const tensor<T> & pp2=Tis_one(cp2)?p2:p2/cp2;
    if (Tswapped){
      a=pp2;
      b=pp1;
    }
    else {
      a=pp1;
      b=pp2;
    }
    // initializes ua to 1 and ub to 0, the coeff of u in ua*a+va*b=a
    tensor<T> ua(T(1),p1.dim), ub(p1.dim),ur(p1.dim);
    tensor<T> b0pow(p1.dim);
    // loop: ddeg <- deg(a)-deg(b), 
    // TDivRem: b0^(ddeg+1)*a = bq+r 
    // hence ur <- ua*b0^(ddeg+1)-q*ub verifies
    // ur*a+vr*b=r
    // a <- b, b <- r/(g*h^ddeg), ua <- ub and ub<- ur/(g*h^ddeg)
    // g <- b0, h <- b0^(m-n) * h / h^ddeg
    for (;;){
      int n=b.lexsorted_degree();
      int m=a.lexsorted_degree();
      if (!n){ // b is cst !=0 hence is the gcd, ub is valid
	break;
      }
      int ddeg=m-n;
      const tensor<T> & b0=Tfirstcoeff(b);
      // b0pow=Tpow(b0,ddeg+1);
      // (a*b0pow).TDivRem1(b,q,r); // division works always
      a.TPseudoDivRem(b,q,r,b0pow);
      // if r is 0 then b is the gcd and ub the coeff
      if (r.coord.empty())
	break;
      // COUT << ua*b0pow << '\n' << q*ub << '\n' ;
      (ua*b0pow).TSub(q*ub,ur); // ur=ua*b0pow-q*ub;
      // COUT << ur << '\n';
#if defined RTOS_THREADX || defined BESTA_OS || defined USTL
      a=b; 
#else
      swap(a,b); // a=b
#endif
      const tensor<T> & temp=Tpow(h,ddeg);
      // now divides r by g*h^(m-n), result is the new b
      r.TDivRem1(g*temp,b,q); // q is not used anymore
#if defined RTOS_THREADX || defined BESTA_OS || defined USTL
      ua=ub;
#else
      swap(ua,ub); // ua=ub
#endif
      ur.TDivRem1(g*temp,ub,q);
      // COUT << (b-ub*p1) << "/" << p2 << '\n';
      // new g=b0 and new h=b0^(m-n)*h/temp
      if (ddeg==1) // the normal case, remainder deg. decreases by 1 each time
	h=b0;
      else // not sure if it's better to keep temp or divide by h^(m-n+1)
	(Tpow(b0,ddeg)*h).TDivRem1(temp,h,q);
      g=b0;
    }
    // ub is valid and b is the gcd, vb=(b-ub*p1)/p2 if not Tswapped
    // vb is stored in ua
    // COUT << ub << '\n';
    if (Tswapped){
      (b-ub*pp2).TDivRem1(pp1,ua,r);
      ua *= cp2; // ua=ua*cp2;
      ub *= cp1; // ub=ub*cp1;
      b *= cp1; b *= cp2; // b=b*cp1*cp2;
    }
    else {
      (b-ub*pp1).TDivRem1(pp2,ua,r);
      ua *= cp1; // ua=ua*cp1;
      ub *= cp2; // ub=ub*cp2;
      b *= cp1; b *= cp2; // b=b*cp1*cp2;
    }
    // final simplifications
    q.coord.clear();
    Tlgcd(b,q); // q=Tlgcd(b);
    Tlgcd(ua,q);
    Tlgcd(ub,q);
    b.TDivRem1(q,d,r,true);  // d=b/Tlgcd
    if (Tswapped){
      ub.TDivRem1(q,v,r,true); // v=ub/Tlgcd
      ua.TDivRem1(q,u,r,true); // u=ua/Tlgcd
    }
    else {
      ub.TDivRem1(q,u,r,true); // u=ub/Tlgcd
      ua.TDivRem1(q,v,r,true); // v=ua/Tlgcd
    }
  }

  // seems interesting in dimension 1 only
  template<class T>
  void TegcdTlgcd(const tensor<T> &p1, const tensor<T> & p2, tensor<T> & u,tensor<T> & v,tensor<T> & d){
    assert(p1.dim==p2.dim);
    // p1=cp1*a p2=cp2*b
    // a*u+b*v=d
    // hence multiplying by cp1*cp2 p1*(u*cp2)+p2*(v*cp1)=d*cp1*cp2
    tensor<T> cp1(Tlgcd(p1)), cp2(Tlgcd(p2)),pp1(p1/cp1),pp2(p2/cp2);
    tensor<T> a(p1.dim),b(p1.dim),q(p1.dim),r(p1.dim);
    bool Tswapped=false;
    if (p1.lexsorted_degree()<p2.lexsorted_degree())
      Tswapped=true;
    // initializes a and b to p1, p2
    if (Tswapped){
      a=pp2;
      b=pp1;
    }
    else {
      a=pp1;
      b=pp2;
    }
    // initializes ua to 1 and ub to 0, the coeff of u in ua*a+va*b=a
    tensor<T> ua(T(1),p1.dim), ub(p1.dim),ur(p1.dim);
    // loop: ddeg <- deg(a)-deg(b), 
    // TDivRem: b0^(ddeg+1)*a = bq+r 
    // hence ur <- ua*b0^(ddeg+1)-q*ub verifies
    // ur*a+vr*b=r
    // divide r and ur by their common Tlgcd
    // a <- b, b <- r/Tlgcd, ua <- ub and ub<- ur/Tlgcd
    for (;;){
      int n=b.lexsorted_degree();
      // int m=a.lexsorted_degree();
      if (!n){ // b is cst !=0 hence is the gcd, ub is valid
	break;
      }
      // int ddeg=m-n;
      tensor<T> b0(Tfirstcoeff(b)),b0pow(b0.dim);
      // b0pow=Tpow(b0,ddeg+1);
      // (a*b0pow).TDivRem1(b,q,r,true); // division works always
      a.TPseudoDivRem(b,q,r,b0pow);
      // if r is 0 then b is the gcd and ub the coeff
      if (r.coord.empty())
	break;
      ur=ua*b0pow-q*ub;
      a=b;
      tensor<T> pgcd(Tlgcd(ur));
      Tlgcd(r,pgcd);
      // now divides r by pgcd, result is the new b
      r.TDivRem1(pgcd,b,q,true); // q is not used anymore
      ua=ub;
      ur.TDivRem1(pgcd,ub,q,true);
    }
    // ub is valid and b is the gcd, vb=(b-ub*p1)/p2 if not Tswapped
    // vb is stored in ua
    // COUT << ub << '\n';
    if (Tswapped){
      q=b-ub*pp2;
      tensor<T> b0(Tfirstcoeff(pp1));
      int m=q.lexsorted_degree();
      int n=pp1.lexsorted_degree();
      tensor<T> b0pow(Tpow(b0,m-n+1));
      b=b*b0pow;
      ub=ub*b0pow;
      q=q*b0pow;
      q.TDivRem1(pp1,ua,r,true);
      ua=ua*cp2;
      ub=ub*cp1;
      b=b*cp1*cp2;
    }
    else {
      // first multiply b and ub by p1_max_deg^# before doing division
      q=b-ub*pp1;
      tensor<T> b0(Tfirstcoeff(pp2));
      int m=q.lexsorted_degree();
      int n=pp2.lexsorted_degree();
      tensor<T> b0pow(Tpow(b0,m-n+1));
      b=b*b0pow;
      ub=ub*b0pow;
      q=q*b0pow;
      q.TDivRem1(pp2,ua,r,true);
      ua=ua*cp1;
      ub=ub*cp2;
      b=b*cp1*cp2;
    }
    // final simplifications
    q=Tlgcd(ua);
    Tlgcd(ub,q);
    Tlgcd(b,q);
    b.TDivRem1(q,d,r,true);  // d=b/Tlgcd
    if (Tswapped){
      ua.TDivRem1(q,u,r,true); // u=ua/Tlgcd
      ub.TDivRem1(q,v,r,true); // v=ub/Tlgcd
    }
    else {
      ub.TDivRem1(q,u,r,true); // u=ub/Tlgcd
      ua.TDivRem1(q,v,r,true); // v=ua/Tlgcd
    }
  }

  // utility for Bézout identity solving
  template<class T>
  void Tegcdtoabcuv(const tensor<T> & a,const tensor<T> &b, const tensor<T> &c, tensor<T> &u,tensor<T> &v, tensor<T> & d, tensor<T> & C){
    tensor<T> d0(Tfirstcoeff(d));
    int m=c.lexsorted_degree();
    int n=d.lexsorted_degree();
    assert(m>=n); // degree of c must be greater than degree of d
    // IMPROVE if n==0 
    C=Tpow(d0,m-n+1);
    tensor<T> coverd(a.dim),temp(a.dim);
    (c*C).TDivRem1(d,coverd,temp);
    assert(temp.coord.empty()); // division of c by d must be exact
    // now multiply a*u+b*v=d by coverd -> a*u*coverd+b*v*coverd=c*d0pow
    u *= coverd; // u=u*coverd;
    v *= coverd; // v=v*coverd;
    m=u.lexsorted_degree();
    n=b.lexsorted_degree();
    if (m<n)
      return;
    // then reduces the degree of u, a*u+b*v=c*C
    d0=Tpow(Tfirstcoeff(b),m-n+1); 
    C *= d0; // C=C*d0;
    // now a*u*d0+b*v*d0=c*C
    (u*d0).TDivRem1(b,temp,u); // replace u*d0 -> temp*b+u
    // a*b + b*(a*temp+v*d0) = c*C
    v=a*temp+v*d0;
    return ;
  }

  // given a,b,c, find u, v and a cst C s.t. C*c=a*u+b*v and degree(u)<degree(b)
  // requires an egcd implementation, for example egcdpsr above
  template<class T>
  void Tabcuv(const tensor<T> & a,const tensor<T> &b, const tensor<T> &c, tensor<T> &u,tensor<T> &v, tensor<T> & C){
    tensor<T> d(a.dim);
    Tegcdpsr(a,b,u,v,d); // a*u+b*v=d
    Tegcdtoabcuv(a,b,c,u,v,d,C);
  }

  template<class T>
  tensor<T> tensor<T>::derivative() const{
    if (coord.empty())
      return *this;
    tensor<T> res(dim);
    if (dim==0)
      return res;
    res.coord.reserve(coord.size());
    typename std::vector< monomial<T> >::const_iterator itend=coord.end();
    T tmp;
    for (typename std::vector< monomial<T> >::const_iterator it=coord.begin();it!=itend;++it){
      index_t i= it->index.iref() ;
      T n(i.front());
      i[0]--;
      tmp = it->value*n;
      if (!is_zero(tmp))
	res.coord.push_back(monomial<T>(tmp,i));
    }
    return res;
  }

  template<class T>
  tensor<T> tensor<T>::integrate() const{
    if (coord.empty())
      return *this;
    tensor<T> res(dim);
    res.coord.reserve(coord.size());
    typename std::vector< monomial<T> >::const_iterator itend=coord.end();
    for (typename std::vector< monomial<T> >::const_iterator it=coord.begin();it!=itend;++it){
      index_t i = it->index.iref();
      T n(i.front()+1);
      i[0]++;
      if (!is_zero(n))
	res.coord.push_back(monomial<T>(rdiv(it->value,n),i));
    }
    return res;
  }

  template<class T>
  void Tfracadd(const tensor<T> & n1, const tensor<T> & d1,const tensor<T> & n2, const tensor<T> & d2, tensor<T> & num, tensor<T> & den){
    // COUT << n1 << "/" << d1 << "+" << n2 << "/" << d2 << "=";
    if (Tis_one(d1)){
      n2.TAdd(n1*d2,num);  //  num=n1*d2+n2;
      den=d2;
      // COUT << num << "/" << den << '\n';
      return;
    }
    if (Tis_one(d2)){
      n1.TAdd(n2*d1,num); // num=n2*d1+n1;
      den=d1;
      // COUT << num << "/" << den << '\n';
      return;
    }
    // n1/d1+n2/d2 with g=gcd(d1,d2), d1=d1g*g, d2=d2g*g is
    // (n1*d2g+n2*d1g)/g * 1/(d1g*d2g)
    tensor<T> d1g(d1),d2g(d2);
    den=simplify(d1g,d2g);
    (n1*d2g).TAdd(n2*d1g,num);
    simplify(num,den);
    den=den*d1g*d2g;

    /* 
       (n1*d2).TAdd(n2*d1,num); //    num=n1*d2+n2*d1;
       den=d1*d2;
       simplify(num,den);
    */
  }

  template <class T>
  void Tfracmul(const tensor<T> & n1, const tensor<T> & d1,const tensor<T> & n2, const tensor<T> & d2, tensor<T> & num, tensor<T> & den){
    // COUT << n1 << "/" << d1 << "*" << n2 << "/" << d2 << "=";
    if (Tis_one(d1)){
      num=n1;
      den=d2;
      simplify(num,den);
      num=num*n2;
      // COUT << num << "/" << den << '\n';
      return;
    }
    if (Tis_one(d2)){
      num=n2;
      den=d1;
      simplify(num,den);
      num=num*n1;
      // COUT << num << "/" << den << '\n';
      return;
    }
    num=n1;
    den=d2;
    simplify(num,den);
    tensor<T> ntemp(n2),dtemp(d1);
    simplify(ntemp,dtemp);
    num=num*ntemp;
    den=den*dtemp;
    // COUT << num << "/" << den << '\n';
  }

  /*
    template<class T>
    std::ostream & operator << (std::ostream & os, const std::vector< facteur<T> > & v){
    std::vector< facteur<T> >::const_iterator itend=v.end();
    for (std::vector< facteur<T> >::const_iterator it=v.begin();;){
    os << *it ;
    ++it;
    if (it==itend)
    break;
    else
    os << "*";
    }
    return(os);
    }
  */

  template<class T>
  T Tproduct(const std::vector< facteur<T> > & v){
    assert (!v.empty());
    typename std::vector< facteur<T> >::const_iterator it=v.begin();
    typename std::vector< facteur<T> >::const_iterator itend=v.end();
    T prod(Tpow(it->fact,it->mult));
    ++it;
    for (;it!=itend;++it)
      prod *= it->mult==1?it->fact:Tpow(it->fact,it->mult); //  prod=prod*Tpow(it->fact,it->mult);
    return prod;
  }

  template<class T>
  T Tproduct(typename std::vector< facteur<T> >::const_iterator it,
	     typename std::vector< facteur<T> >::const_iterator itend){
    assert (it!=itend);
    T prod(Tpow(it->fact,it->mult));
    ++it;
    for (;it!=itend;++it)
      prod *= it->mult==1?it->fact:Tpow(it->fact,it->mult); // prod=prod*Tpow(it->fact,it->mult);
    return prod;
  }

  // square-free factorization of a polynomial
  // T must be a 0 characteristic field
  template<class T>
  std::vector< facteur< tensor<T> > > Tsqff_char0(const tensor<T> &p ){
    tensor<T> y(p.derivative()),w(p);
    tensor<T> c(simplify(w,y));
    // p=p_1*p_2^2*...*p_n^n, c=gcd(p,p')=p_2*...*p_n^(n-1), 
    // w=p/c=pi_i p_i, y=p'/c=sum_{i>=1} ip_i'*pi_{j!=i} p_j
    y.TSub(w.derivative(),y); // y=y-w.derivative(); // y= sum_{i>=2} (i-1) p_i' * pi_{j!=i} p_j
    std::vector< facteur< tensor<T> > > v;
    int k=1; // multiplicity counter
    while(!y.coord.empty()){
      // y=sum_{i>=k+1} (i-k) p_i' * pi_{j!=i, j>=k} p_j
      const tensor<T> & g=simplify(w,y);
      if (!Tis_one(g))
	v.push_back(facteur< tensor<T> >(g,k));
      // this push p_k, now w=pi_{i>=k+1} p_i and 
      // y=sum_{i>=k+1} (i-k) p_i' * pi_{j!=i, j>=k+1} p_j
      y.TSub(w.derivative(),y); // y=y-w.derivative();
      // y=sum_{i>=k+1} (i-(k+1)) p_i' * pi_{j!=i, j>=k+1} p_j
      k++;
    }
    if (!Tis_one(w))
      v.push_back(facteur< tensor<T> >(w,k));
    return std::vector< facteur< tensor<T> > > (v);
  }

  // class for partial fraction decomposition
  template<class T>
  class pf{
  public:
    tensor<T> num;
    tensor<T> fact;
    tensor<T> den;
    int mult; // den=cste*fact^mult
    pf(): num(),fact(),den(),mult(0) {}
    pf(const pf & a) : num(a.num), fact(a.fact), den(a.den), mult(a.mult) {}
    pf(const tensor<T> &n, const tensor<T> & d, const tensor<T> & f,int m) : num(n), fact(f), den(d), mult(m) {};
  };

  template<class T>
  std::ostream & operator << (std::ostream & os, const pf<T> & v){
    os << v.num << "/" << v.den ;
    return os;
  }

#ifdef NSPIRE
  template<class T,class I>
  nio::ios_base<I> & operator << (nio::ios_base<I> & os, const std::vector< pf<T> > & v){
    typename std::vector< pf<T> >::const_iterator itend=v.end();
    for (typename std::vector< pf<T> >::const_iterator it=v.begin();;){
      os << *it ;
      ++it;
      if (it==itend)
	break;
      else
	os << "+";
    }
    return(os);
  }
#else
  template<class T>
  std::ostream & operator << (std::ostream & os, const std::vector< pf<T> > & v){
    typename std::vector< pf<T> >::const_iterator itend=v.end();
    for (typename std::vector< pf<T> >::const_iterator it=v.begin();;){
      os << *it ;
      ++it;
      if (it==itend)
	break;
      else
	os << "+";
    }
    return(os);
  }
#endif

  template<class T>
  tensor<T> TsimplifybyTlgcd(tensor<T>& a,tensor<T> &b){
    const tensor<T> & Tlgcdg=gcd(Tlgcd(a),Tlgcd(b));
    if (!Tis_one(Tlgcdg)){
      a=a/Tlgcdg;
      b=b/Tlgcdg;
    }
    return Tlgcdg;
  }

  // utility for Tpartfrac, see below
  template<class T>
  void Tpartfrac(const tensor<T> & num, const tensor<T> & den, /* const tensor<T> & dendiff ,*/const std::vector< facteur< tensor<T> > > & w , int n, int m, std::vector < pf <T> > & pfdecomp ){
    if (m==n)
      return;
    if (m-n==1){
      tensor<T> nums(num), dens(den);
      TsimplifybyTlgcd(nums,dens);
      pfdecomp.push_back(pf<T>(nums,dens,w[n].fact,w[n].mult));    
      return ;
    }
    typename std::vector< facteur< tensor<T> > >::const_iterator it=w.begin()+n; // &w[n];
    typename std::vector< facteur< tensor<T> > >::const_iterator it_end=w.begin()+m; // &w[m];
    /*
    // check if all factors of degree 1 and mult 1
    for (;it!=itend;++it){
      if (it->mult!=1 || it->fact.lexsorted_degree()!=1)
	break;
    }
    if (it==itend){
      // add rem(num,it->fact)/rem(dendiff,it->fact) / it->fact for all factors
      it=w.begin()+n;
      for (;it!=itend;++it){
	
      }
    }
    */
    // split v in 2 parts, then apply recursively Tpartfrac on each part
    it=w.begin()+n;
    int p=(m+n)/2;
    typename std::vector< facteur< tensor<T> > >::const_iterator it_milieu=w.begin()+p; // &w[p];
    const tensor<T> & fn=Tproduct< tensor<T> >(it,it_milieu);
    const tensor<T> & fm=Tproduct< tensor<T> >(it_milieu,it_end);
    // write C*num=u*fn+v*fm
    tensor<T> C(den.dim),u(den.dim),v(den.dim);
    Tabcuv(fn,fm,num,u,v,C);
    C=C*(den/(fn*fm));
    // num/den= (u*fn+v*fm)/(C*fn*fm)=u/(C*fm)+v/(C*fn)
    Tpartfrac(v,C*fn,w,n,p,pfdecomp);
    Tpartfrac(u,C*fm,w,p,m,pfdecomp);
  }

  // if num/den with den=cste*pi_i d_i^{mult_i}
  // Tpartfrac rewrites num/den as
  // sum_{i} pfi.num/pfi.den with pfi.den = cste* pfi.den^pfi.mult
  // num and den are assumed to be prime together
  template<class T>
  void Tpartfrac(const tensor<T> & num, const tensor<T> & den, const std::vector< facteur< tensor<T> > > & v , std::vector < pf <T> > & pfdecomp, tensor<T> & ipnum, tensor<T> & ipden ){
    int n=int(v.size());
    pfdecomp.reserve(n);
    // compute ip and call Tpartfrac
    tensor<T> rem(num.dim);
    num.TPseudoDivRem(den,ipnum,rem,ipden);
    // ipden*num=den*ipnum+rem hence num/den=ipnum/ipden+rem/(ipden*den)
    const tensor<T> & temp=ipden*den;
    // simplify(rem,temp);
    if (n==1)
      pfdecomp.push_back(pf<T>(rem,temp,v.front().fact,v.front().mult));
    else
      Tpartfrac(rem,temp,/* temp.derivative() ,*/ v,0,n,pfdecomp);
  }

  // reduction of a fraction with multiple poles to single poles by integration
  // by part, use the relation
  // ilaplace(P'/P^(k+1))=laplacevar/k*ilaplace(1/P^k)
  template<class T>
  pf<T> Tlaplace_reduce_pf(const pf<T> & p_cst, tensor<T> & laplacevar ){
    pf<T> p(p_cst);
    assert(p.mult>0);
    if (p.mult==1)
      return p_cst;
    tensor<T> fprime=p.fact.derivative();
    tensor<T> d(fprime.dim),C(fprime.dim),u(fprime.dim),v(fprime.dim);
    Tegcdpsr(p.fact,fprime,u,v,d); // f*u+f'*v=d
    tensor<T> usave(u),vsave(v);
    // int initial_mult=p.mult-1;
    while (p.mult>1){
      Tegcdtoabcuv(p.fact,fprime,p.num,u,v,d,C);
      p.mult--;
      p.den=(p.den/p.fact)*C*T(p.mult);
      p.num=u*T(p.mult)+v.derivative()+v*laplacevar;
      if ( (p.mult % 5)==1) // simplify from time to time
	TsimplifybyTlgcd(p.num,p.den);
      if (p.mult==1)
	break;
      u=usave;
      v=vsave;
    }
    return pf<T>(p);
  }

  // reduction of a fraction with multiple poles to single poles by integration
  // by part, the integrated part is added to intdecomp
  template<class T>
  pf<T> Tintreduce_pf(const pf<T> & p_cst, std::vector< pf<T> > & intdecomp ){
    assert(p_cst.mult>0);
    if (p_cst.mult==1)
      return p_cst;
    pf<T> p(p_cst);
    tensor<T> fprime=p.fact.derivative();
    tensor<T> d(fprime.dim),u(fprime.dim),v(fprime.dim),C(fprime.dim);
    tensor<T> resnum(fprime.dim),resden(T(1),fprime.dim),numtemp(fprime.dim),dentemp(fprime.dim);
    Tegcdpsr(p.fact,fprime,u,v,d); // f*u+f'*v=d
    tensor<T> usave(u),vsave(v);
    int initial_mult=p.mult-1;
    while (p.mult>1){
      Tegcdtoabcuv(p.fact,fprime,p.num,u,v,d,C);
      p.mult--;
      p.den=(p.den/p.fact)*C*T(p.mult);
      p.num=u*T(p.mult)+v.derivative();
      u=-p.den;
      TsimplifybyTlgcd(u,v);
      Tfracadd<T>(resnum,resden,v,u,numtemp,dentemp);
      resnum=numtemp;
      resden=dentemp;
      if ( (p.mult % 5)==1) // simplify from time to time
	TsimplifybyTlgcd(p.num,p.den);
      if (p.mult==1)
	break;
      u=usave;
      v=vsave;
    }
    intdecomp.push_back(pf<T>(resnum,resden,p.fact,initial_mult));
    return pf<T>(p);
  }

  template<class T>
  std::vector<T> makevector(const T & a, const T &b){
    std::vector<T> v;
    v.push_back(a);
    v.push_back(b);
    return v;
  }

  template<class T>
  std::vector<T> makevector(const T & a, const T &b,const T & c){
    std::vector<T> v;
    v.push_back(a);
    v.push_back(b);
    v.push_back(c);
    return v;
  }

  template<class T>
  std::vector<T> makevector(const T & a, const T &b,const T & c,const T & d){
    std::vector<T> v;
    v.push_back(a);
    v.push_back(b);
    v.push_back(c);
    v.push_back(d);
    return v;
  }

  template<class T>
  std::vector<T> merge(const std::vector<T> & a, const std::vector<T> &b){
    std::vector<T> v(a);
    int as=a.size();
    int bs=b.size();
    v.reserve(as+bs);
    typename std::vector<T>::const_iterator it=b.begin(),itend=b.end();
    for (;it!=itend;++it)
      v.push_back(*it);
    return v;
  }

#ifndef NO_NAMESPACE_GIAC
} // namespace giac
#endif // ndef NO_NAMESPACE_GIAC


#endif // ndef _GIAC_POLY_H