From f39f8353dae2542f654c170dc7fb9c97c5f29820 Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Thu, 15 Feb 2001 19:47:57 +0000 Subject: [PATCH] pseries::power_const: fix a critical bug that led to wrong computations. power::series: better handling of divergent basis cases. --- ginac/pseries.cpp | 47 ++++++++++++++++++++++++++++++++++++----------- 1 file changed, 36 insertions(+), 11 deletions(-) diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp index 36fb4429..9adc8fc4 100644 --- a/ginac/pseries.cpp +++ b/ginac/pseries.cpp @@ -755,17 +755,34 @@ ex mul::series(const relational & r, int order, unsigned options) const * @param deg truncation order of series calculation */ ex pseries::power_const(const numeric &p, int deg) const { - int i; + // method: + // 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. const symbol *s = static_cast(var.bp); int ldeg = ldegree(*s); - // Calculate coefficients of powered series + // Compute coefficients of the powered series exvector co; co.reserve(deg); - ex co0; - co.push_back(co0 = power(coeff(*s, ldeg), p)); + co.push_back(power(coeff(*s, ldeg), p)); bool all_sums_zero = true; - for (i=1; i