From 4ee761760b3db8649b8b616256cd7466fe2cd033 Mon Sep 17 00:00:00 2001 From: Jens Vollinga Date: Thu, 18 Sep 2008 20:26:37 +0200 Subject: [PATCH] Fixed various bugs in multivariate factorization. --- check/exam_factor.cpp | 41 ++++- ginac/factor.cpp | 362 +++++++++++++++++++++++++++++------------- 2 files changed, 292 insertions(+), 111 deletions(-) diff --git a/check/exam_factor.cpp b/check/exam_factor.cpp index e3552bfc..8851b7cc 100644 --- a/check/exam_factor.cpp +++ b/check/exam_factor.cpp @@ -118,12 +118,51 @@ static unsigned exam_factor2() e = ex("x+y", syms); result += check_factor(e); - e = ex("x+y", syms); + e = ex("(x^2-y+1)*(x+y)", syms); result += check_factor(e); e = ex("-2*(x+y)*(x-y)", syms); result += check_factor(e); + e = ex("(16+x^2*z^3)*(-17+3*x-5*z)*(2*x+3*z)*(x-y^2-z^3)", syms); + result += check_factor(e); + + e = ex("(x-y*z)*(x-y^2-z^3)*(x+y+z)", syms); + result += check_factor(e); + + e = ex("-(y^2-x+z^3)*x*(x+y+z)", syms); + result += check_factor(e); + + e = ex("-316*(3*x-4*z)*(2*x+3*z)*(x+y)*(-1+x)", syms); + result += check_factor(e); + + e = ex("(x+x^3+z^2)*(3*x-4*z)", syms); + result += check_factor(e); + + e = ex("250*(-3+x)*(4*z-3*x)*(x^3+z^2+x)*x", syms); + result += check_factor(e); + + e = ex("327*(x+z^2+x^3)*(3*x-4*z)*(-7+5*x-x^3)*(1+x+x^2)", syms); + result += check_factor(e); + + e = ex("x-y^2-z^3", syms); + result += check_factor(e); + + e = ex("-390*(7+3*x^4)*(2+x^2)*(x-z^3-y^2)", syms); + result += check_factor(e); + + e = ex("55*(1+x)^2*(3*x-4*z)*(1+x+x^2)*(x+x^3+z^2)", syms); + result += check_factor(e); + + e = ex("x+y*x-1", syms); + result += check_factor(e); + + e = ex("390*(-1+x^6-x)*(7+3*x^4)*(2+x^2)*(y+x)*(-1+y-x^2)*(1+x^2+x)^2", syms); + result += check_factor(e); + + e = ex("310*(y+x)*(-1+y-x^2)", syms); + result += check_factor(e); + return result; } diff --git a/ginac/factor.cpp b/ginac/factor.cpp index 73e7d6c1..84b96eb8 100644 --- a/ginac/factor.cpp +++ b/ginac/factor.cpp @@ -1554,12 +1554,13 @@ vector multivar_diophant(const vector& a_, const ex& x, const ex& c, con vector anew = a; for ( size_t i=0; i Inew = I; Inew.pop_back(); - vector sigma = multivar_diophant(anew, x, cnew, Inew, d, p, k); + sigma = multivar_diophant(anew, x, cnew, Inew, d, p, k); + DCOUTVAR(sigma); ex buf = c; for ( size_t i=0; i multivar_diophant(const vector& a_, const ex& x, const ex& c, con ex e = buf; e = make_modular(e, R); + e = e.expand(); + DCOUTVAR(e); + DCOUTVAR(d); ex monomial = 1; for ( size_t m=1; m<=d; ++m ) { - while ( !e.is_zero() ) { + DCOUTVAR(m); + while ( !e.is_zero() && e.has(xnu) ) { monomial *= (xnu - alphanu); monomial = expand(monomial); + DCOUTVAR(xnu); + DCOUTVAR(alphanu); ex cm = e.diff(ex_to(xnu), m).subs(xnu==alphanu) / factorial(m); + DCOUTVAR(cm); if ( !cm.is_zero() ) { vector delta_s = multivar_diophant(anew, x, cm, Inew, d, p, k); + DCOUTVAR(delta_s); ex buf = e; for ( size_t j=0; j& I, unsigne DCOUTVAR(j); vector U1 = U; ex monomial = 1; + DCOUTVAR(U); for ( size_t m=0; m& I, unsigne U[m] = U[m] - U[m].lcoeff(x) * pow(x,deg) + coef * pow(x,deg); } } + DCOUTVAR(U); ex Uprod = 1; for ( size_t i=0; i& I, unsigne ex e = expand(A[j-1] - Uprod); DCOUTVAR(e); + vector newI; + for ( size_t i=1; i<=j-2; ++i ) { + newI.push_back(I[i-1]); + } + DCOUTVAR(newI); + ex xj = I[j-2].x; int alphaj = I[j-2].evalpoint; size_t deg = A[j-1].degree(xj); @@ -1762,8 +1780,6 @@ ex hensel_multivar(const ex& a, const ex& x, const vector& I, unsigne ex c = dif.subs(xj==alphaj) / factorial(k); DCOUTVAR(c); if ( !c.is_zero() ) { - vector newI = I; - newI.pop_back(); vector deltaU = multivar_diophant(U1, x, c, newI, maxdeg, p, cl_I_to_uint(l)); for ( size_t i=0; i& I, unsigne deltaU[i] *= monomial; U[i] += deltaU[i]; U[i] = make_modular(U[i], R); + U[i] = U[i].expand(); + DCOUTVAR(U[i]); } ex Uprod = 1; for ( size_t i=0; i(e) ) { + if ( is_a(e) || is_a(e) ) { result.append(1); result.append(e); result.append(1); @@ -1862,30 +1882,65 @@ static ex put_factors_into_lst(const ex& e) throw runtime_error("put_factors_into_lst: bad term."); } +#ifdef DEBUGFACTOR +ostream& operator<<(ostream& o, const vector& v) +{ + for ( size_t i=0; i& d) { + DCOUT(checkdivisors); const int k = f.nops()-2; + DCOUTVAR(k); + DCOUTVAR(d.size()); numeric q, r; d[0] = ex_to(f.op(0) * f.op(f.nops()-1)); + if ( d[0] == 1 && k == 1 && abs(f.op(1)) != 1 ) { + DCOUT(false); + DCOUT(END checkdivisors); + return false; + } + DCOUTVAR(d[0]); for ( int i=1; i<=k; ++i ) { - q = ex_to(abs(f.op(i-1))); + DCOUTVAR(i); + DCOUTVAR(abs(f.op(i))); + q = ex_to(abs(f.op(i))); + DCOUTVAR(q); for ( int j=i-1; j>=0; --j ) { r = d[j]; + DCOUTVAR(r); do { r = gcd(r, q); + DCOUTVAR(r); q = q/r; + DCOUTVAR(q); } while ( r != 1 ); if ( q == 1 ) { + DCOUT(true); + DCOUT(END checkdivisors); return true; } } d[i] = q; } + DCOUT(false); + DCOUT(END checkdivisors); return false; } -static void generate_set(const ex& u, const ex& vn, const exset& syms, const ex& f, const numeric& modulus, vector& a, vector& d) +static bool generate_set(const ex& u, const ex& vn, const exset& syms, const ex& f, const numeric& modulus, vector& a, vector& d) { + // computation of d is actually not necessary + DCOUT(generate_set); + DCOUTVAR(u); + DCOUTVAR(vn); + DCOUTVAR(f); + DCOUTVAR(modulus); const ex& x = *syms.begin(); bool trying = true; do { @@ -1895,50 +1950,62 @@ static void generate_set(const ex& u, const ex& vn, const exset& syms, const ex& exset::const_iterator s = syms.begin(); ++s; for ( size_t i=0; i(x))) != 1 ) { continue; } if ( is_a(vn) ) { - d = a; trying = false; } else { + DCOUT(do substitution); lst fnum; lst::const_iterator i = ex_to(f).begin(); fnum.append(*i++); + bool problem = false; while ( i!=ex_to(f).end() ) { ex fs = *i; - s = syms.begin(); - ++s; - for ( size_t j=0; j(fs) ) { + s = syms.begin(); + ++s; + for ( size_t j=0; j& v) -{ - for ( size_t i=0; i(cont) ) { +#ifdef DEBUGFACTOR + return ::factor(cont) * ::factor(pp); +#else return factor(cont) * factor(pp); +#endif } /* factor leading coefficient */ @@ -1973,25 +2044,36 @@ static ex factor_multivariate(const ex& poly, const exset& syms) vnlst = lst(vn); } else { +#ifdef DEBUGFACTOR + ex vnfactors = ::factor(vn); +#else ex vnfactors = factor(vn); +#endif vnlst = put_factors_into_lst(vnfactors); } DCOUTVAR(vnlst); const numeric maxtrials = 3; numeric modulus = (vnlst.nops()-1 > 3) ? vnlst.nops()-1 : 3; + DCOUTVAR(modulus); numeric minimalr = -1; - vector a(syms.size()-1); - vector d(syms.size()-1); + vector a(syms.size()-1, 0); + vector d((vnlst.nops()-1)/2+1, 0); while ( true ) { numeric trialcount = 0; ex u, delta; unsigned int prime; - UniPolyVec uvec; + size_t factor_count; + ex ufac; + ex ufaclst; while ( trialcount < maxtrials ) { - uvec.clear(); - generate_set(pp, vn, syms, vnlst, modulus, a, d); + bool problem = generate_set(pp, vn, syms, vnlst, modulus, a, d); + DCOUTVAR(problem); + if ( problem ) { + ++modulus; + continue; + } DCOUTVAR(a); DCOUTVAR(d); u = pp; @@ -2002,134 +2084,185 @@ static ex factor_multivariate(const ex& poly, const exset& syms) ++s; } delta = u.content(x); + DCOUTVAR(u); // determine proper prime prime = 3; + DCOUTVAR(prime); cl_modint_ring R = find_modint_ring(prime); + DCOUTVAR(u.lcoeff(x)); while ( true ) { if ( irem(ex_to(u.lcoeff(x)), prime) != 0 ) { UniPoly modpoly(R, u, x); UniFactorVec sqrfree_ufv; squarefree(modpoly, sqrfree_ufv); + DCOUTVAR(sqrfree_ufv); if ( sqrfree_ufv.factors.size() == 1 && sqrfree_ufv.factors.front().exp == 1 ) break; } prime = next_prime(prime); + DCOUTVAR(prime); R = find_modint_ring(prime); } - UniPoly umod(R, u, x); - DCOUTVAR(u); - factor_modular(umod, uvec); - DCOUTVAR(uvec); +#ifdef DEBUGFACTOR + ufac = ::factor(u); +#else + ufac = factor(u); +#endif + DCOUTVAR(ufac); + ufaclst = put_factors_into_lst(ufac); + DCOUTVAR(ufaclst); + factor_count = (ufaclst.nops()-1)/2; + DCOUTVAR(factor_count); - if ( uvec.size() == 1 ) { + if ( factor_count <= 1 ) { DCOUTVAR(poly); DCOUT(END factor_multivariate); return poly; } if ( minimalr < 0 ) { - minimalr = uvec.size(); + minimalr = factor_count; } - else if ( minimalr == uvec.size() ) { + else if ( minimalr == factor_count ) { ++trialcount; ++modulus; } - else if ( minimalr > uvec.size() ) { - minimalr = uvec.size(); + else if ( minimalr > factor_count ) { + minimalr = factor_count; trialcount = 0; } DCOUTVAR(trialcount); DCOUTVAR(minimalr); - if ( minimalr == 0 ) { + if ( minimalr <= 1 ) { DCOUTVAR(poly); DCOUT(END factor_multivariate); return poly; } } - vector C; - if ( vnlst.nops() == 1 ) { - C.resize(uvec.size(), 1); + vector ftilde((vnlst.nops()-1)/2+1); + ftilde[0] = ex_to(vnlst.op(0)); + for ( size_t i=1; i(ft); } - else { + DCOUTVAR(ftilde); - vector ftilde((vnlst.nops()-1)/2); - for ( size_t i=0; i used_flag((vnlst.nops()-1)/2+1, false); + vector D(factor_count, 1); + for ( size_t i=0; i<=factor_count; ++i ) { + DCOUTVAR(i); + numeric prefac; + if ( i == 0 ) { + prefac = ex_to(ufaclst.op(0)); + ftilde[0] = ftilde[0] / prefac; + vnlst.let_op(0) = vnlst.op(0) / prefac; + continue; + } + else { + prefac = ex_to(ufaclst.op(2*(i-1)+1).lcoeff(x)); + } + DCOUTVAR(prefac); + for ( size_t j=(vnlst.nops()-1)/2+1; j>0; --j ) { + DCOUTVAR(j); + DCOUTVAR(prefac); + DCOUTVAR(ftilde[j-1]); + if ( abs(ftilde[j-1]) == 1 ) { + used_flag[j-1] = true; + continue; } - ftilde[i] = ex_to(ft); - } - DCOUTVAR(ftilde); - - vector D; - vector fflag(ftilde.size(), false); - for ( size_t i=0; i(ui.lcoeff(x)); - for ( size_t j=0; j 0 ) { + if ( j == 1 ) { + D[i-1] = D[i-1] * pow(vnlst.op(0), count); + } + else { + D[i-1] = D[i-1] * pow(vnlst.op(2*(j-2)+1), count); + } } + else { + ftilde[j-1] = ftilde[j-1] / prefac; + DCOUT(BREAK); + DCOUTVAR(ftilde[j-1]); + break; + } + ++j; } - D.push_back(Di.expand()); } - for ( size_t i=0; i C(factor_count); + DCOUTVAR(C); + DCOUTVAR(delta); + if ( delta == 1 ) { + for ( size_t i=0; i epv; @@ -2140,6 +2273,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms) ep.evalpoint = a[i].to_int(); epv.push_back(ep); } + DCOUTVAR(epv); // calc bound B ex maxcoeff; @@ -2148,22 +2282,30 @@ static ex factor_multivariate(const ex& poly, const exset& syms) } cl_I normmc = ceiling1(the(cln::sqrt(ex_to(maxcoeff).to_cl_N()))); unsigned int maxdegree = 0; - for ( size_t i=0; i maxdegree ) { - maxdegree = uvec[i].degree(); + for ( size_t i=0; i (int)maxdegree ) { + maxdegree = ufaclst[2*i+1].degree(x); } } cl_I B = normmc * expt_pos(cl_I(2), maxdegree); cl_I l = 1; cl_I pl = prime; while ( pl < B ) { - l += 1; + l = l + 1; pl = pl * prime; } - ex res = hensel_multivar(pp, x, epv, prime, l, uvec, C); + UniPolyVec uvec; + cl_modint_ring R = find_modint_ring(expt_pos(cl_I(prime),l)); + for ( size_t i=0; i<(ufaclst.nops()-1)/2; ++i ) { + UniPoly newu(R, ufaclst.op(i*2+1), x); + uvec.push_back(newu); + } + DCOUTVAR(uvec); + + ex res = hensel_multivar(ufaclst.op(0)*pp, x, epv, prime, l, uvec, C); if ( res != lst() ) { - ex result = cont; + ex result = cont * ufaclst.op(0); for ( size_t i=0; i