polynomial.h 6.33 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<std::pair<Monomial, T> >();
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 29
  const std::list<std::pair<Monomial, T> >& terms() const {
    return terms_;
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
30 31
  }

SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
32 33
  void normalize();

34 35
  void SetToZero() {
    terms_.clear();
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
36
  }
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
37

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

42
  const Monomial& LM() const {
43 44
    assert(!terms_.back().second.IsZero());
    return terms_.back().first;
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
45 46
  }

SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
47
  const T& LC() const {
48
    return terms_.back().second;
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
49 50 51
  }

  void set_monic();
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
52 53

  unsigned degree() const {
54
    // TODO(pj): extend for non-degree orderings
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
55 56 57 58 59 60 61
    return LM().degree();
  }

#ifdef BILIN
  unsigned degree_bilin() const;
#endif

62 63 64
  bool topred(const Polynomial &g);
  bool istopred(const Polynomial &g) const;
  T coeff(const Monomial& m) const;
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
65

66 67 68
  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
69 70
  void operator*=(const T&);
  void operator/=(const T&);
SPAENLEHAUER Pierre-Jean's avatar
doc  
SPAENLEHAUER Pierre-Jean committed
71 72

  void print(std::ostream& out) const;
73
  inline bool operator<(const Polynomial<T>& P) const;
74

75
  friend std::ostream& operator<< <T>(std::ostream&, const Polynomial<T>&);
76 77 78 79 80 81 82 83 84 85 86 87 88 89
 private:
  // Adds coefficients of the two leading terms if they share the same
  // monomials. Useful for adding polynomials.
  void CombineLeadingMonomials();

  bool IsEmpty() const {
    return terms_.empty();
  }

  // Removes terms with coefficient 0
  void RemoveZeroes();

  std::list<Term> terms_;

SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
90
};
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
91 92 93

#ifdef BILIN
template<class T>
94
unsigned Polynomial<T>::degree_bilin() const {
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
95
  unsigned R = 0;
96
  for (typename std::list<std::pair<Monomial, T> >::const_iterator it = \
97
      terms_.begin(); it != terms_.end(); ++it)
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
98 99 100 101 102 103 104
    if (it->first.degree_bilin() > R)
      R = it->first.degree_bilin();
  return R;
}
#endif

template<class T>
105 106
void Polynomial<T>::RemoveZeroes() {
  for (auto it = terms_.begin(); it != terms_.end(); ++it) {
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
107
    if (it->second == T(0)) {
108
      it = terms_.erase(it);
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
109 110
      continue;
    }
111 112 113 114 115 116 117
  }
}

template<class T>
void Polynomial<T>::CombineLeadingMonomials() {
  typename std::list<std::pair<Monomial, T> >::iterator it2;
  for (auto it = terms_.begin(); it != terms_.end();) {
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
118 119
    it2 = it;
    ++it;
120
    if (it == terms_.end()) {
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
121 122 123
      break;
    } else if (it->first == it2->first) {
      it->second += it2->second;
124
      terms_.erase(it2);
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
125 126 127 128 129
    }
  }
}

template<class T>
130 131
T Polynomial<T>::coeff(const Monomial& m) const {
  for (typename std::list<std::pair<Monomial, T> >::const_iterator it = \
132
      terms_.begin(); it != terms_.end(); ++it)
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
133 134 135 136 137 138
    if (it->first == m)
      return it->second;
  return T(0);
}

template<class T>
139
void Polynomial<T>::set_monic() {
140 141 142
  RemoveZeroes();
  CombineLeadingMonomials();
  T leading_coeff = terms_.back().second;
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
143 144
  leading_coeff.inv();
  // TODO(pj):  test if 1
145
  for (auto& it : terms_)
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
146 147 148 149
    it.second *= leading_coeff;
}

template<class T>
150
void Polynomial<T>::normalize() {
151 152 153
  terms_.sort();
  RemoveZeroes();
  CombineLeadingMonomials();
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
154 155 156
}

