- int i;
- const symbol *s = static_cast<symbol *>(var.bp);
- int ldeg = ldegree(*s);
-
- // Calculate coefficients of powered series
- exvector co;
- co.reserve(deg);
- ex co0;
- co.push_back(co0 = power(coeff(*s, ldeg), p));
- bool all_sums_zero = true;
- for (i=1; i<deg; i++) {
- ex sum = _ex0();
- for (int j=1; j<=i; j++) {
- ex c = coeff(*s, j + ldeg);
- if (is_order_function(c)) {
- co.push_back(Order(_ex1()));
- break;
- } else
- sum += (p * j - (i - j)) * co[i - j] * c;
- }
- if (!sum.is_zero())
- all_sums_zero = false;
- co.push_back(co0 * sum / numeric(i));
- }
-
- // Construct new series (of non-zero coefficients)
- epvector new_seq;
- bool higher_order = false;
- for (i=0; i<deg; i++) {
- if (!co[i].is_zero())
- new_seq.push_back(expair(co[i], numeric(i) + p * ldeg));
- if (is_order_function(co[i])) {
- higher_order = true;
- break;
- }
- }
- if (!higher_order && !all_sums_zero)
- new_seq.push_back(expair(Order(_ex1()), numeric(deg) + p * ldeg));
- return pseries(var, point, new_seq);
+ // method:
+ // (due to Leonhard Euler)
+ // let A(x) be this series and for the time being let it start with a
+ // constant (later we'll generalize):
+ // A(x) = a_0 + a_1*x + a_2*x^2 + ...
+ // We want to compute
+ // C(x) = A(x)^p
+ // C(x) = c_0 + c_1*x + c_2*x^2 + ...
+ // Taking the derivative on both sides and multiplying with A(x) one
+ // immediately arrives at
+ // C'(x)*A(x) = p*C(x)*A'(x)
+ // Multiplying this out and comparing coefficients we get the recurrence
+ // formula
+ // c_i = (i*p*a_i*c_0 + ((i-1)*p-1)*a_{i-1}*c_1 + ...
+ // ... + (p-(i-1))*a_1*c_{i-1})/(a_0*i)
+ // which can easily be solved given the starting value c_0 = (a_0)^p.
+ // For the more general case where the leading coefficient of A(x) is not
+ // a constant, just consider A2(x) = A(x)*x^m, with some integer m and
+ // repeat the above derivation. The leading power of C2(x) = A2(x)^2 is
+ // then of course x^(p*m) but the recurrence formula still holds.
+
+ if (seq.empty()) {
+ // as a special case, handle the empty (zero) series honoring the
+ // usual power laws such as implemented in power::eval()
+ if (p.real().is_zero())
+ throw std::domain_error("pseries::power_const(): pow(0,I) is undefined");
+ else if (p.real().is_negative())
+ throw pole_error("pseries::power_const(): division by zero",1);
+ else
+ return *this;
+ }
+
+ const int ldeg = ldegree(var);
+ if (!(p*ldeg).is_integer())
+ throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
+
+ // adjust number of coefficients
+ int numcoeff = deg - (p*ldeg).to_int();
+ if (numcoeff <= 0) {
+ epvector epv;
+ epv.reserve(1);
+ epv.push_back(expair(Order(_ex1), deg));
+ return (new pseries(relational(var,point), epv))
+ ->setflag(status_flags::dynallocated);
+ }
+
+ // O(x^n)^(-m) is undefined
+ if (seq.size() == 1 && is_order_function(seq[0].rest) && p.real().is_negative())
+ throw pole_error("pseries::power_const(): division by zero",1);
+
+ // Compute coefficients of the powered series
+ exvector co;
+ co.reserve(numcoeff);
+ co.push_back(power(coeff(var, ldeg), p));
+ for (int i=1; i<numcoeff; ++i) {
+ ex sum = _ex0;
+ for (int j=1; j<=i; ++j) {
+ ex c = coeff(var, j + ldeg);
+ if (is_order_function(c)) {
+ co.push_back(Order(_ex1));
+ break;
+ } else
+ sum += (p * j - (i - j)) * co[i - j] * c;
+ }
+ co.push_back(sum / coeff(var, ldeg) / i);
+ }
+
+ // Construct new series (of non-zero coefficients)
+ epvector new_seq;
+ bool higher_order = false;
+ for (int i=0; i<numcoeff; ++i) {
+ if (!co[i].is_zero())
+ new_seq.push_back(expair(co[i], p * ldeg + i));
+ if (is_order_function(co[i])) {
+ higher_order = true;
+ break;
+ }
+ }
+ if (!higher_order)
+ new_seq.push_back(expair(Order(_ex1), p * ldeg + numcoeff));
+
+ return pseries(relational(var,point), new_seq);
+}
+
+
+/** Return a new pseries object with the powers shifted by deg. */
+pseries pseries::shift_exponents(int deg) const
+{
+ epvector newseq = seq;
+ epvector::iterator i = newseq.begin(), end = newseq.end();
+ while (i != end) {
+ i->coeff += deg;
+ ++i;
+ }
+ return pseries(relational(var, point), newseq);