derivative.cpp
1 #include <symengine/visitor.h>
2 #include <symengine/subs.h>
3 #include <symengine/symengine_casts.h>
4 #include <symengine/derivative.h>
5 
6 namespace SymEngine
7 {
8 
9 extern RCP<const Basic> i2;
10 
11 // Needs create(vec_basic) method to be used.
12 template <typename T>
13 static inline RCP<const Basic> fdiff(const T &self, RCP<const Symbol> x,
14  DiffVisitor &visitor)
15 {
16  RCP<const Basic> diff = zero;
17  RCP<const Basic> ret;
18  bool know_deriv;
19 
20  vec_basic v = self.get_args();
21  vec_basic vdiff(v.size());
22 
23  unsigned count = 0;
24  for (unsigned i = 0; i < v.size(); i++) {
25  vdiff[i] = visitor.apply(v[i]);
26  if (neq(*vdiff[i], *zero)) {
27  count++;
28  }
29  }
30 
31  if (count == 0) {
32  return diff;
33  }
34 
35  for (unsigned i = 0; i < v.size(); i++) {
36  if (eq(*vdiff[i], *zero))
37  continue;
38  know_deriv = fdiff(outArg(ret), self, i);
39  if (know_deriv) {
40  diff = add(diff, mul(ret, vdiff[i]));
41  } else {
42  if (count == 1 and eq(*v[i], *x)) {
43  return Derivative::create(self.rcp_from_this(), {x});
44  }
45  vec_basic new_args = v;
47  stm << (i + 1);
48  new_args[i] = get_dummy(self, "xi_" + stm.str());
49  map_basic_basic m;
50  insert(m, new_args[i], v[i]);
51  diff = add(diff, mul(vdiff[i],
52  make_rcp<const Subs>(
53  Derivative::create(self.create(new_args),
54  {new_args[i]}),
55  m)));
56  }
57  }
58  return diff;
59 }
60 
61 static bool fdiff(const Ptr<RCP<const Basic>> &ret, const Zeta &self,
62  unsigned index)
63 {
64  if (index == 1) {
65  *ret = mul(mul(minus_one, self.get_arg1()),
66  zeta(add(self.get_arg1(), one), self.get_arg2()));
67  return true;
68  } else {
69  return false;
70  }
71 }
72 
73 static bool fdiff(const Ptr<RCP<const Basic>> &ret, const UpperGamma &self,
74  unsigned index)
75 {
76  if (index == 1) {
77  *ret = mul(mul(pow(self.get_arg2(), sub(self.get_arg1(), one)),
78  exp(neg(self.get_arg2()))),
79  minus_one);
80  return true;
81  } else {
82  return false;
83  }
84 }
85 
86 static bool fdiff(const Ptr<RCP<const Basic>> &ret, const LowerGamma &self,
87  unsigned index)
88 {
89  if (index == 1) {
90  *ret = mul(pow(self.get_arg2(), sub(self.get_arg1(), one)),
91  exp(neg(self.get_arg2())));
92  return true;
93  } else {
94  return false;
95  }
96 }
97 
98 static bool fdiff(const Ptr<RCP<const Basic>> &ret, const PolyGamma &self,
99  unsigned index)
100 {
101  if (index == 1) {
102  *ret = polygamma(add(self.get_arg1(), one), self.get_arg2());
103  return true;
104  } else {
105  return false;
106  }
107 }
108 
109 static bool fdiff(const Ptr<RCP<const Basic>> &ret, const Function &self,
110  unsigned index)
111 {
112  // Don't know the derivative, fallback to `Derivative` instances
113  return false;
114 }
115 
116 template <typename P>
117 static inline RCP<const Basic> diff_upolyflint(const P &self, const Symbol &x)
118 {
119  if (self.get_var()->__eq__(x)) {
120  return P::from_container(self.get_var(), self.get_poly().derivative());
121  } else {
122  return P::from_dict(self.get_var(), {{}});
123  }
124 }
125 
126 template <typename Poly, typename Dict>
127 static inline RCP<const Basic> diff_upoly(const Poly &self, const Symbol &x)
128 {
129  if (self.get_var()->__eq__(x)) {
130  Dict d;
131  for (auto it = self.begin(); it != self.end(); ++it) {
132  if (it->first != 0)
133  d[it->first - 1] = it->second * it->first;
134  }
135  return Poly::from_dict(self.get_var(), std::move(d));
136  } else {
137  return Poly::from_dict(self.get_var(), {{}});
138  }
139 }
140 
141 template <typename Container, typename Poly>
142 static RCP<const Basic> diff_mpoly(const MSymEnginePoly<Container, Poly> &self,
143  const RCP<const Symbol> &x)
144 {
145  using Dict = typename Container::dict_type;
146  using Vec = typename Container::vec_type;
147  Dict dict;
148 
149  if (self.get_vars().find(x) != self.get_vars().end()) {
150  auto i = self.get_vars().begin();
151  unsigned int index = 0;
152  while (!(*i)->__eq__(*x)) {
153  i++;
154  index++;
155  } // find the index of the variable we are differentiating WRT.
156  for (auto bucket : self.get_poly().dict_) {
157  if (bucket.first[index] != 0) {
158  Vec v = bucket.first;
159  v[index]--;
160  dict.insert({v, bucket.second * bucket.first[index]});
161  }
162  }
163  vec_basic v;
164  v.insert(v.begin(), self.get_vars().begin(), self.get_vars().end());
165  return Poly::from_dict(v, std::move(dict));
166  } else {
167  vec_basic vs;
168  vs.insert(vs.begin(), self.get_vars().begin(), self.get_vars().end());
169  return Poly::from_dict(vs, {{}});
170  }
171 }
172 
173 #ifndef debug_methods
174 
175 void DiffVisitor::bvisit(const Basic &self)
176 {
177  result_ = Derivative::create(self.rcp_from_this(), {x});
178 }
179 
180 #else
181 // Here we do not have a 'Basic' fallback, but rather must implement all
182 // virtual methods explicitly (if we miss one, the code will not compile).
183 // This is useful to check that we have implemented all methods that we
184 // wanted.
185 
186 #define DIFF0(CLASS) \
187  void DiffVisitor::bvisit(const CLASS &self) \
188  { \
189  result_ = Derivative::create(self.rcp_from_this(), {x}); \
190  }
191 
192 DIFF0(UnivariateSeries)
193 DIFF0(Max)
194 DIFF0(Min)
195 #endif
196 
197 void DiffVisitor::bvisit(const Number &self)
198 {
199  result_ = zero;
200 }
201 
202 void DiffVisitor::bvisit(const Constant &self)
203 {
204  result_ = zero;
205 }
206 
207 void DiffVisitor::bvisit(const Symbol &self)
208 {
209  if (x->get_name() == self.get_name()) {
210  result_ = one;
211  } else {
212  result_ = zero;
213  }
214 }
215 
216 void DiffVisitor::bvisit(const Log &self)
217 {
218  apply(self.get_arg());
219  result_ = mul(div(one, self.get_arg()), result_);
220 }
221 
222 void DiffVisitor::bvisit(const Abs &self)
223 {
224  apply(self.get_arg());
225  if (eq(*result_, *zero)) {
226  result_ = zero;
227  } else {
228  result_ = Derivative::create(self.rcp_from_this(), {x});
229  }
230 }
231 
232 void DiffVisitor::bvisit(const Zeta &self)
233 {
234  result_ = fdiff(self, x, *this);
235 }
236 
237 void DiffVisitor::bvisit(const LowerGamma &self)
238 {
239  result_ = fdiff(self, x, *this);
240 }
241 void DiffVisitor::bvisit(const UpperGamma &self)
242 {
243  result_ = fdiff(self, x, *this);
244 }
245 void DiffVisitor::bvisit(const PolyGamma &self)
246 {
247  result_ = fdiff(self, x, *this);
248 }
249 void DiffVisitor::bvisit(const UnevaluatedExpr &self)
250 {
251  result_ = Derivative::create(self.rcp_from_this(), {x});
252 }
253 void DiffVisitor::bvisit(const TwoArgFunction &self)
254 {
255  result_ = fdiff(self, x, *this);
256 }
257 
258 void DiffVisitor::bvisit(const ASech &self)
259 {
260  apply(self.get_arg());
261  result_ = mul(div(minus_one, mul(sqrt(sub(one, pow(self.get_arg(), i2))),
262  self.get_arg())),
263  result_);
264 }
265 
266 void DiffVisitor::bvisit(const ACoth &self)
267 {
268  apply(self.get_arg());
269  result_ = mul(div(one, sub(one, pow(self.get_arg(), i2))), result_);
270 }
271 
272 void DiffVisitor::bvisit(const ATanh &self)
273 {
274  apply(self.get_arg());
275  result_ = mul(div(one, sub(one, pow(self.get_arg(), i2))), result_);
276 }
277 
278 void DiffVisitor::bvisit(const ACosh &self)
279 {
280  apply(self.get_arg());
281  result_ = mul(div(one, sqrt(sub(pow(self.get_arg(), i2), one))), result_);
282 }
283 
284 void DiffVisitor::bvisit(const ACsch &self)
285 {
286  apply(self.get_arg());
287  result_ = mul(
288  div(minus_one, mul(sqrt(add(one, div(one, pow(self.get_arg(), i2)))),
289  pow(self.get_arg(), i2))),
290  result_);
291 }
292 
293 void DiffVisitor::bvisit(const ASinh &self)
294 {
295  apply(self.get_arg());
296  result_ = mul(div(one, sqrt(add(pow(self.get_arg(), i2), one))), result_);
297 }
298 
299 void DiffVisitor::bvisit(const Coth &self)
300 {
301  apply(self.get_arg());
302  result_ = mul(div(minus_one, pow(sinh(self.get_arg()), i2)), result_);
303 }
304 
305 void DiffVisitor::bvisit(const Tanh &self)
306 {
307  apply(self.get_arg());
308  result_ = mul(sub(one, pow(tanh(self.get_arg()), i2)), result_);
309 }
310 
311 void DiffVisitor::bvisit(const Sech &self)
312 {
313  apply(self.get_arg());
314  result_
315  = mul(mul(mul(minus_one, sech(self.get_arg())), tanh(self.get_arg())),
316  result_);
317 }
318 
319 void DiffVisitor::bvisit(const Cosh &self)
320 {
321  apply(self.get_arg());
322  result_ = mul(sinh(self.get_arg()), result_);
323 }
324 
325 void DiffVisitor::bvisit(const Csch &self)
326 {
327  apply(self.get_arg());
328  result_
329  = mul(mul(mul(minus_one, csch(self.get_arg())), coth(self.get_arg())),
330  result_);
331 }
332 
333 void DiffVisitor::bvisit(const Sinh &self)
334 {
335  apply(self.get_arg());
336  result_ = mul(cosh(self.get_arg()), result_);
337 }
338 
339 void DiffVisitor::bvisit(const Subs &self)
340 {
341  RCP<const Basic> d = zero, t;
342  if (self.get_dict().count(x) == 0) {
343  apply(self.get_arg());
344  d = result_->subs(self.get_dict());
345  }
346  for (const auto &p : self.get_dict()) {
347  apply(p.second);
348  t = result_;
349  if (neq(*t, *zero)) {
350  if (is_a<Symbol>(*p.first)) {
351  d = add(d, mul(t, diff(self.get_arg(),
352  rcp_static_cast<const Symbol>(p.first))
353  ->subs(self.get_dict())));
354  } else {
355  result_ = Derivative::create(self.rcp_from_this(), {x});
356  return;
357  }
358  }
359  }
360  result_ = d;
361 }
362 
363 void DiffVisitor::bvisit(const Derivative &self)
364 {
365  apply(self.get_arg());
366  RCP<const Basic> ret = result_;
367  if (eq(*ret, *zero)) {
368  result_ = zero;
369  }
370  multiset_basic t = self.get_symbols();
371  for (auto &p : t) {
372  // If x is already there in symbols multi-set add x to the symbols
373  // multi-set
374  if (eq(*p, *x)) {
375  t.insert(x);
376  result_ = Derivative::create(self.get_arg(), t);
377  return;
378  }
379  }
380  // Avoid cycles
381  if (is_a<Derivative>(*ret)
382  && eq(*down_cast<const Derivative &>(*ret).get_arg(),
383  *self.get_arg())) {
384  t.insert(x);
385  result_ = Derivative::create(self.get_arg(), t);
386  return;
387  }
388  for (auto &p : t) {
389  ret = diff(ret, rcp_static_cast<const Symbol>(p));
390  }
391  result_ = ret;
392 }
393 
394 static inline RCP<const Symbol> get_dummy(const Basic &b, std::string name)
395 {
396  RCP<const Symbol> s;
397  do {
398  name = "_" + name;
399  s = symbol(name);
400  } while (has_symbol(b, *s));
401  return s;
402 }
403 
404 void DiffVisitor::bvisit(const OneArgFunction &self)
405 {
406  result_ = fdiff(self, x, *this);
407 }
408 
409 void DiffVisitor::bvisit(const MultiArgFunction &self)
410 {
411  result_ = fdiff(self, x, *this);
412 }
413 
414 void DiffVisitor::bvisit(const LambertW &self)
415 {
416  // check http://en.wikipedia.org/wiki/Lambert_W_function#Derivative
417  // for the equation
418  apply(self.get_arg());
419  RCP<const Basic> lambertw_val = lambertw(self.get_arg());
420  result_
421  = mul(div(lambertw_val, mul(self.get_arg(), add(lambertw_val, one))),
422  result_);
423 }
424 
425 void DiffVisitor::bvisit(const Add &self)
426 {
428  RCP<const Number> coef = zero, coef2;
429  RCP<const Basic> t;
430  for (auto &p : self.get_dict()) {
431  apply(p.first);
432  RCP<const Basic> term = result_;
433  if (is_a<Integer>(*term)
434  && down_cast<const Integer &>(*term).is_zero()) {
435  continue;
436  } else if (is_a_Number(*term)) {
437  iaddnum(outArg(coef),
438  mulnum(p.second, rcp_static_cast<const Number>(term)));
439  } else if (is_a<Add>(*term)) {
440  for (auto &q : (down_cast<const Add &>(*term)).get_dict())
441  Add::dict_add_term(d, mulnum(q.second, p.second), q.first);
442  iaddnum(outArg(coef),
443  mulnum(p.second, down_cast<const Add &>(*term).get_coef()));
444  } else {
445  Add::as_coef_term(mul(p.second, term), outArg(coef2), outArg(t));
446  Add::dict_add_term(d, coef2, t);
447  }
448  }
449  result_ = Add::from_dict(coef, std::move(d));
450 }
451 
452 void DiffVisitor::bvisit(const Mul &self)
453 {
454  RCP<const Number> overall_coef = zero;
455  umap_basic_num add_dict;
456  for (auto &p : self.get_dict()) {
457  RCP<const Number> coef = self.get_coef();
458  apply(pow(p.first, p.second));
459  RCP<const Basic> factor = result_;
460  if (is_a<Integer>(*factor)
461  && down_cast<const Integer &>(*factor).is_zero())
462  continue;
463  map_basic_basic d = self.get_dict();
464  d.erase(p.first);
465  if (is_a_Number(*factor)) {
466  imulnum(outArg(coef), rcp_static_cast<const Number>(factor));
467  } else if (is_a<Mul>(*factor)) {
468  RCP<const Mul> tmp = rcp_static_cast<const Mul>(factor);
469  imulnum(outArg(coef), tmp->get_coef());
470  for (auto &q : tmp->get_dict()) {
471  Mul::dict_add_term_new(outArg(coef), d, q.second, q.first);
472  }
473  } else {
474  RCP<const Basic> exp, t;
475  Mul::as_base_exp(factor, outArg(exp), outArg(t));
476  Mul::dict_add_term_new(outArg(coef), d, exp, t);
477  }
478  if (d.size() == 0) {
479  iaddnum(outArg(overall_coef), coef);
480  } else {
481  RCP<const Basic> mul = Mul::from_dict(one, std::move(d));
482  Add::coef_dict_add_term(outArg(overall_coef), add_dict, coef, mul);
483  }
484  }
485  result_ = Add::from_dict(overall_coef, std::move(add_dict));
486 }
487 
488 void DiffVisitor::bvisit(const Pow &self)
489 {
490  if (is_a_Number(*(self.get_exp()))) {
491  apply(self.get_base());
492  result_ = mul(
493  mul(self.get_exp(), pow(self.get_base(), sub(self.get_exp(), one))),
494  result_);
495  } else {
496  apply(mul(self.get_exp(), log(self.get_base())));
497  result_ = mul(self.rcp_from_this(), result_);
498  }
499 }
500 
501 void DiffVisitor::bvisit(const Sin &self)
502 {
503  apply(self.get_arg());
504  result_ = mul(cos(self.get_arg()), result_);
505 }
506 
507 void DiffVisitor::bvisit(const Cos &self)
508 {
509  apply(self.get_arg());
510  result_ = mul(mul(minus_one, sin(self.get_arg())), result_);
511 }
512 
513 void DiffVisitor::bvisit(const Tan &self)
514 {
515  apply(self.get_arg());
516  RCP<const Integer> two = integer(2);
517  result_ = mul(add(one, pow(tan(self.get_arg()), two)), result_);
518 }
519 
520 void DiffVisitor::bvisit(const Cot &self)
521 {
522  apply(self.get_arg());
523  RCP<const Integer> two = integer(2);
524  result_
525  = mul(mul(add(one, pow(cot(self.get_arg()), two)), minus_one), result_);
526 }
527 
528 void DiffVisitor::bvisit(const Csc &self)
529 {
530  apply(self.get_arg());
531  result_ = mul(mul(mul(cot(self.get_arg()), csc(self.get_arg())), minus_one),
532  result_);
533 }
534 
535 void DiffVisitor::bvisit(const Sec &self)
536 {
537  apply(self.get_arg());
538  result_ = mul(mul(tan(self.get_arg()), sec(self.get_arg())), result_);
539 }
540 
541 void DiffVisitor::bvisit(const ASin &self)
542 {
543  apply(self.get_arg());
544  result_ = mul(div(one, sqrt(sub(one, pow(self.get_arg(), i2)))), result_);
545 }
546 
547 void DiffVisitor::bvisit(const ACos &self)
548 {
549  apply(self.get_arg());
550  result_
551  = mul(div(minus_one, sqrt(sub(one, pow(self.get_arg(), i2)))), result_);
552 }
553 
554 void DiffVisitor::bvisit(const ASec &self)
555 {
556  apply(self.get_arg());
557  result_
558  = mul(div(one, mul(pow(self.get_arg(), i2),
559  sqrt(sub(one, div(one, pow(self.get_arg(), i2)))))),
560  result_);
561 }
562 
563 void DiffVisitor::bvisit(const ACsc &self)
564 {
565  apply(self.get_arg());
566  result_ = mul(
567  div(minus_one, mul(pow(self.get_arg(), i2),
568  sqrt(sub(one, div(one, pow(self.get_arg(), i2)))))),
569  result_);
570 }
571 
572 void DiffVisitor::bvisit(const ATan &self)
573 {
574  apply(self.get_arg());
575  result_ = mul(div(one, add(one, pow(self.get_arg(), i2))), result_);
576 }
577 
578 void DiffVisitor::bvisit(const ACot &self)
579 {
580  apply(self.get_arg());
581  result_ = mul(div(minus_one, add(one, pow(self.get_arg(), i2))), result_);
582 }
583 
584 void DiffVisitor::bvisit(const ATan2 &self)
585 {
586  apply(div(self.get_num(), self.get_den()));
587  result_ = mul(div(pow(self.get_den(), i2),
588  add(pow(self.get_den(), i2), pow(self.get_num(), i2))),
589  result_);
590 }
591 
592 void DiffVisitor::bvisit(const Erf &self)
593 {
594  apply(self.get_arg());
595  result_ = mul(
596  div(mul(integer(2), exp(neg(mul(self.get_arg(), self.get_arg())))),
597  sqrt(pi)),
598  result_);
599 }
600 
601 void DiffVisitor::bvisit(const Erfc &self)
602 {
603  apply(self.get_arg());
604  result_ = neg(
605  mul(div(mul(integer(2), exp(neg(mul(self.get_arg(), self.get_arg())))),
606  sqrt(pi)),
607  result_));
608 }
609 
610 void DiffVisitor::bvisit(const Gamma &self)
611 {
612  apply(self.get_arg());
613  result_ = mul(mul(self.rcp_from_this(), polygamma(zero, self.get_arg())),
614  result_);
615 }
616 
617 void DiffVisitor::bvisit(const LogGamma &self)
618 {
619  apply(self.get_arg());
620  result_ = mul(polygamma(zero, self.get_arg()), result_);
621 }
622 
623 void DiffVisitor::bvisit(const UIntPoly &self)
624 {
625  result_ = diff_upoly<UIntPoly, map_uint_mpz>(self, *x);
626 }
627 
628 void DiffVisitor::bvisit(const URatPoly &self)
629 {
630  result_ = diff_upoly<URatPoly, map_uint_mpq>(self, *x);
631 }
632 
633 #ifdef HAVE_SYMENGINE_PIRANHA
634 void DiffVisitor::bvisit(const UIntPolyPiranha &self)
635 {
636  result_ = diff_upoly<UIntPolyPiranha, map_uint_mpz>(self, *x);
637 }
638 void DiffVisitor::bvisit(const URatPolyPiranha &self)
639 {
640  result_ = diff_upoly<URatPolyPiranha, map_uint_mpq>(self, *x);
641 }
642 #endif
643 
644 #ifdef HAVE_SYMENGINE_FLINT
645 void DiffVisitor::bvisit(const UIntPolyFlint &self)
646 {
647  result_ = diff_upolyflint(self, *x);
648 }
649 
650 void DiffVisitor::bvisit(const URatPolyFlint &self)
651 {
652  result_ = diff_upolyflint(self, *x);
653 }
654 #endif
655 void DiffVisitor::bvisit(const MIntPoly &self)
656 {
657  result_ = diff_mpoly(self, x);
658 }
659 
660 void DiffVisitor::bvisit(const MExprPoly &self)
661 {
662  result_ = diff_mpoly(self, x);
663 }
664 
665 void DiffVisitor::bvisit(const UExprPoly &self)
666 {
667  result_ = diff_upoly<UExprPoly, map_int_Expr>(self, *x);
668 }
669 
670 void DiffVisitor::bvisit(const FunctionWrapper &self)
671 {
672  result_ = self.diff_impl(x);
673 }
674 
675 void DiffVisitor::bvisit(const Beta &self)
676 {
677  RCP<const Basic> beta_arg0 = self.get_args()[0];
678  RCP<const Basic> beta_arg1 = self.get_args()[1];
679  apply(beta_arg0);
680  RCP<const Basic> diff_beta_arg0 = result_;
681  apply(beta_arg1);
682  RCP<const Basic> diff_beta_arg1 = result_;
683  result_ = mul(self.rcp_from_this(),
684  add(mul(polygamma(zero, beta_arg0), diff_beta_arg0),
685  sub(mul(polygamma(zero, beta_arg1), diff_beta_arg1),
686  mul(polygamma(zero, add(beta_arg0, beta_arg1)),
687  add(diff_beta_arg0, diff_beta_arg1)))));
688 }
689 
690 void DiffVisitor::bvisit(const Set &self)
691 {
692  throw SymEngineException("Derivative doesn't exist.");
693 }
694 
695 void DiffVisitor::bvisit(const Tuple &self)
696 {
697  throw NotImplementedError("Derivative not implemented");
698 }
699 
700 void DiffVisitor::bvisit(const IdentityMatrix &self)
701 {
702  throw NotImplementedError("Derivative not implemented");
703 }
704 
705 void DiffVisitor::bvisit(const ZeroMatrix &self)
706 {
707  throw NotImplementedError("Derivative not implemented");
708 }
709 
710 void DiffVisitor::bvisit(const MatrixSymbol &self)
711 {
712  throw NotImplementedError("Derivative not implemented");
713 }
714 
715 void DiffVisitor::bvisit(const DiagonalMatrix &self)
716 {
717  throw NotImplementedError("Derivative not implemented");
718 }
719 
720 void DiffVisitor::bvisit(const ImmutableDenseMatrix &self)
721 {
722  throw NotImplementedError("Derivative not implemented");
723 }
724 
725 void DiffVisitor::bvisit(const MatrixAdd &self)
726 {
727  throw NotImplementedError("Derivative not implemented");
728 }
729 
730 void DiffVisitor::bvisit(const HadamardProduct &self)
731 {
732  throw NotImplementedError("Derivative not implemented");
733 }
734 
735 void DiffVisitor::bvisit(const MatrixMul &self)
736 {
737  throw NotImplementedError("Derivative not implemented");
738 }
739 
740 void DiffVisitor::bvisit(const ConjugateMatrix &self)
741 {
742  throw NotImplementedError("Derivative not implemented");
743 }
744 
745 void DiffVisitor::bvisit(const Transpose &self)
746 {
747  throw NotImplementedError("Derivative not implemented");
748 }
749 
750 void DiffVisitor::bvisit(const Trace &self)
751 {
752  throw NotImplementedError("Derivative not implemented");
753 }
754 
755 void DiffVisitor::bvisit(const Boolean &self)
756 {
757  throw SymEngineException("Derivative doesn't exist.");
758 }
759 
760 void DiffVisitor::bvisit(const GaloisField &self)
761 {
762  GaloisFieldDict d;
763  if (self.get_var()->__eq__(*x)) {
764  d = self.get_poly().gf_diff();
765  result_ = GaloisField::from_dict(self.get_var(), std::move(d));
766  } else {
767  result_ = GaloisField::from_dict(self.get_var(), std::move(d));
768  }
769 }
770 
771 void DiffVisitor::bvisit(const Piecewise &self)
772 {
773  PiecewiseVec v = self.get_vec();
774  for (auto &p : v) {
775  apply(p.first);
776  p.first = result_;
777  }
778  result_ = piecewise(std::move(v));
779 }
780 
781 const RCP<const Basic> &DiffVisitor::apply(const Basic &b)
782 {
783  apply(b.rcp_from_this());
784  return result_;
785 }
786 
787 const RCP<const Basic> &DiffVisitor::apply(const RCP<const Basic> &b)
788 {
789  if (not cache) {
790  b->accept(*this);
791  return result_;
792  }
793  auto it = visited.find(b);
794  if (it == visited.end()) {
795  b->accept(*this);
796  insert(visited, b, result_);
797  } else {
798  result_ = it->second;
799  }
800  return result_;
801 }
802 
803 RCP<const Basic> diff(const RCP<const Basic> &arg, const RCP<const Symbol> &x,
804  bool cache)
805 {
806  DiffVisitor v(x, cache);
807  return v.apply(arg);
808 }
809 
810 RCP<const Basic> Basic::diff(const RCP<const Symbol> &x, bool cache) const
811 {
812  return SymEngine::diff(this->rcp_from_this(), x, cache);
813 }
814 
816 // Since SymPy's differentiation makes no sense mathematically, it is
817 // defined separately here for compatibility
818 RCP<const Basic> sdiff(const RCP<const Basic> &arg, const RCP<const Basic> &x,
819  bool cache)
820 {
821  if (is_a<Symbol>(*x)) {
822  return arg->diff(rcp_static_cast<const Symbol>(x), cache);
823  } else {
824  RCP<const Symbol> d = get_dummy(*arg, "x");
825  return ssubs(ssubs(arg, {{x, d}})->diff(d, cache), {{d, x}});
826  }
827 }
828 
829 } // namespace SymEngine
T begin(T... args)
static RCP< const Basic > from_dict(const RCP< const Number > &coef, umap_basic_num &&d)
Create an appropriate instance from dictionary quickly.
Definition: add.cpp:140
static void coef_dict_add_term(const Ptr< RCP< const Number >> &coef, umap_basic_num &d, const RCP< const Number > &c, const RCP< const Basic > &term)
Updates the numerical coefficient and the dictionary.
Definition: add.cpp:261
static void dict_add_term(umap_basic_num &d, const RCP< const Number > &coef, const RCP< const Basic > &t)
Adds a new term to the expression.
Definition: add.cpp:237
static void as_coef_term(const RCP< const Basic > &self, const Ptr< RCP< const Number >> &coef, const Ptr< RCP< const Basic >> &term)
Converts a Basic self into the form of coefficient * term.
Definition: add.cpp:306
RCP< Basic > rcp_from_this()
Get RCP<T> pointer to self (it will cast the pointer to T)
static void as_base_exp(const RCP< const Basic > &self, const Ptr< RCP< const Basic >> &exp, const Ptr< RCP< const Basic >> &base)
Convert to a base and exponent form.
Definition: mul.cpp:320
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
T count(T... args)
T end(T... args)
T erase(T... args)
T find(T... args)
T insert(T... args)
T move(T... args)
Main namespace for SymEngine package.
Definition: add.cpp:19
bool is_a_Number(const Basic &b)
Definition: number.h:130
RCP< const Basic > div(const RCP< const Basic > &a, const RCP< const Basic > &b)
Division.
Definition: mul.cpp:431
std::enable_if< std::is_integral< T >::value, RCP< const Integer > >::type integer(T i)
Definition: integer.h:197
RCP< const Basic > sec(const RCP< const Basic > &arg)
Canonicalize Sec:
Definition: functions.cpp:1202
RCP< const Symbol > symbol(const std::string &name)
inline version to return Symbol
Definition: symbol.h:82
RCP< const Number > mulnum(const RCP< const Number > &self, const RCP< const Number > &other)
Multiply self and other
Definition: number.h:93
RCP< const Basic > polygamma(const RCP< const Basic > &n_, const RCP< const Basic > &x_)
Canonicalize PolyGamma.
Definition: functions.cpp:3397
RCP< const Basic > zeta(const RCP< const Basic > &s, const RCP< const Basic > &a)
Create a new Zeta instance:
Definition: functions.cpp:2800
bool eq(const Basic &a, const Basic &b)
Checks equality for a and b
Definition: basic-inl.h:21
RCP< const Basic > coth(const RCP< const Basic > &arg)
Canonicalize Coth:
Definition: functions.cpp:2333
RCP< const Basic > sech(const RCP< const Basic > &arg)
Canonicalize Sech:
Definition: functions.cpp:2251
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 > tan(const RCP< const Basic > &arg)
Canonicalize Tan:
Definition: functions.cpp:1007
RCP< const Basic > cosh(const RCP< const Basic > &arg)
Canonicalize Cosh:
Definition: functions.cpp:2212
RCP< const Basic > tanh(const RCP< const Basic > &arg)
Canonicalize Tanh:
Definition: functions.cpp:2290
RCP< const Basic > cot(const RCP< const Basic > &arg)
Canonicalize Cot:
Definition: functions.cpp:1073
RCP< const Basic > mul(const RCP< const Basic > &a, const RCP< const Basic > &b)
Multiplication.
Definition: mul.cpp:352
void insert(T1 &m, const T2 &first, const T3 &second)
Definition: dict.h:83
RCP< const Basic > cos(const RCP< const Basic > &arg)
Canonicalize Cos:
Definition: functions.cpp:942
RCP< const Basic > log(const RCP< const Basic > &arg)
Returns the Natural Logarithm from argument arg
Definition: functions.cpp:1774
RCP< const Basic > csc(const RCP< const Basic > &arg)
Canonicalize Csc:
Definition: functions.cpp:1138
int factor(const Ptr< RCP< const Integer >> &f, const Integer &n, double B1)
Definition: ntheory.cpp:371
RCP< const Basic > add(const RCP< const Basic > &a, const RCP< const Basic > &b)
Adds two objects (safely).
Definition: add.cpp:425
RCP< const Basic > csch(const RCP< const Basic > &arg)
Canonicalize Csch:
Definition: functions.cpp:2169
bool neq(const Basic &a, const Basic &b)
Checks inequality for a and b
Definition: basic-inl.h:29
RCP< const Basic > lambertw(const RCP< const Basic > &arg)
Create a new LambertW instance:
Definition: functions.cpp:1846
RCP< const Basic > sinh(const RCP< const Basic > &arg)
Canonicalize Sinh:
Definition: functions.cpp:2127
RCP< const Basic > ssubs(const RCP< const Basic > &x, const map_basic_basic &subs_dict, bool cache=true)
SymPy compatible subs.
Definition: subs.h:559
RCP< const Basic > neg(const RCP< const Basic > &a)
Negation.
Definition: mul.cpp:443
RCP< const Basic > sin(const RCP< const Basic > &arg)
Canonicalize Sin:
Definition: functions.cpp:874
T pow(T... args)
T sqrt(T... args)
T str(T... args)