multigrobner.H

#ifndef _MULTIGROBNER_H_
#define _MULTIGROBNER_H_

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


template <class K>
multiPolynomial<K> division_remainder_with_coefficients(const multiPolynomial<K>& f, 
                                                        const list< multiPolynomial<K> >& F,
                                                        Matrix< Polynomial<K> >& U){
  multiPolynomial<K> h,r;
  Polynomial<K> temp;
  K c;
  multiPolynomial<K> a,b;
  typename list< multiPolynomial<K> >::const_iterator it;
  long i;

  //setup
  h= f;
  r= f;
  r.sets_to_zero(); //for the alphabet...
  //find the right size for U, and set it to 0
  for(it= F.begin(), i=0; it != F.end(); it++, i++);
  U.resize(1,i);
  for(i=0; i< U.size(); i++){
    U[i].sets_to_zero();
    U[i].alphabet= (Alphabet *) f[0].alphabet;
  }

  temp.alphabet= (Alphabet *) f[0].alphabet;
  //good, now let us start
  while( !h.is_zero() ){

    a=h.lm(); 
    for(it= F.begin(), i=0; it != F.end(); it++, i++){
      b= it->lm();
      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);
      temp= a / b;
      temp*= c;
      //temp= Polynomial<K>( a ) * c;

      U[i] += temp;
      h += (*it) * temp * (-1);
    }     
    else{
      r += a * h.coeff_of(a);
      h += a * h.coeff_of(a)*(-1);
    }
    
  }

  //  cout << "U= " << U << endl;

  return r;
}

template <class K>
multiPolynomial<K> division_remainder(const multiPolynomial<K>& f, 
                                      const list< multiPolynomial<K> >& F){
  multiPolynomial<K> h,r;
  Polynomial<K> temp;
  K c;
  multiPolynomial<K> a,b;
  typename list< multiPolynomial<K> >::const_iterator it;
  long i;

  //setup
  h= f;
  r= f;
  r.sets_to_zero(); //for the alphabet...
  temp.alphabet= (Alphabet *) f[0].alphabet;
  //good, now let us start
  while( !h.is_zero() ){

    a=h.lm(); 
    for(it= F.begin(), i=0; it != F.end(); it++, i++){
      b= it->lm();
      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);
      temp= a / b;
      temp*= c;
      //temp= Polynomial<K>( a ) * c;

      h += (*it) * temp * (-1);
    }     
    else{
      r += a * h.coeff_of(a);
      h += a * h.coeff_of(a)*(-1);
    }
    
  }

  return r;
}




template <class K>
class multiIdeal {
private:

  list< Matrix< Polynomial<K> > > coeffs;

  void insert(const multiPolynomial<K>& mP, 
              const Matrix< Polynomial<K> >&coords){
    long i;
    typename list< Matrix< Polynomial<K> > >::iterator it;
    Matrix< Polynomial<K> > M, temp;
    
    //the following is a matrix multiplication really, 
    //except that one of the matrices presents itself 
    //as a list of lines
    M= (*coeffs.begin()) * coords[0];

    for(it= coeffs.begin(), it++, i=1; i< coords.size(); it++, i++){
      M += *it * coords[i];
    }

    coeffs.push_back(M);
    generators.push_back(mP);
  }

  long simplify(){
    long n= coeffs.begin()->ncols();
    typename list< multiPolynomial<K> >::iterator it, itprime, gen_begin;
    typename list< Matrix< Polynomial<K> > >::iterator itcoeff, c_begin;
    long ans= 0;

    gen_begin= generators.begin();
    c_begin= coeffs.begin();
    while(n > 0){
      gen_begin++;
      c_begin++;
      n--;
    }
    //good -- we have moved gen_begin and c_begin
    //to the beginning of the generators added by
    //grobnerize_with_coefficients. The first generators
    //we prefer to keep.
    for(it= gen_begin, itcoeff= c_begin;
        it!= generators.end();){

      for(itprime= generators.begin();
          itprime!= generators.end();
          itprime++)
        if( (itprime != it) && itprime->lm().divides( it->lm() ) )
          break;
      
      if( itprime != generators.end() ){ // found one ?
        it= generators.erase( it );
        itcoeff= coeffs.erase( itcoeff );
        ans++;
      }
      else{
        it++;
        itcoeff++;
      }
    }//for it
    return ans;
  }//simplfy




public:
  list< multiPolynomial<K> > generators;






