From: Richard Kreckel <kreckel@ginac.de>
Date: Fri, 21 Aug 2020 15:47:19 +0000 (+0200)
Subject: Avoid unnecessary expansion in sqrfree_yun().
X-Git-Tag: release_1-8-0~26
X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=e65cbada7281938475aa0cf31a3d56405ac91d4e;ds=inline

Avoid unnecessary expansion in sqrfree_yun().

Don't make the polynomial primitive inside Yun's algorithm. Computing the
primitive part will incur an extra expansion of the polynomial and this
may distroy a pre-factored structure in the other variables. In the multi-
variate case, working with primitive polynomials isn't sufficient anyways
to avoid a final division to discover lost factors.
---

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 6ad61137..f8615499 100644
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -1797,8 +1797,7 @@ ex lcm(const ex &a, const ex &b, bool check_args)
  *  @return   vector of expairs (factor, exponent), sorted by exponents */
 static epvector sqrfree_yun(const ex &a, const symbol &x)
 {
-	// make polynomial primitive w.r.t. x
-	ex w = a.primpart(x);
+	ex w = a;
 	ex z = w.diff(x);
 	ex g = gcd(w, z);
 	if (g.is_zero()) {
@@ -1814,6 +1813,10 @@ static epvector sqrfree_yun(const ex &a, const symbol &x)
 	ex i = 0;  // exponent
 	do {
 		w = quo(w, g, x);
+		if (w.is_zero()) {
+			// hidden zero
+			break;
+		}
 		z = quo(z, g, x) - w.diff(x);
 		i += 1;
 		if (w.is_equal(x)) {
@@ -1829,20 +1832,18 @@ static epvector sqrfree_yun(const ex &a, const symbol &x)
 	} while (!z.is_zero());
 
 	// correct for lost factor
-	// (NB: This factor is not necessarily the unit and content lost when
-	// we made the polynomial primitive since Yun's algorithm only finds
-	// factors up to a unit.)
+	// (being based on GCDs, Yun's algorithm only finds factors up to a unit)
 	const ex lost_factor = quo(a, mul{factors}, x);
 	if (lost_factor.is_equal(_ex1)) {
 		// trivial lost factor
 		return factors;
 	}
-	if (factors[0].coeff.is_equal(1)) {
+	if (!factors.empty() && factors[0].coeff.is_equal(1)) {
 		// multiply factor^1 with lost_factor
 		factors[0].rest *= lost_factor;
 		return factors;
 	}
-	// only factors with higher exponent: prepend lost_factor^1 to the results
+	// no factor^1: prepend lost_factor^1 to the results
 	epvector results = {expair(lost_factor, 1)};
 	std::move(factors.begin(), factors.end(), std::back_inserter(results));
 	return results;