exp_classes.H

#ifndef _EXP_CLASSES_H_
#define _EXP_CLASSES_H_

#include "my_gmp.H"
#include <sstream>
#include <string>
#include "algebras.H"
#include "tables.H"
#include "repring.H"
using namespace std;


template <class K>
class Exponentiator {
public:
  // source and target *********************************
  // ***************************************************
  RepresentationRing *source;
  GradedAlgebra<K> *target;

  // a bunch of things to do with variables ************
  // ***************************************************

  long last_actual_variable;

  long max_dim;

  long max_dim_rel;

  Table< Polynomial<K> > var;

  Tuple< Polynomial<K> > a, b, s, sp;  
  // things to do with the relations *******************
  // ***************************************************

  list< pair< Polynomial<QQ>, Polynomial<QQ> > > equations;
  list< Polynomial<K> > total_relations;
  list< Polynomial<K> >  exp_relations;
  // misc ********************************************
  // *************************************************
  long factor;
  // some universal polynomials

  map< pair<long, long> , Polynomial<K> > univ_computed;
  map< pair<long, long> , Polynomial<K> > lambda_univ_computed;

  //**************************************************
  // methods: destructor *****************************
  // *************************************************

  virtual ~Exponentiator<K>(){}


  // variables, set-up **********************
  // *************************************************


virtual void create_regular_variables(string prefix){
    long i, j, dim;
    QQ d;
    ostringstream autoname;


    // looks for the largest dimension
    dim= 0;
    for(i=1; i <= source->variables_in_use(); i++){
      d= source->epsilon[i];
      if( d > dim)
        dim= ZZtolong(d);
    }

    max_dim= dim;
    var.resize(source->variables_in_use() +1, dim+1);

    //now create the variables and populate var
    for(i=1; i <= source->variables_in_use(); i++){

      dim= ZZtolong( source->epsilon[i] );
      for(j=1; j<= dim; j++){ 

        autoname.str("");
        autoname << prefix << "_" << j
                 << "(" << source->nameof(i) << ")";
        var(i,j)= target->new_variable(autoname.str(), factor*j);
      }
    }

  }


  void create_extra_variables(string prefix){
    long i;
    ostringstream autoname;

    last_actual_variable= target->variables_in_use();

    a.resize(max_dim+1);
    for(i=1; i<= max_dim; i++){
      autoname.str("");
      autoname << "a_" << i;      
      a[i]= target->new_variable(autoname.str(), 1);
    }

    b.resize(max_dim_rel+1);
    for(i=1; i<= max_dim_rel; i++){
      autoname.str("");
      autoname << "b_" << i;      
      b[i]= target->new_variable(autoname.str(), 1);
    }

    s.resize(max_dim+1);
    for(i=1; i<= max_dim; i++){
      autoname.str("");
      autoname << prefix << "_" << i;      
      s[i]= target->new_variable(autoname.str(), i);
    }

    sp.resize(max_dim_rel+1);
    for(i=1; i<= max_dim_rel; i++){
      autoname.str("");
      autoname << prefix << "'_" << i;      
      sp[i]= target->new_variable(autoname.str(), i);
    }
  }

  void kill_extra_variables(){
    long i,n;
    Polynomial<K> zero_poly;

    n= target->variables_in_use() - last_actual_variable;
    for(i=1; i<= n;i++)
      target->kill_top_variable(zero_poly);
  }





  void clean_up_equations(){
    Polynomial<QQ> left, right;
    Polynomial<QQ> temp, dummy;;
    QQ num, den;
    ZZ thegcd;
    long d, max_dim_pow;
    typename list< Polynomial<QQ> >::iterator it;
    typename map< PowProd, QQ >::iterator pow_it;

    left.alphabet= source;
    right.alphabet= source;
    max_dim_pow= 0;

    for(it= source->relations.generators.begin();
        it != source->relations.generators.end();
        it++){
      temp= *it;
      thegcd= 0;
      for(pow_it= temp.coeffs.begin();
          pow_it != temp.coeffs.end();
          pow_it ++){
        // checks out the dimension of the rep
        dummy= Polynomial<QQ>(pow_it->first);
        d= ZZtolong(dummy.evaluate(source->epsilon));
        if(d > max_dim_pow)
          max_dim_pow= d;
        // clears up denominators
        // (there should not be any anyway)
        num= pow_it->second.get_num();
        den= pow_it->second.get_den();
        temp *= den;
        // updates gcp
        thegcd= gcd(thegcd, num);
        // decides if we write to the LHS or RHS
        if(num > 0)
          left.add_term(pow_it->first, num);
        else
          right.add_term(pow_it->first, -num);
      }
      // divides by gcd
      left*= 1/QQ(thegcd);
      right*= 1/QQ(thegcd);

      equations.push_back( make_pair( left, right ) );
      left.sets_to_zero();
      right.sets_to_zero();
    }
    max_dim_rel= max_dim_pow;    
}




