polynomials.H

Go to the documentation of this file.

#ifndef _POLYNOMIALS_H_
#define _POLYNOMIALS_H_

#include "powerprod.H"
#include <map>
#include <iostream>
#include <fstream>
using namespace std;


// a silly little something
template <class K>
K _raise_to_power(K c, long n){
  K ans;

  ans= 1;
  for(long i=1; i<= n;i++)
    ans*= c;

  return ans;
}


template <class K> class Polynomial {
private:

  static int order_override;
  static long variable_priority;

public:
  Alphabet *alphabet;
  map<PowProd, K> coeffs;

  // constructors *******************************************
  //*********************************************************
  Polynomial<K>(){
    alphabet= PowProd::get_default_alphabet();
  }

  Polynomial<K>(const Polynomial<K>& that){
    coeffs= that.coeffs;
    alphabet= that.alphabet;
  }
  Polynomial<K>(const PowProd& pow){
    coeffs[pow]=1;
    alphabet= PowProd::get_default_alphabet();
  }

  // methods to handle the static data, 
  //order_override and variable_priority ************************
  //*************************************************************
  static void set_order_override(int theorder){
    order_override= theorder;
  }

  static void default_order(){
    order_override= 0;
  }

  static int get_order_override(){
    return order_override;
  }

  static void set_variable_priority(long thepriority){
    variable_priority= thepriority;
  }

  static void no_variable_priority(){
    variable_priority= 0;
  }

  static long get_variable_priority(){
    return variable_priority;
  }

  // a few things to do with the polynomial being 0 ******************
  //*******************************************************************
  bool is_zero() const {
    typename map<PowProd,K>::const_iterator it;
    
    it=coeffs.begin();
    return (it == coeffs.end());
  }
  void sets_to_zero() {
    coeffs.clear();
  }

  // some checks for variables showing up *****************************
  //*******************************************************************
  bool involves_variable(long i) const {
    typename map<PowProd,K>::const_iterator it;

    for(it= coeffs.begin(); it!= coeffs.end(); it++)
      if( it->first.power_of(i) > 0 )
        return true;

    return false;
  }

  long nvars() const {
    long ans=0;
    typename map<PowProd,K>::const_iterator it;

    for(it= coeffs.begin(); it!= coeffs.end(); it++)
      if(it->first.nvars() > ans)
        ans= it->first.nvars();

    return ans;
  }



  long has_lonely_variable(Polynomial<K>& replacement){
    typename map<PowProd,K>::const_iterator it;
    long var;

    for(it= coeffs.begin(); it!= coeffs.end(); it++){
      if( it->first.degree()==1 ){
        var= it->first.nvars();
        replacement= *this;
        replacement.add_term( it->first, - it->second);
        if(!replacement.involves_variable(var)){
          replacement*= (-1/it->second);
          return var;
        }
      }
    }

    return 0;
  }



  long has_unitary_variable(Polynomial<K>& replacement){
    typename map<PowProd,K>::const_iterator it;
    long var;

    for(it= coeffs.begin(); it!= coeffs.end(); it++){
      if( it->first.degree()==1 && (it->second == 1 || it->second == -1) ){
        var= it->first.nvars();
        replacement= *this;
        replacement.add_term( it->first, - it->second);
        if(!replacement.involves_variable(var)){
          replacement*= (-1/it->second);
          return var;
        }
      }
    }

    return 0;
  }




  // returning information on the polynomial **************************
  //*******************************************************************

  K coeff_of(const PowProd& pow) const {
    typename map<PowProd,K>::const_iterator it;

    it= coeffs.find(pow);
    if( it != coeffs.end() )
      return it->second;
    else 
      return 0;
  }

  PowProd lp() const {
    typename map<PowProd,K>::const_iterator it, itprime;

    if( order_override == 0){// default: we return the last powprod
      it= coeffs.end(); // this is the leading powprod with respect
      if( it == coeffs.begin() ) // to the order specified by
        return PowProd( 0 ); // PowProd::order.
      else{
        it--;
        return it->first;
      }
    }
    else{ // using an alternative order
    
      itprime= coeffs.begin();
      it= coeffs.begin();
      if( it == coeffs.end() )
        return PowProd( 0 ); // return 1 if poly is zero
      

      for(it++; it != coeffs.end(); it++){
        if( powprod_compare(itprime->first, it->first, order_override, variable_priority) )
          itprime= it;
      }

      return itprime->first;
      }
  }


  Polynomial<K> lt() const {
    PowProd leading;
    Polynomial<K> ans;
    
    leading= lp();
    ans= Polynomial<K>( leading );
    ans.alphabet= this->alphabet;
    ans*= coeff_of( leading );
    return ans;
  }

  K lc() const {
    return coeff_of( lp() );    
  }




  // arithmetic, operators *****************************************
  //*******************************************************************

  void add_term(const PowProd& pow, const K& c){

    if( coeffs.count( pow ) == 0)
      coeffs[pow]=c;
    else
      coeffs[pow]+=c;
    if(coeffs[pow] == 0)
      coeffs.erase(pow);
  }

  Polynomial<K>& operator+=(const Polynomial<K>& that){
    typename map<PowProd,K>::const_iterator it;

    for(it= that.coeffs.begin(); it != that.coeffs.end(); it++)
      add_term(it->first, it->second);      

    return *this;
  }

 Polynomial<K> operator+(const Polynomial<K>& that) const {
   Polynomial<K> ans;

   ans= *this;
   ans+=that;
   return ans;

 }


 
  Polynomial<K>& operator*=(const K& c){
    typename map<PowProd,K>::const_iterator it;

    if( c == 0 )
      sets_to_zero();
    else{
      for(it=coeffs.begin(); it != coeffs.end(); it++)
        coeffs[it->first] *= c;
    }

    return *this;
  }

  friend void polyprod(const Polynomial<K>& P, const Polynomial<K>& Q, Polynomial<K>& R){
    PowProd pow;
    K c;
    typename map<PowProd,K>::const_iterator i_this, i_that;

    R.coeffs.clear();
    for(i_this= P.coeffs.begin(); i_this != P.coeffs.end(); i_this ++)
      for(i_that=Q.coeffs.begin(); i_that != Q.coeffs.end(); i_that++){
        pow= i_this->first;
        pow *= i_that->first;
        c= i_this->second;
        c*= i_that->second;
        R.add_term(pow,c);
      }
    R.alphabet= P.alphabet;
  }

  Polynomial<K> operator*(const Polynomial<K>& that) const {
    Polynomial<K> ans;

    polyprod(*this, that, ans);
    return ans;
  }

  Polynomial<K>& operator*=(const Polynomial<K>& that){
    Polynomial<K> temp;

    polyprod(*this, that, temp);
    coeffs= temp.coeffs;   

    return *this;
  }


  Polynomial<K> operator*(const K& c) const {
    Polynomial<K> ans;
    
    ans= *this;
    ans *= c;
    return ans;
  }
  
  // a couple of things related to substitutions *********************
  //*******************************************************************
  void swap_variables(long i, long j){
    typename map<PowProd,K>::const_iterator it;
    Polynomial<K> temp;
    PowProd pow;
    
    for(it= coeffs.begin(); it != coeffs.end(); it++){
      pow= it->first;
      pow.swap(i,j);
      temp.add_term(pow, it->second);
    }

    coeffs= temp.coeffs;
  }

  void shift_variables(long n){
    typename map<PowProd,K>::const_iterator it;
    Polynomial<K> temp;
    PowProd pow;
    
    for(it= coeffs.begin(); it != coeffs.end(); it++){
      pow= it->first;
      pow.shift(n);
      temp.add_term(pow, it->second);
    }
    
    coeffs= temp.coeffs;
  }



  void substitute_variable(long i, const Polynomial<K>& f) {
    Tuple< Polynomial<K> > computed;
    long max_degree, n;
    typename map<PowProd,K>::const_iterator it;
    Polynomial<K> ans;
    PowProd pow;

    max_degree= 0;
    for(it= coeffs.begin(); it != coeffs.end(); it++){
      n= it->first.power_of(i);
      if(n > max_degree)
        max_degree= n;
    }

    if(max_degree==0)
      return;

    // first set computed[n]=f^n
    computed.resize(max_degree + 1);
    computed[1]= f;
    for(n=2; n <= max_degree; n++)
      computed[n]= f * computed[n-1]; 

    for(it= coeffs.begin(); it != coeffs.end(); it++){
      n= it->first.power_of(i);
      if(n==0)
        ans.add_term(it->first, it->second);
      else{
        pow= it->first;
        pow/=PowProd(i,n);
        ans+= computed[n] * Polynomial<K>( pow ) * it->second;
      }
    }// end of for

    coeffs= ans.coeffs;
  }

  K evaluate(const map<long,K>& augmentation) const {
    K ans, a, b;
    long i;
    typename map<PowProd,K>::const_iterator it;
    
    ans= 0;
    for(it= coeffs.begin(); it != coeffs.end(); it++){
      a= it->second;      
      for(i=1; i<= it->first.nvars(); i++){
        b= augmentation.find(i)->second;
        a*= _raise_to_power(b,it->first.power_of(i));
      }
      ans += a;
    }

    return ans;
  }

  // friends, creating polynomials ***********************************
  //*******************************************************************
  // creates a one-letter polynomial on a variable with specified name
  /*
   * See alphabet.H for an example of code using this.
   * DELETED! becomes ambiguous when several K's are
   * used in the same program.
   */
  /*
  friend Polynomial<K> new_indeterminate(SimpleAlphabet& alpha, const string& name){
    Polynomial<K> ans;

    ans= Polynomial<K>( PowProd( alpha.new_variable(name) ) );
    ans.alphabet= &alpha;

    return ans;
  }
  // returns the constant polynomial, equal to 1, on the alphabet alpha.
  friend Polynomial<K> constant_polynomial(const SimpleAlphabet& alpha){
    Polynomial<K> ans;

    ans= Polynomial<K>( PowProd(0) );
    ans.alphabet= (Alphabet*) &alpha;

    return ans;
  }
  */
  // friends printing **************************************************
  //*******************************************************************

  friend ostream& operator<<(ostream& os, const Polynomial<K>& P){
    
    typename map<PowProd,K>::const_iterator it;
    Alphabet *temp;

    temp= PowProd::get_current_alphabet();
    PowProd::set_alphabet(P.alphabet);
    
    // prints 0 if the polynomial is 0...
    if( P.is_zero() ){
      os << "0";
      return os;
    }
    // else...
    for(it=P.coeffs.begin(); it != P.coeffs.end();){
      if( it-> second != 1)
        os << it->second;
      os << it->first;

      //it->first.debug();

      if( (++it) != P.coeffs.end() )
        os << " + ";
    }
    PowProd::set_alphabet(temp);
    return os;
  }

  friend ofstream& operator<<(ofstream& file, Polynomial<K>& P){
    typename map< PowProd, K >::iterator it;
    K scalar;
    PowProd pow;
    
    
    for(it= P.coeffs.begin(); it!= P.coeffs.end(); it++){
      scalar= it->second;
      pow= it->first;
      file << scalar;
      file << " ";
      file << pow;
    }
    
    file << 0 << " ";
    return file;
  }
  
  friend void operator>>(ifstream& file, Polynomial<K>& P){
    PowProd pow;
    K scalar;
    
    P.coeffs.clear();
    
    file >> scalar;

    while( scalar != 0){
      file >> pow;
      P.add_term(pow, scalar);
      file >> scalar;

    }
    
  }




  // friend computing the S polynomial ********************************
  //*******************************************************************

  friend Polynomial<K> S_polynomial(const Polynomial<K>& f, const Polynomial<K>& g){
    PowProd a,b,c;
    K lambda, mu;
    Polynomial<K> ans, temp;

    a= f.lp(); 
    b= g.lp(); 
    c= lcm(a,b); 
    lambda= 1 / f.lc();
    mu= -1 / g.lc();
    a= c/a; 
    b= c/b; 
    // ans = f*a*lambda + g*b*mu
    ans= f; 
    ans*= Polynomial<K>( a ); 
    ans*= lambda; 
    temp= g; 
    temp*= Polynomial<K>( b ); 
    temp*= mu; 
    ans+= temp; 
    ans.alphabet= f.alphabet;
    
    return ans;
  }
  friend Polynomial<K> S_polynomial_with_coefficients(const Polynomial<K>& f, 
                                                      const Polynomial<K>& g,
                                                      pair< Polynomial<K>, Polynomial<K> >& coeffs){
    PowProd a,b,c;
    K lambda, mu;
    Polynomial<K> ans, temp;

    a= f.lp(); 
    b= g.lp(); 
    c= lcm(a,b); 
    lambda= 1 / f.lc();
    mu= -1 / g.lc();
    a= c/a; 
    b= c/b; 
    // ans = f*a*lambda + g*b*mu
    coeffs.first= Polynomial<K>(a) * lambda;
    coeffs.first.alphabet= f.alphabet;
    coeffs.second= Polynomial<K>(b) * mu;
    coeffs.second.alphabet= f.alphabet;
    ans= f * coeffs.first + g * coeffs.second; 
    ans.alphabet= f.alphabet;
    
    return ans;
 }
  
  
};

template <class K> int Polynomial<K>::order_override=0;
template <class K> long Polynomial<K>::variable_priority=0;


#endif

---