From: Jens Vollinga Date: Mon, 8 Sep 2008 16:32:27 +0000 (+0200) Subject: Added missing code for multivariate factorization. X-Git-Tag: release_1-5-0~67 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=ccec535ede5a414447d69fba95af194df7327bfb Added missing code for multivariate factorization. --- diff --git a/ginac/factor.cpp b/ginac/factor.cpp index 12405dd6..3afaec54 100644 --- a/ginac/factor.cpp +++ b/ginac/factor.cpp @@ -1660,14 +1660,27 @@ vector multivar_diophant(const vector& a_, const ex& x, const ex& c, con return sigma; } +#ifdef DEBUGFACTOR +ostream& operator<<(ostream& o, const vector& v) +{ + for ( size_t i=0; i& I, unsigned int p, const cl_I& l, const UniPolyVec& u, const vector& lcU) { DCOUT(hensel_multivar); DCOUTVAR(a); DCOUTVAR(x); + DCOUTVAR(I); DCOUTVAR(p); DCOUTVAR(l); DCOUTVAR(u); + DCOUTVAR(lcU); const size_t nu = I.size() + 1; const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),l)); @@ -1783,17 +1796,375 @@ ex hensel_multivar(const ex& a, const ex& x, const vector& I, unsigne for ( size_t i=0; i(e) ) { + result.append(e); + DCOUT(END put_factors_into_lst); + DCOUTVAR(result); + return result; + } + if ( is_a(e) ) { + result.append(1); + result.append(e.op(0)); + result.append(e.op(1)); + DCOUT(END put_factors_into_lst); + DCOUTVAR(result); + return result; + } + if ( is_a(e) ) { + result.append(1); + result.append(e); + result.append(1); + DCOUT(END put_factors_into_lst); + DCOUTVAR(result); + return result; + } + if ( is_a(e) ) { + ex nfac = 1; + for ( size_t i=0; i(op) ) { + nfac = op; + } + if ( is_a(op) ) { + result.append(op.op(0)); + result.append(op.op(1)); + } + if ( is_a(op) || is_a(op) ) { + result.append(op); + result.append(1); + } + } + result.prepend(nfac); + DCOUT(END put_factors_into_lst); + DCOUTVAR(result); + return result; + } + throw runtime_error("put_factors_into_lst: bad term."); +} + +static bool checkdivisors(const lst& f, vector& d) +{ + const int k = f.nops()-2; + numeric q, r; + d[0] = ex_to(f.op(0) * f.op(f.nops()-1)); + for ( int i=1; i<=k; ++i ) { + q = ex_to(abs(f.op(i-1))); + for ( int j=i-1; j>=0; --j ) { + r = d[j]; + do { + r = gcd(r, q); + q = q/r; + } while ( r != 1 ); + if ( q == 1 ) { + return true; + } + } + d[i] = q; + } + 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) +{ + const ex& x = *syms.begin(); + bool trying = true; + do { + ex u0 = u; + ex vna = vn; + ex vnatry; + 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 { + lst fnum; + lst::const_iterator i = ex_to(f).begin(); + fnum.append(*i++); + while ( i!=ex_to(f).end() ) { + ex fs = *i; + s = syms.begin(); + ++s; + for ( size_t j=0; j& v) +{ + for ( size_t i=0; i=p.ldegree(x); --i ) { + cont = gcd(cont, p.coeff(x,ex_to(i).to_int())); + if ( cont == 1 ) break; + } + DCOUTVAR(cont); + ex pp = expand(normal(p / cont)); + DCOUTVAR(pp); + if ( !is_a(cont) ) { + return factor(cont) * factor(pp); + } + + /* factor leading coefficient */ + pp = pp.collect(x); + ex vn = p.lcoeff(x); + ex vnlst; + if ( is_a(vn) ) { + vnlst = lst(vn); + } + else { + ex vnfactors = factor(vn); + vnlst = put_factors_into_lst(vnfactors); + } + DCOUTVAR(vnlst); + + const numeric maxtrials = 3; + numeric modulus = (vnlst.nops()-1 > 3) ? vnlst.nops()-1 : 3; + numeric minimalr = -1; + vector a(syms.size()-1); + vector d(syms.size()-1); + + while ( true ) { + numeric trialcount = 0; + ex u, delta; + unsigned int prime; + UniPolyVec uvec; + while ( trialcount < maxtrials ) { + uvec.clear(); + generate_set(pp, vn, syms, vnlst, modulus, a, d); + DCOUTVAR(a); + DCOUTVAR(d); + u = pp; + s = syms.begin(); + ++s; + for ( size_t i=0; i(u.lcoeff(x)), prime) != 0 ) { + UniPoly modpoly(R, u, x); + UniFactorVec sqrfree_ufv; + squarefree(modpoly, sqrfree_ufv); + if ( sqrfree_ufv.factors.size() == 1 && sqrfree_ufv.factors.front().exp == 1 ) break; + } + prime = next_prime(prime); + R = find_modint_ring(prime); + } + + UniPoly umod(R, u, x); + DCOUTVAR(u); + factor_modular(umod, uvec); + DCOUTVAR(uvec); + + if ( uvec.size() == 1 ) { + DCOUTVAR(poly); + DCOUT(END factor_multivariate); + return poly; + } + + if ( minimalr < 0 ) { + minimalr = uvec.size(); + } + else if ( minimalr == uvec.size() ) { + ++trialcount; + ++modulus; + } + else if ( minimalr > uvec.size() ) { + minimalr = uvec.size(); + trialcount = 0; + } + DCOUTVAR(trialcount); + DCOUTVAR(minimalr); + if ( minimalr == 0 ) { + DCOUTVAR(poly); + DCOUT(END factor_multivariate); + return poly; + } + } + + vector C; + if ( vnlst.nops() == 1 ) { + C.resize(uvec.size(), 1); + } + else { + + vector ftilde((vnlst.nops()-1)/2); + for ( size_t i=0; i(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 epv; + s = syms.begin(); + ++s; + for ( size_t i=0; i=u.ldegree(x); --i ) { + maxcoeff += pow(abs(u.coeff(x, i)),2); + } + 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(); + } + } + unsigned int B = cl_I_to_uint(normmc * expt_pos(cl_I(2), maxdegree)); + + ex res = hensel_multivar(poly, x, epv, prime, B, uvec, C); + if ( res != lst() ) { + ex result = cont; + for ( size_t i=0; i 0 ) { int ld = poly.ldegree(x); @@ -1818,8 +2190,9 @@ static ex factor_sqrfree(const ex& poly) } } - // multivariate case not yet implemented! - throw runtime_error("multivariate case not yet implemented!"); + // multivariate case + ex res = factor_multivariate(poly, findsymbols.syms); + return res; } } // anonymous namespace