From: Richard Kreckel <kreckel@ginac.de>
Date: Sun, 6 Aug 2023 17:27:51 +0000 (+0200)
Subject: Fix series of positive powers of polynomials.
X-Git-Tag: release_1-8-7~3
X-Git-Url: https://ginac.de/ginac.git/tutorial/ginac.git?a=commitdiff_plain;h=9435d80ac31c97c2a698928e59cc4489f1fad5e7;p=ginac.git

Fix series of positive powers of polynomials.

Before calling pseries::power_const(), make sure the basis series has
sufficient terms. In particular, do not shrink the order of expansion
- only grow it.

Fix various hairy problems in pseries::power_const() when a polynomial
is raised to a positive integer power.

(The special cases here can actually also be computed by Taylor
expansion but the fixes should be more general.)

Also add test cases.

Thanks to Vitaly Magerya <vmagerya@gmail.com> for reporting this.
---

diff --git a/check/exam_pseries.cpp b/check/exam_pseries.cpp
index f0fead05..671d5751 100644
--- a/check/exam_pseries.cpp
+++ b/check/exam_pseries.cpp
@@ -390,6 +390,24 @@ static unsigned exam_series14()
 	return result;
 }
 
+// Test expansion of powers of polynomials.
+static unsigned exam_series15()
+{
+	unsigned result = 0;
+
+	ex e = pow(x + pow(x,2), 2);
+
+	result += check_series(e, 0, Order(1), 0);
+	result += check_series(e, 0, Order(x), 1);
+	result += check_series(e, 0, Order(pow(x,2)), 2);
+	result += check_series(e, 0, pow(x,2) + Order(pow(x,3)), 3);
+	result += check_series(e, 0, pow(x,2) + 2*pow(x,3) + Order(pow(x,4)), 4);
+	result += check_series(e, 0, pow(x,2) + 2*pow(x,3) + pow(x,4), 5);
+	result += check_series(e, 0, pow(x,2) + 2*pow(x,3) + pow(x,4), 6);
+
+	return result;
+}
+
 unsigned exam_pseries()
 {
 	unsigned result = 0;
@@ -410,6 +428,7 @@ unsigned exam_pseries()
 	result += exam_series12();  cout << '.' << flush;
 	result += exam_series13();  cout << '.' << flush;
 	result += exam_series14();  cout << '.' << flush;
+	result += exam_series15();  cout << '.' << flush;
 	
 	return result;
 }
diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp
index 00188f01..af26a6dd 100644
--- a/ginac/pseries.cpp
+++ b/ginac/pseries.cpp
@@ -992,8 +992,8 @@ ex pseries::power_const(const numeric &p, int deg) const
 	// 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.
+	// repeat the above derivation.  The leading power of C2(x) = A2(x)^p is
+	// then of course a_0^p*x^(p*m) but the recurrence formula still holds.
 	
 	if (seq.empty()) {
 		// as a special case, handle the empty (zero) series honoring the
@@ -1005,16 +1005,27 @@ ex pseries::power_const(const numeric &p, int deg) const
 		else
 			return *this;
 	}
-	
-	const int ldeg = ldegree(var);
-	if (!(p*ldeg).is_integer())
+
+	const int base_ldeg = ldegree(var);
+	if (!(p*base_ldeg).is_integer())
 		throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
+	int new_ldeg = (p*base_ldeg).to_int();
+
+	const int base_deg = degree(var);
+	int new_deg = deg;
+	if (p.is_pos_integer()) {
+		// No need to compute beyond p*base_deg.
+		new_deg = std::min((p*base_deg).to_int(), deg);
+	}
 
 	// adjust number of coefficients
-	int numcoeff = deg - (p*ldeg).to_int();
+	int numcoeff = new_deg - new_ldeg;
+	if (new_deg < deg)
+		numcoeff += 1;
+
 	if (numcoeff <= 0) {
-		epvector epv { expair(Order(_ex1), deg) };
-		return dynallocate<pseries>(relational(var,point), std::move(epv));
+		return dynallocate<pseries>(relational(var, point),
+		                            epvector{{Order(_ex1), deg}});
 	}
 	
 	// O(x^n)^(-m) is undefined
@@ -1024,33 +1035,35 @@ ex pseries::power_const(const numeric &p, int deg) const
 	// Compute coefficients of the powered series
 	exvector co;
 	co.reserve(numcoeff);
-	co.push_back(pow(coeff(var, ldeg), p));
+	co.push_back(pow(coeff(var, base_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);
+			ex c = coeff(var, j + base_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);
+		co.push_back(sum / coeff(var, base_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.emplace_back(expair(co[i], p * ldeg + i));
+		if (!co[i].is_zero()) {
+			new_seq.emplace_back(expair(co[i], new_ldeg + i));
+		}
 		if (is_order_function(co[i])) {
 			higher_order = true;
 			break;
 		}
 	}
-	if (!higher_order)
-		new_seq.emplace_back(expair(Order(_ex1), p * ldeg + numcoeff));
+	if (!higher_order && new_deg == deg) {
+		new_seq.emplace_back(expair{Order(_ex1), new_deg});
+	}
 
 	return pseries(relational(var,point), std::move(new_seq));
 }
@@ -1150,7 +1163,8 @@ ex power::series(const relational & r, int order, unsigned options) const
 
 	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);
+	int extra_terms = (real_ldegree*(1-numexp)).to_int();
+	ex e = basis.series(r, order + std::max(0, extra_terms), options);
 	
 	ex result;
 	try {