/* Copyright (c) 1997-2006
   Ewgenij Gawrilow, Michael Joswig (Technische Universitaet Berlin, Germany)
   http://www.math.tu-berlin.de/polymake,  mailto:polymake@math.tu-berlin.de

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version: http://www.gnu.org/licenses/gpl.txt.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
*/

#ifndef _POLYMAKE_POLYNOMIAL_H
#define _POLYMAKE_POLYNOMIAL_H "$Project: polymake $$Id: Polynomial.h,v 1.3 2006/03/09 16:21:36 gawrilow Exp $"

#include <Rational.h>
#include <Monomial.h>
#include <Vector.h>
#include <GenericVector.h>
#include <SparseMatrix.h>

#include <ext/hash_map>
#include <string>
#include <sstream>

namespace pm {

template <typename Coefficient=Rational, typename Exponent=int>
class Polynomial {

  // no coefficients == 0, unless it's the trivial polynomial
  // no monomial appears twice

protected:
   typedef Monomial<Exponent> monom;
   typedef std_ext::hash_map<monom, Coefficient> term_hash;

   term_hash the_terms;
   mutable int n_variables;  // lazy: -1 <=> not computed
   mutable monom the_lm;
   mutable bool the_lm_set;

public:
   typedef typename term_hash::const_iterator const_iterator;
   typedef typename term_hash::iterator iterator;
   typedef typename iterator::reference reference;
   typedef reference const_reference;
  
   Polynomial()
   {
      // default is the polynomial x=0
      the_terms[ monom() ] = Coefficient(0);
      n_variables = 0;
   }

   void clear()
   {  // the default polynomial
      the_terms.clear();
      the_terms[ monom() ] = Coefficient(0);
      n_variables = 0;
      the_lm.clear();
      the_lm_set = true;
   }

   const_iterator lt() const
   {
      compute_lm();
      return find(the_lm);
   }

   const monom& lm() const
   {
      compute_lm();
      return the_lm;
   }

   Coefficient lc() const
   {
      compute_lm();
      return exists(the_lm);
   }

   int vars() const
   {
      compute_vars();
      return n_variables;
   }
  
   int terms() const { return the_terms.size(); }
   const_iterator begin() const { return the_terms.begin(); }
   const_iterator end() const { return the_terms.end(); }
   int size() const { return the_terms.size(); }
   bool empty() const { return the_terms.empty(); }

   SparseMatrix<Exponent> monomials_as_matrix()  const
   {
      SparseMatrix<Exponent> mon;
      if (n_variables != -1)
	 mon.resize(0,n_variables);

      for (typename Entire<term_hash>::const_iterator m=entire(the_terms);
	   !m.at_end(); ++m) {
	 SparseVector<Exponent> mm = m->first.monomial();
	 if (mm.dim() > mon.cols())
	    mon.resize(mon.rows(),mm.dim());

	 mm.resize(mon.cols());
	 mon /= mm;
      }
    
      n_variables = mon.cols();
      return mon;
   }
  
   Vector<Coefficient> coefficients_as_vector() const
   {
      Vector<Coefficient> coe(the_terms.size());
      typename term_hash::const_iterator t_it=the_terms.begin();
      for (typename Entire< Vector<Coefficient> >::iterator c_it=entire(coe);
	   !c_it.at_end(); ++c_it, ++t_it)
	 *c_it = t_it->second;

      return coe;
   }

   bool trivial() const
   {
#ifdef POLY_DEBUG
      if ( the_terms.size()==1 && the_terms.begin()->second==Coefficient(0) && 
	   (!the_terms.begin()->first.empty() || n_variables>0) )
	 throw std::runtime_error("Polynomial: Wrong representation of trivial polynomial.");
#endif

      return the_terms.size()==1 && the_terms.begin()->second==Coefficient(0);
   }

