-ex pseries::series(const relational & r, int order, bool branchcut) const
-{
- const ex p = r.rhs();
- GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol));
- const symbol *s = static_cast<symbol *>(r.lhs().bp);
-
- 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, branchcut);
+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;