  // methods which do the actual work ***********************
  // ********************************************************

  Polynomial<K> compute_tensor_relation(long max_dim1, long max_dim2){
    long i, j, deg;
    Polynomial<K> P, Q, I, debug;
    Tuple< Polynomial<K> > sigma, sigmap;
    PowProd pow;
    Polynomial<K> ans;
    
    I= target->one();
    
    // compute the elementary symmetric functions
    // in the a_i's and the b_i's.
    sigma.resize(max_dim1 + 1);
    sigmap.resize(max_dim2 + 1);
    P= I;
    Q= I;
    debug= I;
    for(i=1; i<=max_dim1; i++){
      P*= I + a[i];
    }
    for(i=1; i<=max_dim2; i++){
      Q*= I + b[i];
    }

    for(i=1; i<=max_dim1; i++){
      sigma[i]= target->hom_part_of(P,i);
    }
    for(i=1; i<=max_dim2; i++){
      sigmap[i]= target->hom_part_of(Q,i);
    }

    // compute the polynomial that needs to be
    // rewritten as a poly in the sigma_i and sigma_i'

    P= I;

    for(i=1; i<=max_dim1; i++)
      for(j=1; j<=max_dim2; j++)
        P*= I + a[i] + b[j];// done !
    
    
    // do the actual work 
    ans=I; // for the alphabet !!
    ans.sets_to_zero();


    while( !P.is_zero() ){


      pow= P.lp();
      
      // compute a poly in the sigma_i's with the same lp()
     
      Q= I;
      debug= I;
      for(i=1; i<=max_dim1; i++){
        deg= pow.power_of(last_actual_variable + max_dim1 + 1 - i);
        if(i != max_dim1)       
          deg-= pow.power_of(last_actual_variable + max_dim1 - i);
        // now multiply Q by sigma[i] to the power deg
        for(j=1; j<= deg; j++){
          Q *= sigma[i];
          debug *= s[i];
        }
      }
      // now same with sigmaprime
      for(i=1; i<=max_dim2; i++){
        deg= pow.power_of(last_actual_variable + max_dim + max_dim2 + 1 - i);
        if(i != max_dim2)       
          deg-= pow.power_of(last_actual_variable + max_dim + max_dim2 - i);
        // now multiply Q by sigma[i] to the power deg
        for(j=1; j<= deg; j++){
          Q *= sigmap[i];
          debug *= sp[i];
        }
      }
      
      //       cout << "le poly obtenu: " << debug << endl;
      // cout << "qui developpe donne " << Q << endl;
      
      ans+= debug*P.coeff_of(pow);
      P+= Q*(-1)*P.coeff_of(pow);
    }

    return ans;

  }

  void write_lambda_factors(list< Polynomial<K> >& result, long ext, long first, long last){
    list< Polynomial<K> > temp1, temp2;
    typename list< Polynomial<K> >::iterator it;


    // trivial cases: ext too small or too large
    if( (ext < 1) ){
      result.clear();
      result.push_front( target->one() );
      return;
    }
  
    if(  (ext > last - first + 1) ){
      result.clear();
      return;
    }
    
    // else...
    // use recursion:
    write_lambda_factors(result, ext, first + 1, last);
    write_lambda_factors(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 + a[ first ] );
    return;

  }

