From eaaf0b45634f43e67de7dcfe2a2245c7a9300bb1 Mon Sep 17 00:00:00 2001
From: Jens Vollinga <vollinga@thep.physik.uni-mainz.de>
Date: Mon, 16 Feb 2004 22:41:49 +0000
Subject: [PATCH] New series expansion scheme

---
 ginac/pseries.cpp | 43 ++++++++++++++++++++++++++++++++++---------
 1 file changed, 34 insertions(+), 9 deletions(-)

diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp
index dd4f51ab..8330d4db 100644
--- a/ginac/pseries.cpp
+++ b/ginac/pseries.cpp
@@ -520,7 +520,7 @@ ex basic::series(const relational & r, int order, unsigned options) const
 		seq.push_back(expair(coeff, _ex0));
 	
 	int n;
-	for (n=1; n<order; ++n) {
+	for (n=1; n<=order; ++n) {
 		fac = fac.mul(n);
 		// We need to test for zero in order to see if the series terminates.
 		// The problem is that there is no such thing as a perfect test for
@@ -535,9 +535,34 @@ ex basic::series(const relational & r, int order, unsigned options) const
 	}
 	
 	// Higher-order terms, if present
-	deriv = deriv.diff(s);
-	if (!deriv.expand().is_zero())
-		seq.push_back(expair(Order(_ex1), n));
+	int ldeg;
+	try {
+		ldeg = std::abs(deriv.ldegree(s));
+	}
+	catch (runtime_error) {
+		ldeg = 0;
+	}
+	if (ldeg != 0) {
+		// pure polynomial
+		if (!deriv.subs(r, subs_options::no_pattern).is_zero()) {
+			seq.push_back(expair(Order(_ex1), n));
+		} else {
+		seq.push_back(expair(Order(_ex1), n+ldeg-1));
+		}
+	} else {
+		// something more complicated -> loop until next coefficient is found
+		for (;; ++n) {
+			deriv = deriv.diff(s).expand();
+			if (deriv.is_zero()) {
+				break;
+			}
+			if (!deriv.subs(r, subs_options::no_pattern).is_zero()) {
+				seq.push_back(expair(Order(_ex1), n));
+				break;
+			}
+		}
+	}
+	
 	return pseries(r, seq);
 }
 
@@ -774,7 +799,7 @@ ex mul::series(const relational & r, int order, unsigned options) const
 			acc = ex_to<pseries>(acc.mul_series(ex_to<pseries>(op)));
 	}
 
-	if (negldegree && (seq.size() > 1)) {
+	if (seq.size() > 1) {
 
 		// re-calculation of series
 
@@ -857,7 +882,7 @@ ex pseries::power_const(const numeric &p, int deg) const
 	co.reserve(deg);
 	co.push_back(power(coeff(var, ldeg), p));
 	bool all_sums_zero = true;
-	for (int i=1; i<deg; ++i) {
+	for (int i=1; i<=deg; ++i) {
 		ex sum = _ex0;
 		for (int j=1; j<=i; ++j) {
 			ex c = coeff(var, j + ldeg);
@@ -875,7 +900,7 @@ ex pseries::power_const(const numeric &p, int deg) const
 	// Construct new series (of non-zero coefficients)
 	epvector new_seq;
 	bool higher_order = false;
-	for (int i=0; i<deg; ++i) {
+	for (int i=0; i<=deg; ++i) {
 		if (!co[i].is_zero())
 			new_seq.push_back(expair(co[i], p * ldeg + i));
 		if (is_order_function(co[i])) {
@@ -884,7 +909,7 @@ ex pseries::power_const(const numeric &p, int deg) const
 		}
 	}
 	if (!higher_order && !all_sums_zero)
-		new_seq.push_back(expair(Order(_ex1), p * ldeg + deg));
+		new_seq.push_back(expair(Order(_ex1), p * ldeg + deg + 1));
 
 	return pseries(relational(var,point), new_seq);
 }
@@ -923,7 +948,7 @@ ex power::series(const relational & r, int order, unsigned options) const
 	// Is the expression of type something^(-int)?
 	if (!must_expand_basis && !exponent.info(info_flags::negint))
 		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())
 		return basic::series(r, order, options);
-- 
2.50.0