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