   std::string pretty_print() const
   {
      std::ostringstream str;
      bool first = true;
      for (typename Entire<term_hash>::const_iterator t=entire(the_terms);
	   !t.at_end(); ++t) {
	 if (first)
	    first = false;
	 else
	    str << "  +  ";

	 if (t->second != Coefficient(1) || t->first.empty())
	    str << t->second;
	 if (t->second != Coefficient(1) && !t->first.empty())
	    str << " * ";
	 if (!t->first.empty())
	    str << t->first.pretty_print();
      }

      return str.str();
   }
  
protected:
   Polynomial(const term_hash& t, const int vars=-1, const monom& lm=monom(), const bool lm_set=false)
   {
      the_terms = t;
      n_variables = vars;
      the_lm = lm;
      the_lm_set = lm_set;
   }
  
   const term_hash& give_terms() const { return the_terms; }

   void compute_lm() const
   {
      if (trivial())
	 throw std::runtime_error("Polynomial: Zero polynomial has no leading term.");

      if (!the_lm_set) {
	 const_iterator it=the_terms.begin();
	 the_lm = it->first;  ++it;
	 for ( ; it!=the_terms.end(); ++it)
	    if (it->first>the_lm)
	       the_lm = it->first;
      }

      the_lm_set = true;
   }

   void compute_vars() const
   {
      if (n_variables==-1) {
	 n_variables = 0;
	 for (typename Entire<term_hash>::const_iterator it=entire(the_terms);
	      !it.at_end(); ++it)
	    if (n_variables < it->first.dim())
	       n_variables = it->first.dim();
      }
   }

#ifdef POLY_DEBUG
public:
#endif

   bool lt_set() const { return the_lm_set; }

   bool vars_set() const { return n_variables!=-1; }

#ifdef POLY_DEBUG
protected:
#endif

   void add_term(const monom& m, const Coefficient& c)
   {
      if (c!=0) {
	 typename term_hash::iterator find = the_terms.find(m);

	 if (find!=the_terms.end())  // monom exists allready
	    if (-c==find->second)
	       if (the_terms.size()==1)
		  clear();  // create the trivial polynomial
	       else {
		  if (the_lm_set && the_lm==m)
		     the_lm_set = false;
		  if (m.dim()==n_variables)
		     n_variables = -1;
		  the_terms.erase(m);
	       }
	    else
	       find->second += c;
      
	 else {  // new monomial
	    if (trivial()) {
	       the_terms.clear();
	       the_terms[m] = c;
	       n_variables = m.dim();
	       the_lm = m;
	       the_lm_set = true;

	    } else {
	       if (the_lm_set && m>the_lm)
		  the_lm = m;

	       if (n_variables!=-1 && m.dim()>n_variables)
		  n_variables = m.dim();

	       the_terms[m] = c;
	    }
	 }
      }
   }
  
   template <typename _Vector>
   void add_term(const GenericVector<_Vector,Exponent>& m, const Coefficient& c)
   {
      add_term(monom(m),c);
   }

public:
   template <typename SequenceType1, typename SequenceType2>
   Polynomial(const SequenceType1& mon, const SequenceType2& coe)
   {    
      the_terms[ monom() ] = Coefficient(0);
      n_variables = 0;

      DIMENSION_CHECK(size_t(coe.size()) != size_t(mon.size()), "Polynomial: dimension mismatch");
    
      typename SequenceType1::const_iterator m=mon.begin();
      for (typename Entire<SequenceType2>::const_iterator c=entire(coe);
	   !c.at_end(); ++m,++c)
	 add_term(*m,*c);
   }

   template <typename Iterator>
   explicit Polynomial(Iterator b, typename enable_if<void, check_iterator_feature<Iterator,end_sensitive>::answer>::type* =0)
   {
      the_terms[ monom() ] = Coefficient(0);
      n_variables = 0;

      for (; !b.at_end(); ++b)
	 add_term(b->first,b->second);
   }

   // returns the coefficient of the monomial m or 0 iff m does not exists
   Coefficient exists(const monom& m) const
   {
      const_iterator find = the_terms.find(m);
      if (find==the_terms.end())
	 return Coefficient(0);

      return find->second;
   }

