derive.cpp

#include "derive.H"
using namespace std;


void DerivationAnalyzer::setup(GradedAlgebra<F_2> *thesource){
  long i, n;
  list< Polynomial<F_2> > pols, total;
  list< Polynomial<F_2> >::iterator it;

  source= thesource;

  n= 1;
  //write the generators for *source as a module
  //over the subalgebra of squares
  for(i=1; i<= source->variables_in_use(); i++){
    write_gens( pols, i, 1, source->variables_in_use() );
    for(it= pols.begin(); it!= pols.end(); it++, n++)
      index_of[ it->lp() ]= n; 
    total.splice( total.end(), pols);
  }
  //turns the list into a tuple
  mod_gens.resize(n);
  mod_gens[0]= source->one();
  for(i= 1, it= total.begin(); i< n; i++, it++)
    mod_gens[i]= *it;

  d.resize(n);
  d[0].resize(n);
  d[0].set_alphabet( source );
  d[0].sets_to_zero();

  cout << "generators: " << mod_gens << endl;
  //now construct the subalgebra of squares
  //together with the inclusion homomorphism
  Tuple< Polynomial<F_2> > x= source->get_variables();
  sq_inclusion.source= &squares;
  sq_inclusion.target= source;
  cs_inclusion.source= &constants;
  cs_inclusion.target= source;
  for(i=1; i<= source->variables_in_use(); i++){
    squares.new_variable("s" + source->nameof(i), 
                         2*source->hom_degrees[i] );
    constants.new_variable("s" + source->nameof(i),
                           2*source->hom_degrees[i] );
    sq_inclusion.set_image(i, x[i]*x[i] );
    cs_inclusion.set_image(i, x[i]*x[i] );
  }

  cs_inclusion.populate_associated_algebra();

  //  cout << squares << endl;
  //cout << sq_inclusion << endl;

}

void DerivationAnalyzer::write_gens(list< Polynomial<F_2> >& result, 
                long ext, long first, long last) {
  list< Polynomial<F_2> > temp1, temp2;
  list< Polynomial<F_2> >::iterator it;


  // trivial cases: ext too small or too large
  if( (ext < 1) ){
    result.clear();
    result.push_front( source->one() );
    return;
  }
  
  if( (ext > last - first + 1) ){
    result.clear();
    return;
  }
    
  // else...
  // use recursion:
  write_gens(result, ext, first + 1, last);
  write_gens(temp2, ext - 1, first + 1, last);
  // now add the 'first' variable to the terms
  // in temp2
  for(it= temp2.begin(); it != temp2.end(); it++)
    result.push_back( *it * Polynomial<F_2>(PowProd( first )) );
  return;

}

multiPolynomial<F_2> DerivationAnalyzer::coords_of(PowProd P)

{
  PowProd tail(0);
  PowProd half;
  long i;
  multiPolynomial<F_2> ans;
  Polynomial<F_2> poly;

  for(i= 1; i<= P.nvars(); i++)
    if( P.power_of(i) % 2 != 0 )        // odd ?
      tail*= PowProd(i);

  P /= tail;
  half.powers.resize( P.powers.size() );
  for(i=0; i< P.powers.size(); i++)
    half.powers[i]= P.powers[i] / 2;



  //now simply format the answer
  ans.resize( mod_gens.size() );
  ans.set_alphabet( &squares );
  poly= Polynomial<F_2>( half );
  poly.alphabet= &squares;
  ans[ index_of[tail] ]= poly;
  return ans;
}

multiPolynomial<F_2> DerivationAnalyzer::coords_of(const Polynomial<F_2>& P) {
  map< PowProd, F_2 >::const_iterator it;
  multiPolynomial<F_2> ans;

  ans.resize( mod_gens.size() );
  ans.set_alphabet( &squares );

  for(it= P.coeffs.begin();
      it != P.coeffs.end();
      it++)
    ans+= this->coords_of( it->first );
  
  return ans;
}

