homfinder.cpp

#include "homfinder.H"
using namespace std;

HomFinder::HomFinder(){
}

HomFinder::HomFinder(GradedAlgebra<F_2> *thesource, GradedAlgebra<F_2> *thetarget){
  source= thesource;
  target= thetarget;
}

void HomFinder::read_solutions_from(ifstream& file)
{
  long n;
  AffineHomomorphism<F_2> phi(source, target);
  
  solutions.clear();
  file >> n;
  //  cout << n << " homomorphisms to read" << endl;
  solution_count= n;
  
  while(n>0){
    file >> phi;
    solutions.push_back(phi);
    n--;
  }
  
} 



void HomFinder::write_solutions_to(ofstream& file)
{
  list< AffineHomomorphism<F_2> >::const_iterator ith;
  AffineHomomorphism<F_2> phi(source, target);
  
  file << solution_count << " ";
  
  for(ith= solutions.begin();
      ith!= solutions.end();
      ith++){
    phi= *ith;
    file << phi;
  }
}


ZZ HomFinder::populate_GL(long n){
  ZZ total, browse, bit, count;
  long i;
  Matrix<F_2> M(n, n);
  
  GL.clear();
  count= 0;
  
  // set total= 2^(n*n)
  total= 1;
  for(i=1; i <= n*n; i++)
    total*= 2;

  for(browse= 0; browse < total; browse++){

    bit= 1;
    for(i=0; i< n*n; i++){
      if( (browse & bit) != 0 )
        M[i]=1;
      else
        M[i]=0;
      bit *= 2;
    }

    if(M.is_invertible()){
      count++;  
      GL.push_back(M);
    }

  }

  //  cout << count << " matrices to consider" << endl;
  return count;
}

ZZ HomFinder::filter_GL()
{
  list< Matrix<F_2> >::iterator M, N;
  Matrix<F_2> mat;
  ZZ count=0;
  AffineHomomorphism<F_2> phi(&basic_source, target);
  

  for(M= GL.begin(); M!= GL.end();){
    
    set_hom_from_matrix(phi, basic_variables, *M);

    if(phi.is_well_defined()){
      M++;
      count++;
    }
    else
      M= GL.erase(M);
    
  }// matrix M

  //now perform the second reduction
  long m=1;
  
   for(M= GL.begin(); M!= GL.end(); m++, M++)
     for(N=M, N++; N!=GL.end();){
       mat= *N * M->inverse();
       phi.make_identity(target);
       set_hom_from_matrix(phi, target_basic_variables, mat);
       if( phi.is_well_defined() ){
//       cout << "wooooooot ! killing a matrix !" << endl;
//       cout << "automorphism was: " << phi << endl;
         //cout << "it was a clone of matrix # " << m << endl;
         N= GL.erase(N);
         count--;
       }
       else
         N++;
     } 
   if(verbose)
     cout << "*** Step 1:" << endl
          << count << " possible choices, after reduction"
          << endl;
  return count;
}



void HomFinder::init(){
  long n=0;
  long d, i;  
  Tuple< Polynomial<F_2> > pows;

  
  // populate higher_variables, basic_variables,
  // monomials, polvars, and compute n
  for(i=1; i <= source->variables_in_use(); i++){
    d= source->hom_degrees.find(i)->second;
    if( d > 1){
      higher_variables.push_back(i);
      if( !source->relations.variable_shows_up(i) )
        polvars.push_back(i);
    }
    else{
      basic_variables.push_back(i);
      n++;
    }
    
    if( monomials.count( d ) == 0){
      target->list_monomials(d, pows);      
      monomials[d]= pows;
    }
    
  }

  //as the name says...
  populate_GL(n);
  //now the same with target_basic_variables
  for(i=1; i <= target->variables_in_use(); i++){
    d= target->hom_degrees.find(i)->second;
    if( d == 1)
      target_basic_variables.push_back(i);
  }


  list< long >::iterator hvar;

  //populates basic_source
  basic_source= *source;
  basic_source.relations.generators.clear(); 
  list< Polynomial<F_2> >::const_iterator it;
  bool add;  
  for(it= source->relations.generators.begin();
      it!= source->relations.generators.end();
      it++){
    add= true;
    for(hvar= higher_variables.begin();
        hvar!= higher_variables.end();
        hvar++)
      if( it->involves_variable(*hvar) ){
        add= false;
        break;
      }
    if( add )
      basic_source.add_relation( *it );
    else
      higher_relations.push_back( *it );
  }

  // keep only the compatible matrices
  filter_GL();
  //sort the highest_relation according to their weight
  sort();
}



