+ex power::series(const relational & r, int order, unsigned options) const
+{
+ // If basis is already a series, just power it
+ if (is_exactly_a<pseries>(basis))
+ return ex_to<pseries>(basis).power_const(ex_to<numeric>(exponent), order);
+
+ // Basis is not a series, may there be a singularity?
+ bool must_expand_basis = false;
+ try {
+ basis.subs(r, subs_options::no_pattern);
+ } catch (pole_error) {
+ must_expand_basis = true;
+ }
+
+ // Is the expression of type something^(-int)?
+ if (!must_expand_basis && !exponent.info(info_flags::negint)
+ && (!is_a<add>(basis) || !is_a<numeric>(exponent)))
+ return basic::series(r, order, options);
+
+ // Is the expression of type 0^something?
+ if (!must_expand_basis && !basis.subs(r, subs_options::no_pattern).is_zero()
+ && (!is_a<add>(basis) || !is_a<numeric>(exponent)))
+ return basic::series(r, order, options);
+
+ // Singularity encountered, is the basis equal to (var - point)?
+ if (basis.is_equal(r.lhs() - r.rhs())) {
+ epvector new_seq;
+ if (ex_to<numeric>(exponent).to_int() < order)
+ new_seq.push_back(expair(_ex1, exponent));
+ else
+ new_seq.push_back(expair(Order(_ex1), exponent));
+ return pseries(r, new_seq);
+ }
+
+ // No, expand basis into series
+
+ numeric numexp;
+ if (is_a<numeric>(exponent)) {
+ numexp = ex_to<numeric>(exponent);
+ } else {
+ numexp = 0;
+ }
+ const ex& sym = r.lhs();
+ // find existing minimal degree
+ int real_ldegree = basis.expand().ldegree(sym-r.rhs());
+ if (real_ldegree == 0) {
+ int orderloop = 0;
+ do {
+ orderloop++;
+ real_ldegree = basis.series(r, orderloop, options).ldegree(sym);
+ } while (real_ldegree == orderloop);
+ }
+
+ if (!(real_ldegree*numexp).is_integer())
+ throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
+ ex e = basis.series(r, (order + real_ldegree*(1-numexp)).to_int(), options);
+
+ ex result;
+ try {
+ result = ex_to<pseries>(e).power_const(numexp, order);
+ } catch (pole_error) {
+ epvector ser;
+ ser.push_back(expair(Order(_ex1), order));
+ result = pseries(r, ser);
+ }
+
+ return result;
+}
+
+
+/** Re-expansion of a pseries object. */
+ex pseries::series(const relational & r, int order, unsigned options) const
+{
+ const ex p = r.rhs();
+ GINAC_ASSERT(is_a<symbol>(r.lhs()));
+ const symbol &s = ex_to<symbol>(r.lhs());
+
+ if (var.is_equal(s) && point.is_equal(p)) {
+ if (order > degree(s))
+ return *this;
+ else {
+ epvector new_seq;
+ epvector::const_iterator it = seq.begin(), itend = seq.end();
+ while (it != itend) {
+ int o = ex_to<numeric>(it->coeff).to_int();
+ if (o >= order) {
+ new_seq.push_back(expair(Order(_ex1), o));
+ break;
+ }
+ new_seq.push_back(*it);
+ ++it;
+ }
+ return pseries(r, new_seq);
+ }
+ } else
+ return convert_to_poly().series(r, order, options);
+}
+
+ex integral::series(const relational & r, int order, unsigned options) const
+{
+ if (x.subs(r) != x)
+ throw std::logic_error("Cannot series expand wrt dummy variable");
+
+ // Expanding integrant with r substituted taken in boundaries.
+ ex fseries = f.series(r, order, options);
+ epvector fexpansion;
+ fexpansion.reserve(fseries.nops());
+ for (size_t i=0; i<fseries.nops(); ++i) {
+ ex currcoeff = ex_to<pseries>(fseries).coeffop(i);
+ currcoeff = (currcoeff == Order(_ex1))
+ ? currcoeff
+ : integral(x, a.subs(r), b.subs(r), currcoeff);
+ if (currcoeff != 0)
+ fexpansion.push_back(
+ expair(currcoeff, ex_to<pseries>(fseries).exponop(i)));
+ }
+
+ // Expanding lower boundary
+ ex result = (new pseries(r, fexpansion))->setflag(status_flags::dynallocated);
+ ex aseries = (a-a.subs(r)).series(r, order, options);
+ fseries = f.series(x == (a.subs(r)), order, options);
+ for (size_t i=0; i<fseries.nops(); ++i) {
+ ex currcoeff = ex_to<pseries>(fseries).coeffop(i);
+ if (is_order_function(currcoeff))
+ break;
+ ex currexpon = ex_to<pseries>(fseries).exponop(i);
+ int orderforf = order-ex_to<numeric>(currexpon).to_int()-1;
+ currcoeff = currcoeff.series(r, orderforf);
+ ex term = ex_to<pseries>(aseries).power_const(ex_to<numeric>(currexpon+1),order);
+ term = ex_to<pseries>(term).mul_const(ex_to<numeric>(-1/(currexpon+1)));
+ term = ex_to<pseries>(term).mul_series(ex_to<pseries>(currcoeff));
+ result = ex_to<pseries>(result).add_series(ex_to<pseries>(term));
+ }
+
+ // Expanding upper boundary
+ ex bseries = (b-b.subs(r)).series(r, order, options);
+ fseries = f.series(x == (b.subs(r)), order, options);
+ for (size_t i=0; i<fseries.nops(); ++i) {
+ ex currcoeff = ex_to<pseries>(fseries).coeffop(i);
+ if (is_order_function(currcoeff))
+ break;
+ ex currexpon = ex_to<pseries>(fseries).exponop(i);
+ int orderforf = order-ex_to<numeric>(currexpon).to_int()-1;
+ currcoeff = currcoeff.series(r, orderforf);
+ ex term = ex_to<pseries>(bseries).power_const(ex_to<numeric>(currexpon+1),order);
+ term = ex_to<pseries>(term).mul_const(ex_to<numeric>(1/(currexpon+1)));
+ term = ex_to<pseries>(term).mul_series(ex_to<pseries>(currcoeff));
+ result = ex_to<pseries>(result).add_series(ex_to<pseries>(term));
+ }
+
+ return result;