multipoly.H

#ifndef _MULTIPOLY_H_
#define _MULTIPOLY_H_

#include "polynomials.H"
#include "tables.H"
using namespace std;


template <class K>
class multiPolynomial {
public:
  // note that a Matrix is a Tuple, so we can
  // use coords[i] instead of coords(0,i) for
  // example.
  Matrix< Polynomial<K> > coords;
  multiPolynomial(){
    coords.set_sep('|');
  }
  
  void resize(long n) {

    coords.resize(1,n);
  }

  long size() const {

    return coords.size();
  }

  void set_alphabet(const Alphabet *alpha) {
    long i;

    for(i= 0; i< coords.size(); i++)
      coords[i].alphabet= (Alphabet *) alpha;
  }


  //is zero or not ?
  bool is_zero() const {
    long i;

    for(i=0; i< coords.size(); i++)
      if(!coords[i].is_zero())
        return false;

    return true;
  }

  void sets_to_zero() {
    long i;

    for(i=0; i< coords.size(); i++)
      coords[i].sets_to_zero();
  }




  multiPolynomial<K> lm() const {
    long i;
    PowProd max(0), other;
    long j= -1;
    multiPolynomial<K> ans;

    ans.resize( coords.size() );
    ans.set_alphabet((Alphabet *) coords[0].alphabet );
    
    for(i= 0; i< coords.size(); i++){
      if( !coords[i].is_zero() ){
          other= coords[i].lp();
          if( max < other || (j==-1 )){
            max= other;
            j= i;
          }
        }
    }

    if(j != -1)
      ans[j].add_term(max, 1);
    
    return ans;
  }

  PowProd lp() const {
    long i;
    PowProd max(0), other;
    long j= -1;
    
    for(i= 0; i< coords.size(); i++){
      if( !coords[i].is_zero() ){
          other= coords[i].lp();
          if( max < other || (j==-1 )){
            max= other;
            j= i;
          }
        }
    }

    return max;

  }

  K lc() const {
    long i;
    PowProd max(0), other;
    long j= -1;
    
    for(i= 0; i< coords.size(); i++){
      if( !coords[i].is_zero() ){
          other= coords[i].lp();
          if( max < other || (j==-1 )){
            max= other;
            j= i;
          }
        }
    }
    if( j== -1 )                // was all 0 ?
      return 0;

    return coords[j].coeff_of(max);
   }




  multiPolynomial<K> lm_POT() const {
    long i;
    multiPolynomial<K> ans;

    ans.coords= coords; //this resizes and sets the alphabets
    ans.sets_to_zero();

    for(i=0; i < coords.size(); i++)
      if( !coords[i].is_zero() )
        break;

    if(i< coords.size()) //found a nonzero poly ?
      ans[i].add_term(coords[i].lp(), 1);
    
    return ans;
  }

  PowProd lp_POT() const {
    long i;

    for(i=0; i< coords.size(); i++)
      if( !coords[i].is_zero() )
        return coords[i].lp();

    return PowProd(0);
  }


  K lc_POT() const {
    long i;

    for(i=0; i< coords.size(); i++)
      if( !coords[i].is_zero() )
        return coords[i].lc();

    return 1;
  }


  K coeff_of(const multiPolynomial<K>& a) const {
    long i;
      
    for(i=0; i< a.size(); i++)
      if(!a[i].is_zero())
        break;
    if( i == a.size() ) //zero ?
      return 0;
    else
      return coords[i].coeff_of( a[i].lp() );
  }



  //operators
  Polynomial<K>& operator[](long n) {
    return coords(0,n);
  }

  const Polynomial<K>& operator[](long n) const {
    return coords(0,n);
  }

  multiPolynomial<K>& operator+=(const multiPolynomial<K>& that){

    coords+= that.coords;
    return *this;
  }

  multiPolynomial<K> operator+(const multiPolynomial<K>& that) const {
    multiPolynomial<K> ans;
    
    ans= *this;
    ans+= that;
    return ans;
  }


  multiPolynomial<K>& operator*=(const Polynomial<K>& P) {

    coords*= P;
    return *this;
  }

  multiPolynomial<K> operator*(const Polynomial<K>& P) const {
    multiPolynomial ans;

    ans= *this;
    ans*= P;
    return ans;
  }