  void grobnerize_with_coefficients(Matrix< Polynomial<K> >& M){
    list< pair< multiPolynomial<K>*, multiPolynomial<K>* > > couples;
    pair< multiPolynomial<K>*, multiPolynomial<K>* > fg;
    pair< Polynomial<K>, Polynomial<K> > S_coeffs;
    multiPolynomial<K> Spoly, h;
    Polynomial<K> I;
    typename list< multiPolynomial<K> >::iterator it, itprime;
    Matrix< Polynomial<K> > U, tuple;
    //n= nb of original generators, 
    //m= nb of generators added
    long i, j, n, m;
    long S_count= 0;
    


    if(generators.begin() == generators.end())
      return;
    //I= one
    I= Polynomial<K>( PowProd(0) );
    I.alphabet= generators.begin()->coords[0].alphabet;
    
    //write the beginning of the matrix, as a list of 1-line
    //matrices. It starts with the identity.

    n= 0;
    tuple.resize(1, generators.size() );
    for(i=0; i< (long) generators.size(); i++){
      tuple[i].sets_to_zero();
      tuple[i].alphabet= I.alphabet;
    }
    for(it= generators.begin(); it!= generators.end(); it++){
      tuple[n]= I;
      if(n>0)
        tuple[n-1].sets_to_zero();
      n++;      
      coeffs.push_back( tuple );
    }
    //n is now set; set m to zero
    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_back(   make_pair( &(*it) , &(*itprime) )     );
    
    
    while( !couples.empty() ){
      S_count++;
      
      // we pick the couples alternatively from the beginning 
      // and the end of the list -- this in order to avoid
      // the worst-case scenario when one of the two choices
      // gives rise to an exponentially growing computation.
      // On one example this has reduced a computation time
      // from 15 minutes to 5 seconds.

      if(S_count % 2 == 0){
        fg= couples.front(); // take (f,g) from the stack
        couples.pop_front();
      }
      else{
        fg= couples.back();
        couples.pop_back();
      }


      Spoly= S_polynomial_with_coefficients( *fg.first , *fg.second, S_coeffs ); 
      //computes S(f,g)
      

      if( S_count % 1000 == 0 )
        cout << S_count << " S polynomials computed, " 
             << 100*(S_count)*2/((n+m)*(n+m)) 
             << "% done"  << endl;
      
      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 we now multiply by -1:
      U*= I*(-1);

      
      if( !h.is_zero() ){
        
        //build (coords of h)
        //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;
        U[j]+= S_coeffs.second;
        //add it
        //generators.push_back(h); //old !
        insert(h, U); 
        //add one to m
        m++;

        //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_back(  make_pair(  &(*it), &(*itprime)  )  )  ;
        
      }
      
    }// end of while loop
    
    cout << m << " polynomials added" << endl;  

    cout << simplify() << " removed" << endl;

    M.get_data_from_matrix_list(coeffs);
    coeffs.clear();
//     if( m==0 ) { //nothing added?
//       M.get_data_from_matrix_list(coeffs);
//       return;
//     }
//     //else M would not be square
//     coeffs.begin()->resize(1,n+m); //this remedies it!
//     //M will be square of size n + m now:
//     M.get_data_from_matrix_list( coeffs );
//     //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< multiPolynomial<K> >::const_iterator it, itprime;
    list< Matrix< Polynomial<K> > > answer;
    Matrix< Polynomial<K> > tuple;
    pair< Polynomial<K>, Polynomial<K> > S_coeffs;
    multiPolynomial<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_matrix_list(answer);
  }
  //   everything below is commented out -- a bunch of methods
  // from the Ideal class which I may or may not adapt to the
  // multiIdeal case.
  

//   /// turns a minimal grobner basis into a reduced + minimal grobner basis
//   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
//   }
//   /// checks if a variable shows up among the ideal generators
//   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

---