grobner.H

#ifndef _GROBNER_H_
#define _GROBNER_H_

#include "polynomials.H"
#include "matrix.H"
#include <list>
#include <stack>
#include <utility>

/*
 * In one variable for example, this is the remainder in the Euclidean
 * long division (of f by the one element in the list).
 */
template <class K>
Polynomial<K> division_remainder(const Polynomial<K>& f, const list< Polynomial<K> >& F){
  Polynomial<K> h,r,temp;
  K c;
  PowProd a,b;
  typename list< Polynomial<K> >::const_iterator it;
  
  h= f;
  
  while( !h.is_zero() ){

    a=h.lp(); 
    for(it= F.begin(); it != F.end(); it++){
      b= it->lp();
      if(b.divides(a))
        break;
    }
    if(it != F.end()){ // lp(h) was divisible by some lp(f_i) ?

      c= h.coeff_of(a);
      c/= it->coeff_of(b);
      c *= -1;
      a /= b;
      temp= Polynomial<K>( a );
      temp *= c;
      temp *= (*it);
      h += temp;
    }     
    else{

      temp= h.lt();      
      r += temp;
      temp *= -1;
      h += temp;
    }
    
  }
  r.alphabet= (Alphabet *) f.alphabet;
  return r;
}

template <class K>
Polynomial<K> division_remainder_with_coefficients(const Polynomial<K>& f, 
                                                   const list< Polynomial<K> >& F,
                                                   Tuple< Polynomial<K> >& U){
  Polynomial<K> h,r,temp;
  K c;
  PowProd a,b;
  typename list< Polynomial<K> >::const_iterator it;
  long i;
  
  h= f;
  //find the right size for U, and set it to 0
  for(it= F.begin(), i=0; it != F.end(); it++, i++);
  U.resize(i);
  for(i=0; i< U.size(); i++){
    U[i].sets_to_zero();
    U[i].alphabet= (Alphabet *) f.alphabet;
  }
  
  while( !h.is_zero() ){

    a=h.lp(); 
    for(it= F.begin(), i=0; it != F.end(); it++, i++){
      b= it->lp();
      if(b.divides(a))
        break;
    }
    if(it != F.end()){ // lp(h) was divisible by some lp(f_i) ?

      c= h.coeff_of(a);
      c/= it->coeff_of(b);
      a /= b;
      temp= Polynomial<K>( a ) * c;
      U[i] += temp;
      temp *= (*it) * (-1);
      h += temp;
    }     
    else{

      temp= h.lt();      
      r += temp;
      temp *= -1;
      h += temp;
    }
    
  }
  r.alphabet= (Alphabet *) f.alphabet;
  return r;
}



template <class K>
class Ideal {
private:
  list< Polynomial<K> > G, R, P;
  list< pair< Polynomial<K>*, Polynomial<K>* > > B;
  
  void filter_B(Polynomial<K> *ptr)
  {
    typename list< pair< Polynomial<K>*, Polynomial<K>* > >::iterator it;
    
    for(it= B.begin(); it!= B.end();)
      if( (it->first == ptr) || (it->second == ptr) )
        it= B.erase(it);
      else
        it++;
    
  }
  
  void reduce_all()
  {
    Polynomial<K> h;
    typename list< Polynomial<K> >::iterator it;
    
      
    while( !R.empty() ){

      h= R.front();
      R.pop_front();
      if( !P.empty() ){ 
        // P empty only happens the 1st time this method is called
        //(from grobnerize()). In this case G is empty also.
        it= P.begin();
        G.splice(G.end(), P); // G= G union P; P emptied
        h= division_remainder(h, G);
        P.splice(P.begin(), G, it, G.end() ); // G and P back to normal
      }
      

      if( !h.is_zero() ){

                

        for(it= G.begin(); it!= G.end();){
         
          if( h.lp().divides( it->lp() ) ){
            filter_B( &(*it) );
            R.push_back( *it );
            it= G.erase(it);
          }
          else
            it++;
        }
        

        for(it= P.begin(); it!= P.end();){
          
          
          if( h.lp().divides( it->lp() ) ){
            R.push_back( *it );
            it= P.erase(it);
          }
          else
            it++;
        }
        


        P.push_back(h);
        
        
      }//if h not zero
      
    }//while
    
      
  }//reduce_all
/*
 * This adds to B all the pairs of the form
 * (g, p) with g in G union P    (g != p)
 *             p in P 
 * 
 */
  void populate_B()
  {
    typename list< Polynomial<K> >::iterator itG, itP, book;
    
    
    for(itG= G.begin(); itG!= G.end(); itG++)
      for(itP= P.begin(); itP!= P.end(); itP++)
        B.push_back( make_pair( &(*itG), &(*itP) ) );
    
    for(itG= P.begin(); itG!= P.end(); itG++)
      for(itP= itG, itP++; itP!= P.end(); itP++)
        B.push_back( make_pair( &(*itG), &(*itP) ) );
          
  }
  