template<class T>
157
bool Polynomial<T>::topred(const Polynomial &g) {
158
  if (g.IsEmpty() || IsEmpty())
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
159 160 161
    return false;
  if (!LM().isDivisibleBy(g.LM()))
    return false;
162
  Polynomial reductor = g*(LM()/g.LM());
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
163 164 165 166 167 168 169
  reductor.set_monic();
  reductor *= LC();
  operator-=(reductor);
  return true;
}

template<class T>
170
bool Polynomial<T>::istopred(const Polynomial<T> &g) const {
171
  if (g.IsEmpty() || IsEmpty()) return false;
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
172 173 174 175
  return LM().isDivisibleBy(g.LM());
}

template<class T>
176 177
Polynomial<T> Polynomial<T>::operator*(const Monomial &m) const {
  Polynomial R;
178
  R.terms_.insert(R.terms_.begin(), this->terms_.begin(), this->terms_.end());
179
  for (typename std::list<std::pair<Monomial, T> >::iterator it = \
180
      R.terms_.begin(); it != R.terms_.end(); ++it)
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
181 182 183 184 185
    (it->first) *= m;
  return R;
}

template<class T>
186
void Polynomial<T>::operator+=(const Polynomial<T> &p2) {
187 188 189
  terms_.merge(p2.terms_);
  RemoveZeroes();
  CombineLeadingMonomials();
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
190 191 192
}

template<class T>
193
void Polynomial<T>::operator*=(const T &e) {
194 195
  for (typename std::list<std::pair<Monomial, T> >::iterator it = terms_.begin(); \
      it != terms_.end(); ++it)
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
196 197 198 199
    it->second *= e;
}

template<class T>
200
void Polynomial<T>::operator/=(const T& e) {
201 202
  for (typename std::list<std::pair<Monomial, T> >::iterator it = terms_.begin(); \
      it != terms_.end(); ++it)
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
203 204 205 206
    it->second /= e;
}

template <class T>
207
void Polynomial<T>::operator-=(const Polynomial<T>& p2) {
208
  typename std::list<std::pair<Monomial, T> >::iterator it = terms_.begin();
209
  typename std::list<std::pair<Monomial, T> >::const_iterator it2 = \
210
       p2.terms_.begin();
211
  typename std::list<std::pair<Monomial, T> >::const_iterator it2end = \
212
       p2.terms_.end();
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
213

214 215 216 217
  while (it2 != p2.terms_.end()) {
    if (it == terms_.end()) {
      terms_.insert(it, it2, it2end);
      while (it != terms_.end()) {
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
218 219 220 221 222 223 224
        it->second.neg(); ++it;
      }
      break;
    } else if (it->first < it2->first) {
      ++it;
    } else if (it->first == it2->first) {
      if (it->second == it2->second) {
225
        it = terms_.erase(it);
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
226 227 228 229 230 231
        ++it2;
      } else {
        it->second -= it2->second;
        ++it2;
      }
    } else {
232
      it = terms_.insert(it, *it2);
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
233 234 235 236 237 238 239
      it->second.neg();
      ++it2;
    }
  }
}

template <class T>
240
void Polynomial<T>::print(std::ostream& o) const {
241
  for (auto& it : std::as_const(terms_)) {
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
242 243 244 245 246 247 248 249
    o << it.second << " ";
    it.first.print(o);
    o << std::endl;
  }
  o << ";";
}

template <class T>
250
std::ostream& operator<<(std::ostream& o, const Polynomial<T>& p) {
251
  if (p.IsEmpty()) {
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
252 253 254
    o << "0";
    return o;
  }
255 256
  typename std::list<std::pair<Monomial, T> >::const_iterator it = p.terms_.begin();
  if (it == p.terms_.end())
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
257 258 259 260
    return o;
  if (it->second != T(0))
    o << it->second << it->first;
  ++it;
261
  for (; it != p.terms_.end() ; ++it)
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
262 263 264 265 266 267
    if (it->second != T(0))
      o << "+" << it->second << it->first;
  return o;
}

template <class T>
268
inline bool Polynomial<T>::operator<(const Polynomial<T>& P) const {
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
269 270
  return LM() < P.LM();
}
271
}  // namespace tinygb
SPAENLEHAUER Pierre-Jean's avatar
SPAENLEHAUER Pierre-Jean committed
272

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