void DerivationAnalyzer::constant_subalgebra(string prefix, 
                                             long offset){
  multiIdeal<F_2> module;
  long i, j;
  Matrix< Polynomial<F_2> > M, T;
  Polynomial<F_2> P, Q, replacement;
  list< Polynomial<F_2> >::iterator it;
  list< Polynomial<F_2> > solutions;
  list< Polynomial<F_2> > homo_solutions;
  ostringstream autoname;

  //add the d(generators) to the multiIdeal
  for(i=1; i< d.size(); i++) {
    module.generators.push_back( d[i] );
    //the syzygy methods below have an oddity:
    //they do not treat 0 elements well, in fact
    //they are ignored (or things amount to this).
    //we need to add manually the corresponding
    //solutions:
    if( d[i].is_zero() ) {
      solutions.push_back( mod_gens[i] );
    }
    //also note that d[0]= 0, but one() is
    //automatically part of the algebra of
    //solutions
  }
  //done, now add the relations of *source
  
  for(it= source->relations.generators.begin();
      it != source->relations.generators.end();
      it++)
    for(i= 0; i< mod_gens.size(); i++){
      module.generators.push_back( this->coords_of(mod_gens[i]*(*it) ) );
    }

  cout << "computing the kernel of the derivation..." << endl;

  //  module.grobnerize();
  //cout << "done !" << endl;
  
  module.grobnerize_with_coefficients(T);

  cout << "for information: T has "
       << T.nlines() << " lines " << endl;

  module.syzygy(M);
  cout << "performing large matrix multiplication... (" 
       << M.nlines() << " lines)" << endl;
  M*= T;
  //now resize M to keep the first columns
  M.resize( M.nlines(), mod_gens.size()-1 );


  //now forming the algebra of constants
  cout << "finding a minimal presentation" << endl;

  Q.alphabet= source;
  

  cout << "nb of generators: " << M.nlines() << endl;

  for(i=0; i< M.nlines(); i++){
    
    Q.sets_to_zero();
    for(j=1; j < mod_gens.size(); j++)
      Q+= sq_inclusion.of(M(i, j-1)) * mod_gens[j];
    //Q is computed, now reduce it
    Q= source->reduced_form(Q);
    solutions.push_back(Q);
  }

  //split_and_sort_relations...
  i= 0;

  while( !solutions.empty() ) {
    i++;

    for(it= solutions.begin();
        it!= solutions.end();) {
      Q= source->hom_part_of( *it, i);
      if( !Q.is_zero() ) {
        homo_solutions.push_back( Q );
        *it += Q;
        //cout << "  added " << Q << endl;
      }      
      if( it->is_zero() )
        it= solutions.erase( it );
      else
        it++;
    }

  }

  i= 0;

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

    if( !cs_inclusion.maps_onto(*it) ){
        
      cout << "adding " << *it << endl;
    
      //we need to find a name for the variable
      autoname.str("");
      i++;
      autoname << prefix << "_" << offset + i;
      constants.new_variable( autoname.str(),
                              source->hom_degree_of(*it) );
      cs_inclusion.set_image( constants.variables_in_use(), *it);
      cs_inclusion.populate_associated_algebra();
    }
      
  }
 
  constants.relations= cs_inclusion.kernel;
  //cout << constants << endl << cs_inclusion << endl;

  

}























// MilnorAnalyzer code /////////////////////////

void MilnorAnalyzer::start(UnstableAlgebra *theA) {
  long i, k;
  Polynomial<F_2> P, R;
  char letter= 'a';
  ostringstream autoname;
  Tuple< Polynomial<F_2> > x;

  A= theA;

  //the first thing to do is compute the kernel
  //of Sq^1, which is different from the other
  //derivations for technical reasons.
  cout << "dealing with Sq^1" << endl;
  k= 0;
  //get things ready; in particular this
  //makes a copy of A into *source, and
  //'inclusion' is from *source to A
  setup_for_Sq1();
  //now compute the kernel
  autoname << letter;
  constant_subalgebra(autoname.str(), 0);
  //'inclusion' must be from 'constants' to A,
  // both as the final answer and within done():
  inclusion *= cs_inclusion;
  inclusion.populate_associated_algebra();

  //now the higher milnor derivations

  while( !done() ) {
    k++;
    letter++;
    cout << endl 
         << "dealing with Q_" << k << "..." << endl;
    //the following swaps things around:
    //constants is copied to *source(== alg),
    //squares is updated... also, 'inclusion'
    //is again from *source to A :
    re_setup(k);
    //now compute the kernel    
    autoname.str("");
    autoname << letter;
    constant_subalgebra(autoname.str(), 0);
    inclusion *= cs_inclusion;
    inclusion.populate_associated_algebra();
  }

  cout << "wooooot !! done !!" << endl;

}


