solve.cpp
1 #include <symengine/solve.h>
2 #include <symengine/polys/basic_conversions.h>
3 #include <symengine/logic.h>
4 #include <symengine/mul.h>
5 #include <symengine/as_real_imag.cpp>
6 
7 namespace SymEngine
8 {
9 
10 RCP<const Set> solve_poly_linear(const vec_basic &coeffs,
11  const RCP<const Set> &domain)
12 {
13  if (coeffs.size() != 2) {
14  throw SymEngineException("Expected a polynomial of degree 1. Try with "
15  "solve() or solve_poly()");
16  }
17  auto root = neg(div(coeffs[0], coeffs[1]));
18  return set_intersection({domain, finiteset({root})});
19 }
20 
21 RCP<const Set> solve_poly_quadratic(const vec_basic &coeffs,
22  const RCP<const Set> &domain)
23 {
24  if (coeffs.size() != 3) {
25  throw SymEngineException("Expected a polynomial of degree 2. Try with "
26  "solve() or solve_poly()");
27  }
28 
29  auto a = coeffs[2];
30  auto b = div(coeffs[1], a), c = div(coeffs[0], a);
31  RCP<const Basic> root1, root2;
32  if (eq(*c, *zero)) {
33  root1 = neg(b);
34  root2 = zero;
35  } else if (eq(*b, *zero)) {
36  root1 = sqrt(neg(c));
37  root2 = neg(root1);
38  } else {
39  auto discriminant = sub(mul(b, b), mul(integer(4), c));
40  auto lterm = div(neg(b), integer(2));
41  auto rterm = div(sqrt(discriminant), integer(2));
42  root1 = add(lterm, rterm);
43  root2 = sub(lterm, rterm);
44  }
45  return set_intersection({domain, finiteset({root1, root2})});
46 }
47 
48 RCP<const Set> solve_poly_cubic(const vec_basic &coeffs,
49  const RCP<const Set> &domain)
50 {
51  if (coeffs.size() != 4) {
52  throw SymEngineException("Expected a polynomial of degree 3. Try with "
53  "solve() or solve_poly()");
54  }
55 
56  auto a = coeffs[3];
57  auto b = div(coeffs[2], a), c = div(coeffs[1], a), d = div(coeffs[0], a);
58 
59  // ref :
60  // https://en.wikipedia.org/wiki/Cubic_function#General_solution_to_the_cubic_equation_with_real_coefficients
61  auto i2 = integer(2), i3 = integer(3), i4 = integer(4), i9 = integer(9),
62  i27 = integer(27);
63 
64  RCP<const Basic> root1, root2, root3;
65  if (eq(*d, *zero)) {
66  root1 = zero;
67  auto fset = solve_poly_quadratic({c, b, one}, domain);
68  SYMENGINE_ASSERT(is_a<FiniteSet>(*fset));
69  auto cont = down_cast<const FiniteSet &>(*fset).get_container();
70  if (cont.size() == 2) {
71  root2 = *cont.begin();
72  root3 = *std::next(cont.begin());
73  } else {
74  root2 = root3 = *cont.begin();
75  }
76  } else {
77  auto delta0 = sub(mul(b, b), mul(i3, c));
78  auto delta1
79  = add(sub(mul(pow(b, i3), i2), mul({i9, b, c})), mul(i27, d));
80  auto delta = div(sub(mul(i4, pow(delta0, i3)), pow(delta1, i2)), i27);
81  if (eq(*delta, *zero)) {
82  if (eq(*delta0, *zero)) {
83  root1 = root2 = root3 = div(neg(b), i3);
84  } else {
85  root1 = root2
86  = div(sub(mul(i9, d), mul(b, c)), mul(i2, delta0));
87  root3 = div(sub(mul({i4, b, c}), add(mul(d, i9), pow(b, i3))),
88  delta0);
89  }
90  } else {
91  auto temp = sqrt(mul(neg(i27), delta));
92  auto Cexpr = div(add(delta1, temp), i2);
93  if (eq(*Cexpr, *zero)) {
94  Cexpr = div(sub(delta1, temp), i2);
95  }
96  auto C = pow(Cexpr, div(one, i3));
97  root1 = neg(div(add(b, add(C, div(delta0, C))), i3));
98  auto coef = div(mul(I, sqrt(i3)), i2);
99  temp = neg(div(one, i2));
100  auto cbrt1 = add(temp, coef);
101  auto cbrt2 = sub(temp, coef);
102  root2 = neg(div(
103  add(b, add(mul(cbrt1, C), div(delta0, mul(cbrt1, C)))), i3));
104  root3 = neg(div(
105  add(b, add(mul(cbrt2, C), div(delta0, mul(cbrt2, C)))), i3));
106  }
107  }
108  return set_intersection({domain, finiteset({root1, root2, root3})});
109 }
110 
111 RCP<const Set> solve_poly_quartic(const vec_basic &coeffs,
112  const RCP<const Set> &domain)
113 {
114  if (coeffs.size() != 5) {
115  throw SymEngineException("Expected a polynomial of degree 4. Try with "
116  "solve() or solve_poly()");
117  }
118 
119  auto i2 = integer(2), i3 = integer(3), i4 = integer(4), i8 = integer(8),
120  i16 = integer(16), i64 = integer(64), i256 = integer(256);
121 
122  // ref : http://mathforum.org/dr.math/faq/faq.cubic.equations.html
123  auto lc = coeffs[4];
124  auto a = div(coeffs[3], lc), b = div(coeffs[2], lc), c = div(coeffs[1], lc),
125  d = div(coeffs[0], lc);
126  set_basic roots;
127 
128  if (eq(*d, *zero)) {
129  vec_basic newcoeffs(4);
130  newcoeffs[0] = c, newcoeffs[1] = b, newcoeffs[2] = a,
131  newcoeffs[3] = one;
132  auto rcubic = solve_poly_cubic(newcoeffs, domain);
133  SYMENGINE_ASSERT(is_a<FiniteSet>(*rcubic));
134  roots = down_cast<const FiniteSet &>(*rcubic).get_container();
135  roots.insert(zero);
136  } else {
137  // substitute x = y-a/4 to get equation of the form y**4 + e*y**2 + f*y
138  // + g = 0
139  auto sqa = mul(a, a);
140  auto cba = mul(sqa, a);
141  auto aby4 = div(a, i4);
142  auto e = sub(b, div(mul(i3, sqa), i8));
143  auto ff = sub(add(c, div(cba, i8)), div(mul(a, b), i2));
144  auto g = sub(add(d, div(mul(sqa, b), i16)),
145  add(div(mul(a, c), i4), div(mul({i3, cba, a}), i256)));
146 
147  // two special cases
148  if (eq(*g, *zero)) {
149  vec_basic newcoeffs(4);
150  newcoeffs[0] = ff, newcoeffs[1] = e, newcoeffs[2] = zero,
151  newcoeffs[3] = one;
152  auto rcubic = solve_poly_cubic(newcoeffs, domain);
153  SYMENGINE_ASSERT(is_a<FiniteSet>(*rcubic));
154  auto rtemp = down_cast<const FiniteSet &>(*rcubic).get_container();
155  SYMENGINE_ASSERT(rtemp.size() > 0 and rtemp.size() <= 3);
156  for (auto &r : rtemp) {
157  roots.insert(sub(r, aby4));
158  }
159  roots.insert(neg(aby4));
160  } else if (eq(*ff, *zero)) {
161  vec_basic newcoeffs(3);
162  newcoeffs[0] = g, newcoeffs[1] = e, newcoeffs[2] = one;
163  auto rquad = solve_poly_quadratic(newcoeffs, domain);
164  SYMENGINE_ASSERT(is_a<FiniteSet>(*rquad));
165  auto rtemp = down_cast<const FiniteSet &>(*rquad).get_container();
166  SYMENGINE_ASSERT(rtemp.size() > 0 and rtemp.size() <= 2);
167  for (auto &r : rtemp) {
168  auto sqrtr = sqrt(r);
169  roots.insert(sub(sqrtr, aby4));
170  roots.insert(sub(neg(sqrtr), aby4));
171  }
172  } else {
173  // Leonhard Euler's method
174  vec_basic newcoeffs(4);
175  newcoeffs[0] = neg(div(mul(ff, ff), i64)),
176  newcoeffs[1] = div(sub(mul(e, e), mul(i4, g)), i16),
177  newcoeffs[2] = div(e, i2);
178  newcoeffs[3] = one;
179 
180  auto rcubic = solve_poly_cubic(newcoeffs);
181  SYMENGINE_ASSERT(is_a<FiniteSet>(*rcubic));
182  roots = down_cast<const FiniteSet &>(*rcubic).get_container();
183  SYMENGINE_ASSERT(roots.size() > 0 and roots.size() <= 3);
184  auto p = sqrt(*roots.begin());
185  auto q = p;
186  if (roots.size() > 1) {
187  q = sqrt(*std::next(roots.begin()));
188  }
189  auto r = div(neg(ff), mul({i8, p, q}));
190  roots.clear();
191  roots.insert(add({p, q, r, neg(aby4)}));
192  roots.insert(add({p, neg(q), neg(r), neg(aby4)}));
193  roots.insert(add({neg(p), q, neg(r), neg(aby4)}));
194  roots.insert(add({neg(p), neg(q), r, neg(aby4)}));
195  }
196  }
197  return set_intersection({domain, finiteset(roots)});
198 }
199 
200 RCP<const Set> solve_poly_heuristics(const vec_basic &coeffs,
201  const RCP<const Set> &domain)
202 {
203  auto degree = coeffs.size() - 1;
204  switch (degree) {
205  case 0: {
206  if (eq(*coeffs[0], *zero)) {
207  return domain;
208  } else {
209  return emptyset();
210  }
211  }
212  case 1:
213  return solve_poly_linear(coeffs, domain);
214  case 2:
215  return solve_poly_quadratic(coeffs, domain);
216  case 3:
217  return solve_poly_cubic(coeffs, domain);
218  case 4:
219  return solve_poly_quartic(coeffs, domain);
220  default:
221  throw SymEngineException(
222  "expected a polynomial of order between 0 to 4");
223  }
224 }
225 
226 inline RCP<const Basic> get_coeff_basic(const integer_class &i)
227 {
228  return integer(i);
229 }
230 
231 inline RCP<const Basic> get_coeff_basic(const Expression &i)
232 {
233  return i.get_basic();
234 }
235 
236 template <typename Poly>
237 inline vec_basic extract_coeffs(const RCP<const Poly> &f)
238 {
239  int degree = f->get_degree();
240  vec_basic coeffs;
241  for (int i = 0; i <= degree; i++)
242  coeffs.push_back(get_coeff_basic(f->get_coeff(i)));
243  return coeffs;
244 }
245 
246 RCP<const Set> solve_poly(const RCP<const Basic> &f,
247  const RCP<const Symbol> &sym,
248  const RCP<const Set> &domain)
249 {
250 
251 #if defined(HAVE_SYMENGINE_FLINT) && __FLINT_RELEASE > 20502
252  try {
253  auto poly = from_basic<UIntPolyFlint>(f, sym);
254  auto fac = factors(*poly);
255  set_set solns;
256  for (const auto &elem : fac) {
257  auto uip = UIntPoly::from_poly(*elem.first);
258  auto degree = uip->get_poly().degree();
259  if (degree <= 4) {
260  solns.insert(
261  solve_poly_heuristics(extract_coeffs(uip), domain));
262  } else {
263  solns.insert(
264  conditionset(sym, logical_and({Eq(uip->as_symbolic(), zero),
265  domain->contains(sym)})));
266  }
267  }
268  return SymEngine::set_union(solns);
269  } catch (SymEngineException &x) {
270  // Try next
271  }
272 #endif
273  RCP<const Basic> gen = rcp_static_cast<const Basic>(sym);
274  auto uexp = from_basic<UExprPoly>(f, gen);
275  auto degree = uexp->get_degree();
276  if (degree <= 4) {
277  return solve_poly_heuristics(extract_coeffs(uexp), domain);
278  } else {
279  return conditionset(sym,
280  logical_and({Eq(f, zero), domain->contains(sym)}));
281  }
282 }
283 
284 RCP<const Set> solve_rational(const RCP<const Basic> &f,
285  const RCP<const Symbol> &sym,
286  const RCP<const Set> &domain)
287 {
288  RCP<const Basic> num, den;
289  as_numer_denom(f, outArg(num), outArg(den));
290  if (has_symbol(*den, *sym)) {
291  auto numsoln = solve(num, sym, domain);
292  auto densoln = solve(den, sym, domain);
293  return set_complement(numsoln, densoln);
294  }
295  return solve_poly(num, sym, domain);
296 }
297 
298 /* Helper Visitors for solve_trig */
299 
301  : public BaseVisitor<IsALinearArgTrigVisitor, LocalStopVisitor>
302 {
303 protected:
304  Ptr<const Symbol> x_;
305  bool is_;
306 
307 public:
308  IsALinearArgTrigVisitor(Ptr<const Symbol> x) : x_(x) {}
309 
310  bool apply(const Basic &b)
311  {
312  stop_ = false;
313  is_ = true;
314  preorder_traversal_local_stop(b, *this);
315  return is_;
316  }
317 
318  bool apply(const RCP<const Basic> &b)
319  {
320  return apply(*b);
321  }
322 
323  void bvisit(const Basic &x)
324  {
325  local_stop_ = false;
326  }
327 
328  void bvisit(const Symbol &x)
329  {
330  if (x_->__eq__(x)) {
331  is_ = 0;
332  stop_ = true;
333  }
334  }
335 
336  template <typename T,
337  typename
338  = enable_if_t<std::is_base_of<TrigFunction, T>::value
340  void bvisit(const T &x)
341  {
342  is_ = (from_basic<UExprPoly>(x.get_args()[0], (*x_).rcp_from_this())
343  ->get_degree()
344  <= 1);
345  if (not is_)
346  stop_ = true;
347  local_stop_ = true;
348  }
349 };
350 
351 bool is_a_LinearArgTrigEquation(const Basic &b, const Symbol &x)
352 {
353  IsALinearArgTrigVisitor v(ptrFromRef(x));
354  return v.apply(b);
355 }
356 
357 class InvertComplexVisitor : public BaseVisitor<InvertComplexVisitor>
358 {
359 protected:
360  RCP<const Set> result_;
361  RCP<const Set> gY_;
362  RCP<const Dummy> nD_;
363  RCP<const Symbol> sym_;
364  RCP<const Set> domain_;
365 
366 public:
367  InvertComplexVisitor(RCP<const Set> gY, RCP<const Dummy> nD,
368  RCP<const Symbol> sym, RCP<const Set> domain)
369  : gY_(gY), nD_(nD), sym_(sym), domain_(domain)
370  {
371  }
372 
373  void bvisit(const Basic &x)
374  {
375  result_ = gY_;
376  }
377 
378  void bvisit(const Add &x)
379  {
380  vec_basic f1X, f2X;
381  for (auto &elem : x.get_args()) {
382  if (has_symbol(*elem, *sym_)) {
383  f1X.push_back(elem);
384  } else {
385  f2X.push_back(elem);
386  }
387  }
388  auto depX = add(f1X), indepX = add(f2X);
389  if (not eq(*indepX, *zero)) {
390  gY_ = imageset(nD_, sub(nD_, indepX), gY_);
391  result_ = apply(*depX);
392  } else {
393  result_ = gY_;
394  }
395  }
396 
397  void bvisit(const Mul &x)
398  {
399  vec_basic f1X, f2X;
400  for (auto &elem : x.get_args()) {
401  if (has_symbol(*elem, *sym_)) {
402  f1X.push_back(elem);
403  } else {
404  f2X.push_back(elem);
405  }
406  }
407  auto depX = mul(f1X), indepX = mul(f2X);
408  if (not eq(*indepX, *one)) {
409  if (eq(*indepX, *NegInf) or eq(*indepX, *Inf)
410  or eq(*indepX, *ComplexInf)) {
411  result_ = emptyset();
412  } else {
413  gY_ = imageset(nD_, div(nD_, indepX), gY_);
414  result_ = apply(*depX);
415  }
416  } else {
417  result_ = gY_;
418  }
419  }
420 
421  void bvisit(const Pow &x)
422  {
423  if (eq(*x.get_base(), *E) and is_a<FiniteSet>(*gY_)) {
424  set_set inv;
425  for (const auto &elem :
426  down_cast<const FiniteSet &>(*gY_).get_container()) {
427  if (eq(*elem, *zero))
428  continue;
429  RCP<const Basic> re, im;
430  as_real_imag(elem, outArg(re), outArg(im));
431  auto logabs = log(add(mul(re, re), mul(im, im)));
432  auto logarg = atan2(im, re);
433  inv.insert(imageset(
434  nD_,
435  add(mul(add(mul({integer(2), nD_, pi}), logarg), I),
436  div(logabs, integer(2))),
437  interval(NegInf, Inf, true,
438  true))); // TODO : replace interval(-oo,oo) with
439  // Set of Integers once Class for Range is implemented.
440  }
441  gY_ = set_union(inv);
442  apply(*x.get_exp());
443  return;
444  }
445  result_ = gY_;
446  }
447 
448  RCP<const Set> apply(const Basic &b)
449  {
450  result_ = gY_;
451  b.accept(*this);
452  return set_intersection({domain_, result_});
453  }
454 };
455 
456 RCP<const Set> invertComplex(const RCP<const Basic> &fX,
457  const RCP<const Set> &gY,
458  const RCP<const Symbol> &sym,
459  const RCP<const Dummy> &nD,
460  const RCP<const Set> &domain)
461 {
462  InvertComplexVisitor v(gY, nD, sym, domain);
463  return v.apply(*fX);
464 }
465 
466 RCP<const Set> solve_trig(const RCP<const Basic> &f,
467  const RCP<const Symbol> &sym,
468  const RCP<const Set> &domain)
469 {
470  // TODO : first simplify f using `fu`.
471  auto exp_f = rewrite_as_exp(f);
472  RCP<const Basic> num, den;
473  as_numer_denom(exp_f, outArg(num), outArg(den));
474 
475  auto xD = dummy("x");
476  map_basic_basic d;
477  auto temp = exp(mul(I, sym));
478  d[temp] = xD;
479  num = expand(num), den = expand(den);
480  num = num->subs(d);
481  den = den->subs(d);
482 
483  if (has_symbol(*num, *sym) or has_symbol(*den, *sym)) {
484  return conditionset(sym, logical_and({Eq(f, zero)}));
485  }
486 
487  auto soln = set_complement(solve(num, xD), solve(den, xD));
488  if (eq(*soln, *emptyset()))
489  return emptyset();
490  else if (is_a<FiniteSet>(*soln)) {
491  set_set res;
492  auto nD
493  = dummy("n"); // use the same dummy for finding every solution set.
494  auto n = symbol(
495  "n"); // replaces the above dummy in final set of solutions.
496  map_basic_basic d;
497  d[nD] = n;
498  for (const auto &elem :
499  down_cast<const FiniteSet &>(*soln).get_container()) {
500  res.insert(
501  invertComplex(exp(mul(I, sym)), finiteset({elem}), sym, nD));
502  }
503  auto ans = set_union(res)->subs(d);
504  if (not is_a_Set(*ans))
505  throw SymEngineException("Expected an object of type Set");
506  return set_intersection({rcp_static_cast<const Set>(ans), domain});
507  }
508  return conditionset(sym, logical_and({Eq(f, zero), domain->contains(sym)}));
509 }
510 
511 RCP<const Set> solve(const RCP<const Basic> &f, const RCP<const Symbol> &sym,
512  const RCP<const Set> &domain)
513 {
514  if (eq(*f, *boolTrue))
515  return domain;
516  if (eq(*f, *boolFalse))
517  return emptyset();
518  if (is_a<Equality>(*f)) {
519  return solve(sub(down_cast<const Relational &>(*f).get_arg1(),
520  down_cast<const Relational &>(*f).get_arg2()),
521  sym, domain);
522  } else if (is_a<Unequality>(*f)) {
523  auto soln = solve(sub(down_cast<const Relational &>(*f).get_arg1(),
524  down_cast<const Relational &>(*f).get_arg2()),
525  sym, domain);
526  return set_complement(domain, soln);
527  } else if (is_a_Relational(*f)) {
528  // Solving inequalities is not implemented yet.
529  return conditionset(sym, logical_and({rcp_static_cast<const Boolean>(f),
530  domain->contains(sym)}));
531  }
532 
533  if (is_a_Number(*f)) {
534  if (eq(*f, *zero)) {
535  return domain;
536  } else {
537  return emptyset();
538  }
539  }
540 
541  if (not has_symbol(*f, *sym))
542  return emptyset();
543 
544  if (is_a_LinearArgTrigEquation(*f, *sym)) {
545  return solve_trig(f, sym, domain);
546  }
547 
548  if (is_a<Mul>(*f)) {
549  auto args = f->get_args();
550  set_set solns;
551  for (auto &a : args) {
552  solns.insert(solve(a, sym, domain));
553  }
554  return SymEngine::set_union(solns);
555  }
556 
557  return solve_rational(f, sym, domain);
558 }
559 
560 vec_basic linsolve_helper(const DenseMatrix &A, const DenseMatrix &b)
561 {
562  DenseMatrix res(A.nrows(), 1);
563  fraction_free_gauss_jordan_solve(A, b, res);
564  vec_basic fs;
565  for (unsigned i = 0; i < res.nrows(); i++) {
566  fs.push_back(res.get(i, 0));
567  }
568  return fs;
569 }
570 
571 vec_basic linsolve(const DenseMatrix &system, const vec_sym &syms)
572 {
573  DenseMatrix A(system.nrows(), system.ncols() - 1), b(system.nrows(), 1);
574  system.submatrix(A, 0, 0, system.nrows() - 1, system.ncols() - 2);
575  system.submatrix(b, 0, system.ncols() - 1, system.nrows() - 1,
576  system.ncols() - 1);
577  return linsolve_helper(A, b);
578 }
579 
580 vec_basic linsolve(const vec_basic &system, const vec_sym &syms)
581 {
582  auto mat = linear_eqns_to_matrix(system, syms);
583  DenseMatrix A = mat.first, b = mat.second;
584  return linsolve_helper(A, b);
585 }
586 
587 set_basic get_set_from_vec(const vec_sym &syms)
588 {
589  set_basic sb;
590  for (auto &s : syms)
591  sb.insert(s);
592  return sb;
593 }
594 
596 linear_eqns_to_matrix(const vec_basic &equations, const vec_sym &syms)
597 {
598  auto size = numeric_cast<unsigned int>(syms.size());
599  DenseMatrix A(numeric_cast<unsigned int>(equations.size()), size);
600  zeros(A);
601  vec_basic bvec;
602 
603  int row = 0;
604  auto gens = get_set_from_vec(syms);
605  umap_basic_uint index_of_sym;
606  for (unsigned int i = 0; i < size; i++) {
607  index_of_sym[syms[i]] = i;
608  }
609  for (const auto &eqn : equations) {
610  auto neqn = eqn;
611  if (is_a<Equality>(*eqn)) {
612  neqn = sub(down_cast<const Equality &>(*eqn).get_arg2(),
613  down_cast<const Equality &>(*eqn).get_arg1());
614  }
615 
616  auto mpoly = from_basic<MExprPoly>(neqn, gens);
617  RCP<const Basic> rem = zero;
618  for (const auto &p : mpoly->get_poly().dict_) {
619  RCP<const Basic> res = (p.second.get_basic());
620  int whichvar = 0, non_zero = 0;
621  RCP<const Basic> cursim;
622  for (auto &sym : gens) {
623  if (0 != p.first[whichvar]) {
624  non_zero++;
625  cursim = sym;
626  if (p.first[whichvar] != 1 or non_zero == 2) {
627  throw SymEngineException("Expected a linear equation.");
628  }
629  }
630  whichvar++;
631  }
632  if (not non_zero) {
633  rem = res;
634  } else {
635  A.set(row, index_of_sym[cursim], res);
636  }
637  }
638  bvec.push_back(neg(rem));
639  ++row;
640  }
641  return std::make_pair(
642  A, DenseMatrix(numeric_cast<unsigned int>(equations.size()), 1, bvec));
643 }
644 } // namespace SymEngine
The base class for representing addition in symbolic expressions.
Definition: add.h:27
vec_basic get_args() const override
Returns the arguments of the Add.
Definition: add.cpp:397
The lowest unit of symbolic representation.
Definition: basic.h:97
vec_basic get_args() const override
Returns the list of arguments.
Definition: mul.cpp:507
RCP< const Basic > get_exp() const
Definition: pow.h:42
RCP< const Basic > get_base() const
Definition: pow.h:37
T insert(T... args)
T make_pair(T... args)
Main namespace for SymEngine package.
Definition: add.cpp:19
bool is_a_Number(const Basic &b)
Definition: number.h:130
RCP< const Set > interval(const RCP< const Number > &start, const RCP< const Number > &end, const bool left_open=false, const bool right_open=false)
Definition: sets.h:611
RCP< const Basic > div(const RCP< const Basic > &a, const RCP< const Basic > &b)
Division.
Definition: mul.cpp:431
RCP< const Dummy > dummy()
inline version to return Dummy
Definition: symbol.h:88
std::enable_if< std::is_integral< T >::value, RCP< const Integer > >::type integer(T i)
Definition: integer.h:197
RCP< const Symbol > symbol(const std::string &name)
inline version to return Symbol
Definition: symbol.h:82
bool eq(const Basic &a, const Basic &b)
Checks equality for a and b
Definition: basic-inl.h:21
RCP< const EmptySet > emptyset()
Definition: sets.h:590
RCP< const Basic > atan2(const RCP< const Basic > &num, const RCP< const Basic > &den)
Canonicalize ATan2:
Definition: functions.cpp:1614
RCP< const Basic > sub(const RCP< const Basic > &a, const RCP< const Basic > &b)
Substracts b from a.
Definition: add.cpp:495
RCP< const Basic > exp(const RCP< const Basic > &x)
Returns the natural exponential function E**x = pow(E, x)
Definition: pow.cpp:271
RCP< const Basic > mul(const RCP< const Basic > &a, const RCP< const Basic > &b)
Multiplication.
Definition: mul.cpp:352
RCP< const Set > conditionset(const RCP< const Basic > &sym, const RCP< const Boolean > &condition)
Definition: sets.cpp:1781
RCP< const Basic > log(const RCP< const Basic > &arg)
Returns the Natural Logarithm from argument arg
Definition: functions.cpp:1774
RCP< const Basic > add(const RCP< const Basic > &a, const RCP< const Basic > &b)
Adds two objects (safely).
Definition: add.cpp:425
RCP< const Boolean > Eq(const RCP< const Basic > &lhs)
Returns the canonicalized Equality object from a single argument.
Definition: logic.cpp:653
RCP< const Set > finiteset(const set_basic &container)
Definition: sets.h:602
RCP< const Basic > expand(const RCP< const Basic > &self, bool deep=true)
Expands self
Definition: expand.cpp:369
RCP< const Basic > neg(const RCP< const Basic > &a)
Negation.
Definition: mul.cpp:443
T next(T... args)
T pow(T... args)
T push_back(T... args)
T set_intersection(T... args)
T set_union(T... args)
T sqrt(T... args)
T system(T... args)