  Polynomial<K> compute_lambda_relation(long ext, long dim){
    long i, j, deg;
    Polynomial<K> P, Q, I, debug;
    Tuple< Polynomial<K> > sigma;
    PowProd pow;
    Polynomial<K> ans;
    list< Polynomial<K> > splits;
    typename list< Polynomial<K> >::iterator it;

    I= target->one();

    // compute the elementary symmetric functions
    // in the a_i's.
    sigma.resize(dim + 1);
    P= I;
    
    for(i=1; i<=dim; i++){
      P*= I + a[i];
    }

    for(i=1; i<= dim; i++){
      sigma[i]= target->hom_part_of(P,i);
    }


    // compute the polynomial that needs to be
    // rewritten as a poly in the sigma_i

    write_lambda_factors(splits, ext, 1, dim);
    // now just multiply out the factors
    
    for(P= I, it= splits.begin(); it != splits.end(); it++)
      P*= *it;


    // do the actual work 
    ans=I; // for the alphabet !!
    ans.sets_to_zero();
    
    
    while( !P.is_zero() ){
      
      
      pow= P.lp();
      
      // compute a poly in the sigma_i's with the same lp()
      
      Q= I;
      debug= I;
      for(i=1; i<=dim; i++){
        deg= pow.power_of(last_actual_variable + dim + 1 - i);
        if(i != dim)    
          deg-= pow.power_of(last_actual_variable + dim - i);
        // now multiply Q by sigma[i] to the power deg
        for(j=1; j<= deg; j++){
          Q *= sigma[i];
          debug *= s[i];
        }
      }
        
      ans+= debug*P.coeff_of(pow);
      P+= Q*(-1)*P.coeff_of(pow);
    }
   
    return ans;
  }




  void exp_classes_after_tensoring(const Polynomial<K>& exp,
                                   long dim,
                                   const PowProd& pow,
                                   Polynomial<K>& result) {

    long i, j, d;
    PowProd new_pow;
    Polynomial<K> universal;
   

    if(pow == PowProd(0)){ // trivial case
      result= exp;
      return;
    }
    // else...
    // let's look at the first variable in pow 
    // with positive power;
    for(i=1; i<= pow.nvars(); i++)
      if( pow.power_of(i) > 0 )
        break; //done !
    d= ZZtolong(source->epsilon[i]);
    // get the universal polynomial for a tensor prod
    // already computed?
    if( univ_computed.count( make_pair(d,dim) )==0 )
      univ_computed[ make_pair(d,dim) ] = compute_tensor_relation(d, dim);
    
    universal=  univ_computed[ make_pair(d,dim) ] ;
    // do the substitutions:
    // s[j] gets replaced by var(i,j)
    for(j=1; j<= d; j++)
      universal.substitute_variable(last_actual_variable + max_dim
                                 + max_dim_rel +j, var(i,j) );
    // sp[j] gets replaced by hom part of(exp)
    for(j=1; j<= dim; j++)
      universal.substitute_variable(last_actual_variable + 2*max_dim
                                    + max_dim_rel + j, 
                                    target->hom_part_of(exp,j*factor));
    
    // good ! now we have computed the exp classes of the given
    // rep tensored with the rep given by variable i in *source.
    // now use recursion:
    new_pow= pow/PowProd(i);
    exp_classes_after_tensoring(universal, d*dim, new_pow, result);
    return;
  }


  Polynomial<K> exp_classes_of_ext_power(long thevar, long ext){
    long dim, j;
    Polynomial<K> universal;

    dim= ZZtolong(source->epsilon[thevar]);

    // get the universal polynomial for an exterior power
    // already computed?
    if( lambda_univ_computed.count( make_pair(ext,dim) )==0 )
      lambda_univ_computed[ make_pair(ext,dim) ] = compute_lambda_relation(ext, dim);
    
    universal=  lambda_univ_computed[ make_pair(ext,dim) ] ;
    // do the substitutions:
    // s[j] gets replaced by var(thevar,j)
    for(j=1; j<= dim; j++)
      universal.substitute_variable(last_actual_variable + max_dim
                                 + max_dim_rel +j, var(thevar,j) );
    return universal;
  }


