[GiNaC-devel] [PATCH 4/4] [BUGFIX] pgcd(), chinrem_gcd(): use appropriate definition of the degree.

[Alexei Sheplyakov]

Effect: fixes (rare) GCD miscalculation.

pgcd() treats polynomials Z_p[x_0, ..., x_n] as R[x_0, ..., x_{n - 1}], where
R = Z_p[x_n]. Therefore one should use correct definition of the degree
(i.e. compute it w.r.t. x_0, ..., x_{n-1}, but NOT w.r.t. x_n!).

One should use appropriate definition of degree (that is, w.r.t.  x_0, ..., x_n)
in chinrem_gcd() too.

Thanks to Aless Lasaruk for a bug report.
---
 check/pgcd_infinite_loop.cpp |    8 ++++++--
 ginac/polynomial/mgcd.cpp    |   12 ++++++++----
 ginac/polynomial/pgcd.cpp    |   10 ++++++----
 3 files changed, 20 insertions(+), 10 deletions(-)

diff --git a/check/pgcd_infinite_loop.cpp b/check/pgcd_infinite_loop.cpp
index 013b746..ae39f1d 100644
--- a/check/pgcd_infinite_loop.cpp
+++ b/check/pgcd_infinite_loop.cpp
@@ -26,12 +26,16 @@ static const string srep("\
 
 int main(int argc, char** argv)
 {
-	cout << "Checking for more pgcd() bugs (infinite loop) ... " << flush;
+	cout << "Checking for more pgcd() bugs (infinite loop, miscalculation) ... " << flush;
 	parser the_parser;
 	ex e = the_parser(srep);
 	const symbol x = ex_to<symbol>(the_parser.get_syms()["x"]);
 	ex g = gcd(e, e.diff(x));
+	ex should_be = the_parser(string("u^4*z^2"));
+	if (!(g-should_be).expand().is_zero()) {
+		cout << "GCD was miscomputed. " << flush;
+		return 1;
+	}
 	cout << "not found. " << flush;
-	cout << g << endl << flush;
 	return 0;
 }
diff --git a/ginac/polynomial/mgcd.cpp b/ginac/polynomial/mgcd.cpp
index bf30b5c..24da934 100644
--- a/ginac/polynomial/mgcd.cpp
+++ b/ginac/polynomial/mgcd.cpp
@@ -27,6 +27,7 @@
 #include "primes_factory.h"
 #include "divide_in_z_p.h"
 #include "poly_cra.h"
+#include <numeric> // std::accumulate
 
 #include <cln/integer.h>
 
@@ -58,12 +59,15 @@ ex chinrem_gcd(const ex& A_, const ex& B_, const exvector& vars)
 	const cln::cl_I g_lc = cln::gcd(a_lc, b_lc);
 
 	const ex& x(vars.back());
-	int n = std::min(A.degree(x), B.degree(x));
+	exp_vector_t n = std::min(degree_vector(A, vars), degree_vector(B, vars));
+	const int nTot = std::accumulate(n.begin(), n.end(), 0);
 	const cln::cl_I A_max_coeff = to_cl_I(A.max_coefficient()); 
 	const cln::cl_I B_max_coeff = to_cl_I(B.max_coefficient());
-	const cln::cl_I lcoeff_limit = (cln::cl_I(1) << n)*cln::abs(g_lc)*
+
+	const cln::cl_I lcoeff_limit = (cln::cl_I(1) << nTot)*cln::abs(g_lc)*
 		std::min(A_max_coeff, B_max_coeff);
 
+
 	cln::cl_I q = 0;
 	ex H = 0;
 
@@ -83,8 +87,8 @@ ex chinrem_gcd(const ex& A_, const ex& B_, const exvector& vars)
 		const cln::cl_I Cp_lc = integer_lcoeff(Cp, vars);
 		const cln::cl_I nlc = smod(recip(Cp_lc, p)*g_lcp, p);
 		Cp = (Cp*numeric(nlc)).expand().smod(pnum);
-		int cp_deg = Cp.degree(x);
-		if (cp_deg == 0)
+		exp_vector_t cp_deg = degree_vector(Cp, vars);
+		if (zerop(cp_deg))
 			return numeric(g_lc);
 		if (zerop(q)) {
 			H = Cp;
diff --git a/ginac/polynomial/pgcd.cpp b/ginac/polynomial/pgcd.cpp
index fac3735..a33b02d 100644
--- a/ginac/polynomial/pgcd.cpp
+++ b/ginac/polynomial/pgcd.cpp
@@ -80,7 +80,8 @@ ex pgcd(const ex& A, const ex& B, const exvector& vars, const long p)
 	const ex lc_gcd = euclid_gcd(AL, BL, mainvar, p);
 
 	// The estimate of degree of the gcd of Ab and Bb
-	int gcd_deg = std::min(degree(Aprim, mainvar), degree(Bprim, mainvar));
+	exp_vector_t gcd_deg = std::min(degree_vector(Aprim, restvars),
+					degree_vector(Bprim, restvars));
 	eval_point_finder::value_type b;
 
 	eval_point_finder find_eval_point(p);
@@ -105,8 +106,11 @@ ex pgcd(const ex& A, const ex& B, const exvector& vars, const long p)
 		const cln::cl_I correct_lc = smod(lcb_gcd*recip(Cblc, p), p);
 		Cb = (Cb*numeric(correct_lc)).smod(pn);
 
+		const exp_vector_t img_gcd_deg = degree_vector(Cb, restvars);
+		// Test for relatively prime polynomials
+		if (zerop(img_gcd_deg))
+			return cont_gcd;
 		// Test for unlucky homomorphisms
-		const int img_gcd_deg = Cb.degree(restvars.back());
 		if (img_gcd_deg < gcd_deg) {
 			// The degree decreased, previous homomorphisms were
 			// bad, so we have to start it all over.
@@ -121,8 +125,6 @@ ex pgcd(const ex& A, const ex& B, const exvector& vars, const long p)
 			// evaluation point is bad. Skip it.
 			continue;
 		}
-		if (img_gcd_deg == 0)
-			return cont_gcd;
 
 		// Image has the same degree as the previous one
 		// (or at least not higher than the limit)
-- 
1.7.1