void HomFinder::set_hom_from_matrix(AffineHomomorphism<F_2>& phi,
                                    const list< long >& basics,
                                    const Matrix<F_2>& M)
{
  long n= M.nlines();
  long i,j;
  list< long >::const_iterator bvars;
  Polynomial<F_2> P= target->one();//alphabet...

  for(bvars= basics.begin(), j=0;
      bvars!= basics.end();
      bvars++, j++){
    P.sets_to_zero();
    for(i=0; i<n; i++)
      if( M(i,j)==1 )
        P+= monomials[1][i];
    phi.set_image(*bvars, P);
  }
}


void HomFinder::set_f_from_matrix(const Matrix<F_2>& M){
  long n= M.nlines();
  long i,j;
  list< long >::iterator bvars;
  Polynomial<F_2> P= target->one();//alphabet...

  for(bvars= basic_variables.begin(), j=0;
      bvars!= basic_variables.end();
      bvars++, j++){
    P.sets_to_zero();
    for(i=0; i<n; i++)
      if( M(i,j)==1 )
        P+= monomials[1][i];
    f.set_image(*bvars, P);
  }
}


void HomFinder::start(){
  list< Matrix<F_2> >::iterator M;
  ZZ code;
  long m=0;

  f.target= target;
  f.source= source;
  
  solution_count= 0;
  

  for(M= GL.begin(); M!= GL.end(); M++){
    m++;
  
    set_f_from_matrix(*M);
    
    //    cout << f << endl;
    
    branch(f, higher_relations.begin(), 0);
    



  }// matrix M

}//done!


void HomFinder::set_f_from_integer(ZZ code){
  ZZ bit=1;
  list< long >::iterator hvars;
  long deg, i;
  Polynomial<F_2> P=target->one(); //for the alphabet !


  for(hvars= higher_variables.begin();
      hvars!= higher_variables.end();
      hvars++){
    P.sets_to_zero();
    deg= source->hom_degrees.find( *hvars )->second;
    for(i=0; i< monomials[deg].size(); i++){
      if( (code & bit) != 0 )
        P+= monomials[deg][i];
      bit*=2;
    }
    f.set_image(*hvars, P);
  }// for hvars
}

// receives:  phi: *source -> *target
// a morphism, well-defined so far.

void HomFinder::branch(const AffineHomomorphism<F_2>& phi,
                       list< Polynomial<F_2> >::iterator next,
                       long depth){

  list< long > newvars;
  long i;
  AffineHomomorphism<F_2> local_phi= phi;
  string tab;
  // a tabulation, for display purposes
  for(i=0; i<depth; i++)
    tab+="    "; 



  if( next == higher_relations.end() ){
    //cout << "homomorphism found ! " <<endl;
    
    solutions.push_back(phi);
    solution_count++;
    cout << "homomorphism found, total "
         << solution_count << endl;
           
    return;
  }
  // else, 

  // DEBUGGING PURPOSES
  /* char stop;
     cout << tab << "depth " << depth << endl
     << tab << "next equation: " << *next << endl
     << phi;
     cin >> stop;
  */





  list< Polynomial<F_2> >::iterator next_next= next;
  next_next++;

  //we look at *next and see what variables it involves
  for(i=1; i<= source->variables_in_use(); i++)
    if( next->involves_variable(i) && !phi.has_image_set(i) )
      newvars.push_back(i);
 
  // compute "possibilities"
  list< long >::iterator nv;
  ZZ poss;
  ZZ dim=0;
  long d;
  for(nv= newvars.begin();
      nv!= newvars.end();
      nv++){
    d= source->hom_degrees.find(*nv)->second;
    dim+= monomials[d].size();
  }
  // now set possibilities = 2^dim 
  poss=1;
  
  for(i=0; i< dim; i++){
    poss*=2;
  }

  //now try them all !
  ZZ code;
  ZZ bit;
  long deg;
  Polynomial<F_2> P=target->one(); //for the alphabet !

  for(code=0; code < poss; code++){

    bit= 1;

    for(nv= newvars.begin();
        nv!= newvars.end();
        nv++){
      P.sets_to_zero();
      deg= source->hom_degrees.find( *nv )->second;
      for(i=0; i< monomials[deg].size(); i++){
        if( (code & bit) != 0 )
          P+= monomials[deg][i];
        bit*=2;
      }
      local_phi.set_image(*nv, P);
    }// for nv
    // now local_phi is set.
    if( target->reduced_form( local_phi.of(*next) ).is_zero() ){
      cout << tab << "(depth " << depth
           << " code " << code << " out of "
           << poss <<") match found" << endl;
      
      branch(local_phi, next_next, depth + 1);
    }
  }//for code

  return;
}