  void compute_exp_relations(){
    Polynomial<K> templeft, tempright, P, Q;
    long i, j, k, n;
    typename list< pair< Polynomial<QQ>, Polynomial<QQ> > >::iterator it;
    typename map< PowProd, QQ>::iterator it_pow;
   
    //we go through the multiplicative equations one by one

    if(target->verbose)
      cout << endl 
           <<"{\\bf Translating the relations}"
           << endl << endl;

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

      if(target->very_verbose)
        cout << "we examine equation $" << it->first
             << " = " << it->second 
             << "$" << endl << endl;

      templeft= target->one();
      tempright= target->one();
      //we go through the powprods on the left
      for(it_pow= it->first.coeffs.begin();
          it_pow!= it->first.coeffs.end();
          it_pow++){
        // we compute in P the exp classes of the monomial
        exp_classes_after_tensoring(target->one(), 1, it_pow->first, P);
        // we put this to the appropriate power
        n= ZZtolong(it_pow->second);
        Q= P;
        for(i=2; i<= n; i++)
          P*= Q;
        //we multiply templeft by P
        templeft*= P;
      }
      //we go through the powprods on the right
      for(it_pow= it->second.coeffs.begin();
          it_pow!= it->second.coeffs.end();
          it_pow++){
        // we compute in P the exp classes of the monomial
        exp_classes_after_tensoring(target->one(), 1, it_pow->first, P);
        // we put this to the appropriate power
        n= ZZtolong(it_pow->second);
        Q= P;
        for(i=2; i<= n; i++)
          P*= Q;
        //we multiply tempright by P
        tempright*= P;
      }
      // we add the equation templeft-tempright = 0
      templeft+= tempright*(-1);
      total_relations.push_back(templeft);

      if(target->very_verbose)
        cout << "translated into $" << templeft
             << " = 0" << "$" << endl << endl;

    }// end of main for loop

    // next, we examine the lambda relations

    for(i= 1; i <= source->variables_in_use(); i++)
      for(j= 2; j <= source->epsilon[i]; j++){
        // we are now looking at lambda^j( x_i )
      if(target->very_verbose)
        cout << "we examine equation $\\lambda^" << j
             << "(" << source->nameof(i) << ") = "
             << source->lambda_operations[i][j] 
             << "$" << endl << endl;


        templeft= exp_classes_of_ext_power(i, j);
        tempright= target->one();
        // we explore the pows in the expression 
        // for lambda^j(x_i) given in *source   
        for(it_pow= source->lambda_operations[i][j].coeffs.begin();
            it_pow != source->lambda_operations[i][j].coeffs.end();
            it_pow++){
        exp_classes_after_tensoring(target->one(), 1, it_pow->first, P);
        n= ZZtolong( it_pow->second );
        if( n < 0 ) // then we multiply templeft by P^-n
          for(k=1; k <= -n; k++)
            templeft*= P;
        else        // we multiply tempright by P^n
          for(k=1; k <= n; k++)
            tempright*= P;
        }
      // we add the equation templeft-tempright = 0
      templeft+= tempright*(-1);
      total_relations.push_back(templeft);

      if(target->very_verbose)
        cout << "translated into $" << templeft
             << " = 0$" << endl << endl;
      }



  }

  void split_and_sort(){
    long i;
    Polynomial<K> P;
    typename list< Polynomial<K> >::iterator total_it;
    
    // now we split the equations in total_relations
    // into homogeneous parts, sort them according
    // to their degrees, and add them to exp_relations

    if(target->verbose)
      cout << endl 
           << "{\\bf Now splitting into homogeneous parts}" 
           << endl << endl;

    // now we split the relations we have just
    // obtained into homogeneous parts, sorted by degree
    i= 0;
    while( !total_relations.empty() ){
      i++;

      if(target->very_verbose)
        cout << "adding the relations in degree " << i 
             << endl << endl;

      for(total_it= total_relations.begin();
          total_it != total_relations.end();){
        P= target->hom_part_of(*total_it, i);
        if( !P.is_zero() ){
          if( target->very_verbose)
            cout << "adding $" << P << " = 0$" 
                 << endl << endl;
          exp_relations.push_back(P);
          *total_it += P*(-1);
        }
        if( total_it->is_zero() )
          total_it= total_relations.erase(total_it);
        else
          total_it++;
      }
    }// end of WHILE



  }



  void add_relations(){

    target->relations.generators= exp_relations;
    exp_relations.clear();
    total_relations.clear();
  }


  void start(RepresentationRing *thesource,
             GradedAlgebra<K> *thetarget,
             string prefix,
             long thefactor){

    source= thesource;
    target= thetarget;
    factor= thefactor;

    this->create_regular_variables(prefix);
    this->clean_up_equations();
    this->create_extra_variables(prefix);
    this->compute_exp_relations();
    this->split_and_sort();
    this->kill_extra_variables();
    this->add_relations();
  }

};


#endif

---