Loading...
Searching...
No Matches
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>
13#include <memory>
14
15#ifdef HAVE_SYMENGINE_FLINT
16#include <symengine/flint_wrapper.h>
17using fz_t = SymEngine::fmpz_wrapper;
18using 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
25namespace SymEngine
26{
27// misc methods
28
29#if SYMENGINE_INTEGER_CLASS == SYMENGINE_GMPXX \
30 || SYMENGINE_INTEGER_CLASS == SYMENGINE_GMP
31#ifdef HAVE_SYMENGINE_FLINT
32inline 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}
38inline 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
47inline 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}
53inline 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
63inline 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}
69inline 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
80inline integer_class to_mp_class(const piranha::integer &x)
81{
82 return integer_class(x.get_mpz_view());
83}
84inline 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
92inline integer_class to_mp_class(const integer_class &i)
93{
94 return i;
95}
96
97inline rational_class to_mp_class(const rational_class &i)
98{
99 return i;
100}
101
102// dict wrapper
103template <typename Key, typename Value, typename Wrapper>
105{
106public:
108 typedef Key key_type;
109
110public:
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
351umap_basic_num _find_gens_poly(const RCP<const Basic> &x);
352
353template <typename Container, typename Poly>
354class UPolyBase : public Basic
355{
356private:
357 RCP<const Basic> var_;
358 Container poly_;
359
360public:
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
406template <typename Cont, typename Poly>
407class UExprPolyBase : public UPolyBase<Cont, Poly>
408{
409public:
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.
457template <typename Container, typename Poly, typename Cf>
458class UNonExprPoly : public UPolyBase<Container, Poly>
459{
460public:
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
500template <typename Container, typename Poly>
501class UIntPolyBase : public UNonExprPoly<Container, Poly, integer_class>
502{
503public:
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
541template <typename Container, typename Poly>
542class URatPolyBase : public UNonExprPoly<Container, Poly, rational_class>
543{
544public:
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) {
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
583template <typename T, typename Int>
585{
586protected:
587 RCP<const T> ptr_;
588 long i_;
589
590public:
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
617template <typename T, typename Int>
619{
620public:
621 ContainerForIter(RCP<const T> ptr, long 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
643template <typename T, typename Int>
645{
646public:
647 ContainerRevIter(RCP<const T> ptr, long 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
665template <typename P>
667 static const bool value
669};
670
671template <typename Poly>
672RCP<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
682template <typename Poly>
683RCP<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
690template <typename Poly>
691RCP<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
701template <typename Poly>
702RCP<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
712template <typename Poly>
713RCP<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:97
static RCP< const Basic > from_dict(const RCP< const Number > &coef, map_basic_basic &&d)
Create a Mul from a dict.
Definition: mul.cpp:115
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
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:352
RCP< const Basic > add(const RCP< const Basic > &a, const RCP< const Basic > &b)
Adds two objects (safely).
Definition: add.cpp:425
std::enable_if< std::is_integral< T >::value, RCP< constInteger > >::type integer(T i)
Definition: integer.h:200
T push_back(T... args)
T rbegin(T... args)
T size(T... args)