void MilnorAnalyzer::setup_for_Sq1(){
  long i, n;
  list< Polynomial<F_2> > pols, total;
  list< Polynomial<F_2> >::iterator it;
  Tuple< Polynomial<F_2> > x;
  Polynomial<F_2> P, R;
  
  //first we make a copy of *source
  //and we setup 'inclusion' to be the inclusion
  //of this copy into *A
  x= A->get_variables();
  inclusion.source= &alg;
  inclusion.target= A;
  for(i=1; i<= A->variables_in_use(); i++) {
    alg.new_variable(A->nameof(i),
                     A->hom_degrees[i]);
    inclusion.set_image(i, x[i]); 
  }
  alg.relations= A->relations;
  alg.fix_alphabets();
  //or, we could have done something like
  //alg= (cast) A;
  //and inclusion.make_identity() or
  //whatever (there would be details anyway)
  source= &alg;//done
  //note that A still points at the original
  //unstable algebra
  //now see what variables in *source satisfy
  //sq^1 = 0, and put them at the beginning
  n_cs_vars= 0;
  for(i=1; i<= source->variables_in_use(); i++) {
    P= A->Sq(1, x[i]);
    if( A->reduced_form(P).is_zero() ) {
      n_cs_vars++;
      alg.swap_variables_order(i, n_cs_vars);
      inclusion.swap(i, n_cs_vars);
    }
  }//the following will be useful:
  inclusion.populate_associated_algebra();
  //ok, done
  //now add to squares and constants 
  //the variables for which sq^1=0 and
  //the squares of the others
  x= source->get_variables();
  sq_inclusion.source= &squares;
  sq_inclusion.target= source;
  cs_inclusion.source= &constants;
  cs_inclusion.target= source;
  for(i=1; i<= n_cs_vars; i++) { //those with sq1=0
    squares.new_variable(source->nameof(i), 
                         source->hom_degrees[i] );
    constants.new_variable(source->nameof(i),
                           source->hom_degrees[i] );
    sq_inclusion.set_image(i, x[i] );
    cs_inclusion.set_image(i, x[i] );

  }//now the others:
  for(i= n_cs_vars+1; i<=source->variables_in_use(); i++){
      squares.new_variable("s" + source->nameof(i), 
                           2*source->hom_degrees[i] );
      constants.new_variable("s" + source->nameof(i),
                             2*source->hom_degrees[i] );
      sq_inclusion.set_image(i, x[i]*x[i] );
      cs_inclusion.set_image(i, x[i]*x[i] );
  }
  //cs_inclusion is supposed to have its associated
  //algebra already computed:
  cs_inclusion.populate_associated_algebra();
  //now we setup newvar_squared: for sq1 the name
  //'newvar' is not appropriate, btw! in this context
  //the "new" variables are those with sq1 != 0.
  //the purpose of newvar_squared, anyway, remains
  //to express the square of variable i in terms
  //of elements in the algebra 'squares': here for
  //sq1 it is just the identity, by definition.
  newvar_squared.clear();
  x= squares.get_variables();
  for(i= n_cs_vars+1; i<= squares.variables_in_use(); i++)
    newvar_squared.set_image(i, x[i]);
  //good, now write the generators for *source as a module
  //over the subalgebra of squares
  n= 1;
  for(i=1;
      i<= source->variables_in_use() - n_cs_vars; 
      i++){
    write_gens( pols, i, n_cs_vars+1, 
                source->variables_in_use() );
    for(it= pols.begin(); it!= pols.end(); it++, n++)
      index_of[ it->lp() ]= n; 
    total.splice( total.end(), pols);
  }
  //turns the list into a tuple
  mod_gens.resize(n);
  mod_gens[0]= source->one();
  for(i= 1, it= total.begin(); i< n; i++, it++)
    mod_gens[i]= *it;
  cout << "initial generators: " << mod_gens << endl;
  //finally set d= Sq^1
  d.resize(n);
  d[0].resize(n);
  d[0].set_alphabet( source );
  d[0].sets_to_zero();
  for(i=1; i< mod_gens.size(); i++) {
    P= A->Sq(1, inclusion.of(mod_gens[i]));
    inclusion.maps_onto(P, R);
    //this has worked since
    //inclusion.populate_associated_algebra()
    //has been called
    d[i]= this->coords_of(R);
  }
  //that's all !


}