// void HomFinder::old_define_polynomial_variables()
// {
//   list< AffineHomomorphism<F_2> >::iterator it;
//   list< AffineHomomorphism<F_2> > old_solutions;
//   list< long > polvars;
//   list< long >::iterator itpol;
//   long i, d, kept=0;
//   ZZ dim, poss, code, bit, precise_poss, precise_code, precise_bit;
//   GradedAlgebra<F_2> alg;
//   map< long, Tuple< Polynomial<F_2> > > choices;
//   map< long, Tuple< Polynomial<F_2> > > precise_choices;
//   AffineHomomorphism<F_2> hom;
//   Tuple< Polynomial<F_2> > pows;
//   Polynomial<F_2> P= target->one(), Q;
  
//   //first we copy the solutions
//   old_solutions.clear();
//   old_solutions.splice( old_solutions.begin(), solutions );
//   //we look for the polynomial variables
//   polvars.clear();
//   it= old_solutions.begin();

//   for(i=1; i<= source->variables_in_use(); i++)
//     if( !it->has_image_set(i) )
//       polvars.push_back(i);
//   //ok now we compute the possible values
//   precise_choices.clear();
//   dim= 0;
  
//   for(itpol= polvars.begin(); itpol!= polvars.end(); itpol++){
      
//       d= source->hom_degrees.find( *itpol )->second;//d= degree of variable *itpol
//       if( precise_choices.count( d ) == 0 ){
//      target->list_monomials(d, pows);      
//      precise_choices[d]= pows;
//       }
//       dim+= precise_choices[d].size();
//   }
//   //set poss= 2**dim
//   precise_poss= 1;
//   while( dim > 0){
//     precise_poss*= 2;
//     dim--;
//   }
   
//   cout << "total possibilities: " << precise_poss << endl;

//   //ok! now we go through the solutions one by one and
//   //see each time what possibilities we have for the values
//   //of the polynomial variables *up to the ideal image* of
//   //the other variables. Much fewer possibilities to check;
//   //and when we find one that works, we check *all* possibilities.
 
//   for(it= old_solutions.begin();
//       it!= old_solutions.end();
//       it++){

//     hom= *it;

//     //let's create an algebra that's a quotient of target
//     //by the image generated by the image of *it
//     alg= *target;


//     for(i=1; i<= source->variables_in_use(); i++)
//       if( it->has_image_set(i) )
//      alg.add_relation( it->of_generator(i) );
    
//     alg.fix_alphabets();
//     alg.find_grobner_basis();
    
//     // now let's see what the possible choices are
//     choices.clear();
//     dim= 0;
    
//     for(itpol= polvars.begin(); itpol!= polvars.end(); itpol++){
      
//       d= source->hom_degrees.find( *itpol )->second;//d= degree of variable *itpol
//       if( choices.count( d ) == 0 ){
//      alg.list_monomials(d, pows);      
//      choices[d]= pows;
//       }
//       dim+= choices[d].size();
//     }
//     //set poss= 2**dim
//     poss= 1;
//     while( dim > 0){
//       poss*= 2;
//       dim--;
//     }
    
//     //now try them!
    
//     for(code= 0; code < poss; code++){
      
//       bit= 1;    
//       for(itpol= polvars.begin(); itpol!= polvars.end(); itpol++){
//      d= source->hom_degrees.find( *itpol )->second;//d= degree of variable *itpol
//      P.sets_to_zero();
        
