From 63f9b0510fecd2b50f75df37c9525bc215dc2578 Mon Sep 17 00:00:00 2001
From: Richard Kreckel <kreckel@ginac.de>
Date: Wed, 3 Feb 2016 08:13:58 +0100
Subject: [PATCH] [bugfix] fix sqrfree(poly) for zero polynomials in disguise.

Yun's algorithm now handles polynomials which would become zero after
expanding.
---
 check/exam_paranoia.cpp | 47 +++++++++++++++++++++++++++++++++++++++++
 ginac/normal.cpp        | 14 ++++++++----
 2 files changed, 57 insertions(+), 4 deletions(-)

diff --git a/check/exam_paranoia.cpp b/check/exam_paranoia.cpp
index 2882cff0..2b9430d0 100644
--- a/check/exam_paranoia.cpp
+++ b/check/exam_paranoia.cpp
@@ -578,6 +578,52 @@ static unsigned exam_paranoia22()
 	return 0;
 }
 
+// Bug in sqrfree_yun (fixed 2016-02-02).
+static unsigned exam_paranoia23()
+{
+	unsigned result = 0;
+	symbol x("x");
+	ex e;
+
+	e = (x-1)*(x+1) - x*x + 1;  // an unexpanded 0...
+	try {
+		ex f = sqrfree(e);
+		if (!f.is_zero()) {
+			clog << "sqrfree(" << e << ") returns " << f << " instead of 0\n";
+			++result;
+		}
+	} catch (const exception &err) {
+		clog << "sqrfree(" << e << ") throws " << err.what() << endl;
+		++result;
+	}
+
+	e = pow(x-1,3) - expand(pow(x-1,3));  // ...still after differentiating...
+	try {
+		ex f = sqrfree(e);
+		if (!f.is_zero()) {
+			clog << "sqrfree(" << e << ") returns " << f << " instead of 0\n";
+			++result;
+		}
+	} catch (const exception &err) {
+		clog << "sqrfree(" << e << ") throws " << err.what() << endl;
+		++result;
+	}
+
+	e = pow(x-1,4) - expand(pow(x-1,4));  // ...and after differentiating twice.
+	try {
+		ex f = sqrfree(e);
+		if (!f.is_zero()) {
+			clog << "sqrfree(" << e << ") returns " << f << " instead of 0\n";
+			++result;
+		}
+	} catch (const exception &err) {
+		clog << "sqrfree(" << e << ") throws " << err.what() << endl;
+		++result;
+	}
+
+	return result;
+}
+
 unsigned exam_paranoia()
 {
 	unsigned result = 0;
@@ -607,6 +653,7 @@ unsigned exam_paranoia()
 	result += is_polynomial_false_positive(); cout << '.' << flush;
 	result += exam_paranoia21();  cout << '.' << flush;
 	result += exam_paranoia22();  cout << '.' << flush;
+	result += exam_paranoia23();  cout << '.' << flush;
 	
 	return result;
 }
diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 25c3c80c..89fe532a 100644
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -1790,7 +1790,7 @@ ex lcm(const ex &a, const ex &b, bool check_args)
  *  Yun's algorithm.  Used internally by sqrfree().
  *
  *  @param a  multivariate polynomial over Z[X], treated here as univariate
- *            polynomial in x.
+ *            polynomial in x (needs not be expanded).
  *  @param x  variable to factor in
  *  @return   vector of factors sorted in ascending degree */
 static exvector sqrfree_yun(const ex &a, const symbol &x)
@@ -1799,6 +1799,9 @@ static exvector sqrfree_yun(const ex &a, const symbol &x)
 	ex w = a;
 	ex z = w.diff(x);
 	ex g = gcd(w, z);
+	if (g.is_zero()) {
+		return res;
+	}
 	if (g.is_equal(_ex1)) {
 		res.push_back(a);
 		return res;
@@ -1806,6 +1809,9 @@ static exvector sqrfree_yun(const ex &a, const symbol &x)
 	ex y;
 	do {
 		w = quo(w, g, x);
+		if (w.is_zero()) {
+			return res;
+		}
 		y = quo(z, g, x);
 		z = y - w.diff(x);
 		g = gcd(w, z);
@@ -1817,7 +1823,7 @@ static exvector sqrfree_yun(const ex &a, const symbol &x)
 
 /** Compute a square-free factorization of a multivariate polynomial in Q[X].
  *
- *  @param a  multivariate polynomial over Q[X]
+ *  @param a  multivariate polynomial over Q[X] (needs not be expanded)
  *  @param l  lst of variables to factor in, may be left empty for autodetection
  *  @return   a square-free factorization of \p a.
  *
@@ -1852,8 +1858,8 @@ static exvector sqrfree_yun(const ex &a, const symbol &x)
  */
 ex sqrfree(const ex &a, const lst &l)
 {
-	if (is_exactly_a<numeric>(a) ||     // algorithm does not trap a==0
-	    is_a<symbol>(a))        // shortcut
+	if (is_exactly_a<numeric>(a) ||
+	    is_a<symbol>(a))        // shortcuts
 		return a;
 
 	// If no lst of variables to factorize in was specified we have to
-- 
2.44.0