  void new_basis()
  {
    list< Polynomial<K> > H;
    typename list< Polynomial<K> >::iterator itG, itH;
    Polynomial<K> poly;
    
  
    populate_B();
    G.splice(G.end(), P);// G= G union P; P emptied
    H= G;
    
    for(itG= G.begin(), itH= H.begin(); itG!= G.end(); itG++){
      itH= H.erase(itH);
      poly= division_remainder( *itG, H);
      H.push_front(*itG);
      *itG= poly;
    }

      
  }
  
  

public:
  list< Polynomial<K> > generators;

  void grobnerize()
  {
    pair< Polynomial<K>*, Polynomial<K>* > fg;
    Polynomial<K> Spoly, h;
    typename list< Polynomial<K> >::iterator it;
    

    //init
    R.clear();
    R.splice( R.begin(), generators );
    P.clear();
    G.clear();
    B.clear();
    reduce_all();
    new_basis();
    
    while( !B.empty() ){

      fg=B.front();
      B.pop_front();
      
      if( fg.first->lp().prime_to( fg.second->lp() ) )
        h.sets_to_zero();
      else{
        Spoly= S_polynomial( *fg.first , *fg.second ); //computes S(f,g)
        h= division_remainder(Spoly, G); //S(f,g) -> h by G
      }
      
      if( !h.is_zero() ){
        P.clear();
        P.push_back( h );
        R.clear();
        
        for(it= G.begin(); it!= G.end();)
          if(h.lp().divides( it->lp() ) ) {
            R.push_back( *it );
            filter_B( &(*it) );
            it= G.erase(it);
          }
          else
            it++;
          
        reduce_all();
        new_basis();
      }// if h non zero
    }//main while loop
    
    //finally:
    generators.splice( generators.begin(), G);
    
  }
  








