-ex power::series(const relational & r, int order) const
-{
- ex e;
- if (!is_ex_exactly_of_type(basis, pseries)) {
- // Basis is not a series, may there be a singulary?
- if (!exponent.info(info_flags::negint))
- return basic::series(r, order);
-
- // Expression is of type something^(-int), check for singularity
- if (!basis.subs(r).is_zero())
- return basic::series(r, order);
-
- // Singularity encountered, expand basis into series
- e = basis.series(r, order);
- } else {
- // Basis is a series
- e = basis;
- }
-
- // Power e
- return ex_to_pseries(e).power_const(ex_to_numeric(exponent), order);
+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;
+ }
+
+ bool exponent_is_regular = true;
+ try {
+ exponent.subs(r, subs_options::no_pattern);
+ } catch (pole_error) {
+ exponent_is_regular = false;
+ }
+
+ if (!exponent_is_regular) {
+ ex l = exponent*log(basis);
+ // this == exp(l);
+ ex le = l.series(r, order, options);
+ // Note: expanding exp(l) won't help, since that will attempt
+ // Taylor expansion, and fail (because exponent is "singular")
+ // Still l itself might be expanded in Taylor series.
+ // Examples:
+ // sin(x)/x*log(cos(x))
+ // 1/x*log(1 + x)
+ return exp(le).series(r, order, options);
+ // Note: if l happens to have a Laurent expansion (with
+ // negative powers of (var - point)), expanding exp(le)
+ // will barf (which is The Right Thing).
+ }
+
+ // 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
+ ex eb = basis.expand();
+ int real_ldegree = 0;
+ if (eb.info(info_flags::rational_function))
+ real_ldegree = eb.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;