   template <typename _Vector>
   Coefficient exists(const GenericVector<_Vector,Exponent>& m) const
   {
      return exists(monom(m));
   }

   const_iterator find(const monom& m) const
   {
      return the_terms.find(m);
   }

   template <typename _Vector>
   const_iterator find(const GenericVector<_Vector,Exponent>& m) const 
   {
      return the_terms.find(monom(m));
   }

   Polynomial operator* (const Coefficient& a) const
   {
      term_hash prod = the_terms;
      for (typename Entire<term_hash>::iterator t=entire(prod); !t.at_end(); ++t)
	 t->second *= a;

      return Polynomial(prod,n_variables,the_lm,the_lm_set);
   }

   friend Polynomial operator* (const Coefficient& a, const Polynomial& p)
   {
      term_hash prod = p.give_terms();
      for (typename Entire<term_hash>::iterator t=entire(prod); !t.at_end(); ++t)
	 t->second = a*t->second;

      const int p_vars= p.vars_set() ? p.vars() : -1;
      const monom p_lm= p.lt_set() ? p.lm() : monom();

      return Polynomial(prod,p_vars,p_lm,p.lt_set());
   }

   Polynomial  operator* (const Polynomial& p) const
   {
      Polynomial prod;
      for (typename Entire<term_hash>::const_iterator t1=entire(the_terms);
	   !t1.at_end(); ++t1)
	 for (typename Entire<Polynomial>::const_iterator t2=entire(p);
	      !t2.at_end(); ++t2)
	    prod.add_term(t1->first*t2->first, t1->second*t2->second);

      return prod;
   }

   Polynomial operator+ (const std::pair<monom,Coefficient>& t)
   {
      Polynomial sum = *this;
      sum.add_term(t.first,t.second);
      return sum;
   }

   Polynomial& operator+= (const std::pair<monom,Coefficient>& t)
   {
      add_term(t.first,t.second);
      return *this;
   }

   Polynomial operator+ (const Polynomial& p) const
   {
      Polynomial sum=p;
      for (typename Entire<term_hash>::const_iterator t=entire(the_terms);
	   !t.at_end(); ++t)
	 sum.add_term(t->first, t->second);
    
      return sum;
   }

   Polynomial operator- () const
   {
      return operator*(Coefficient(-1));
   }

   Polynomial operator- (const std::pair<monom,Coefficient>& t)
   {
      Polynomial sum = *this;
      sum.add_term(t.first,-t.second);
      return sum;
   }

   Polynomial& operator-= (const std::pair<monom,Coefficient>& t)
   {
      add_term(t.first,-t.second);
      return *this;
   }

   Polynomial operator- (const Polynomial& p) const
   {
      return operator+( p.operator-() );
   }
  
   bool operator== (const Polynomial& p) const
   {
      if (terms()!=p.terms())
	 return false;

      for (typename Entire<Polynomial>::const_iterator t=entire(p);
	   !t.at_end(); ++t)
	 if (exists(t->first) != t->second)
	    return false;

      return true;
   }

   // FIXME
   // all operators with single term
   // this* operator*= (const Coefficient& a)
   // this* operator*= (const Polynomial& p)
   // this* operator+= (const Polynomial& p)
   // this* operator-= (const Polynomial& p)
   // IO operator
   // efficiency check
   // operator>
   // operator<

#if _POLYMAKE_ICC_ONLY
   // the hash_map implementation by Intel is bidirectional,
   // but this makes absolutely no sense here
   typedef forward_iterator_tag category;
#endif
};

#if _POLYMAKE_ICC_ONLY
template <typename Coefficient, typename Exponent>
struct container_has_category< Polynomial<Coefficient,Exponent> > : True { };
#endif

} // end namespace pm

namespace polymake {
   using pm::Polynomial;
}

#endif // _POLYMAKE_POLYNOMIAL_H


// Local Variables:
// mode:C++
// c-basis-offset:3
// End:
