Hello again, On Fri, May 14, 2010 at 10:21:09AM +0200, Aless Lasaruk wrote:
i would like to report a bug in polynomial factorization.
I want to factorize the following multivariate polynomial:
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
GiNaC crashes with the following exception:
malloc: *** mmap(size=3774734336) failed (error code=12) *** error: can't allocate region *** set a breakpoint in malloc_error_break to debug terminate called after throwing an instance of 'std::bad_alloc' what(): St9bad_alloc
Could you please check if the patch below fixes the problem? (@developers: I've already sent this to ginac-devel) Best regards, Alexei 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..ae39f1d --- /dev/null +++ b/check/pgcd_infinite_loop.cpp @@ -0,0 +1,41 @@ +#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, 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; + return 0; +} diff --git a/ginac/polynomial/collect_vargs.cpp b/ginac/polynomial/collect_vargs.cpp index 89c3b8b..c4368fb 100644 --- a/ginac/polynomial/collect_vargs.cpp +++ b/ginac/polynomial/collect_vargs.cpp @@ -148,8 +148,14 @@ ex_collect_to_ex(const ex_collect_t& ec, const GiNaC::exvector& vars) exvector tv; tv.reserve(vars.size() + 1); for (std::size_t j = 0; j < vars.size(); ++j) { - if (ec[i].first[j] != 0) - tv.push_back(power(vars[j], ec[i].first[j])); + const exp_vector_t& exp_vector(ec[i].first); + + bug_on(exp_vector.size() != vars.size(), + "expected " << vars.size() << " variables, " + "expression has " << exp_vector.size() << " instead"); + + if (exp_vector[j] != 0) + tv.push_back(power(vars[j], exp_vector[j])); } tv.push_back(ec[i].second); ex tmp = (new mul(tv))->setflag(status_flags::dynallocated); @@ -171,6 +177,18 @@ ex lcoeff_wrt(ex e, const exvector& x) return ec.rbegin()->second; } +exp_vector_t degree_vector(ex e, const exvector& vars) +{ + e = e.expand(); + exp_vector_t dvec(vars.size()); + for (std::size_t i = vars.size(); i-- != 0; ) { + const int deg_i = e.degree(vars[i]); + e = e.coeff(vars[i], deg_i); + dvec[i] = deg_i; + } + return dvec; +} + cln::cl_I integer_lcoeff(const ex& e, const exvector& vars) { ex_collect_t ec; diff --git a/ginac/polynomial/collect_vargs.h b/ginac/polynomial/collect_vargs.h index 44c3d72..a927c3f 100644 --- a/ginac/polynomial/collect_vargs.h +++ b/ginac/polynomial/collect_vargs.h @@ -28,10 +28,34 @@ #include <cln/integer.h> #include <utility> // for std::pair #include <vector> +#include <algorithm> // std::lexicographical_compare namespace GiNaC { typedef std::vector<int> exp_vector_t; + +static inline bool operator<(const exp_vector_t& v1, const exp_vector_t& v2) +{ + return std::lexicographical_compare(v1.rbegin(), v1.rend(), + v2.rbegin(), v2.rend()); +} + +static inline bool operator>(const exp_vector_t& v1, const exp_vector_t& v2) +{ + if (v1 == v2) + return false; + return !(v1 < v2); +} + +static inline bool zerop(const exp_vector_t& v) +{ + for (exp_vector_t::const_reverse_iterator i = v.rbegin(); i != v.rend(); ++i) { + if (*i != 0) + return false; + } + return true; +} + typedef std::vector<std::pair<exp_vector_t, ex> > ex_collect_t; extern void @@ -46,6 +70,13 @@ ex_collect_to_ex(const ex_collect_t& ec, const exvector& x); */ extern ex lcoeff_wrt(ex e, const exvector& x); + +/** + * Degree vector of a leading term of a multivariate polynomial. + * (generalization of degree(expr, var)) + */ +extern exp_vector_t degree_vector(ex e, const exvector& vars); + /** * Leading coefficient c \in R (where R = Z or Z_p) of a multivariate * polynomial e \in R[x_0, \ldots, x_n] 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 0de0f70..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) @@ -134,8 +136,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 diff --git a/ginac/polynomial/primpart_content.cpp b/ginac/polynomial/primpart_content.cpp index f703ee7..694f29d 100644 --- a/ginac/polynomial/primpart_content.cpp +++ b/ginac/polynomial/primpart_content.cpp @@ -62,7 +62,7 @@ void primpart_content(ex& pp, ex& c, ex e, const exvector& vars, // p_1(x_n) p_2(x_0, \ldots, x_{n-1}) c = ec.rbegin()->second; ec.rbegin()->second = ex1; - pp = ex_collect_to_ex(ec, vars).expand().smod(numeric(p)); + pp = ex_collect_to_ex(ec, rest_vars).expand().smod(numeric(p)); return; }