//      for(i=0; i < choices[d].size(); i++){
//        if( (code & bit) > 0 )
//          P+= choices[d][i];
//        bit*= 2;
//      }
        
//      //now P is set
//      it->set_image(*itpol, P);
                
//       }//for itpol
//       // now *it is fully defined !
//       if( it->is_finite() ){
//      //ok we've got one
//      cout << "found one!" << endl;
//      //we start again, but with the 'precise' version
//      for(precise_code= 0; precise_code < precise_poss; precise_code++){
          
//        cout << "precise code: " << precise_code << endl;

//        precise_bit= 1;    
//        for(itpol= polvars.begin(); itpol!= polvars.end(); itpol++){
//          d= source->hom_degrees.find( *itpol )->second;//d= degree of variable *itpol
//          P.sets_to_zero();
            
//          for(i=0; i < precise_choices[d].size(); i++){
//            if( (precise_code & precise_bit) > 0 )
//              P+= precise_choices[d][i];
//            precise_bit*= 2;
//          }
        
//          //now P is set
//          //is P equal to our previous choice for
//          //it->image_of(*itpol) modulo ... ?
//          Q= P + it->of_generator(*itpol);
//          if( !alg.reduced_form(Q).is_zero() ) //if not
//            break;                             //we break
//          else
//            hom.set_image(*itpol, P);
//        }//for itpol
//        if( itpol == polvars.end() ){ //if we did not break
//          kept++;
//          // it->compute_kernel();
//          solutions.push_back(hom);
//          cout << "found a precise solution" << endl
//               << "precise code " << precise_code
//               << " out of " << precise_poss
//               << " kept " << kept << endl;
//        }//if we did not break
//      }//for precise code
//       }//if we had something finite
      
//     }//for code
//   }//for it
  
//   cout << "kept: " << kept << endl;
//   solution_count= kept;
  
// }








void HomFinder::define_polynomial_variables()
{
  list< AffineHomomorphism<F_2> >::iterator it;
  list< AffineHomomorphism<F_2> > old_solutions;
  //list< long > polvars; //now made global
  list< long >::iterator itpol;
  long i, d, kept=0;
  ZZ dim, poss, code, bit;
  GradedAlgebra<F_2> alg;
  map< long, Tuple< Polynomial<F_2> > > choices;
  Tuple< Polynomial<F_2> > pows;
  Polynomial<F_2> P= target->one();

  old_solutions.clear();
  old_solutions.splice( old_solutions.begin(), solutions );
  
  for(it= old_solutions.begin();
      it!= old_solutions.end();
      it++){
    

    //let's find out which variables are "polynomial"
    //also, create an algebra that's a quotient of target
    //by the image generated by the image of *it
    
    alg= *target;

                        
    for(i=1; i<= source->variables_in_use(); i++)
      if( it->has_image_set(i) )
        alg.add_relation( it->of_generator(i) );
    
    alg.fix_alphabets();
    alg.find_grobner_basis();
    
    // now let's see what the possible choices are
    choices.clear();
    dim= 0;
    
    for(itpol= polvars.begin(); itpol!= polvars.end(); itpol++){
      
      d= source->hom_degrees.find( *itpol )->second;//d= degree of variable *itpol
      if( choices.count( d ) == 0){
        alg.list_monomials(d, pows);      
        choices[d]= pows;
      }
      dim+= choices[d].size();
      
            
    }
    //set poss= 2**dim
    poss= 1;
    while( dim > 0){
      poss*= 2;
      dim--;
    }
    
    //now try them!
    
    for(code= 0; code < poss; code++){
      bit= 1;    
      for(itpol= polvars.begin(); itpol!= polvars.end(); itpol++){
        d= source->hom_degrees.find( *itpol )->second;//d= degree of variable *itpol
        P.sets_to_zero();
        
        for(i=0; i < choices[d].size(); i++){
          if( (code & bit) > 0 )
            P+= choices[d][i];
          bit*= 2;
        }
        
        //now P is set
        it->set_image(*itpol, P);
        
        
      }//for itpol
      // now *it is fully defined !
      if( it->is_finite() ){
        kept++;
        //it->compute_kernel();
        solutions.push_back(*it);
      }
      
    }//for code
  }//for it

  if(verbose)
    cout << kept << " possibilities kept after Steps 1,2,3 and Test 1"
         << endl;
  solution_count= kept;
  
}