//get ready for Q_k
void MilnorAnalyzer::re_setup(long k) {
  long i;
  Tuple< Polynomial<F_2> > x;
  Polynomial<F_2> P, R;
  GradedAlgebra<F_2> empty;


  //first thing first: see what variables in
  // 'constants' satisfy Q_k(var)==0 and move
  //them to the beginning
  n_cs_vars= 0;

  for(i= n_cs_vars + 1; 
      i <= constants.variables_in_use(); i++){
    P= A->Q(k, inclusion.of_generator(i) );
    if( A->reduced_form(P).is_zero() ) {
      // cout << "simplifying!! works ?????" 
      //           << " P= " << P << " k= " << k
      //           << endl;
      n_cs_vars++;
      constants.swap_variables_order(i, n_cs_vars); 
      cs_inclusion.swap(i, n_cs_vars);
      inclusion.swap(i, n_cs_vars);
    }
  }
  //swapping the variables unfortunately
  //destroys the grobner basis in the algebra
  //associated to cs_inclusion. It is unavoidable...
  cs_inclusion.populate_associated_algebra(); //see (*)
  inclusion.populate_associated_algebra();    //see (**)

  //good now add them to 'squares'
  squares= empty;
  for(i= 1; i <= n_cs_vars; i++) 
    squares.new_variable( constants.nameof(i),
                          constants.hom_degrees[i] );
  cout << "n_cs_vars= " << n_cs_vars << endl;
  cout << "new algebra of squares: "
       << squares << endl;

  //now we populate newvar_squared.
  newvar_squared.clear();
  x= constants.get_variables();

  for(i= squares.variables_in_use() + 1;
      i <= constants.variables_in_use();
      i++) {

    P= cs_inclusion.of( x[i]*x[i] );
    cs_inclusion.maps_onto(P, R); // (*)
    //now R maps onto P but it is reduced
    //in the associated algebra to cs_inclusion
    //and therefore it is really in the image of
    //squares inside constants
    R.alphabet= &squares;
    newvar_squared.set_image(i, R);
  }
  //constants becomes the main algebra
  alg= constants;
  alg.fix_alphabets();
  source= &alg;
  //temporarily, we view 'inclusion' as
  //defined on *source:
  inclusion.source= source;
  //the subalgebra of squares is a subalgebra
  //of the newly created alg
  sq_inclusion.target= source;
  sq_inclusion.clear();
  x= source->get_variables();

  for(i=1; i<= squares.variables_in_use(); i++)
    sq_inclusion.set_image(i, x[i]);
  
  //set the constants= the squares, for the moment
  constants= squares;
  constants.fix_alphabets();
  cs_inclusion= sq_inclusion;
  cs_inclusion.source= &constants;
  cs_inclusion.fix_alphabets();
  cs_inclusion.populate_associated_algebra(); //again...
  // all the swapping around is done !!
  //now rewrite the generators
  index_of.clear();
  mod_gens.resize(0);

  long n;
  list< Polynomial<F_2> > pols, total;
  list< Polynomial<F_2> >::iterator it;

  n= 1;
  //write the generators for *source as a module
  //over the subalgebra of squares
  for(i=1;
      i<= source->variables_in_use() - n_cs_vars; 
      i++){
    write_gens( pols, i, n_cs_vars + 1, source->variables_in_use() );
    for(it= pols.begin(); it!= pols.end(); it++, n++)
      index_of[ it->lp() ]= n; 
    total.splice( total.end(), pols);
  }
  //turns the list into a tuple
  mod_gens.resize(n);
  mod_gens[0]= source->one();
  for(i= 1, it= total.begin(); i< n; i++, it++)
    mod_gens[i]= *it;

  cout << "generators: " << mod_gens << endl;
  //finally setup d=Q_k:
  d.resize(n);
  d[0].resize(n);
  d[0].set_alphabet( source );
  d[0].sets_to_zero();
  for(i=1; i< mod_gens.size(); i++) {
    P= A->Q(k, inclusion.of( mod_gens[i] ));
    inclusion.maps_onto(P, R); // (**)
    d[i]= this->coords_of(R);
  }
  //that's all !
}



multiPolynomial<F_2> MilnorAnalyzer::coords_of(PowProd P) {
  PowProd tail(0);
  PowProd half, begin;
  long i;
  multiPolynomial<F_2> ans;
  Polynomial<F_2> poly;


  begin= P;

  if( P.nvars() > n_cs_vars ) { // P involves 'new' variables ?

    for(i= n_cs_vars - 1; i >= 0; i--)
      if( begin.powers[i] != 0 )
        break;
    begin.powers.resize(i+1);
    
  }


  P/= begin;

 
  
  for(i= 1; i<= P.nvars(); i++)
    if( P.power_of(i) % 2 != 0 )        // odd ?
      tail*= PowProd(i);
  
 
  P /= tail;

  half.powers.resize( P.powers.size() );
  for(i=0; i< P.powers.size(); i++)
    half.powers[i]= P.powers[i] / 2;

  poly= newvar_squared.of( Polynomial<F_2>(half) );
  poly*= Polynomial<F_2>( begin );
  
  //  half*= begin;

  //now simply format the answer
  ans.resize( mod_gens.size() );
  ans.set_alphabet( &squares );

  poly.alphabet= &squares;
  ans[ index_of[tail] ]= poly;
  return ans;
}

multiPolynomial<F_2> MilnorAnalyzer::coords_of(const Polynomial<F_2>& P) {

  return DerivationAnalyzer::coords_of(P);
}

bool MilnorAnalyzer::done() {
  
  constants.even_subalgebra(ans, ans_inclusion);
  ans_inclusion= inclusion * ans_inclusion;
  ans_inclusion.populate_associated_algebra();

  if( A->has_closure_of( ans_inclusion ) ) {
      ans.relations= ans_inclusion.kernel;
      return true;
    }
  else
    return false;
}

---