polynomial.h

#ifndef POLYNOMIAL_H_
#define POLYNOMIAL_H_

#include <list>
#include <utility>
#include "./param.h"
#include "./monomial.h"

namespace tinygb {

template<class T>
class Polynomial;

template<class T>
std::ostream& operator<<(std::ostream &o, const Polynomial<T> &p);

template<class T>
class Polynomial {
 public:
  typedef std::pair<Monomial, T> Term;

  Polynomial() {
      terms_ = std::list<Term>();
    }
  explicit Polynomial(const std::list<Term>& terms_)
    : terms_(terms_) {}

  const std::list<Term>& terms() const {
    return terms_;
  }

  // Put the list of terms in standard form. Depends on the monomial ordering.
  // The list is ordered from smallest to largest
  // TODO(pj): leading monomial should be first
  void Normalize();

  void SetToZero() {
    terms_.clear();
  }

  void insert(Monomial m, T e) {
    terms_.push_back(std::pair<Monomial, T>(m, e));
  }

  // Assumes that the polynomial's list of terms is normalized
  const Monomial& LeadingMonomial() const {
    assert(!terms_.back().second.is_zero());
    return terms_.back().first;
  }

  const T& LeadingCoefficient() const {
    return terms_.back().second;
  }

  void DivideByLeadingCoefficient();

  unsigned Degree() const {
    // TODO(pj): Extend for non-degree orderings.
    return LeadingMonomial().Degree();
  }

#ifdef BILIN
  unsigned Degree_bilin() const;
#endif

  // Reduces leading monomial by the polynomial given.
  // Returns true if leading monomial is divisible by the leading monomial of
  // the argument. When this function returns false, then nothing has been
  // modified.
  // Assumes that terms_ is normalized.
  bool TopReductionBy(const Polynomial &);

  // Returns true if the leading monomial is divisible by the leading monomial
  // of the argument.
  // Assumes that terms_ is normalized.
  bool IsTopReducibleBy(const Polynomial &) const;

  Polynomial<T> operator*(const Monomial&) const;
  void operator+=(const Polynomial<T>&);
  void operator-=(const Polynomial<T>&);
  void operator*=(const T&);
  void operator/=(const T&);
  void Print(std::ostream& out) const;
  inline bool operator<(const Polynomial<T>& P) const;

  friend std::ostream& operator<< <T>(std::ostream&, const Polynomial<T>&);

 private:
  // Adds coefficients of consecutive terms if they share the same monomials.
  // Useful for adding polynomials. Do not use if terms_ is not sorted.
  void CombineMonomials();
  bool IsEmpty() const {
    return terms_.empty();
  }
  // Removes terms with coefficient 0
  void RemoveZeroTerms();

  std::list<Term> terms_;

};

#ifdef BILIN
template<class T>
unsigned Polynomial<T>::Degree_bilin() const {
  unsigned R = 0;
  for (auto& it = std::as_const(terms_))
    if (it->first.Degree_bilin() > R)
      R = it->first.Degree_bilin();
  return R;
}
#endif

template<class T>
void Polynomial<T>::RemoveZeroTerms() {
  for (auto it = terms_.begin(); it != terms_.end(); ++it) {
    while (it != terms_.end() && it->second == T(0))
      it = terms_.erase(it);
  }
}

template<class T>
void Polynomial<T>::CombineMonomials() {
  typename std::list<std::pair<Monomial, T> >::iterator it2;
  for (auto it1 = terms_.begin(); it1 != terms_.end();) {
    it2 = it1;
    ++it1;
    if (it1 == terms_.end()) {
      break;
    } else if (it1->first == it2->first) {
      it1->second += it2->second;
      terms_.erase(it2);
    }
  }
}

template<class T>
void Polynomial<T>::DivideByLeadingCoefficient() {
  T leading_coefficient = terms_.back().second;
  assert(leading_coefficient != T(0));
  leading_coefficient.inv();
  if (leading_coefficient == T(1)) return;
  for (auto& it : terms_)
    it.second *= leading_coefficient;
}

template<class T>
void Polynomial<T>::Normalize() {
  terms_.sort();
  CombineMonomials();
  RemoveZeroTerms();
}

template<class T>
bool Polynomial<T>::TopReductionBy(const Polynomial &g) {
  assert(!g.IsEmpty());
  if (IsEmpty()) return false;
  if (!LeadingMonomial().IsDivisibleBy(g.LeadingMonomial())) return false;
  Polynomial reductor = g*(LeadingMonomial()/g.LeadingMonomial());
  reductor.DivideByLeadingCoefficient();
  reductor *= LeadingCoefficient();
  operator-=(reductor);
  return true;
}

template<class T>
bool Polynomial<T>::IsTopReducibleBy(const Polynomial<T> &g) const {
  assert(!g.IsEmpty());
  if (IsEmpty()) return false;
  return LeadingMonomial().IsDivisibleBy(g.LeadingMonomial());
}

template<class T>
Polynomial<T> Polynomial<T>::operator*(const Monomial &m) const {
  Polynomial R;
  R.terms_.insert(R.terms_.begin(), this->terms_.begin(), this->terms_.end());
  for (auto& it : R.terms_)
    it.first *= m;
  return R;
}

template<class T>
void Polynomial<T>::operator+=(const Polynomial<T> &P) {
  std::list<Term> tmp_terms = P.terms_;
  terms_.merge(tmp_terms);
  CombineMonomials();
  RemoveZeroTerms();
}

template<class T>
void Polynomial<T>::operator*=(const T &e) {
  for (auto& it : terms_)
    it.second *= e;
}

template<class T>
void Polynomial<T>::operator/=(const T& e) {
  for (auto& it : terms_)
    it.second /= e;
}

template <class T>
void Polynomial<T>::operator-=(const Polynomial<T>& P) {
  std::list<Term> tmp_terms = P.terms_;
  for (auto& it : tmp_terms)
    it.second.neg();
  terms_.merge(tmp_terms);
  CombineMonomials();
  RemoveZeroTerms();
}

template <class T>
void Polynomial<T>::Print(std::ostream& o) const {
  for (auto& it : std::as_const(terms_)) {
    o << it.second << " ";
    it.first.Print(o);
    o << std::endl;
  }
  o << ";";
}

template <class T>
std::ostream& operator<<(std::ostream& o, const Polynomial<T>& P) {
  if (P.IsEmpty()) {
    o << "0";
    return o;
  }
  auto it = std::as_const(P.terms_).begin();
  if (it == P.terms_.end())
    return o;
  if (it->second != T(0))
    o << it->second << it->first;
  ++it;
  for (; it != P.terms_.end() ; ++it)
    if (it->second != T(0))
      o << "+" << it->second << it->first;
  return o;
}

template <class T>
inline bool Polynomial<T>::operator<(const Polynomial<T>& P) const {
  return LeadingMonomial() < P.LeadingMonomial();
}
}  // namespace tinygb

#endif  // POLYNOMIAL_H_