2 #include <symengine/subs.h>
3 #include <symengine/symengine_casts.h>
9 extern RCP<const Basic> i2;
13 static inline RCP<const Basic> fdiff(
const T &
self, RCP<const Symbol> x,
16 RCP<const Basic> diff = zero;
20 vec_basic v =
self.get_args();
21 vec_basic vdiff(v.size());
24 for (
unsigned i = 0; i < v.size(); i++) {
25 vdiff[i] = visitor.apply(v[i]);
26 if (
neq(*vdiff[i], *zero)) {
35 for (
unsigned i = 0; i < v.size(); i++) {
36 if (
eq(*vdiff[i], *zero))
38 know_deriv = fdiff(outArg(ret),
self, i);
40 diff =
add(diff,
mul(ret, vdiff[i]));
42 if (count == 1 and
eq(*v[i], *x)) {
43 return Derivative::create(
self.rcp_from_this(), {x});
45 vec_basic new_args = v;
48 new_args[i] = get_dummy(
self,
"xi_" + stm.
str());
50 insert(m, new_args[i], v[i]);
51 diff =
add(diff,
mul(vdiff[i],
53 Derivative::create(
self.create(new_args),
61 static bool fdiff(
const Ptr<RCP<const Basic>> &ret,
const Zeta &
self,
65 *ret =
mul(
mul(minus_one,
self.get_arg1()),
66 zeta(
add(
self.get_arg1(), one),
self.get_arg2()));
73 static bool fdiff(
const Ptr<RCP<const Basic>> &ret,
const UpperGamma &
self,
77 *ret =
mul(
mul(
pow(
self.get_arg2(),
sub(
self.get_arg1(), one)),
78 exp(
neg(
self.get_arg2()))),
86 static bool fdiff(
const Ptr<RCP<const Basic>> &ret,
const LowerGamma &
self,
90 *ret =
mul(
pow(
self.get_arg2(),
sub(
self.get_arg1(), one)),
91 exp(
neg(
self.get_arg2())));
98 static bool fdiff(
const Ptr<RCP<const Basic>> &ret,
const PolyGamma &
self,
102 *ret =
polygamma(
add(
self.get_arg1(), one),
self.get_arg2());
109 static bool fdiff(
const Ptr<RCP<const Basic>> &ret,
const Function &
self,
116 template <
typename P>
117 static inline RCP<const Basic> diff_upolyflint(
const P &
self,
const Symbol &x)
119 if (
self.get_var()->__eq__(x)) {
120 return P::from_container(
self.get_var(),
self.get_poly().derivative());
122 return P::from_dict(
self.get_var(), {{}});
126 template <
typename Poly,
typename Dict>
127 static inline RCP<const Basic> diff_upoly(
const Poly &
self,
const Symbol &x)
129 if (
self.get_var()->__eq__(x)) {
131 for (
auto it =
self.
begin(); it !=
self.end(); ++it) {
133 d[it->first - 1] = it->second * it->first;
135 return Poly::from_dict(
self.get_var(),
std::move(d));
137 return Poly::from_dict(
self.get_var(), {{}});
141 template <
typename Container,
typename Poly>
142 static RCP<const Basic> diff_mpoly(
const MSymEnginePoly<Container, Poly> &
self,
143 const RCP<const Symbol> &x)
145 using Dict =
typename Container::dict_type;
146 using Vec =
typename Container::vec_type;
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)) {
156 for (
auto bucket :
self.get_poly().dict_) {
157 if (bucket.first[index] != 0) {
158 Vec v = bucket.first;
160 dict.insert({v, bucket.second * bucket.first[index]});
164 v.
insert(v.begin(),
self.get_vars().begin(),
self.get_vars().end());
165 return Poly::from_dict(v,
std::move(dict));
168 vs.
insert(vs.begin(),
self.get_vars().begin(),
self.get_vars().end());
169 return Poly::from_dict(vs, {{}});
173 #ifndef debug_methods
175 void DiffVisitor::bvisit(
const Basic &
self)
177 result_ = Derivative::create(
self.rcp_from_this(), {x});
186 #define DIFF0(CLASS) \
187 void DiffVisitor::bvisit(const CLASS &self) \
189 result_ = Derivative::create(self.rcp_from_this(), {x}); \
192 DIFF0(UnivariateSeries)
197 void DiffVisitor::bvisit(
const Number &
self)
202 void DiffVisitor::bvisit(
const Constant &
self)
207 void DiffVisitor::bvisit(
const Symbol &
self)
209 if (x->get_name() ==
self.get_name()) {
216 void DiffVisitor::bvisit(
const Log &
self)
218 apply(
self.get_arg());
219 result_ =
mul(
div(one,
self.get_arg()), result_);
222 void DiffVisitor::bvisit(
const Abs &
self)
224 apply(
self.get_arg());
225 if (
eq(*result_, *zero)) {
228 result_ = Derivative::create(
self.rcp_from_this(), {x});
232 void DiffVisitor::bvisit(
const Zeta &
self)
234 result_ = fdiff(
self, x, *
this);
237 void DiffVisitor::bvisit(
const LowerGamma &
self)
239 result_ = fdiff(
self, x, *
this);
241 void DiffVisitor::bvisit(
const UpperGamma &
self)
243 result_ = fdiff(
self, x, *
this);
245 void DiffVisitor::bvisit(
const PolyGamma &
self)
247 result_ = fdiff(
self, x, *
this);
249 void DiffVisitor::bvisit(
const UnevaluatedExpr &
self)
251 result_ = Derivative::create(
self.rcp_from_this(), {x});
253 void DiffVisitor::bvisit(
const TwoArgFunction &
self)
255 result_ = fdiff(
self, x, *
this);
258 void DiffVisitor::bvisit(
const ASech &
self)
260 apply(
self.get_arg());
266 void DiffVisitor::bvisit(
const ACoth &
self)
268 apply(
self.get_arg());
269 result_ =
mul(
div(one,
sub(one,
pow(
self.get_arg(), i2))), result_);
272 void DiffVisitor::bvisit(
const ATanh &
self)
274 apply(
self.get_arg());
275 result_ =
mul(
div(one,
sub(one,
pow(
self.get_arg(), i2))), result_);
278 void DiffVisitor::bvisit(
const ACosh &
self)
280 apply(
self.get_arg());
284 void DiffVisitor::bvisit(
const ACsch &
self)
286 apply(
self.get_arg());
289 pow(
self.get_arg(), i2))),
293 void DiffVisitor::bvisit(
const ASinh &
self)
295 apply(
self.get_arg());
299 void DiffVisitor::bvisit(
const Coth &
self)
301 apply(
self.get_arg());
302 result_ =
mul(
div(minus_one,
pow(
sinh(
self.get_arg()), i2)), result_);
305 void DiffVisitor::bvisit(
const Tanh &
self)
307 apply(
self.get_arg());
308 result_ =
mul(
sub(one,
pow(
tanh(
self.get_arg()), i2)), result_);
311 void DiffVisitor::bvisit(
const Sech &
self)
313 apply(
self.get_arg());
319 void DiffVisitor::bvisit(
const Cosh &
self)
321 apply(
self.get_arg());
322 result_ =
mul(
sinh(
self.get_arg()), result_);
325 void DiffVisitor::bvisit(
const Csch &
self)
327 apply(
self.get_arg());
333 void DiffVisitor::bvisit(
const Sinh &
self)
335 apply(
self.get_arg());
336 result_ =
mul(
cosh(
self.get_arg()), result_);
339 void DiffVisitor::bvisit(
const Subs &
self)
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());
346 for (
const auto &p :
self.get_dict()) {
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())));
355 result_ = Derivative::create(
self.rcp_from_this(), {x});
363 void DiffVisitor::bvisit(
const Derivative &
self)
365 apply(
self.get_arg());
366 RCP<const Basic> ret = result_;
367 if (
eq(*ret, *zero)) {
370 multiset_basic t =
self.get_symbols();
376 result_ = Derivative::create(
self.get_arg(), t);
381 if (is_a<Derivative>(*ret)
382 &&
eq(*down_cast<const Derivative &>(*ret).get_arg(),
385 result_ = Derivative::create(
self.get_arg(), t);
389 ret = diff(ret, rcp_static_cast<const Symbol>(p));
394 static inline RCP<const Symbol> get_dummy(
const Basic &b,
std::string name)
400 }
while (has_symbol(b, *s));
404 void DiffVisitor::bvisit(
const OneArgFunction &
self)
406 result_ = fdiff(
self, x, *
this);
409 void DiffVisitor::bvisit(
const MultiArgFunction &
self)
411 result_ = fdiff(
self, x, *
this);
414 void DiffVisitor::bvisit(
const LambertW &
self)
418 apply(
self.get_arg());
419 RCP<const Basic> lambertw_val =
lambertw(
self.get_arg());
421 =
mul(
div(lambertw_val,
mul(
self.get_arg(),
add(lambertw_val, one))),
425 void DiffVisitor::bvisit(
const Add &
self)
428 RCP<const Number> coef = zero, coef2;
430 for (
auto &p :
self.get_dict()) {
432 RCP<const Basic> term = result_;
433 if (is_a<Integer>(*term)
434 && down_cast<const Integer &>(*term).is_zero()) {
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())
442 iaddnum(outArg(coef),
443 mulnum(p.second, down_cast<const Add &>(*term).get_coef()));
452 void DiffVisitor::bvisit(
const Mul &
self)
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())
463 map_basic_basic d =
self.get_dict();
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);
474 RCP<const Basic>
exp, t;
476 Mul::dict_add_term_new(outArg(coef), d,
exp, t);
479 iaddnum(outArg(overall_coef), coef);
488 void DiffVisitor::bvisit(
const Pow &
self)
491 apply(
self.get_base());
493 mul(
self.get_exp(),
pow(
self.get_base(),
sub(
self.get_exp(), one))),
496 apply(
mul(
self.get_exp(),
log(
self.get_base())));
497 result_ =
mul(
self.rcp_from_this(), result_);
501 void DiffVisitor::bvisit(
const Sin &
self)
503 apply(
self.get_arg());
504 result_ =
mul(
cos(
self.get_arg()), result_);
507 void DiffVisitor::bvisit(
const Cos &
self)
509 apply(
self.get_arg());
510 result_ =
mul(
mul(minus_one,
sin(
self.get_arg())), result_);
513 void DiffVisitor::bvisit(
const Tan &
self)
515 apply(
self.get_arg());
516 RCP<const Integer> two =
integer(2);
517 result_ =
mul(
add(one,
pow(
tan(
self.get_arg()), two)), result_);
520 void DiffVisitor::bvisit(
const Cot &
self)
522 apply(
self.get_arg());
523 RCP<const Integer> two =
integer(2);
525 =
mul(
mul(
add(one,
pow(
cot(
self.get_arg()), two)), minus_one), result_);
528 void DiffVisitor::bvisit(
const Csc &
self)
530 apply(
self.get_arg());
531 result_ =
mul(
mul(
mul(
cot(
self.get_arg()),
csc(
self.get_arg())), minus_one),
535 void DiffVisitor::bvisit(
const Sec &
self)
537 apply(
self.get_arg());
538 result_ =
mul(
mul(
tan(
self.get_arg()),
sec(
self.get_arg())), result_);
541 void DiffVisitor::bvisit(
const ASin &
self)
543 apply(
self.get_arg());
547 void DiffVisitor::bvisit(
const ACos &
self)
549 apply(
self.get_arg());
554 void DiffVisitor::bvisit(
const ASec &
self)
556 apply(
self.get_arg());
563 void DiffVisitor::bvisit(
const ACsc &
self)
565 apply(
self.get_arg());
567 div(minus_one,
mul(
pow(
self.get_arg(), i2),
572 void DiffVisitor::bvisit(
const ATan &
self)
574 apply(
self.get_arg());
575 result_ =
mul(
div(one,
add(one,
pow(
self.get_arg(), i2))), result_);
578 void DiffVisitor::bvisit(
const ACot &
self)
580 apply(
self.get_arg());
581 result_ =
mul(
div(minus_one,
add(one,
pow(
self.get_arg(), i2))), result_);
584 void DiffVisitor::bvisit(
const ATan2 &
self)
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))),
592 void DiffVisitor::bvisit(
const Erf &
self)
594 apply(
self.get_arg());
601 void DiffVisitor::bvisit(
const Erfc &
self)
603 apply(
self.get_arg());
610 void DiffVisitor::bvisit(
const Gamma &
self)
612 apply(
self.get_arg());
613 result_ =
mul(
mul(
self.rcp_from_this(),
polygamma(zero,
self.get_arg())),
617 void DiffVisitor::bvisit(
const LogGamma &
self)
619 apply(
self.get_arg());
620 result_ =
mul(
polygamma(zero,
self.get_arg()), result_);
623 void DiffVisitor::bvisit(
const UIntPoly &
self)
625 result_ = diff_upoly<UIntPoly, map_uint_mpz>(
self, *x);
628 void DiffVisitor::bvisit(
const URatPoly &
self)
630 result_ = diff_upoly<URatPoly, map_uint_mpq>(
self, *x);
633 #ifdef HAVE_SYMENGINE_PIRANHA
634 void DiffVisitor::bvisit(
const UIntPolyPiranha &
self)
636 result_ = diff_upoly<UIntPolyPiranha, map_uint_mpz>(
self, *x);
638 void DiffVisitor::bvisit(
const URatPolyPiranha &
self)
640 result_ = diff_upoly<URatPolyPiranha, map_uint_mpq>(
self, *x);
644 #ifdef HAVE_SYMENGINE_FLINT
645 void DiffVisitor::bvisit(
const UIntPolyFlint &
self)
647 result_ = diff_upolyflint(
self, *x);
650 void DiffVisitor::bvisit(
const URatPolyFlint &
self)
652 result_ = diff_upolyflint(
self, *x);
655 void DiffVisitor::bvisit(
const MIntPoly &
self)
657 result_ = diff_mpoly(
self, x);
660 void DiffVisitor::bvisit(
const MExprPoly &
self)
662 result_ = diff_mpoly(
self, x);
665 void DiffVisitor::bvisit(
const UExprPoly &
self)
667 result_ = diff_upoly<UExprPoly, map_int_Expr>(
self, *x);
670 void DiffVisitor::bvisit(
const FunctionWrapper &
self)
672 result_ =
self.diff_impl(x);
675 void DiffVisitor::bvisit(
const Beta &
self)
677 RCP<const Basic> beta_arg0 =
self.get_args()[0];
678 RCP<const Basic> beta_arg1 =
self.get_args()[1];
680 RCP<const Basic> diff_beta_arg0 = result_;
682 RCP<const Basic> diff_beta_arg1 = result_;
683 result_ =
mul(
self.rcp_from_this(),
687 add(diff_beta_arg0, diff_beta_arg1)))));
690 void DiffVisitor::bvisit(
const Set &
self)
692 throw SymEngineException(
"Derivative doesn't exist.");
695 void DiffVisitor::bvisit(
const Tuple &
self)
697 throw NotImplementedError(
"Derivative not implemented");
700 void DiffVisitor::bvisit(
const IdentityMatrix &
self)
702 throw NotImplementedError(
"Derivative not implemented");
705 void DiffVisitor::bvisit(
const ZeroMatrix &
self)
707 throw NotImplementedError(
"Derivative not implemented");
710 void DiffVisitor::bvisit(
const MatrixSymbol &
self)
712 throw NotImplementedError(
"Derivative not implemented");
715 void DiffVisitor::bvisit(
const DiagonalMatrix &
self)
717 throw NotImplementedError(
"Derivative not implemented");
720 void DiffVisitor::bvisit(
const ImmutableDenseMatrix &
self)
722 throw NotImplementedError(
"Derivative not implemented");
725 void DiffVisitor::bvisit(
const MatrixAdd &
self)
727 throw NotImplementedError(
"Derivative not implemented");
730 void DiffVisitor::bvisit(
const HadamardProduct &
self)
732 throw NotImplementedError(
"Derivative not implemented");
735 void DiffVisitor::bvisit(
const MatrixMul &
self)
737 throw NotImplementedError(
"Derivative not implemented");
740 void DiffVisitor::bvisit(
const ConjugateMatrix &
self)
742 throw NotImplementedError(
"Derivative not implemented");
745 void DiffVisitor::bvisit(
const Transpose &
self)
747 throw NotImplementedError(
"Derivative not implemented");
750 void DiffVisitor::bvisit(
const Trace &
self)
752 throw NotImplementedError(
"Derivative not implemented");
755 void DiffVisitor::bvisit(
const Boolean &
self)
757 throw SymEngineException(
"Derivative doesn't exist.");
760 void DiffVisitor::bvisit(
const GaloisField &
self)
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));
767 result_ = GaloisField::from_dict(
self.get_var(),
std::move(d));
771 void DiffVisitor::bvisit(
const Piecewise &
self)
773 PiecewiseVec v =
self.get_vec();
781 const RCP<const Basic> &DiffVisitor::apply(
const Basic &b)
783 apply(b.rcp_from_this());
787 const RCP<const Basic> &DiffVisitor::apply(
const RCP<const Basic> &b)
793 auto it = visited.
find(b);
794 if (it == visited.
end()) {
796 insert(visited, b, result_);
798 result_ = it->second;
803 RCP<const Basic> diff(
const RCP<const Basic> &arg,
const RCP<const Symbol> &x,
810 RCP<const Basic> Basic::diff(
const RCP<const Symbol> &x,
bool cache)
const
818 RCP<const Basic> sdiff(
const RCP<const Basic> &arg,
const RCP<const Basic> &x,
821 if (is_a<Symbol>(*x)) {
822 return arg->diff(rcp_static_cast<const Symbol>(x), cache);
824 RCP<const Symbol> d = get_dummy(*arg,
"x");
825 return ssubs(
ssubs(arg, {{x, d}})->diff(d, cache), {{d, x}});
static RCP< const Basic > from_dict(const RCP< const Number > &coef, umap_basic_num &&d)
Create an appropriate instance from dictionary quickly.
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.
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.
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.
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.
static RCP< const Basic > from_dict(const RCP< const Number > &coef, map_basic_basic &&d)
Create a Mul from a dict.
Main namespace for SymEngine package.
bool is_a_Number(const Basic &b)
RCP< const Basic > div(const RCP< const Basic > &a, const RCP< const Basic > &b)
Division.
std::enable_if< std::is_integral< T >::value, RCP< const Integer > >::type integer(T i)
RCP< const Basic > sec(const RCP< const Basic > &arg)
Canonicalize Sec:
RCP< const Symbol > symbol(const std::string &name)
inline version to return Symbol
RCP< const Number > mulnum(const RCP< const Number > &self, const RCP< const Number > &other)
Multiply self and other
RCP< const Basic > polygamma(const RCP< const Basic > &n_, const RCP< const Basic > &x_)
Canonicalize PolyGamma.
RCP< const Basic > zeta(const RCP< const Basic > &s, const RCP< const Basic > &a)
Create a new Zeta instance:
bool eq(const Basic &a, const Basic &b)
Checks equality for a and b
RCP< const Basic > coth(const RCP< const Basic > &arg)
Canonicalize Coth:
RCP< const Basic > sech(const RCP< const Basic > &arg)
Canonicalize Sech:
RCP< const Basic > sub(const RCP< const Basic > &a, const RCP< const Basic > &b)
Substracts b from a.
RCP< const Basic > exp(const RCP< const Basic > &x)
Returns the natural exponential function E**x = pow(E, x)
RCP< const Basic > tan(const RCP< const Basic > &arg)
Canonicalize Tan:
RCP< const Basic > cosh(const RCP< const Basic > &arg)
Canonicalize Cosh:
RCP< const Basic > tanh(const RCP< const Basic > &arg)
Canonicalize Tanh:
RCP< const Basic > cot(const RCP< const Basic > &arg)
Canonicalize Cot:
RCP< const Basic > mul(const RCP< const Basic > &a, const RCP< const Basic > &b)
Multiplication.
void insert(T1 &m, const T2 &first, const T3 &second)
RCP< const Basic > cos(const RCP< const Basic > &arg)
Canonicalize Cos:
RCP< const Basic > log(const RCP< const Basic > &arg)
Returns the Natural Logarithm from argument arg
RCP< const Basic > csc(const RCP< const Basic > &arg)
Canonicalize Csc:
int factor(const Ptr< RCP< const Integer >> &f, const Integer &n, double B1)
RCP< const Basic > add(const RCP< const Basic > &a, const RCP< const Basic > &b)
Adds two objects (safely).
RCP< const Basic > csch(const RCP< const Basic > &arg)
Canonicalize Csch:
bool neq(const Basic &a, const Basic &b)
Checks inequality for a and b
RCP< const Basic > lambertw(const RCP< const Basic > &arg)
Create a new LambertW instance:
RCP< const Basic > sinh(const RCP< const Basic > &arg)
Canonicalize Sinh:
RCP< const Basic > ssubs(const RCP< const Basic > &x, const map_basic_basic &subs_dict, bool cache=true)
SymPy compatible subs.
RCP< const Basic > neg(const RCP< const Basic > &a)
Negation.
RCP< const Basic > sin(const RCP< const Basic > &arg)
Canonicalize Sin: