polynomial.h 5.94 KB
Newer Older
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
1 2
#ifndef POLYNOMIAL_H_
#define POLYNOMIAL_H_
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
3 4

#include <list>
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
5 6 7
#include <utility>
#include "./param.h"
#include "./monomial.h"
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
8

9 10
namespace tinygb {

SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
11
template<class T>
12
class Polynomial;
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
13 14

template<class T>
15
std::ostream& operator<<(std::ostream &o, const Polynomial<T> &p);
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
16 17

template<class T>
18
class Polynomial {
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
19
 public:
20 21
  typedef std::pair<Monomial, T> Term;

22
  Polynomial() {
23
      terms_ = std::list<Term>();
SPAENLEHAUER Pierre-Jean's avatar
doc  
SPAENLEHAUER Pierre-Jean committed
24
    }
25 26
  explicit Polynomial(const std::list<Term>& terms_)
    : terms_(terms_) {}
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
27

28
  const std::list<Term>& terms() const {
29
    return terms_;
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
30 31
  }

32 33 34 35
  // 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();
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
36

37 38
  void SetToZero() {
    terms_.clear();
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
39
  }
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
40

41
  void insert(Monomial m, T e) {
42
    terms_.push_back(std::pair<Monomial, T>(m, e));
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
43
  }
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
44

45 46 47
  // Assumes that the polynomial's list of terms is normalized
  const Monomial& LeadingMonomial() const {
    assert(!terms_.back().second.is_zero());
48
    return terms_.back().first;
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
49 50
  }

51
  const T& LeadingCoefficient() const {
52
    return terms_.back().second;
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
53 54
  }

55
  void DivideByLeadingCoefficient();
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
56

57 58 59
  unsigned Degree() const {
    // TODO(pj): Extend for non-degree orderings.
    return LeadingMonomial().Degree();
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
60 61 62
  }

#ifdef BILIN
63
  unsigned Degree_bilin() const;
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
64 65
#endif

66 67 68 69 70 71 72 73 74 75 76
  // 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;
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
77

78 79 80
  Polynomial<T> operator*(const Monomial&) const;
  void operator+=(const Polynomial<T>&);
  void operator-=(const Polynomial<T>&);
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
81 82
  void operator*=(const T&);
  void operator/=(const T&);
83
  void Print(std::ostream& out) const;
84
  inline bool operator<(const Polynomial<T>& P) const;
85

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

88 89 90 91
 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();
92 93 94 95
  bool IsEmpty() const {
    return terms_.empty();
  }
  // Removes terms with coefficient 0
96
  void RemoveZeroTerms();
97 98 99

  std::list<Term> terms_;

SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
100
};
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
101 102 103

#ifdef BILIN
template<class T>
104
unsigned Polynomial<T>::Degree_bilin() const {
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
105
  unsigned R = 0;
106 107 108
  for (auto& it = std::as_const(terms_))
    if (it->first.Degree_bilin() > R)
      R = it->first.Degree_bilin();
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
109 110 111 112 113
  return R;
}
#endif

template<class T>
114
void Polynomial<T>::RemoveZeroTerms() {
115
  for (auto it = terms_.begin(); it != terms_.end(); ++it) {
116
    while (it != terms_.end() && it->second == T(0))
117 118 119 120 121
      it = terms_.erase(it);
  }
}

template<class T>
122
void Polynomial<T>::CombineMonomials() {
123
  typename std::list<std::pair<Monomial, T> >::iterator it2;
124 125 126 127
  for (auto it1 = terms_.begin(); it1 != terms_.end();) {
    it2 = it1;
    ++it1;
    if (it1 == terms_.end()) {
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
128
      break;
129 130
    } else if (it1->first == it2->first) {
      it1->second += it2->second;
131
      terms_.erase(it2);
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
132 133 134 135 136
    }
  }
}

template<class T>
137 138 139 140 141
void Polynomial<T>::DivideByLeadingCoefficient() {
  T leading_coefficient = terms_.back().second;
  assert(leading_coefficient != T(0));
  leading_coefficient.inv();
  if (leading_coefficient == T(1)) return;
142
  for (auto& it : terms_)
143
    it.second *= leading_coefficient;
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
144 145 146
}

template<class T>
147
void Polynomial<T>::Normalize() {
148
  terms_.sort();
149 150
  CombineMonomials();
  RemoveZeroTerms();
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
151 152 153
}

template<class T>
154 155 156 157 158 159 160
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();
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
161 162 163 164 165
  operator-=(reductor);
  return true;
}

template<class T>
166 167 168 169
bool Polynomial<T>::IsTopReducibleBy(const Polynomial<T> &g) const {
  assert(!g.IsEmpty());
  if (IsEmpty()) return false;
  return LeadingMonomial().IsDivisibleBy(g.LeadingMonomial());
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
170 171 172
}

template<class T>
173 174
Polynomial<T> Polynomial<T>::operator*(const Monomial &m) const {
  Polynomial R;
175
  R.terms_.insert(R.terms_.begin(), this->terms_.begin(), this->terms_.end());
176 177
  for (auto& it : R.terms_)
    it.first *= m;
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
178 179 180 181
  return R;
}

template<class T>
182 183 184 185 186
void Polynomial<T>::operator+=(const Polynomial<T> &P) {
  std::list<Term> tmp_terms = P.terms_;
  terms_.merge(tmp_terms);
  CombineMonomials();
  RemoveZeroTerms();
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
187 188 189
}

template<class T>
190
void Polynomial<T>::operator*=(const T &e) {
191 192
  for (auto& it : terms_)
    it.second *= e;
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
193 194 195
}

template<class T>
196
void Polynomial<T>::operator/=(const T& e) {
197 198
  for (auto& it : terms_)
    it.second /= e;
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
199 200 201
}

template <class T>
202 203 204 205 206 207 208
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();
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
209 210 211
}

template <class T>
212
void Polynomial<T>::Print(std::ostream& o) const {
213
  for (auto& it : std::as_const(terms_)) {
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
214
    o << it.second << " ";
215
    it.first.Print(o);
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
216 217 218 219 220 221
    o << std::endl;
  }
  o << ";";
}

template <class T>
222 223
std::ostream& operator<<(std::ostream& o, const Polynomial<T>& P) {
  if (P.IsEmpty()) {
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
224 225 226
    o << "0";
    return o;
  }
227 228
  auto it = std::as_const(P.terms_).begin();
  if (it == P.terms_.end())
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
229 230 231 232
    return o;
  if (it->second != T(0))
    o << it->second << it->first;
  ++it;
233
  for (; it != P.terms_.end() ; ++it)
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
234 235 236 237 238 239
    if (it->second != T(0))
      o << "+" << it->second << it->first;
  return o;
}

template <class T>
240
inline bool Polynomial<T>::operator<(const Polynomial<T>& P) const {
241
  return LeadingMonomial() < P.LeadingMonomial();
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
242
}
243
}  // namespace tinygb
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
244

SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
245
#endif  // POLYNOMIAL_H_