uintpoly_piranha.h
Go to the documentation of this file.
1 
5 #ifndef SYMENGINE_UINTPOLY_PIRANHA_H
6 #define SYMENGINE_UINTPOLY_PIRANHA_H
7 
8 #include <symengine/polys/upolybase.h>
9 #include <symengine/expression.h>
10 #include <symengine/dict.h>
11 #include <memory>
12 
13 #ifdef HAVE_SYMENGINE_PIRANHA
14 #include <piranha/monomial.hpp>
15 #include <piranha/polynomial.hpp>
16 #include <piranha/mp_rational.hpp>
17 #include <piranha/mp_integer.hpp>
18 #include <piranha/math.hpp>
19 #include <piranha/type_traits.hpp>
20 
21 #if SYMENGINE_INTEGER_CLASS != SYMENGINE_PIRANHA
22 namespace piranha
23 {
24 
25 // overloading pow for pirahna::math::evaluate
26 namespace math
27 {
28 using namespace SymEngine;
29 template <typename U, typename V>
30 struct pow_impl<V, U,
31  enable_if_t<std::is_integral<U>::value
32  and (std::is_same<V, integer_class>::value
33  or std::is_same<V, rational_class>::value)>> {
34  template <typename T2>
35  V operator()(const V &r, const T2 &x) const
36  {
37  V res;
38  mp_pow_ui(res, r, x);
39  return res;
40  }
41 };
42 
43 template <>
44 struct gcd_impl<SymEngine::integer_class, SymEngine::integer_class> {
45  SymEngine::integer_class operator()(const SymEngine::integer_class &r,
46  const SymEngine::integer_class &x) const
47  {
48  SymEngine::integer_class res;
49  mp_gcd(res, r, x);
50  return res;
51  }
52 };
53 
54 template <>
55 struct divexact_impl<SymEngine::integer_class> {
56  void operator()(SymEngine::integer_class &r,
57  const SymEngine::integer_class &x,
58  const SymEngine::integer_class &y) const
59  {
60  SymEngine::integer_class rem;
61  mp_tdiv_qr(r, rem, x, y);
62  if (rem != SymEngine::integer_class(0)) {
63  piranha_throw(inexact_division);
64  }
65  }
66 };
67 template <>
68 struct divexact_impl<SymEngine::rational_class> {
69  void operator()(SymEngine::rational_class &r,
70  const SymEngine::rational_class &x,
71  const SymEngine::rational_class &y) const
72  {
73  r = x / y;
74  }
75 };
76 } // namespace math
77 
78 template <>
79 struct has_exact_ring_operations<SymEngine::integer_class> {
80  static const bool value = true;
81 };
82 template <>
83 struct has_exact_ring_operations<SymEngine::rational_class> {
84  static const bool value = true;
85 };
86 } // namespace piranha
87 #endif
88 
89 // need definition for piranha::rational too
90 namespace piranha
91 {
92 namespace math
93 {
94 template <>
95 struct gcd_impl<SymEngine::rational_class, SymEngine::rational_class> {
96  SymEngine::rational_class
97  operator()(const SymEngine::rational_class &r,
98  const SymEngine::rational_class &x) const
99  {
100  return SymEngine::rational_class(1);
101  }
102 };
103 } // namespace math
104 } // namespace piranha
105 
106 namespace SymEngine
107 {
108 using pmonomial = piranha::monomial<unsigned int>;
109 using pintpoly = piranha::polynomial<integer_class, pmonomial>;
110 using pratpoly = piranha::polynomial<rational_class, pmonomial>;
111 
112 template <typename Cf, typename Container>
113 class PiranhaForIter
114 {
115 public:
116  typedef decltype(std::declval<Container &>()._container().begin()) ptr_type;
117 
118 private:
119  ptr_type ptr_;
120 
121 public:
122  PiranhaForIter(ptr_type ptr) : ptr_{ptr} {}
123 
124  bool operator==(const PiranhaForIter &rhs)
125  {
126  return (ptr_ == rhs.ptr_);
127  }
128 
129  bool operator!=(const PiranhaForIter &rhs)
130  {
131  return (ptr_ != rhs.ptr_);
132  }
133 
134  PiranhaForIter operator++()
135  {
136  ptr_++;
137  return *this;
138  }
139 
141  {
142  return std::make_pair(*(ptr_->m_key.begin()), ptr_->m_cf);
143  }
144 
146  {
147  return std::make_shared<std::pair<unsigned int, const Cf &>>(
148  *(ptr_->m_key.begin()), ptr_->m_cf);
149  }
150 };
151 
152 template <typename Container, template <typename X, typename Y> class BaseType,
153  typename Poly>
154 class UPiranhaPoly : public BaseType<Container, Poly>
155 {
156 public:
157  using Cf = typename BaseType<Container, Poly>::coef_type;
158  using term = typename Container::term_type;
159 
160  UPiranhaPoly(const RCP<const Basic> &var, Container &&dict)
161  : BaseType<Container, Poly>(var, std::move(dict))
162  {
163  }
164 
165  int compare(const Basic &o) const
166  {
167  SYMENGINE_ASSERT(is_a<Poly>(o))
168  const Poly &s = down_cast<const Poly &>(o);
169  int cmp = this->get_var()->compare(*s.get_var());
170  if (cmp != 0)
171  return cmp;
172  if (this->get_poly() == s.get_poly())
173  return 0;
174  return (this->get_poly().hash() < s.get_poly().hash()) ? -1 : 1;
175  }
176 
177  static Container container_from_dict(const RCP<const Basic> &var,
179  {
180  Container p;
181  piranha::symbol_set ss({{piranha::symbol(detail::poly_print(var))}});
182  p.set_symbol_set(ss);
183  for (auto &it : d)
184  if (it.second != 0)
185  p.insert(term(it.second, pmonomial{it.first}));
186 
187  return p;
188  }
189 
190  static RCP<const Poly> from_vec(const RCP<const Basic> &var,
191  const std::vector<Cf> &v)
192  {
193  Container p;
194  piranha::symbol_set ss({{piranha::symbol(detail::poly_print(var))}});
195  p.set_symbol_set(ss);
196  for (unsigned int i = 0; i < v.size(); i++) {
197  if (v[i] != 0) {
198  p.insert(term(v[i], pmonomial{i}));
199  }
200  }
201  return make_rcp<const Poly>(var, std::move(p));
202  }
203 
204  template <typename FromPoly>
205  static enable_if_t<is_a_UPoly<FromPoly>::value, RCP<const Poly>>
206  from_poly(const FromPoly &f)
207  {
208  Container p;
209  piranha::symbol_set ss(
210  {{piranha::symbol(detail::poly_print(f.get_var()))}});
211  p.set_symbol_set(ss);
212  for (auto it = f.begin(); it != f.end(); ++it)
213  p.insert(term(it->second, pmonomial{it->first}));
214  return Poly::from_container(f.get_var(), std::move(p));
215  }
216 
217  Cf eval(const Cf &x) const
218  {
220  = {{detail::poly_print(this->get_var()), x}};
221  return piranha::math::evaluate<Cf, Container>(this->get_poly(), t);
222  }
223 
224  Cf get_coeff(unsigned int x) const
225  {
226  return this->get_poly().find_cf(pmonomial{x});
227  }
228 
229  const Cf &get_coeff_ref(unsigned int x) const
230  {
231  static Cf pzero(0);
232 
233  term temp = term(0, pmonomial{x});
234  auto it = this->get_poly()._container().find(temp);
235  if (it == this->get_poly()._container().end())
236  return pzero;
237  return it->m_cf;
238  }
239 
240  int size() const
241  {
242  if (this->get_poly().size() == 0)
243  return 0;
244  return this->get_degree() + 1;
245  }
246 
247  // begin() and end() are unordered
248  // obegin() and oend() are ordered, from highest degree to lowest
249  typedef PiranhaForIter<Cf, Container> iterator;
250  typedef ContainerRevIter<Poly, const Cf &> r_iterator;
251  iterator begin() const
252  {
253  return iterator(this->get_poly()._container().begin());
254  }
255  iterator end() const
256  {
257  return iterator(this->get_poly()._container().end());
258  }
259  r_iterator obegin() const
260  {
261  return r_iterator(this->template rcp_from_this_cast<Poly>(),
262  (long)size() - 1);
263  }
264  r_iterator oend() const
265  {
266  return r_iterator(this->template rcp_from_this_cast<Poly>(), -1);
267  }
268 };
269 
270 class UIntPolyPiranha
271  : public UPiranhaPoly<pintpoly, UIntPolyBase, UIntPolyPiranha>
272 {
273 public:
274  IMPLEMENT_TYPEID(SYMENGINE_UINTPOLYPIRANHA)
276  UIntPolyPiranha(const RCP<const Basic> &var, pintpoly &&dict);
278  hash_t __hash__() const;
279 
280 }; // UIntPolyPiranha
281 
282 class URatPolyPiranha
283  : public UPiranhaPoly<pratpoly, URatPolyBase, URatPolyPiranha>
284 {
285 public:
286  IMPLEMENT_TYPEID(SYMENGINE_URATPOLYPIRANHA)
288  URatPolyPiranha(const RCP<const Basic> &var, pratpoly &&dict);
290  hash_t __hash__() const;
291 };
292 
293 inline RCP<const UIntPolyPiranha> gcd_upoly(const UIntPolyPiranha &a,
294  const UIntPolyPiranha &b)
295 {
296  if (!(a.get_var()->__eq__(*b.get_var())))
297  throw SymEngineException("Error: variables must agree.");
298 
299  pintpoly gcdx(std::get<0>(pintpoly::gcd(a.get_poly(), b.get_poly())));
300  // following the convention, that leading coefficient should be positive
301  if (gcdx.find_cf(pmonomial{gcdx.degree()}) < 0)
302  piranha::math::negate(gcdx);
303  return make_rcp<const UIntPolyPiranha>(a.get_var(), std::move(gcdx));
304 }
305 
306 inline RCP<const URatPolyPiranha> gcd_upoly(const URatPolyPiranha &a,
307  const URatPolyPiranha &b)
308 {
309  if (!(a.get_var()->__eq__(*b.get_var())))
310  throw SymEngineException("Error: variables must agree.");
311 
312  pratpoly gcdx(std::get<0>(pratpoly::gcd(a.get_poly(), b.get_poly())));
313  // following the convention, that polynomial should be monic
314  gcdx *= (1 / gcdx.find_cf(pmonomial{gcdx.degree()}));
315  return make_rcp<const URatPolyPiranha>(a.get_var(), std::move(gcdx));
316 }
317 
318 inline RCP<const UIntPolyPiranha> lcm_upoly(const UIntPolyPiranha &a,
319  const UIntPolyPiranha &b)
320 {
321  if (!(a.get_var()->__eq__(*b.get_var())))
322  throw SymEngineException("Error: variables must agree.");
323 
324  pintpoly lcmx(std::get<0>(pintpoly::gcd(a.get_poly(), b.get_poly())));
325  lcmx = (a.get_poly() * b.get_poly()) / lcmx;
326  if (lcmx.find_cf(pmonomial{lcmx.degree()}) < 0)
327  piranha::math::negate(lcmx);
328  return make_rcp<const UIntPolyPiranha>(a.get_var(), std::move(lcmx));
329 }
330 
331 inline RCP<const URatPolyPiranha> lcm_upoly(const URatPolyPiranha &a,
332  const URatPolyPiranha &b)
333 {
334  if (!(a.get_var()->__eq__(*b.get_var())))
335  throw SymEngineException("Error: variables must agree.");
336 
337  pratpoly lcmx(std::get<0>(pratpoly::gcd(a.get_poly(), b.get_poly())));
338  lcmx = (a.get_poly() * b.get_poly()) / lcmx;
339  lcmx *= (1 / lcmx.find_cf(pmonomial{lcmx.degree()}));
340  return make_rcp<const URatPolyPiranha>(a.get_var(), std::move(lcmx));
341 }
342 
343 template <typename Container, template <typename X, typename Y> class BaseType,
344  typename Poly>
345 RCP<const Poly> pow_upoly(const UPiranhaPoly<Container, BaseType, Poly> &a,
346  unsigned int p)
347 {
348  return make_rcp<const Poly>(a.get_var(),
349  std::move(piranha::math::pow(a.get_poly(), p)));
350 }
351 
352 template <typename Container, template <typename X, typename Y> class BaseType,
353  typename Poly>
354 bool divides_upoly(const UPiranhaPoly<Container, BaseType, Poly> &a,
355  const Poly &b, const Ptr<RCP<const Poly>> &res)
356 {
357  if (!(a.get_var()->__eq__(*b.get_var())))
358  throw SymEngineException("Error: variables must agree.");
359 
360  try {
361  Container z;
362  piranha::math::divexact(z, b.get_poly(), a.get_poly());
363  *res = Poly::from_container(a.get_var(), std::move(z));
364  return true;
365  } catch (const piranha::math::inexact_division &) {
366  return false;
367  }
368 }
369 } // namespace SymEngine
370 
371 #endif // HAVE_SYMENGINE_PIRANHA
372 
373 #endif // SYMENGINE_UINTPOLY_PIRANHA_H
#define IMPLEMENT_TYPEID(SYMENGINE_ID)
Inline members and functions.
Definition: basic.h:340
T begin(T... args)
The lowest unit of symbolic representation.
Definition: basic.h:95
T end(T... args)
T make_pair(T... args)
T move(T... args)
Main namespace for SymEngine package.
Definition: add.cpp:19
RCP< const Symbol > symbol(const std::string &name)
inline version to return Symbol
Definition: symbol.h:82
STL namespace.
T operator!=(T... args)
T size(T... args)