6 #ifndef SYMENGINE_SERIES_H
7 #define SYMENGINE_SERIES_H
12 #include <symengine/symengine_exception.h>
20 virtual RCP<const Basic> as_basic()
const = 0;
22 virtual RCP<const Basic> get_coeff(
int)
const = 0;
23 virtual unsigned get_degree()
const = 0;
27 template <
typename Poly,
typename Coeff,
typename Series>
33 const unsigned degree_;
37 : p_(
std::move(p)), var_(var), degree_(degree)
40 inline unsigned get_degree()
const override
50 inline const Poly &get_poly()
const
87 return (is_a<Series>(o) and var_ == down_cast<const Series &>(o).var_
88 and p_ == down_cast<const Series &>(o).p_
89 and degree_ == down_cast<const Series &>(o).degree_);
92 RCP<const Number>
add(
const Number &other)
const override
94 if (is_a<Series>(other)) {
95 const Series &o = down_cast<const Series &>(other);
96 auto deg =
std::min(degree_, o.degree_);
98 throw NotImplementedError(
99 "Multivariate Series not implemented");
101 return make_rcp<Series>(Poly(p_ + o.p_), var_, deg);
102 }
else if (other.get_type_code() < Series::type_code_id) {
103 Poly p = Series::series(other.
rcp_from_this(), var_, degree_)->p_;
104 return make_rcp<Series>(Poly(p_ + p), var_, degree_);
106 return other.
add(*
this);
110 RCP<const Number>
mul(
const Number &other)
const override
112 if (is_a<Series>(other)) {
113 const Series &o = down_cast<const Series &>(other);
114 auto deg =
std::min(degree_, o.degree_);
115 if (var_ != o.var_) {
116 throw NotImplementedError(
117 "Multivariate Series not implemented");
119 return make_rcp<Series>(Series::mul(p_, o.p_, deg), var_, deg);
120 }
else if (other.get_type_code() < Series::type_code_id) {
121 Poly p = Series::series(other.
rcp_from_this(), var_, degree_)->p_;
122 return make_rcp<Series>(Series::mul(p_, p, degree_), var_, degree_);
124 return other.
mul(*
this);
128 RCP<const Number>
pow(
const Number &other)
const override
132 if (is_a<Series>(other)) {
133 const Series &o = down_cast<const Series &>(other);
135 if (var_ != o.var_) {
136 throw NotImplementedError(
137 "Multivariate Series not implemented");
140 }
else if (is_a<Integer>(other)) {
145 down_cast<const Integer &>(other).
neg()->as_int())),
147 p = Series::series_invert(p, Series::var(var_), deg);
148 return make_rcp<Series>(p, var_, deg);
152 numeric_cast<int>((down_cast<const Integer &>(other).as_int())),
154 return make_rcp<Series>(p, var_, deg);
155 }
else if (other.get_type_code() < Series::type_code_id) {
156 p = Series::series(other.
rcp_from_this(), var_, degree_)->p_;
158 return other.rpow(*
this);
160 p = Series::series_exp(
161 Poly(p * Series::series_log(p_, Series::var(var_), deg)),
162 Series::var(var_), deg);
163 return make_rcp<Series>(p, var_, deg);
166 RCP<const Number> rpow(
const Number &other)
const override
168 if (other.get_type_code() < Series::type_code_id) {
169 Poly p = Series::series(other.
rcp_from_this(), var_, degree_)->p_;
170 p = Series::series_exp(
171 Poly(p_ * Series::series_log(p, Series::var(var_), degree_)),
172 Series::var(var_), degree_);
173 return make_rcp<Series>(p, var_, degree_);
175 throw SymEngineException(
"Unknown type");
182 if (not steps.
empty()) {
183 if (*(steps.
rbegin()) == prec)
189 unsigned int tprec = prec;
191 tprec = 2 + tprec / 2;
199 static inline Poly series_invert(
const Poly &s,
const Poly &var,
203 throw DivisionByZeroError(
204 "Series::series_invert: Division By Zero");
207 const int ldeg = Series::ldegree(s);
208 const Coeff co = Series::find_cf(s, var, ldeg);
209 Poly p(1 / co), ss = s;
211 ss = s * Series::pow(var, -ldeg, prec);
213 auto steps = step_list(prec);
214 for (
const auto step : steps) {
215 p = Series::mul(2 - Series::mul(p, ss, step), p, step);
218 return p * Series::pow(var, -ldeg, prec);
224 static inline Poly series_reverse(
const Poly &s,
const Poly &var,
227 const Coeff co = Series::find_cf(s, var, 0);
229 throw SymEngineException(
"reversion of series with constant term");
230 const Coeff a = Series::find_cf(s, var, 1);
232 throw SymEngineException(
233 "reversion of series with zero term of degree one");
236 for (
unsigned int i = 2; i < prec; i++) {
237 Poly sp = Series::subs(s, var, r, i + 1);
238 r -= Series::pow(var, i, i + 1) * Series::find_cf(sp, var, i) / a;
243 static inline Poly series_nthroot(
const Poly &s,
int n,
const Poly &var,
251 return Series::series_invert(s, var, prec);
253 const int ldeg = Series::ldegree(s);
255 throw NotImplementedError(
"Puiseux series not implemented.");
259 ss = s * Series::pow(var, -ldeg, prec);
261 Coeff ct = Series::find_cf(ss, var, 0);
268 Coeff ctroot = Series::root(ct, n);
269 Poly res_p(1), sn = ss / ct;
270 auto steps = step_list(prec);
271 for (
const auto step : steps) {
272 Poly t = Series::mul(Series::pow(res_p, n + 1, step), sn, step);
273 res_p += (res_p - t) / n;
276 res_p *= Series::pow(var, ldeg / n, prec);
279 return res_p / ctroot;
281 return Series::series_invert(res_p, var, prec) * ctroot;
294 Poly monom(var), vsquare(var * var);
295 for (
unsigned int i = 1; i < prec; i += 2,
sign *= -1) {
296 res_p += monom * (Coeff(
sign) / Coeff(i));
301 const Coeff c(Series::find_cf(s, var, 0));
302 const Poly p(Series::pow(s, 2, prec - 1) + 1);
303 res_p = Series::mul(Series::diff(s, var),
304 Series::series_invert(p, var, prec - 1), prec - 1);
308 return Series::integrate(res_p, var);
310 return Series::integrate(res_p, var) + Series::atan(c);
314 static inline Poly series_tan(
const Poly &s,
const Poly &var,
317 Poly res_p(0), ss = s;
318 const Coeff c(Series::find_cf(s, var, 0));
336 auto steps = step_list(prec);
337 for (
const auto step : steps) {
338 Poly t = Series::pow(res_p, 2, step) + 1;
339 res_p += Series::mul(ss - Series::series_atan(res_p, var, step), t,
347 res_p + Series::tan(c),
348 Series::series_invert(1 + res_p * (-Series::tan(c)), var, prec),
353 static inline Poly series_cot(
const Poly &s,
const Poly &var,
356 return Series::series_invert(Series::series_tan(s, var, prec), var,
360 static inline Poly _series_sin(
const Poly &s,
unsigned int prec)
362 Poly res_p(0), monom(s);
363 Poly ssquare = Series::mul(s, s, prec);
365 for (
unsigned int i = 0; i < prec / 2; i++) {
366 const int j = 2 * i + 1;
370 res_p += Series::mul(monom, Poly(prod), prec);
371 monom = Series::mul(monom, ssquare, prec);
376 static inline Poly series_sin(
const Poly &s,
const Poly &var,
379 const Coeff c(Series::find_cf(s, var, 0));
381 const Poly t = s - c;
382 return Series::cos(c) * _series_sin(t, prec)
383 + Series::sin(c) * _series_cos(t, prec);
385 return _series_sin(s, prec);
406 static inline Poly series_csc(
const Poly &s,
const Poly &var,
409 return Series::series_invert(Series::series_sin(s, var, prec), var,
413 static inline Poly series_asin(
const Poly &s,
const Poly &var,
416 const Coeff c(Series::find_cf(s, var, 0));
419 const Poly t(1 - Series::pow(s, 2, prec - 1));
420 const Poly res_p(Series::integrate(
421 Series::diff(s, var) * Series::series_nthroot(t, -2, var, prec - 1),
425 return res_p + Series::asin(c);
431 static inline Poly series_acos(
const Poly &s,
const Poly &var,
434 const Coeff c(Series::find_cf(s, var, 0));
435 return Series::acos(c) - series_asin(s - c, var, prec);
438 static inline Poly _series_cos(
const Poly &s,
unsigned int prec)
441 Poly ssquare = Series::mul(s, s, prec);
444 for (
unsigned int i = 1; i <= prec / 2; i++) {
449 res_p += Series::mul(monom, Poly(prod), prec);
450 monom = Series::mul(monom, ssquare, prec);
455 static inline Poly series_cos(
const Poly &s,
const Poly &var,
458 const Coeff c(Series::find_cf(s, var, 0));
460 const Poly t = s - c;
461 return Series::cos(c) * _series_cos(t, prec)
462 - Series::sin(c) * _series_sin(t, prec);
464 return _series_cos(s, prec);
486 static inline Poly series_sec(
const Poly &s,
const Poly &var,
489 return Series::series_invert(Series::series_cos(s, var, prec), var,
493 static inline Poly
series_log(
const Poly &s,
const Poly &var,
502 for (
unsigned int i = 1; i < prec; i++) {
503 res_p += monom * Coeff(((i % 2) == 0) ? -1 : 1) / Coeff(i);
509 const Coeff c(Series::find_cf(s, var, 0));
510 res_p = Series::mul(Series::diff(s, var),
511 Series::series_invert(s, var, prec), prec - 1);
512 res_p = Series::integrate(res_p, var);
515 res_p += Series::log(c);
520 static inline Poly
series_exp(
const Poly &s,
const Poly &var,
531 for (
unsigned int i = 1; i < prec; i++) {
533 res_p += monom * coef;
539 const Coeff c(Series::find_cf(s, var, 0));
545 auto steps = step_list(prec);
546 for (
const auto step : steps) {
547 res_p = Series::mul(res_p, t - Series::series_log(res_p, var, step),
551 return res_p * Series::exp(c);
557 static inline Poly series_lambertw(
const Poly &s,
const Poly &var,
560 if (Series::find_cf(s, var, 0) != 0)
561 throw NotImplementedError(
"lambertw(const) not Implemented");
565 auto steps = step_list(prec);
566 for (
const auto step : steps) {
567 const Poly e(Series::series_exp(p1, var, step));
568 const Poly p2(Series::mul(e, p1, step) - s);
569 const Poly p3(Series::series_invert(Series::mul(e, (p1 + 1), step),
571 p1 -= Series::mul(p2, p3, step);
576 static inline Poly series_sinh(
const Poly &s,
const Poly &var,
579 const Coeff c(Series::find_cf(s, var, 0));
580 const Poly p1(Series::series_exp(s - c, var, prec));
581 const Poly p2(Series::series_invert(p1, var, prec));
584 return (p1 - p2) / 2;
586 return Series::cosh(c) * (p1 - p2) / 2
587 + Series::sinh(c) * (p1 + p2) / 2;
591 static inline Poly series_cosh(
const Poly &s,
const Poly &var,
594 const Coeff c(Series::find_cf(s, var, 0));
595 const Poly p1(Series::series_exp(s - c, var, prec));
596 const Poly p2(Series::series_invert(p1, var, prec));
599 return (p1 + p2) / 2;
601 return Series::cosh(c) * (p1 + p2) / 2
602 + Series::sinh(c) * (p1 - p2) / 2;
606 static inline Poly series_atanh(
const Poly &s,
const Poly &var,
609 const Coeff c(Series::find_cf(s, var, 0));
610 const Poly p(1 - Series::pow(s, 2, prec - 1));
611 const Poly res_p(Series::mul(Series::diff(s, var),
612 Series::series_invert(p, var, prec - 1),
616 return Series::integrate(res_p, var);
618 return Series::integrate(res_p, var) + Series::atanh(c);
622 static inline Poly series_asinh(
const Poly &s,
const Poly &var,
625 const Coeff c(Series::find_cf(s, var, 0));
627 const Poly p(Series::series_nthroot(Series::pow(s, 2, prec - 1) + 1, 2,
629 const Poly res_p(Series::diff(s, var)
630 * Series::series_invert(p, var, prec - 1));
633 return Series::integrate(res_p, var);
635 return Series::integrate(res_p, var) + Series::asinh(c);
639 static inline Poly series_tanh(
const Poly &s,
const Poly &var,
642 const Coeff c(Series::find_cf(s, var, 0));
648 auto steps = step_list(prec);
649 for (
const auto step : steps) {
650 const Poly p(s_ - Series::series_atanh(res_p, var, step));
651 res_p += Series::mul(-p, Series::pow(res_p, 2, step) - 1, step);
654 return (res_p + Series::tanh(c))
655 * Series::series_invert(1 + Series::tanh(c) * res_p, var,
662 static inline Coeff sin(
const Coeff &c)
664 throw NotImplementedError(
"sin(const) not implemented");
666 static inline Coeff cos(
const Coeff &c)
668 throw NotImplementedError(
"cos(const) not implemented");
670 static inline Coeff tan(
const Coeff &c)
672 throw NotImplementedError(
"tan(const) not implemented");
674 static inline Coeff asin(
const Coeff &c)
676 throw NotImplementedError(
"asin(const) not implemented");
678 static inline Coeff acos(
const Coeff &c)
680 throw NotImplementedError(
"acos(const) not implemented");
682 static inline Coeff atan(
const Coeff &c)
684 throw NotImplementedError(
"atan(const) not implemented");
686 static inline Coeff sinh(
const Coeff &c)
688 throw NotImplementedError(
"sinh(const) not implemented");
690 static inline Coeff cosh(
const Coeff &c)
692 throw NotImplementedError(
"cosh(const) not implemented");
694 static inline Coeff tanh(
const Coeff &c)
696 throw NotImplementedError(
"tanh(const) not implemented");
698 static inline Coeff asinh(
const Coeff &c)
700 throw NotImplementedError(
"asinh(const) not implemented");
702 static inline Coeff atanh(
const Coeff &c)
704 throw NotImplementedError(
"atanh(const) not implemented");
706 static inline Coeff exp(
const Coeff &c)
708 throw NotImplementedError(
"exp(const) not implemented");
710 static inline Coeff log(
const Coeff &c)
712 throw NotImplementedError(
"log(const) not implemented");
716 RCP<const SeriesCoeffInterface> series(
const RCP<const Basic> &ex,
717 const RCP<const Symbol> &var,
720 RCP<const SeriesCoeffInterface> series_invfunc(
const RCP<const Basic> &ex,
721 const RCP<const Symbol> &var,
The lowest unit of symbolic representation.
RCP< T > rcp_from_this()
Get RCP<T> pointer to self (it will cast the pointer to T)
virtual RCP< const Number > mul(const Number &other) const =0
Multiplication.
virtual RCP< const Number > add(const Number &other) const =0
Addition.
virtual bool is_negative() const =0
bool is_one() const override
RCP< const Number > mul(const Number &other) const override
Multiplication.
RCP< const Number > pow(const Number &other) const override
Power.
static Poly series_atan(const Poly &s, const Poly &var, unsigned int prec)
bool is_minus_one() const override
RCP< const Number > add(const Number &other) const override
Addition.
bool is_negative() const override
bool __eq__(const Basic &o) const override
Test equality.
bool is_positive() const override
static Poly series_exp(const Poly &s, const Poly &var, unsigned int prec)
static Poly series_log(const Poly &s, const Poly &var, unsigned int prec)
bool is_zero() const override
bool is_complex() const override
Main namespace for SymEngine package.
RCP< const Basic > sign(const RCP< const Basic > &arg)
Canonicalize Sign.
RCP< const Basic > neg(const RCP< const Basic > &a)
Negation.