  multiPolynomial<K>& operator*=(const K& c) {
    long i;

    for(i=0; i< coords.size(); i++)
      coords[i] *= c;

    return *this;
  }

  multiPolynomial<K> operator*(const K& c) const {
    multiPolynomial ans;

    ans= *this;
    ans*= c;
    return ans;
  }


  bool divides(const multiPolynomial<K>& that) const {
    long i;

    if(this->size() != that.size())
      return false;
    for(i=0; i< coords.size(); i++)
      if(!coords[i].is_zero())
        break;
    if(i== coords.size() || that.coords[i].is_zero())
      return false;
    return coords[i].lp().divides( that.coords[i].lp() );
  }

  Polynomial<K> operator/(const multiPolynomial<K>& that) const {
    long i;
    PowProd temp;
    Polynomial<K> dummy;
    dummy.sets_to_zero();

    if(this->size() != that.size())
      return dummy;
    for(i=0; i< coords.size(); i++)
      if(!coords[i].is_zero())
        break;
    if(i== coords.size() || that.coords[i].is_zero())
      return dummy;
  
    temp= coords[i].lp();
    temp /= that.coords[i].lp();
    dummy= Polynomial<K>( temp );
    dummy.alphabet= coords[i].alphabet;
 
    return dummy;
  }

  friend ostream& operator<<(ostream& os, const multiPolynomial<K>& M){

    os << M.coords;
    return os;
  }

  friend multiPolynomial<K> S_polynomial_with_coefficients(const multiPolynomial<K>& f, 
                                                           const multiPolynomial<K>& g,
                                                           pair< Polynomial<K>, Polynomial<K> >& coeffs){
    PowProd a,b,c;
    K lambda, mu;
    multiPolynomial<K> A, B, ans;
    long i, j;

    A= f.lm();
    B= g.lm();
    //do A and B have their non zero entries
    //in the same position ?
    for(i= 0; i < A.size(); i++)
      if( !A[i].is_zero() )
        break;
    for(j= 0; j < B.size(); j++)
      if( !B[j].is_zero() )
        break;
    
    if( (i!=j) || (i==A.size()) || (j==B.size()) ){
      coeffs.first.sets_to_zero();
      coeffs.second.sets_to_zero();
    }
    else{
      a= A.lp();
      b= B.lp();
      c= lcm(a,b);
      lambda= 1 / f.lc();
      mu= -1 / g.lc();
      a= c/a; 
      b= c/b; 
      coeffs.first= Polynomial<K>(a) * lambda;    
      coeffs.second= Polynomial<K>(b) * mu;
    }

       
    // ans = f*a*lambda + g*b*mu
    coeffs.first.alphabet= f[0].alphabet;
    coeffs.second.alphabet= f[0].alphabet;
    ans= f * coeffs.first + g * coeffs.second; 
        
    return ans;
 }

  friend multiPolynomial<K> S_polynomial(const multiPolynomial<K>& f, 
                                         const multiPolynomial<K>& g){
    pair< Polynomial<K>, Polynomial<K> > coeffs;
    PowProd a,b,c;
    K lambda, mu;
    multiPolynomial<K> A, B, ans;
    long i, j;

    A= f.lm();
    B= g.lm();
    //do A and B have their non zero entries
    //in the same position ?
    for(i= 0; i < A.size(); i++)
      if( !A[i].is_zero() )
        break;
    for(j= 0; j < B.size(); j++)
      if( !B[j].is_zero() )
        break;
    
    if( (i!=j) || (i==A.size()) || (j==B.size()) ){
      coeffs.first.sets_to_zero();
      coeffs.second.sets_to_zero();
    }
    else{
      a= A.lp();
      b= B.lp();
      c= lcm(a,b);
      lambda= 1 / f.lc();
      mu= -1 / g.lc();
      a= c/a; 
      b= c/b; 
      coeffs.first= Polynomial<K>(a) * lambda;    
      coeffs.second= Polynomial<K>(b) * mu;
    }

       
    // ans = f*a*lambda + g*b*mu
    coeffs.first.alphabet= f[0].alphabet;
    coeffs.second.alphabet= f[0].alphabet;
    ans= f * coeffs.first + g * coeffs.second; 
        
    return ans;
                                         

 }



};

#endif

---