void HomFinder::define_polynomial_variables_brutally()
{
  list< AffineHomomorphism<F_2> >::iterator it;
  list< AffineHomomorphism<F_2> > old_solutions;
  //list< long > polvars; //now made global
  list< long >::iterator itpol;
  long i, d, kept=0;
  ZZ dim, poss, code, bit;
  GradedAlgebra<F_2> alg;
  map< long, Tuple< Polynomial<F_2> > > choices;
  Tuple< Polynomial<F_2> > pows;
  Polynomial<F_2> P= target->one();

  old_solutions.clear();
  old_solutions.splice( old_solutions.begin(), solutions );

  alg= *target;
  alg.fix_alphabets();
  
  for(it= old_solutions.begin();
      it!= old_solutions.end();
      it++){
    

    // now let's see what the possible choices are
    choices.clear();
    dim= 0;
    
    for(itpol= polvars.begin(); itpol!= polvars.end(); itpol++){
      
      d= source->hom_degrees.find( *itpol )->second;//d= degree of variable *itpol
      if( choices.count( d ) == 0){
        alg.list_monomials(d, pows);      
        choices[d]= pows;
      }
      dim+= choices[d].size();
      
            
    }
    //set poss= 2**dim
    poss= 1;
    while( dim > 0){
      poss*= 2;
      dim--;
    }
    
    //now try them!
    
    for(code= 0; code < poss; code++){
      bit= 1;    
      for(itpol= polvars.begin(); itpol!= polvars.end(); itpol++){
        d= source->hom_degrees.find( *itpol )->second;//d= degree of variable *itpol
        P.sets_to_zero();
        
        for(i=0; i < choices[d].size(); i++){
          if( (code & bit) > 0 )
            P+= choices[d][i];
          bit*= 2;
        }
        
        //now P is set
        it->set_image(*itpol, P);
        
        
      }//for itpol
      // now *it is fully defined !
      if( it->is_finite() ){
        kept++;
        //it->compute_kernel();
        solutions.push_back(*it);
      }
      
    }//for code
  }//for it

  if(verbose)
    cout << kept << " possibilities kept after Steps 1 & 2 "
         << " (all choices for polynomial variables included) and Test 1"
         << endl;
  solution_count= kept;
  
}



void HomFinder::identify()
{
  list< AffineHomomorphism<F_2> >::iterator it, itp;
  long cit=0, citp=0;
      
  for(it= solutions.begin();
      it!= solutions.end();
      it++){
    
    cit++;
    citp= 0;
    
    for(itp= it, itp++;
        itp!= solutions.end();){
      
      citp++;
      if(very_verbose)
        cout << "comparing " << cit << " and " << citp << endl;
      
        
      // are f and g essentially the same ?
      
      
      if( it->is_essentially_contained_in(*itp) &&
          itp->is_essentially_contained_in(*it) &&
          it->factors_through(itp->kernel) &&
          itp->factors_through(it->kernel) )
        {
          itp= solutions.erase( itp );
          solution_count--;
        }
      else
        itp++;
    }//for itp 
  }//for it
  
}


long HomFinder::weight(const Polynomial<F_2>& relation, const AffineHomomorphism<F_2>& phi)
{
  list< long >::iterator it;
  list< long > newvars;
  
  newvars.clear();
    
  for(it= higher_variables.begin(); it!= higher_variables.end(); it++)
    if( relation.involves_variable(*it) && !phi.has_image_set(*it) )
      newvars.push_back(*it);

  // compute "possibilities"
  list< long >::iterator nv;
  long dim=0;
  long d;
  for(nv= newvars.begin();
      nv!= newvars.end();
      nv++){
    d= source->hom_degrees.find(*nv)->second;
    dim+= monomials[d].size();
  }

  return dim;
}

