From 2b0ad5c381dc081cc4066b0c5f939f5169ad9ab3 Mon Sep 17 00:00:00 2001 From: Jens Vollinga Date: Thu, 13 Nov 2008 11:13:05 +0100 Subject: [PATCH 1/1] Fixed lots of bugs in factor_multivariate(). factor_univariate() takes now the content of the polynomial into accout when choosing a prime. It also returns the chosen prime. This supports the factor_multivariate() code. Polynomials are expanded before calling factor_sqrfree(). This avoids problems with some input polynomials. Added static qualifiers to all hidden functions. --- ginac/factor.cpp | 427 +++++++++++++++++++++-------------------------- 1 file changed, 187 insertions(+), 240 deletions(-) diff --git a/ginac/factor.cpp b/ginac/factor.cpp index ecf537b0..fb738972 100644 --- a/ginac/factor.cpp +++ b/ginac/factor.cpp @@ -365,10 +365,12 @@ static void umodpoly_from_ex(umodpoly& ump, const ex& e, const ex& x, const cl_m canonicalize(ump); } +#ifdef DEBUGFACTOR static void umodpoly_from_ex(umodpoly& ump, const ex& e, const ex& x, const cl_I& modulus) { umodpoly_from_ex(ump, e, x, find_modint_ring(modulus)); } +#endif static ex upoly_to_ex(const upoly& a, const ex& x) { @@ -1337,7 +1339,7 @@ struct ModFactors upvec factors; }; -static ex factor_univariate(const ex& poly, const ex& x) +static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime) { ex unit, cont, prim_ex; poly.unitcontprim(x, unit, cont, prim_ex); @@ -1345,18 +1347,19 @@ static ex factor_univariate(const ex& poly, const ex& x) upoly_from_ex(prim, prim_ex, x); // determine proper prime and minimize number of modular factors - unsigned int p = 3, lastp = 3; + prime = 3; + unsigned int lastp = prime; cl_modint_ring R; unsigned int trials = 0; unsigned int minfactors = 0; - cl_I lc = lcoeff(prim); + cl_I lc = lcoeff(prim) * the(ex_to(cont).to_cl_N()); upvec factors; while ( trials < 2 ) { umodpoly modpoly; while ( true ) { - p = next_prime(p); - if ( !zerop(rem(lc, p)) ) { - R = find_modint_ring(p); + prime = next_prime(prime); + if ( !zerop(rem(lc, prime)) ) { + R = find_modint_ring(prime); umodpoly_from_upoly(modpoly, prim, R); if ( squarefree(modpoly) ) break; } @@ -1373,15 +1376,15 @@ static ex factor_univariate(const ex& poly, const ex& x) if ( minfactors == 0 || trialfactors.size() < minfactors ) { factors = trialfactors; minfactors = trialfactors.size(); - lastp = p; + lastp = prime; trials = 1; } else { ++trials; } } - p = lastp; - R = find_modint_ring(p); + prime = lastp; + R = find_modint_ring(prime); // lift all factor combinations stack tocheck; @@ -1395,7 +1398,7 @@ static ex factor_univariate(const ex& poly, const ex& x) const size_t n = tocheck.top().factors.size(); factor_partition part(tocheck.top().factors); while ( true ) { - hensel_univar(tocheck.top().poly, p, part.left(), part.right(), f1, f2); + hensel_univar(tocheck.top().poly, prime, part.left(), part.right(), f1, f2); if ( !f1.empty() ) { if ( part.size_left() == 1 ) { if ( part.size_right() == 1 ) { @@ -1462,6 +1465,12 @@ static ex factor_univariate(const ex& poly, const ex& x) return unit * cont * result; } +static inline ex factor_univariate(const ex& poly, const ex& x) +{ + unsigned int prime; + return factor_univariate(poly, x, prime); +} + struct EvalPoint { ex x; @@ -1469,9 +1478,9 @@ struct EvalPoint }; // forward declaration -vector multivar_diophant(const vector& a_, const ex& x, const ex& c, const vector& I, unsigned int d, unsigned int p, unsigned int k); +static vector multivar_diophant(const vector& a_, const ex& x, const ex& c, const vector& I, unsigned int d, unsigned int p, unsigned int k); -upvec multiterm_eea_lift(const upvec& a, const ex& x, unsigned int p, unsigned int k) +static upvec multiterm_eea_lift(const upvec& a, const ex& x, unsigned int p, unsigned int k) { const size_t r = a.size(); cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k)); @@ -1502,7 +1511,7 @@ upvec multiterm_eea_lift(const upvec& a, const ex& x, unsigned int p, unsigned i /** * Assert: a not empty. */ -void change_modulus(const cl_modint_ring& R, umodpoly& a) +static void change_modulus(const cl_modint_ring& R, umodpoly& a) { if ( a.empty() ) return; cl_modint_ring oldR = a[0].ring(); @@ -1513,7 +1522,7 @@ void change_modulus(const cl_modint_ring& R, umodpoly& a) canonicalize(a); } -void eea_lift(const umodpoly& a, const umodpoly& b, const ex& x, unsigned int p, unsigned int k, umodpoly& s_, umodpoly& t_) +static void eea_lift(const umodpoly& a, const umodpoly& b, const ex& x, unsigned int p, unsigned int k, umodpoly& s_, umodpoly& t_) { cl_modint_ring R = find_modint_ring(p); umodpoly amod = a; @@ -1556,7 +1565,7 @@ void eea_lift(const umodpoly& a, const umodpoly& b, const ex& x, unsigned int p, s_ = s; t_ = t; } -upvec univar_diophant(const upvec& a, const ex& x, unsigned int m, unsigned int p, unsigned int k) +static upvec univar_diophant(const upvec& a, const ex& x, unsigned int m, unsigned int p, unsigned int k) { cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k)); @@ -1616,7 +1625,8 @@ static ex make_modular(const ex& e, const cl_modint_ring& R) return map(e.expand()); } -vector multivar_diophant(const vector& a_, const ex& x, const ex& c, const vector& I, unsigned int d, unsigned int p, unsigned int k) +static vector multivar_diophant(const vector& a_, const ex& x, const ex& c, const vector& I, + unsigned int d, unsigned int p, unsigned int k) { vector a = a_; @@ -1727,7 +1737,7 @@ ostream& operator<<(ostream& o, const vector& v) } #endif // def DEBUGFACTOR -ex hensel_multivar(const ex& a, const ex& x, const vector& I, unsigned int p, const cl_I& l, const upvec& u, const vector& lcU) +static ex hensel_multivar(const ex& a, const ex& x, const vector& I, unsigned int p, const cl_I& l, const upvec& u, const vector& lcU) { const size_t nu = I.size() + 1; const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),l)); @@ -1834,13 +1844,11 @@ static ex put_factors_into_lst(const ex& e) if ( is_a(e) ) { result.append(1); result.append(e.op(0)); - result.append(e.op(1)); return result; } if ( is_a(e) || is_a(e) ) { result.append(1); result.append(e); - result.append(1); return result; } if ( is_a(e) ) { @@ -1852,11 +1860,9 @@ static ex put_factors_into_lst(const ex& e) } 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); @@ -1875,15 +1881,18 @@ ostream& operator<<(ostream& o, const vector& v) } #endif // def DEBUGFACTOR -static bool checkdivisors(const lst& f, vector& d) +/** Checks whether in a set of numbers each has a unique prime factor. + * + * @param[in] f list of numbers to check + * @return true: if number set is bad, false: otherwise + */ +static bool checkdivisors(const lst& f) { - const int k = f.nops()-2; + const int k = f.nops(); 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 ) { - return false; - } - for ( int i=1; i<=k; ++i ) { + vector d(k); + d[0] = ex_to(abs(f.op(0))); + for ( int i=1; i(abs(f.op(i))); for ( int j=i-1; j>=0; --j ) { r = d[j]; @@ -1900,13 +1909,30 @@ static bool checkdivisors(const lst& f, vector& d) return false; } -static bool generate_set(const ex& u, const ex& vn, const exset& syms, const ex& f, const numeric& modulus, vector& a, vector& d) +/** Generates a set of evaluation points for a multivariate polynomial. + * The set fulfills the following conditions: + * 1. lcoeff(evaluated_polynomial) does not vanish + * 2. factors of lcoeff(evaluated_polynomial) have each a unique prime factor + * 3. evaluated_polynomial is square free + * See [W1] for more details. + * + * @param[in] u multivariate polynomial to be factored + * @param[in] vn leading coefficient of u in x (x==first symbol in syms) + * @param[in] syms set of symbols that appear in u + * @param[in] f lst containing the factors of the leading coefficient vn + * @param[in,out] modulus integer modulus for random number generation (i.e. |a_i| < modulus) + * @param[out] u0 returns the evaluated (univariate) polynomial + * @param[out] a returns the valid evaluation points. must have initial size equal + * number of symbols-1 before calling generate_set + */ +static void generate_set(const ex& u, const ex& vn, const exset& syms, const lst& f, + numeric& modulus, ex& u0, vector& a) { - // computation of d is actually not necessary const ex& x = *syms.begin(); - bool trying = true; - do { - ex u0 = u; + while ( true ) { + ++modulus; + /* generate a set of integers ... */ + u0 = u; ex vna = vn; ex vnatry; exset::const_iterator s = syms.begin(); @@ -1915,71 +1941,64 @@ static bool generate_set(const ex& u, const ex& vn, const exset& syms, const ex& do { a[i] = mod(numeric(rand()), 2*modulus) - modulus; vnatry = vna.subs(*s == a[i]); + /* ... for which the leading coefficient doesn't vanish ... */ } while ( vnatry == 0 ); vna = vnatry; u0 = u0.subs(*s == a[i]); ++s; } - if ( gcd(u0,u0.diff(ex_to(x))) != 1 ) { + /* ... for which u0 is square free ... */ + ex g = gcd(u0, u0.diff(ex_to(x))); + if ( !is_a(g) ) { continue; } - if ( is_a(vn) ) { - trying = false; - } - else { - 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; - if ( !is_a(fs) ) { + if ( !is_a(vn) ) { + /* ... and for which the evaluated factors have each an unique prime factor */ + lst fnum = f; + fnum.let_op(0) = fnum.op(0) * u0.content(x); + for ( size_t i=1; i(fnum.op(i)) ) { s = syms.begin(); ++s; - for ( size_t j=0; j=p.ldegree(x); --i ) { - cont = gcd(cont, p.coeff(x,ex_to(i).to_int())); + for ( int i=p.degree(x)-1; i>=p.ldegree(x); --i ) { + cont = gcd(cont, p.coeff(x,i)); if ( cont == 1 ) break; } ex pp = expand(normal(p / cont)); if ( !is_a(cont) ) { - return factor(cont) * factor(pp); + return factor_sqrfree(cont) * factor_sqrfree(pp); } /* factor leading coefficient */ - pp = pp.collect(x); - ex vn = pp.lcoeff(x); - pp = pp.expand(); + ex vn = pp.collect(x).lcoeff(x); ex vnlst; if ( is_a(vn) ) { vnlst = lst(vn); @@ -1989,200 +2008,131 @@ static ex factor_multivariate(const ex& poly, const exset& syms) vnlst = put_factors_into_lst(vnfactors); } - const numeric maxtrials = 3; - numeric modulus = (vnlst.nops()-1 > 3) ? vnlst.nops()-1 : 3; - numeric minimalr = -1; + const unsigned int maxtrials = 3; + numeric modulus = (vnlst.nops() > 3) ? vnlst.nops() : 3; vector a(syms.size()-1, 0); - vector d((vnlst.nops()-1)/2+1, 0); + /* try now to factorize until we are successful */ while ( true ) { - numeric trialcount = 0; + + unsigned int trialcount = 0; + unsigned int prime; + int factor_count = 0; + int min_factor_count = -1; ex u, delta; - unsigned int prime = 3; - size_t factor_count = 0; - ex ufac; - ex ufaclst; + ex ufac, ufaclst; + + /* try several evaluation points to reduce the number of modular factors */ while ( trialcount < maxtrials ) { - bool problem = generate_set(pp, vn, syms, vnlst, modulus, a, d); - if ( problem ) { - ++modulus; - continue; - } - u = pp; - s = syms.begin(); - ++s; - for ( size_t i=0; i(u.lcoeff(x)), prime) != 0 ) { - umodpoly modpoly; - umodpoly_from_ex(modpoly, u, x, R); - if ( squarefree(modpoly) ) break; - } - prime = next_prime(prime); - R = find_modint_ring(prime); - } - ufac = factor(u); + /* generate a set of valid evaluation points */ + generate_set(pp, vn, syms, ex_to(vnlst), modulus, u, a); + + ufac = factor_univariate(u, x, prime); ufaclst = put_factors_into_lst(ufac); - factor_count = (ufaclst.nops()-1)/2; - - // veto factorization for which gcd(u_i, u_j) != 1 for all i,j - upvec tryu; - for ( size_t i=0; i<(ufaclst.nops()-1)/2; ++i ) { - umodpoly newu; - umodpoly_from_ex(newu, ufaclst.op(i*2+1), x, R); - tryu.push_back(newu); - } - bool veto = false; - for ( size_t i=0; i factor_count ) { - minimalr = factor_count; + else if ( min_factor_count > factor_count ) { + /* new minimum, reset trial counter */ + min_factor_count = factor_count; trialcount = 0; } - if ( minimalr <= 1 ) { - return poly; - } } - vector ftilde((vnlst.nops()-1)/2+1); - ftilde[0] = ex_to(vnlst.op(0)); - for ( size_t i=1; i C(factor_count); + if ( is_a(vn) ) { + for ( size_t i=1; i(ft); } + else { + vector ftilde(vnlst.nops()-1); + for ( size_t i=0; i(ft); + } - vector used_flag((vnlst.nops()-1)/2+1, false); - vector D(factor_count, 1); - for ( size_t i=0; i<=factor_count; ++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; + vector used_flag(ftilde.size(), false); + vector D(factor_count, 1); + if ( delta == 1 ) { + for ( int i=0; i(ufaclst.op(i+1).lcoeff(x)); + for ( int j=ftilde.size()-1; j>=0; --j ) { + int count = 0; + while ( irem(prefac, ftilde[j]) == 0 ) { + prefac = iquo(prefac, ftilde[j]); + ++count; + } + if ( count ) { + used_flag[j] = true; + D[i] = D[i] * pow(vnlst.op(j+1), count); + } + } + C[i] = D[i] * prefac; + } } else { - prefac = ex_to(ufaclst.op(2*(i-1)+1).lcoeff(x)); - } - for ( size_t j=(vnlst.nops()-1)/2+1; j>0; --j ) { - if ( abs(ftilde[j-1]) == 1 ) { - used_flag[j-1] = true; - continue; - } - numeric g = gcd(prefac, ftilde[j-1]); - if ( g != 1 ) { - prefac = prefac / g; - numeric count = abs(iquo(g, ftilde[j-1])); - used_flag[j-1] = true; - if ( i > 0 ) { - if ( j == 1 ) { - D[i-1] = D[i-1] * pow(vnlst.op(0), count); + for ( int i=0; i(ufaclst.op(i+1).lcoeff(x)); + for ( int j=ftilde.size()-1; j>=0; --j ) { + int count = 0; + while ( irem(prefac, ftilde[j]) == 0 ) { + prefac = iquo(prefac, ftilde[j]); + ++count; } - else { - D[i-1] = D[i-1] * pow(vnlst.op(2*(j-2)+1), count); + while ( irem(ex_to(delta)*prefac, ftilde[j]) == 0 ) { + numeric g = gcd(prefac, ex_to(ftilde[j])); + prefac = iquo(prefac, g); + delta = delta / (ftilde[j]/g); + ufaclst.let_op(i+1) = ufaclst.op(i+1) * (ftilde[j]/g); + ++count; + } + if ( count ) { + used_flag[j] = true; + D[i] = D[i] * pow(vnlst.op(j+1), count); } } - else { - ftilde[j-1] = ftilde[j-1] / prefac; - break; - } - ++j; + C[i] = D[i] * prefac; } } - } - - bool some_factor_unused = false; - for ( size_t i=0; i C(factor_count); - if ( delta == 1 ) { - for ( size_t i=0; i epv; s = syms.begin(); @@ -2195,9 +2145,9 @@ static ex factor_multivariate(const ex& poly, const exset& syms) // calc bound p^l int maxdeg = 0; - for ( size_t i=0; i maxdeg ) { - maxdeg = ufaclst[2*i+1].degree(x); + for ( int i=1; i<=factor_count; ++i ) { + if ( ufaclst.op(i).degree(x) > maxdeg ) { + maxdeg = ufaclst[i].degree(x); } } cl_I B = 2*calc_bound(u, x, maxdeg); @@ -2207,18 +2157,15 @@ static ex factor_multivariate(const ex& poly, const exset& syms) l = l + 1; pl = pl * prime; } - - upvec 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 ) { - umodpoly newu; - umodpoly_from_ex(newu, ufaclst.op(i*2+1), x, R); - uvec.push_back(newu); + upvec modfactors(ufaclst.nops()-1); + for ( size_t i=1; i(sfpoly) ) { -- 2.44.0