+ // 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;
+ }