From: Jens Vollinga Date: Mon, 10 Nov 2008 12:38:11 +0000 (+0100) Subject: Added modular square free factorization. X-Git-Tag: release_1-5-0~26 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=edc92b7a463993da62357fb4afad053e8c6d0771 Added modular square free factorization. Completed distinct degree factorization. Univariate polynomial factorization uses now upoly. Merged class Partition and function split into class factor_partition. --- diff --git a/ginac/factor.cpp b/ginac/factor.cpp index f3b48fdb..204010be 100644 --- a/ginac/factor.cpp +++ b/ginac/factor.cpp @@ -72,11 +72,29 @@ namespace { typedef vector mvec; #ifdef DEBUGFACTOR +ostream& operator<<(ostream& o, const vector& v) +{ + vector::const_iterator i = v.begin(), end = v.end(); + while ( i != end ) { + o << *i++ << " "; + } + return o; +} +ostream& operator<<(ostream& o, const vector& v) +{ + vector::const_iterator i = v.begin(), end = v.end(); + while ( i != end ) { + o << *i << "[" << i-v.begin() << "]" << " "; + ++i; + } + return o; +} ostream& operator<<(ostream& o, const vector& v) { vector::const_iterator i = v.begin(), end = v.end(); while ( i != end ) { - o << *i++ << " "; + o << *i << "[" << i-v.begin() << "]" << " "; + ++i; } return o; } @@ -84,7 +102,8 @@ ostream& operator<<(ostream& o, const vector< vector >& v) { vector< vector >::const_iterator i = v.begin(), end = v.end(); while ( i != end ) { - o << *i++ << endl; + o << i-v.begin() << ": " << *i << endl; + ++i; } return o; } @@ -161,29 +180,16 @@ canonicalize(T& p, const typename T::size_type hint = std::numeric_limitszero()); -// cl_MI norm = recip(a[0]); -// umodpoly an = a; -// for ( size_t i=0; ione(); -// for ( size_t i=1; izero(); + int deg = degree(a); + a.resize(degree(a)*q+1, zero); + for ( int i=deg; i>0; --i ) { + a[i*q] = a[i]; + a[i] = zero; + } } template @@ -410,6 +416,20 @@ static upoly umodpoly_to_upoly(const umodpoly& a) return e; } +static umodpoly umodpoly_to_umodpoly(const umodpoly& a, const cl_modint_ring& R, unsigned int m) +{ + umodpoly e; + if ( a.empty() ) return e; + cl_modint_ring oldR = a[0].ring(); + size_t sa = a.size(); + e.resize(sa+m, R->zero()); + for ( int i=0; icanonhom(oldR->retract(a[i])); + } + canonicalize(e); + return e; +} + /** Divides all coefficients of the polynomial a by the integer x. * All coefficients are supposed to be divisible by x. If they are not, the * the cast will raise an exception. @@ -877,170 +897,153 @@ static void berlekamp(const umodpoly& a, upvec& upv) } } -static void rem_xq(int q, const umodpoly& b, umodpoly& c) +static void expt_1_over_p(const umodpoly& a, unsigned int prime, umodpoly& ap) { - cl_modint_ring R = b[0].ring(); - - int n = degree(b); - if ( n > q ) { - c.resize(q+1, R->zero()); - c[q] = R->one(); - return; + size_t newdeg = degree(a)/prime; + ap.resize(newdeg+1); + ap[0] = a[0]; + for ( size_t i=1; i<=newdeg; ++i ) { + ap[i] = a[i*prime]; } +} - c.clear(); - c.resize(n+1, R->zero()); - int k = q-n; - c[n] = R->one(); - - int ofs = 0; - do { - cl_MI qk = div(c[n-ofs], b[n]); - if ( !zerop(qk) ) { - for ( int i=1; i<=n; ++i ) { - c[n-i+ofs] = c[n-i] - qk * b[n-i]; +static void modsqrfree(const umodpoly& a, upvec& factors, vector& mult) +{ + const unsigned int prime = cl_I_to_uint(a[0].ring()->modulus); + int i = 1; + umodpoly b; + deriv(a, b); + if ( b.size() ) { + umodpoly c; + gcd(a, b, c); + umodpoly w; + div(a, c, w); + while ( unequal_one(w) ) { + umodpoly y; + gcd(w, c, y); + umodpoly z; + div(w, y, z); + factors.push_back(z); + mult.push_back(i); + ++i; + w = y; + umodpoly buf; + div(c, y, buf); + c = buf; + } + if ( unequal_one(c) ) { + umodpoly cp; + expt_1_over_p(c, prime, cp); + size_t previ = mult.size(); + modsqrfree(cp, factors, mult); + for ( size_t i=previ; i& degrees, upvec& ddfactors) { umodpoly a = a_; cl_modint_ring R = a[0].ring(); int q = cl_I_to_int(R->modulus); - int n = degree(a); - size_t nhalf = n/2; + int nhalf = degree(a)/2; - size_t i = 1; - umodpoly w(1, R->one()); + int i = 1; + umodpoly w(2); + w[0] = R->zero(); + w[1] = R->one(); umodpoly x = w; - upvec ai; - + bool nontrivial = false; while ( i <= nhalf ) { - expt_pos(w, q, w); - rem(w, a, w); - + expt_pos(w, q); umodpoly buf; - gcd(a, w-x, buf); - ai.push_back(buf); - - if ( unequal_one(ai.back()) ) { - div(a, ai.back(), a); - rem(w, a, w); + rem(w, a, buf); + w = buf; + umodpoly wx = w - x; + gcd(a, wx, buf); + if ( unequal_one(buf) ) { + degrees.push_back(i); + ddfactors.push_back(buf); + } + if ( unequal_one(buf) ) { + umodpoly buf2; + div(a, buf, buf2); + a = buf2; + nhalf = degree(a)/2; + rem(w, a, buf); + w = buf; } - ++i; } - - result = ai; + if ( unequal_one(a) ) { + degrees.push_back(degree(a)); + ddfactors.push_back(a); + } } -static void same_degree_factor(const umodpoly& a, upvec& result) +static void same_degree_factor(const umodpoly& a, upvec& upv) { cl_modint_ring R = a[0].ring(); int deg = degree(a); - upvec buf; - distinct_degree_factor(a, buf); - int degsum = 0; - - for ( size_t i=0; i degrees; + upvec ddfactors; + distinct_degree_factor(a, degrees, ddfactors); - if ( degsum < deg ) { - result.push_back(a); + for ( size_t i=0; imodulus); - int n = degree(a); - - cl_N pm = 0.3; - int l = cl_I_to_int(ceiling1(the(expt(n, pm)))); - upvec h(l+1); - umodpoly qk(1, R->one()); - h[0] = qk; - for ( int i=1; i<=l; ++i ) { - expt_pos(h[i-1], q, qk); - rem(qk, a, h[i]); - } - - int m = std::ceil(((double)n)/2/l); - upvec H(m); - int ql = std::pow(q, l); - H[0] = h[l]; - for ( int i=1; i mult; + modsqrfree(p, factors, mult); - upvec I(m); - umodpoly one(1, R->one()); - for ( int i=0; i0; --j ) { + upv.insert(upv.end(), upvbuf.begin(), upvbuf.end()); } - rem(I[i], a, I[i]); - } - - upvec F(m, one); - umodpoly f = a; - for ( int i=0; i=0; --j ) { - umodpoly g; - gcd(f, H[i]-h[j], g); - result[l*(i+1)-j-1] = g; - div(f, g, f); +#else + for ( size_t i=0; i0; --j ) { + upv.push_back(factors[i]); + } } } -} - -static void cantor_zassenhaus(const umodpoly& a, upvec& result) -{ -} - -static void factor_modular(const umodpoly& p, upvec& upv) -{ - //same_degree_factor(p, upv); - berlekamp(p, upv); - return; +#endif } /** Calculates polynomials s and t such that a*s+b*t==1. @@ -1101,14 +1104,13 @@ static upoly replace_lc(const upoly& poly, const cl_I& lc) return r; } -static ex hensel_univar(const ex& a_, const ex& x, unsigned int p, const umodpoly& u1_, const umodpoly& w1_, const ex& gamma_ = 0) +static void hensel_univar(const upoly& a_, unsigned int p, const umodpoly& u1_, const umodpoly& w1_, upoly& u, upoly& w) { - upoly a; - upoly_from_ex(a, a_, x); + upoly a = a_; const cl_modint_ring& R = u1_[0].ring(); // calc bound B - cl_R maxcoeff; + cl_R maxcoeff = 0; for ( int i=degree(a); i>=0; --i ) { maxcoeff = maxcoeff + square(abs(a[i])); } @@ -1118,18 +1120,13 @@ static ex hensel_univar(const ex& a_, const ex& x, unsigned int p, const umodpol // step 1 cl_I alpha = lcoeff(a); - cl_I gamma = the(ex_to(gamma_).to_cl_N()); - if ( gamma == 0 ) { - gamma = alpha; - } - cl_I gamma_ui = abs(gamma); - a = a * gamma; + a = a * alpha; umodpoly nu1 = u1_; normalize_in_field(nu1); umodpoly nw1 = w1_; normalize_in_field(nw1); upoly phi; - phi = umodpoly_to_upoly(nu1) * gamma; + phi = umodpoly_to_upoly(nu1) * alpha; umodpoly u1; umodpoly_from_upoly(u1, phi, R); phi = umodpoly_to_upoly(nw1) * alpha; @@ -1142,14 +1139,20 @@ static ex hensel_univar(const ex& a_, const ex& x, unsigned int p, const umodpol exteuclid(u1, w1, s, t); // step 3 - upoly u = replace_lc(umodpoly_to_upoly(u1), gamma); - upoly w = replace_lc(umodpoly_to_upoly(w1), alpha); + u = replace_lc(umodpoly_to_upoly(u1), alpha); + w = replace_lc(umodpoly_to_upoly(w1), alpha); upoly e = a - u * w; cl_I modulus = p; - const cl_I maxmodulus = 2*B*gamma_ui; + const cl_I maxmodulus = 2*B*abs(alpha); // step 4 while ( !e.empty() && modulus < maxmodulus ) { + // ad-hoc divisablity check + for ( size_t k=0; kone()); + split(); } int operator[](size_t i) const { return k[i]; } size_t size() const { return n; } size_t size_first() const { return n-sum; } size_t size_second() const { return sum; } #ifdef DEBUGFACTOR - void get() const - { - for ( size_t i=0; i 0; + if ( sum > 0 ) { + split(); + return true; + } + else { + return false; + } } ++k[i]; ++sum; } return false; } + void split() + { + left = one; + right = one; + for ( size_t i=0; i k; }; -static void split(const upvec& factors, const Partition& part, umodpoly& a, umodpoly& b) -{ - umodpoly one(1, factors.front()[0].ring()->one()); - a = one; - b = one; - for ( size_t i=0; i(prim.lcoeff(x)); + cl_I lc = lcoeff(prim); upvec factors; while ( trials < 2 ) { + umodpoly modpoly; while ( true ) { p = next_prime(p); - if ( irem(lcoeff, p) != 0 ) { + if ( !zerop(rem(lc, p)) ) { R = find_modint_ring(p); - umodpoly modpoly; - umodpoly_from_ex(modpoly, prim, x, R); + umodpoly_from_upoly(modpoly, prim, R); if ( squarefree(modpoly) ) break; } } // do modular factorization - umodpoly modpoly; - umodpoly_from_ex(modpoly, prim, x, R); upvec trialfactors; factor_modular(modpoly, trialfactors); if ( trialfactors.size() <= 1 ) { @@ -1319,7 +1334,6 @@ static ex factor_univariate(const ex& poly, const ex& x) } p = lastp; R = find_modint_ring(p); - cl_univpoly_modint_ring UPR = find_univpoly_ring(R); // lift all factor combinations stack tocheck; @@ -1327,24 +1341,22 @@ static ex factor_univariate(const ex& poly, const ex& x) mf.poly = prim; mf.factors = factors; tocheck.push(mf); + upoly f1, f2; ex result = 1; while ( tocheck.size() ) { const size_t n = tocheck.top().factors.size(); - Partition part(n); + factor_partition part(tocheck.top().factors); while ( true ) { - umodpoly a, b; - split(tocheck.top().factors, part, a, b); - - ex answer = hensel_univar(tocheck.top().poly, x, p, a, b); - if ( answer != lst() ) { + hensel_univar(tocheck.top().poly, p, part.left, part.right, f1, f2); + if ( !f1.empty() ) { if ( part.size_first() == 1 ) { if ( part.size_second() == 1 ) { - result *= answer.op(0) * answer.op(1); + result *= upoly_to_ex(f1, x) * upoly_to_ex(f2, x); tocheck.pop(); break; } - result *= answer.op(0); - tocheck.top().poly = answer.op(1); + result *= upoly_to_ex(f1, x); + tocheck.top().poly = f2; for ( size_t i=0; i 2 ) { upvec s = multiterm_eea_lift(a, x, p, k); for ( size_t j=0; j