void HomFinder::sort()
{
  long minimal_weight, w;
  list< Polynomial<F_2> >::iterator it, itmini;
  list< Polynomial<F_2> > old_gens;
  AffineHomomorphism<F_2> phi(source, target);
  Polynomial<F_2> zeropol= target->one();
  list< long >::iterator hvars;
  
  old_gens.splice(old_gens.begin(), higher_relations );
  zeropol.sets_to_zero();
    
  while( !old_gens.empty() ){
    minimal_weight= -1;
    // find the relation of minimal weight
    for(it= old_gens.begin();
        it!= old_gens.end();
        it++){
      w= weight(*it, phi);
      if( minimal_weight < 0 || minimal_weight > w){
        minimal_weight= w;
        itmini= it;
      }
    }

    //    cout << "minimum: " << minimal_weight << " for relation "
    // << *itmini << endl;
        
    higher_relations.push_back( *itmini );
    //update phi
    for(hvars= higher_variables.begin();
        hvars!= higher_variables.end();
        hvars++)
      if( itmini->involves_variable(*hvars) )
        phi.set_image(*hvars, zeropol);
    //done
    old_gens.erase(itmini);
        
  }//while nonempty
  
   
}

// TO DO: REMOVE SOME GROBNER BASIS COMPUTATIONS !!!!!
void HomFinder::test_steenrod(){
  list< AffineHomomorphism<F_2> >::iterator it;
  list< Polynomial<F_2> >::iterator itker;
  long n=0;
  UnstableAlgebra test;
  UnstableAlgebra *ptr;

  ptr= (UnstableAlgebra *) source;

  for(it= solutions.begin(); it!=solutions.end();){

    test= *ptr;
    test.fix_alphabets();

    n++;
    if(verbose)
      cout << "testing homomorphism number " << n << endl;

    for(itker= it->kernel.generators.begin(); 
        itker!= it->kernel.generators.end();
        itker++){
      //      cout << *itker << endl;
      test.find_grobner_basis();
      if( !test.reduced_form( *itker ).is_zero() )
        test.add_relation( *itker );
      }
    

    test.find_grobner_basis();
    if( test.is_saturated() ){
      if(verbose)
        cout << "steenrod test passed" << endl;
      it++;
    }
    else
      {      
        if(verbose)
          cout << "rejected" << endl;
        it= solutions.erase( it );
        solution_count--;
      }

  }//for it

}


void HomFinder::restrict_to_elemab(ifstream& res_file) {
  //check which kernels are compatible with restrictions
  list< AffineHomomorphism<F_2> >::iterator it;
  list< Polynomial<F_2> >::iterator ker;
  long nb_elemab, rank, i;
  Tuple< AbstractHomomorphism<F_2> > res;
  Polynomial<F_2> P;
  
  res_file >> nb_elemab;
  res.resize(nb_elemab);
  for(i=0; i< nb_elemab; i++){
    res_file >> rank;//unused so far !!
    res_file >> res[i];
  }
  
  for(it= solutions.begin();
      it!= solutions.end();){

    if(verbose) 
      cout << "trying one solution..." << endl;

        //we try each elemab
    for(i=0; i< nb_elemab; i++){
      if(very_verbose)
        cout << "trying subgroup number " << i << endl;
          
      for(ker= it->kernel.generators.begin();
          ker!= it->kernel.generators.end();
          ker++){
            
        P= res[i].of( *ker );
            
        if(!P.is_zero()){
          if(verbose)
            cout << "impossible !" << endl
                 << *ker << " restricts to " << P << endl;
          break; // out of for ker
        }
            
      }//for ker
      if( ker != it->kernel.generators.end() )//did we break ?
        break; // out of for i
    }//for i
    if(i<nb_elemab){ // contradiction found ?
      it= solutions.erase(it);
      solution_count--;
    }
    else
      it++;
  }//for it
      
 
}


void HomFinder::get_equivalence_class(const list< AffineHomomorphism<F_2> >& homs,
                                      const AffineHomomorphism<F_2>& f){
  list< AffineHomomorphism<F_2> >::const_iterator it;

  solution_count= 0;
  solutions.clear();
  
  for(it= homs.begin();
      it!= homs.end();
      it++){
        
    if( it->is_essentially_contained_in(f) &&
        f.is_essentially_contained_in(*it) &&
        it->factors_through(f.kernel) &&
        f.factors_through(it->kernel) ) {

      solutions.push_back(*it);
      solution_count++;
    }

  }//for it

//   if(verbose)
//     cout << "there were " << solution_count
//       << " homomorphisms in the equivalence class"
//       << endl;
}



