[GiNaC-devel] Factorization Bug
Alexei Sheplyakov
alexei.sheplyakov at gmail.com
Mon May 17 08:24:19 CEST 2010
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;
}
More information about the GiNaC-devel
mailing list