upolybase.h
1 
5 #ifndef SYMENGINE_UINT_BASE_H
6 #define SYMENGINE_UINT_BASE_H
7 
8 #include <symengine/basic.h>
9 #include <symengine/pow.h>
10 #include <symengine/add.h>
11 #include <symengine/rational.h>
12 #include <symengine/expression.h>
13 #include <memory>
14 
15 #ifdef HAVE_SYMENGINE_FLINT
16 #include <symengine/flint_wrapper.h>
17 using fz_t = SymEngine::fmpz_wrapper;
18 using fq_t = SymEngine::fmpq_wrapper;
19 #endif
20 #ifdef HAVE_SYMENGINE_PIRANHA
21 #include <piranha/mp_integer.hpp>
22 #include <piranha/mp_rational.hpp>
23 #endif
24 
25 namespace SymEngine
26 {
27 // misc methods
28 
29 #if SYMENGINE_INTEGER_CLASS == SYMENGINE_GMPXX \
30  || SYMENGINE_INTEGER_CLASS == SYMENGINE_GMP
31 #ifdef HAVE_SYMENGINE_FLINT
32 inline integer_class to_mp_class(const fz_t &i)
33 {
34  integer_class x;
35  fmpz_get_mpz(x.get_mpz_t(), i.get_fmpz_t());
36  return x;
37 }
38 inline rational_class to_mp_class(const fq_t &i)
39 {
40  rational_class x;
41  fmpq_get_mpq(x.get_mpq_t(), i.get_fmpq_t());
42  return x;
43 }
44 #endif
45 
46 #ifdef HAVE_SYMENGINE_PIRANHA
47 inline integer_class to_mp_class(const piranha::integer &i)
48 {
49  integer_class x;
50  mpz_set(x.get_mpz_t(), i.get_mpz_view());
51  return x;
52 }
53 inline rational_class to_mp_class(const piranha::rational &i)
54 {
55  rational_class x;
56  mpq_set(x.get_mpq_t(), i.get_mpq_view());
57  return x;
58 }
59 #endif
60 
61 #elif SYMENGINE_INTEGER_CLASS == SYMENGINE_PIRANHA
62 #ifdef HAVE_SYMENGINE_FLINT
63 inline integer_class to_mp_class(const fz_t &i)
64 {
65  integer_class x;
66  fmpz_get_mpz(get_mpz_t(x), i.get_fmpz_t());
67  return x;
68 }
69 inline rational_class to_mp_class(const fq_t &i)
70 {
71  rational_class s;
72  fmpz_get_mpz(get_mpz_t(s._num()), i.get_num().get_fmpz_t());
73  fmpz_get_mpz(get_mpz_t(s._den()), i.get_den().get_fmpz_t());
74  return s;
75 }
76 #endif
77 
78 #elif SYMENGINE_INTEGER_CLASS == SYMENGINE_FLINT
79 #ifdef HAVE_SYMENGINE_PIRANHA
80 inline integer_class to_mp_class(const piranha::integer &x)
81 {
82  return integer_class(x.get_mpz_view());
83 }
84 inline rational_class to_mp_class(const piranha::rational &x)
85 {
86  return rational_class(x.get_mpq_view());
87 }
88 #endif
89 
90 #endif
91 
92 inline integer_class to_mp_class(const integer_class &i)
93 {
94  return i;
95 }
96 
97 inline rational_class to_mp_class(const rational_class &i)
98 {
99  return i;
100 }
101 
102 // dict wrapper
103 template <typename Key, typename Value, typename Wrapper>
105 {
106 public:
107  std::map<Key, Value> dict_;
108  typedef Key key_type;
109 
110 public:
111  ODictWrapper() SYMENGINE_NOEXCEPT {}
112  ~ODictWrapper() SYMENGINE_NOEXCEPT {}
113 
114  ODictWrapper(const int &i)
115  {
116  if (i != 0)
117  dict_ = {{0, Value(i)}};
118  }
119 
121  {
122  for (auto &iter : p) {
123  if (iter.second != Value(0))
124  dict_[iter.first] = iter.second;
125  }
126  }
127 
129  {
130  for (auto &iter : p) {
131  if (iter.second != Value(0)) {
132  auto erase = iter;
133  iter++;
134  p.erase(erase);
135  } else {
136  iter++;
137  }
138  }
139  dict_ = p;
140  }
141 
142  ODictWrapper(const Value &p)
143  {
144  if (p != Value(0))
145  dict_[0] = p;
146  }
147 
149  {
150  dict_[1] = Value(1);
151  }
152 
153  static Wrapper from_vec(const std::vector<Value> &v)
154  {
155  Wrapper x;
156  x.dict_ = {};
157  for (unsigned int i = 0; i < v.size(); i++) {
158  if (v[i] != Value(0)) {
159  x.dict_[i] = v[i];
160  }
161  }
162  return x;
163  }
164 
165  Wrapper &operator=(Wrapper &&other) SYMENGINE_NOEXCEPT
166  {
167  if (this != &other)
168  dict_ = std::move(other.dict_);
169  return static_cast<Wrapper &>(*this);
170  }
171 
172  friend Wrapper operator+(const Wrapper &a, const Wrapper &b)
173  {
174  Wrapper c = a;
175  c += b;
176  return c;
177  }
178 
179  Wrapper &operator+=(const Wrapper &other)
180  {
181  for (auto &iter : other.dict_) {
182  auto t = dict_.lower_bound(iter.first);
183  if (t != dict_.end() and t->first == iter.first) {
184  t->second += iter.second;
185  if (t->second == 0) {
186  dict_.erase(t);
187  }
188  } else {
189  dict_.insert(t, {iter.first, iter.second});
190  }
191  }
192  return static_cast<Wrapper &>(*this);
193  }
194 
195  friend Wrapper operator-(const Wrapper &a, const Wrapper &b)
196  {
197  Wrapper c = a;
198  c -= b;
199  return c;
200  }
201 
202  Wrapper operator-() const
203  {
204  ODictWrapper c = *this;
205  for (auto &iter : c.dict_)
206  iter.second *= -1;
207  return static_cast<Wrapper &>(c);
208  }
209 
210  Wrapper &operator-=(const Wrapper &other)
211  {
212  for (auto &iter : other.dict_) {
213  auto t = dict_.lower_bound(iter.first);
214  if (t != dict_.end() and t->first == iter.first) {
215  t->second -= iter.second;
216  if (t->second == 0) {
217  dict_.erase(t);
218  }
219  } else {
220  dict_.insert(t, {iter.first, -iter.second});
221  }
222  }
223  return static_cast<Wrapper &>(*this);
224  }
225 
226  static Wrapper mul(const Wrapper &a, const Wrapper &b)
227  {
228  if (a.get_dict().empty())
229  return a;
230  if (b.get_dict().empty())
231  return b;
232 
233  Wrapper p;
234  for (const auto &i1 : a.dict_)
235  for (const auto &i2 : b.dict_)
236  p.dict_[i1.first + i2.first] += i1.second * i2.second;
237 
238  for (auto it = p.dict_.cbegin(); it != p.dict_.cend();) {
239  if (it->second == 0) {
240  p.dict_.erase(it++);
241  } else {
242  ++it;
243  }
244  }
245  return p;
246  }
247 
248  static Wrapper pow(const Wrapper &a, unsigned int p)
249  {
250  Wrapper tmp = a, res(1);
251 
252  while (p != 1) {
253  if (p % 2 == 0) {
254  tmp = tmp * tmp;
255  } else {
256  res = res * tmp;
257  tmp = tmp * tmp;
258  }
259  p >>= 1;
260  }
261 
262  return (res * tmp);
263  }
264 
265  template <typename FromPoly>
266  static Wrapper from_poly(const FromPoly &p)
267  {
268  Wrapper t;
269  for (auto it = p.begin(); it != p.end(); ++it)
270  t.dict_[it->first] = it->second;
271  return t;
272  }
273 
274  friend Wrapper operator*(const Wrapper &a, const Wrapper &b)
275  {
276  return Wrapper::mul(a, b);
277  }
278 
279  Wrapper &operator*=(const Wrapper &other)
280  {
281  if (dict_.empty())
282  return static_cast<Wrapper &>(*this);
283 
284  if (other.dict_.empty()) {
285  dict_.clear();
286  return static_cast<Wrapper &>(*this);
287  }
288 
289  // ! other is a just constant term
290  if (other.dict_.size() == 1
291  and other.dict_.find(0) != other.dict_.end()) {
292  auto t = other.dict_.begin();
293  for (auto &i1 : dict_)
294  i1.second *= t->second;
295  return static_cast<Wrapper &>(*this);
296  }
297 
298  Wrapper res = Wrapper::mul(static_cast<Wrapper &>(*this), other);
299  res.dict_.swap(this->dict_);
300  return static_cast<Wrapper &>(*this);
301  }
302 
303  friend bool operator==(const Wrapper &a, const Wrapper &b)
304  {
305  return a.dict_ == b.dict_;
306  }
307 
308  bool operator!=(const Wrapper &other) const
309  {
310  return not(*this == other);
311  }
312 
313  const std::map<Key, Value> &get_dict() const
314  {
315  return dict_;
316  }
317 
318  size_t size() const
319  {
320  return dict_.size();
321  }
322 
323  bool empty() const
324  {
325  return dict_.empty();
326  }
327 
328  Key degree() const
329  {
330  if (dict_.empty())
331  return Key(0);
332  return dict_.rbegin()->first;
333  }
334 
335  Value get_coeff(Key x) const
336  {
337  auto ite = dict_.find(x);
338  if (ite != dict_.end())
339  return ite->second;
340  return Value(0);
341  }
342 
343  Value get_lc() const
344  {
345  if (dict_.empty())
346  return Value(0);
347  return dict_.rbegin()->second;
348  }
349 };
350 
351 umap_basic_num _find_gens_poly(const RCP<const Basic> &x);
352 
353 template <typename Container, typename Poly>
354 class UPolyBase : public Basic
355 {
356 private:
357  RCP<const Basic> var_;
358  Container poly_;
359 
360 public:
361  UPolyBase(const RCP<const Basic> &var, Container &&container)
362  : var_{var}, poly_{container}
363  {
364  }
365 
366  typedef Container container_type;
367 
369  virtual int compare(const Basic &o) const = 0;
370  virtual hash_t __hash__() const = 0;
371 
372  // return `degree` + 1. `0` returned for zero poly.
373  virtual int size() const = 0;
374 
376  inline bool __eq__(const Basic &o) const
377  {
378  if (is_a<Poly>(o))
379  return eq(*var_, *(down_cast<const Poly &>(o).var_))
380  and poly_ == down_cast<const Poly &>(o).poly_;
381  return false;
382  }
383 
384  inline const RCP<const Basic> &get_var() const
385  {
386  return var_;
387  }
388 
389  inline const Container &get_poly() const
390  {
391  return poly_;
392  }
393 
394  inline vec_basic get_args() const
395  {
396  return {};
397  }
398 
399  static RCP<const Poly> from_container(const RCP<const Basic> &var,
400  Container &&d)
401  {
402  return make_rcp<const Poly>(var, std::move(d));
403  }
404 };
405 
406 template <typename Cont, typename Poly>
407 class UExprPolyBase : public UPolyBase<Cont, Poly>
408 {
409 public:
410  typedef Expression coef_type;
411 
412  UExprPolyBase(const RCP<const Basic> &var, Cont &&container)
413  : UPolyBase<Cont, Poly>(var, std::move(container))
414  {
415  }
416 
417  inline int get_degree() const
418  {
419  return this->get_poly().degree();
420  }
421 
422  static RCP<const Poly> from_dict(const RCP<const Basic> &var,
424  {
425  return Poly::from_container(
426  var, Poly::container_from_dict(var, std::move(d)));
427  }
428 
429  RCP<const Basic> as_symbolic() const
430  {
431  auto it = (down_cast<const Poly &>(*this)).begin();
432  auto end = (down_cast<const Poly &>(*this)).end();
433 
434  vec_basic args;
435  for (; it != end; ++it) {
436  if (it->first == 0)
437  args.push_back(it->second.get_basic());
438  else if (it->first == 1) {
439  if (it->second == Expression(1))
440  args.push_back(this->get_var());
441  else
442  args.push_back(
443  mul(it->second.get_basic(), this->get_var()));
444  } else if (it->second == 1)
445  args.push_back(pow(this->get_var(), integer(it->first)));
446  else
447  args.push_back(mul(it->second.get_basic(),
448  pow(this->get_var(), integer(it->first))));
449  }
450  if (this->get_poly().empty())
451  args.push_back(zero);
452  return SymEngine::add(args);
453  }
454 };
455 // super class for all non-expr polys, all methods which are
456 // common for all non-expr polys go here eg. degree, eval etc.
457 template <typename Container, typename Poly, typename Cf>
458 class UNonExprPoly : public UPolyBase<Container, Poly>
459 {
460 public:
461  typedef Cf coef_type;
462 
463  UNonExprPoly(const RCP<const Basic> &var, Container &&container)
464  : UPolyBase<Container, Poly>(var, std::move(container))
465  {
466  }
467 
468  // return coefficient of degree 'i'
469  virtual Cf get_coeff(unsigned int i) const = 0;
470  // return value of poly when ealudated at `x`
471  virtual Cf eval(const Cf &x) const = 0;
472 
473  std::vector<Cf> multieval(const std::vector<Cf> &v) const
474  {
475  // this is not the optimal algorithm
476  std::vector<Cf> res(v.size());
477  for (unsigned int i = 0; i < v.size(); ++i)
478  res[i] = eval(v[i]);
479  return res;
480  }
481 
482  inline int get_degree() const
483  {
484  return numeric_cast<int>(this->get_poly().degree());
485  }
486 
487  Cf get_lc() const
488  {
489  return get_coeff(get_degree());
490  }
491 
492  static RCP<const Poly> from_dict(const RCP<const Basic> &var,
494  {
495  return Poly::from_container(
496  var, Poly::container_from_dict(var, std::move(d)));
497  }
498 };
499 
500 template <typename Container, typename Poly>
501 class UIntPolyBase : public UNonExprPoly<Container, Poly, integer_class>
502 {
503 public:
504  UIntPolyBase(const RCP<const Basic> &var, Container &&container)
506  std::move(container))
507  {
508  }
509 
510  RCP<const Basic> as_symbolic() const
511  {
512  auto it = (down_cast<const Poly &>(*this)).begin();
513  auto end = (down_cast<const Poly &>(*this)).end();
514 
515  vec_basic args;
516  for (; it != end; ++it) {
517  integer_class m = it->second;
518 
519  if (it->first == 0) {
520  args.push_back(integer(m));
521  } else if (it->first == 1) {
522  if (m == 1) {
523  args.push_back(this->get_var());
524  } else {
525  args.push_back(
526  Mul::from_dict(integer(m), {{this->get_var(), one}}));
527  }
528  } else {
529  if (m == 1) {
530  args.push_back(pow(this->get_var(), integer(it->first)));
531  } else {
533  integer(m), {{this->get_var(), integer(it->first)}}));
534  }
535  }
536  }
537  return SymEngine::add(args);
538  }
539 };
540 
541 template <typename Container, typename Poly>
542 class URatPolyBase : public UNonExprPoly<Container, Poly, rational_class>
543 {
544 public:
545  URatPolyBase(const RCP<const Basic> &var, Container &&container)
547  std::move(container))
548  {
549  }
550 
551  RCP<const Basic> as_symbolic() const
552  {
553  auto it = (down_cast<const Poly &>(*this)).begin();
554  auto end = (down_cast<const Poly &>(*this)).end();
555 
556  vec_basic args;
557  for (; it != end; ++it) {
558  rational_class m = it->second;
559 
560  if (it->first == 0) {
561  args.push_back(Rational::from_mpq(m));
562  } else if (it->first == 1) {
563  if (m == 1) {
564  args.push_back(this->get_var());
565  } else {
567  {{this->get_var(), one}}));
568  }
569  } else {
570  if (m == 1) {
571  args.push_back(pow(this->get_var(), integer(it->first)));
572  } else {
575  {{this->get_var(), integer(it->first)}}));
576  }
577  }
578  }
579  return SymEngine::add(args);
580  }
581 };
582 
583 template <typename T, typename Int>
585 {
586 protected:
587  RCP<const T> ptr_;
588  long i_;
589 
590 public:
591  ContainerBaseIter(RCP<const T> ptr, long x) : ptr_{ptr}, i_{x} {}
592 
593  friend bool operator==(const ContainerBaseIter &lhs,
594  const ContainerBaseIter &rhs)
595  {
596  return (lhs.ptr_ == rhs.ptr_) and (lhs.i_ == rhs.i_);
597  }
598 
599  bool operator!=(const ContainerBaseIter &rhs)
600  {
601  return not(*this == rhs);
602  }
603 
604  std::pair<long, Int> operator*()
605  {
606  return std::make_pair(i_, ptr_->get_coeff_ref(i_));
607  }
608 
610  {
611  return std::make_shared<std::pair<unsigned, Int>>(
612  numeric_cast<unsigned>(i_),
613  ptr_->get_coeff_ref(numeric_cast<unsigned>(i_)));
614  }
615 };
616 
617 template <typename T, typename Int>
618 class ContainerForIter : public ContainerBaseIter<T, Int>
619 {
620 public:
621  ContainerForIter(RCP<const T> ptr, long x)
622  : ContainerBaseIter<T, Int>(ptr, x)
623  {
624  if (this->ptr_->get_coeff_ref(numeric_cast<unsigned>(this->i_)) == 0
625  and this->i_ < this->ptr_->size()) {
626  ++(*this);
627  }
628  }
629 
630  ContainerForIter operator++()
631  {
632  this->i_++;
633  while (this->i_ < this->ptr_->size()) {
634  if (this->ptr_->get_coeff_ref(numeric_cast<unsigned>(this->i_))
635  != 0)
636  break;
637  this->i_++;
638  }
639  return *this;
640  }
641 };
642 
643 template <typename T, typename Int>
644 class ContainerRevIter : public ContainerBaseIter<T, Int>
645 {
646 public:
647  ContainerRevIter(RCP<const T> ptr, long x)
648  : ContainerBaseIter<T, Int>(ptr, x)
649  {
650  }
651 
652  ContainerRevIter operator++()
653  {
654  this->i_--;
655  while (this->i_ >= 0) {
656  if (this->ptr_->get_coeff_ref(numeric_cast<unsigned>(this->i_))
657  != 0)
658  break;
659  this->i_--;
660  }
661  return *this;
662  }
663 };
664 
665 template <typename P>
666 struct is_a_UPoly {
667  static const bool value
669 };
670 
671 template <typename Poly>
672 RCP<const Poly> add_upoly(const Poly &a, const Poly &b)
673 {
674  if (!(a.get_var()->__eq__(*b.get_var())))
675  throw SymEngineException("Error: variables must agree.");
676 
677  auto dict = a.get_poly();
678  dict += b.get_poly();
679  return Poly::from_container(a.get_var(), std::move(dict));
680 }
681 
682 template <typename Poly>
683 RCP<const Poly> neg_upoly(const Poly &a)
684 {
685  auto dict = a.get_poly();
686  dict = -dict;
687  return Poly::from_container(a.get_var(), std::move(dict));
688 }
689 
690 template <typename Poly>
691 RCP<const Poly> sub_upoly(const Poly &a, const Poly &b)
692 {
693  if (!(a.get_var()->__eq__(*b.get_var())))
694  throw SymEngineException("Error: variables must agree.");
695 
696  auto dict = a.get_poly();
697  dict -= b.get_poly();
698  return Poly::from_container(a.get_var(), std::move(dict));
699 }
700 
701 template <typename Poly>
702 RCP<const Poly> mul_upoly(const Poly &a, const Poly &b)
703 {
704  if (!(a.get_var()->__eq__(*b.get_var())))
705  throw SymEngineException("Error: variables must agree.");
706 
707  auto dict = a.get_poly();
708  dict *= b.get_poly();
709  return Poly::from_container(a.get_var(), std::move(dict));
710 }
711 
712 template <typename Poly>
713 RCP<const Poly> quo_upoly(const Poly &a, const Poly &b)
714 {
715  if (!(a.get_var()->__eq__(*b.get_var())))
716  throw SymEngineException("Error: variables must agree.");
717 
718  auto dict = a.get_poly();
719  dict /= b.get_poly();
720  return Poly::from_dict(a.get_var(), std::move(dict));
721 }
722 } // namespace SymEngine
723 
724 #endif // SYMENGINE_UINT_BASE_H
Classes and functions relating to the binary operation of addition.
The base class for SymEngine.
The lowest unit of symbolic representation.
Definition: basic.h:95
static RCP< const Basic > from_dict(const RCP< const Number > &coef, map_basic_basic &&d)
Create a Mul from a dict.
Definition: mul.cpp:116
static RCP< const Number > from_mpq(const rational_class &i)
Definition: rational.cpp:23
virtual hash_t __hash__() const =0
bool __eq__(const Basic &o) const
Definition: upolybase.h:376
virtual int compare(const Basic &o) const =0
vec_basic get_args() const
Returns the list of arguments.
Definition: upolybase.h:394
T clear(T... args)
T empty(T... args)
T end(T... args)
T erase(T... args)
T find(T... args)
T insert(T... args)
T lower_bound(T... args)
T make_pair(T... args)
T move(T... args)
Main namespace for SymEngine package.
Definition: add.cpp:19
std::enable_if< std::is_integral< T >::value, RCP< const Integer > >::type integer(T i)
Definition: integer.h:200
bool eq(const Basic &a, const Basic &b)
Checks equality for a and b
Definition: basic-inl.h:21
RCP< const Basic > mul(const RCP< const Basic > &a, const RCP< const Basic > &b)
Multiplication.
Definition: mul.cpp:347
RCP< const Basic > add(const RCP< const Basic > &a, const RCP< const Basic > &b)
Adds two objects (safely).
Definition: add.cpp:425
RCP< const Number > rational(long n, long d)
convenience creator from two longs
Definition: rational.h:328
T push_back(T... args)
T rbegin(T... args)
T size(T... args)