bool HomFinder::polvars_are_involved_in(const AffineHomomorphism<F_2>& f){
  list< long >::iterator itpol;

  for(itpol= polvars.begin();
      itpol!= polvars.end();
      itpol++)
    if(f.kernel.variable_shows_up( *itpol ))
      return true;
  return false;
}

long HomFinder::extend(AffineHomomorphism<F_2>& f){
  long ans=0;  

  // construct B/<f(A)>
  GradedAlgebra<F_2> quotient= *target;
  long i;
  for(i=1; i<= source->variables_in_use(); i++)
    quotient.add_relation(f.of_generator(i));

  quotient.fix_alphabets();
  quotient.find_grobner_basis();
  //done
  
  //now add the missing variables to A
  Tuple< Polynomial<F_2> > Bvars= target->get_variables();
  map< long, Polynomial<F_2> > newvars;

  for(i=1; i<= quotient.variables_in_use(); i++)
    if( !quotient.reduced_form(Bvars[i]).is_zero() ){
      ans++;
      source->new_variable(quotient.nameof(i), quotient.hom_degrees[i]);
      newvars[ source->variables_in_use() ]= Bvars[i];
    }
  //CAREFUL ! at this point the steenrod operations are not defined
  //for the new elements just added in A !!

  // now extend f
  map< long, Polynomial<F_2> >::iterator it_new;

  for(it_new= newvars.begin();
      it_new != newvars.end();
      it_new++)
    f.set_image( it_new->first, it_new->second );
  
  return ans;
}

long HomFinder::extend_all(){
  long ans=0;  
  AffineHomomorphism<F_2> f= *solutions.begin();

  // construct B/<f(A)>
  GradedAlgebra<F_2> quotient= *target;
  long i;
  for(i=1; i<= source->variables_in_use(); i++)
    quotient.add_relation(f.of_generator(i));

  quotient.fix_alphabets();
  quotient.find_grobner_basis();
  //done
  
  //now add the missing variables to A
  Tuple< Polynomial<F_2> > Bvars= target->get_variables();
  map< long, Polynomial<F_2> > newvars;

  for(i=1; i<= quotient.variables_in_use(); i++)
    if( !quotient.reduced_form(Bvars[i]).is_zero() ){
      ans++;
      source->new_variable(quotient.nameof(i), quotient.hom_degrees[i]);
      newvars[ source->variables_in_use() ]= Bvars[i];
    }
  //CAREFUL ! at this point the steenrod operations are not defined
  //for the new elements just added in A !!

  // now extend all the solutions
  list< AffineHomomorphism<F_2> >::iterator it;
  map< long, Polynomial<F_2> >::iterator it_new;

  for(it= solutions.begin();
      it != solutions.end();
      it++){

    for(it_new= newvars.begin();
        it_new != newvars.end();
        it_new++)
      it->set_image( it_new->first, it_new->second );
    if( ans > 0 )
      it->compute_kernel(); //this is new in extend_all(), compared to
                            //extend()
  }
  
  return ans;
}



bool HomFinder::polynomial_test(){
  bool ans;
  long n;
  list< AffineHomomorphism<F_2> >::iterator it;
  Polynomial<F_2> zero= source->one();
  zero.sets_to_zero();
  list< Polynomial<F_2> >::iterator debug;

  for(it= solutions.begin();
      it != solutions.end();
      it++){

    n= extend( *it );
    it->compute_kernel();
    ans= polvars_are_involved_in( *it );
    // we put A back into its original state, whether the test has
    // failed or not, just to be clean
    while(n > 0){
      source->kill_top_variable(zero);
      n--;
    }
    //ok done
    if( ans )
      return false;
    // else, we restore *it also, so it's all ready to proceed
    it->remove_extra_data();
  }//for it

  return true;
}


void HomFinder::compute_all_kernels(){
  list< AffineHomomorphism<F_2> >::iterator it;

  for(it= solutions.begin();
      it != solutions.end();
      it++){
    it->compute_kernel();
  }
}

---