  void grobnerize_with_coefficients(Matrix< Polynomial<K> >& M){
    stack< pair< Polynomial<K>*, Polynomial<K>* > > couples;
    pair< Polynomial<K>*, Polynomial<K>* > fg;
    pair< Polynomial<K>, Polynomial<K> > S_coeffs;
    Polynomial<K> Spoly, h, I;
    typename list< Polynomial<K> >::iterator it, itprime;
    Tuple< Polynomial<K> > U, tuple;
    list< Tuple< Polynomial<K> > > coeffs;
    long i, j, n, m;
    
    if(generators.begin() == generators.end())
      return;
    //I= one
    I= Polynomial<K>( PowProd(0) );
    I.alphabet= generators.begin()->alphabet;
    
    //write the beginning of the matrix, as a list of tuples
    //we write it as -identity, and at the very end we will
    //multiply everything by -1 again. This is because of
    //division_remainder_with_coefficients which returns 
    //the opposite of the coeffs we need.
    n= 0;
    for(it= generators.begin(); it!= generators.end(); it++){
      n++;
      tuple= Tuple< Polynomial<K> >(n);
      tuple[n-1]= I*(-1);
      coeffs.push_back( tuple );
    }

    m= 0;
    // populates the stack with the addresses of all pairs
    // of polynomials from the list
    
    for(it= generators.begin(); it !=generators.end(); it++)
      for(itprime= it, itprime++; itprime != generators.end(); itprime++)
        couples.push(   make_pair( &(*it) , &(*itprime) )     );
    
    
    while( !couples.empty() ){
      
      fg= couples.top(); // take (f,g) from the stack
      couples.pop();

      if( fg.first->lp().prime_to( fg.second->lp() ) )
        h.sets_to_zero();
      else{
        Spoly= S_polynomial_with_coefficients( *fg.first , *fg.second, S_coeffs ); 
        //computes S(f,g)
        h= division_remainder_with_coefficients(Spoly, generators, U); 
        //S(f,g) -> h by F, so S(f,g)=h + sum U[i]F[i]
        //so h = S(f,g) - sum U[i]*F[i], hence all the (-1) signs
      }
      
      if( !h.is_zero() ){
        //add it
        generators.push_back(h);
        //updates the stack of couples
        itprime= generators.end();
        itprime--; //itprime points at the last element, ie h
        for(it= generators.begin(); it != itprime; it++)
          couples.push(  make_pair(  &(*it), &(*itprime)  )  )  ;
        //build -(coords of h). In the end we will multiply by (-1)
        //set i and j...
        for(i=0, it=generators.begin(); it!= generators.end(); i++, it++)
          if( &(*it) == fg.first )
            break;
        for(j=0, it=generators.begin(); it!= generators.end(); j++, it++)
          if( &(*it) == fg.second )
            break;
        U[i]+= S_coeffs.first * (-1);
        U[j]+= S_coeffs.second * (-1);
        coeffs.push_back( U );
        m++;
      }
      
    }// end of while loop
    
  
    if( m==0 ) { //nothing added?
      M.get_data_from_tuple_list(coeffs);
      M*= I*(-1);
      return;
    }
    //else M would not be square
    coeffs.begin()->resize(n+m); //this remedies it!
    //M will be square of size n + m now:
    M.get_data_from_tuple_list( coeffs );
    M*= I*(-1); //finally ! signs corrected !
    //now compute M**m, it will have m columns of 0
    //on the right
    Matrix< Polynomial<K> > N;
    N= M;
    i=m-1;
    while(i>0){
      M*= N;
      i--;
    }
    //get rid of the columns of 0
    M.resize(n+m, n);
    //tada !!
  }
  

  void syzygy(Matrix< Polynomial<K> >& M) const {
    typename list< Polynomial<K> >::const_iterator it, itprime;
    list< Tuple< Polynomial<K> > > answer;
    Tuple< Polynomial<K> > tuple;
    pair< Polynomial<K>, Polynomial<K> > S_coeffs;
    Polynomial<K> S;
    long i, j;

    for(it= generators.begin(), i=0; it!=generators.end(); it++, i++)
      for(itprime= it, itprime++, j=i+1; 
          itprime!= generators.end(); 
          itprime++, j++){
        S= S_polynomial_with_coefficients(*it, *itprime, S_coeffs);
        division_remainder_with_coefficients(S, generators, tuple);
        tuple[i] += S_coeffs.first * (-1);
        tuple[j] += S_coeffs.second * (-1);
        answer.push_back(tuple);
      }

    M.get_data_from_tuple_list(answer);
  }




  void minimalize(){
    typename list< Polynomial<K> >::iterator it, itprime;
    PowProd a,b;
    K c;
    
    for(it= generators.begin(); it!= generators.end(); it++){
      itprime= generators.begin();
      while(itprime != generators.end()){
        
        if(it == itprime)
          itprime++;
        else{
          a=it->lp();
          b=itprime->lp();
          if(a.divides(b))
            itprime= generators.erase(itprime);
          else
            itprime++;
        }//end of else
        
      }//end of while
    }//end of for
    
    for(it= generators.begin(); it!= generators.end(); it++){
      c= it->lc();
      c= 1/c;
      *it *= c; // we divide *it by its leading coefficient
    }

  }

  void reduce(){
    typename list< Polynomial<K> >::iterator it;
    Polynomial<K> f;
    
    it= generators.begin();
    while( it != generators.end() ){
      f= *it;
      it= generators.erase(it);
      f= division_remainder(f,generators);
      generators.push_front(f);
    }//end of while
  }
  bool variable_shows_up(long i) const {
    typename list< Polynomial<K> >::const_iterator it;

    for(it= generators.begin(); it!= generators.end(); it++)
      if( it->involves_variable(i) )
        return true;

    return false;
  }


  friend ostream& operator<<(ostream& os, const Ideal<K>& J){
    typename list< Polynomial<K> >::const_iterator it;

    for(it= J.generators.begin(); it != J.generators.end(); it++)
      os << *it << endl;
    return os;
  }


};

#endif

---