[GiNaC-devel] [PATCH 1/4] [BUGFIX] pgcd(): avoid infinite loop if the first GCD candidate coincides with GCD

[Alexei Sheplyakov]

136	if (H_lcoeff.is_equal(lc_gcd)) {
137		if ((Hprev-H).expand().smod(pn).is_zero()) // (*)
138			continue;

The check (*) can result in infinite loop if the (very) first GCD candidate
gives the correct GCD: any subsequent iterations give the same result, so
Hprev and H will be equal (hence the infinite loop). That check is not very
useful (AFAIR I was trying to avoid extra division checks), so let's remove it.

Thanks to Aless Lasaruk for a bug report.
---
 check/Makefile.am            |    4 ++++
 check/pgcd_infinite_loop.cpp |   37 +++++++++++++++++++++++++++++++++++++
 ginac/polynomial/pgcd.cpp    |    2 --
 3 files changed, 41 insertions(+), 2 deletions(-)
 create mode 100644 check/pgcd_infinite_loop.cpp

diff --git a/check/Makefile.am b/check/Makefile.am
index 3751613..7aaed07 100644
--- a/check/Makefile.am
+++ b/check/Makefile.am
@@ -31,6 +31,7 @@ EXAMS = exam_paranoia \
 	exam_mod_gcd \
 	bugme_chinrem_gcd \
 	pgcd_relatively_prime_bug \
+	pgcd_infinite_loop \
 	exam_cra
 
 TIMES = time_dennyfliegner \
@@ -259,6 +260,9 @@ bugme_chinrem_gcd_LDADD = ../ginac/libginac.la
 pgcd_relatively_prime_bug_SOURCES = pgcd_relatively_prime_bug.cpp
 pgcd_relatively_prime_bug_LDADD = ../ginac/libginac.la
 
+pgcd_infinite_loop_SOURCES =  pgcd_infinite_loop.cpp
+pgcd_infinite_loop_LDADD =  ../ginac/libginac.la
+
 AM_CPPFLAGS = -I$(srcdir)/../ginac -I../ginac -DIN_GINAC
 
 CLEANFILES = exam.gar
diff --git a/check/pgcd_infinite_loop.cpp b/check/pgcd_infinite_loop.cpp
new file mode 100644
index 0000000..013b746
--- /dev/null
+++ b/check/pgcd_infinite_loop.cpp
@@ -0,0 +1,37 @@
+#include <iostream>
+#include <string>
+#include "ginac.h"
+using namespace GiNaC;
+using namespace std;
+
+static const string srep("\
+792*z^8*w^4*x^3*y^4*u^7 + 24*z^4*w^4*x^2*y^3*u^4	\
++ 264*z^8*w^3*x^2*y^7*u^5 + 198*z^4*w^5*x^5*y*u^6	\ 
++ 110*z^2*w^3*x^5*y^4*u^6 - 120*z^8*w*x^4*u^6		\ 
+- 480*z^5*w*x^4*y^6*u^8 - 720*z^7*x^3*y^3*u^7		\
++ 165*z^4*w^2*x^4*y*u^5 + 450*z^8*w^6*x^2*y*u^8		\
++ 40*z^2*w^3*x^3*y^3*u^6 - 288*z^7*w^2*x^3*y^6*u^6	\
++ 250*z^6*w^4*x^2*y^4*u^8 + 576*z^7*w^7*x^2*y^4*u^8	\
+- 80*z^6*w^2*x^5*y^3*u^7 - 144*z^8*w^4*x^5*u^7		\
++ 120*z^4*w*x^2*y^6*u^6 + 320*z^5*w^5*x^2*y^7*u^8	\
++ 192*z^7*w^6*x*y^7*u^6 - 12*z^4*w^3*x^3*y^5*u^6	\
+- 36*z^4*w^4*x^4*y^2*u^8 + 72*z^4*w^5*x^3*u^6		\
+- 20*z^2*w^2*x^4*y^5*u^8 + 660*z^8*w*x^2*y^4*u^6	\
++ 66*z^4*w^4*x^4*y^4*u^4 + 440*z^6*w^2*x^3*y^7*u^7	\
+- 30*z^4*w*x^3*y^2*u^7 - 48*z^8*w^3*x^4*y^3*u^5		\
++ 72*z^6*w^2*x*y^6*u^4 - 864*z^7*w^3*x^4*y^3*u^8	\
++ 480*z^7*w^4*x*y^4*u^7 + 60*z^4*w^2*x^2*u^5		\
++ 375*z^8*w^3*x*y*u^7 + 150*z^8*w^5*x*y^4*u^6		\
++ 180*z^6*x*y^3*u^5 + 216*z^6*w^3*x^2*y^3*u^6");
+
+int main(int argc, char** argv)
+{
+	cout << "Checking for more pgcd() bugs (infinite loop) ... " << 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));
+	cout << "not found. " << flush;
+	cout << g << endl << flush;
+	return 0;
+}
diff --git a/ginac/polynomial/pgcd.cpp b/ginac/polynomial/pgcd.cpp
index 0de0f70..fac3735 100644
--- a/ginac/polynomial/pgcd.cpp
+++ b/ginac/polynomial/pgcd.cpp
@@ -134,8 +134,6 @@ ex pgcd(const ex& A, const ex& B, const exvector& vars, const long p)
 		const ex H_lcoeff = lcoeff_wrt(H, restvars);
 
 		if (H_lcoeff.is_equal(lc_gcd)) {
-			if ((Hprev-H).expand().smod(pn).is_zero())
-				continue;
 			ex C /* primitive part of H */, contH /* dummy */;
 			primpart_content(C, contH, H, vars, p);
 			// Normalize GCD so that leading coefficient is 1
-- 
1.7.1