X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=ginac%2Ffactor.cpp;h=24f828cb793e9ad9d11e5104a4904d64345cd993;hb=b23a2dfdcd478b00a11c6c2184fbe197a3ac5bfb;hp=5ea0520d268f513c9f5d16832bac4a117ce988c2;hpb=931ece4838fcef7b9ec71761d224e2e1cd11a89d;p=ginac.git diff --git a/ginac/factor.cpp b/ginac/factor.cpp index 5ea0520d..24f828cb 100644 --- a/ginac/factor.cpp +++ b/ginac/factor.cpp @@ -1,16 +1,39 @@ /** @file factor.cpp * - * Polynomial factorization code (implementation). + * Polynomial factorization (implementation). + * + * The interface function factor() at the end of this file is defined in the + * GiNaC namespace. All other utility functions and classes are defined in an + * additional anonymous namespace. + * + * Factorization starts by doing a square free factorization and making the + * coefficients integer. Then, depending on the number of free variables it + * proceeds either in dedicated univariate or multivariate factorization code. + * + * Univariate factorization does a modular factorization via Berlekamp's + * algorithm and distinct degree factorization. Hensel lifting is used at the + * end. + * + * Multivariate factorization uses the univariate factorization (applying a + * evaluation homomorphism first) and Hensel lifting raises the answer to the + * multivariate domain. The Hensel lifting code is completely distinct from the + * code used by the univariate factorization. * * Algorithms used can be found in - * [W1] An Improved Multivariate Polynomial Factoring Algorithm, - * P.S.Wang, Mathematics of Computation, Vol. 32, No. 144 (1978) 1215--1231. + * [Wan] An Improved Multivariate Polynomial Factoring Algorithm, + * P.S.Wang, + * Mathematics of Computation, Vol. 32, No. 144 (1978) 1215--1231. * [GCL] Algorithms for Computer Algebra, - * K.O.Geddes, S.R.Czapor, G.Labahn, Springer Verlag, 1992. + * K.O.Geddes, S.R.Czapor, G.Labahn, + * Springer Verlag, 1992. + * [Mig] Some Useful Bounds, + * M.Mignotte, + * In "Computer Algebra, Symbolic and Algebraic Computation" (B.Buchberger et al., eds.), + * pp. 259-263, Springer-Verlag, New York, 1982. */ /* - * GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2011 Johannes Gutenberg University Mainz, Germany * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -27,13 +50,8 @@ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ -#define DEBUGFACTOR +//#define DEBUGFACTOR -#ifdef DEBUGFACTOR -#include -#include -using namespace GiNaC; -#else #include "factor.h" #include "ex.h" @@ -46,70 +64,82 @@ using namespace GiNaC; #include "mul.h" #include "normal.h" #include "add.h" -#endif #include #include #include #include #include +#ifdef DEBUGFACTOR +#include +#endif using namespace std; #include using namespace cln; -#ifdef DEBUGFACTOR -namespace Factor { -#else namespace GiNaC { -#endif #ifdef DEBUGFACTOR #define DCOUT(str) cout << #str << endl #define DCOUTVAR(var) cout << #var << ": " << var << endl #define DCOUT2(str,var) cout << #str << ": " << var << endl -#else -#define DCOUT(str) -#define DCOUTVAR(var) -#define DCOUT2(str,var) -#endif - -// forward declaration -ex factor(const ex& poly, unsigned options); - -// anonymous namespace to hide all utility functions -namespace { - -typedef vector mvec; - -#ifdef DEBUGFACTOR -ostream& operator<<(ostream& o, const mvec& v) +ostream& operator<<(ostream& o, const vector& v) { - mvec::const_iterator i = v.begin(), end = v.end(); + vector::const_iterator i = v.begin(), end = v.end(); while ( i != end ) { o << *i++ << " "; } return o; } -#endif // def DEBUGFACTOR - -#ifdef DEBUGFACTOR -ostream& operator<<(ostream& o, const vector& v) +static ostream& operator<<(ostream& o, const vector& v) { - vector::const_iterator i = v.begin(), end = v.end(); + vector::const_iterator i = v.begin(), end = v.end(); while ( i != end ) { - o << *i++ << endl; + o << *i << "[" << i-v.begin() << "]" << " "; + ++i; } return o; } +static 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) +{ + for ( size_t i=0; i >& v) +{ + vector< vector >::const_iterator i = v.begin(), end = v.end(); + while ( i != end ) { + o << i-v.begin() << ": " << *i << endl; + ++i; + } + return o; +} +#else +#define DCOUT(str) +#define DCOUTVAR(var) +#define DCOUT2(str,var) #endif // def DEBUGFACTOR +// anonymous namespace to hide all utility functions +namespace { + //////////////////////////////////////////////////////////////////////////////// // modular univariate polynomial code -//typedef cl_UP_MI umod; typedef std::vector umodpoly; -//typedef vector umodvec; +typedef std::vector upoly; typedef vector upvec; // COPY FROM UPOLY.HPP @@ -176,37 +206,49 @@ 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 struct enable_if +{ + typedef T type; +}; + +template struct enable_if { /* empty */ }; + +template struct uvar_poly_p +{ + static const bool value = false; +}; + +template<> struct uvar_poly_p +{ + static const bool value = true; +}; + +template<> struct uvar_poly_p +{ + static const bool value = true; +}; + +template +// Don't define this for anything but univariate polynomials. +static typename enable_if::value, T>::type +operator+(const T& a, const T& b) { int sa = a.size(); int sb = b.size(); if ( sa >= sb ) { - umodpoly r(sa); + T r(sa); int i = 0; for ( ; i +// Don't define this for anything but univariate polynomials. Otherwise +// overload resolution might fail (this actually happens when compiling +// GiNaC with g++ 3.4). +static typename enable_if::value, T>::type +operator-(const T& a, const T& b) { int sa = a.size(); int sb = b.size(); if ( sa >= sb ) { - umodpoly r(sa); + T r(sa); int i = 0; for ( ; i degree(a) || (i-j) > degree(b) ) continue; + c[i] = c[i] + a[j] * b[i-j]; + } + } + canonicalize(c); + return c; +} + static umodpoly operator*(const umodpoly& a, const umodpoly& b) { umodpoly c; @@ -270,7 +334,7 @@ static umodpoly operator*(const umodpoly& a, const umodpoly& b) c.resize(n+1, a[0].ring()->zero()); for ( int i=0 ; i<=n; ++i ) { for ( int j=0 ; j<=i; ++j ) { - if ( j > degree(a) || (i-j) > degree(b) ) continue; // TODO optimize! + if ( j > degree(a) || (i-j) > degree(b) ) continue; c[i] = c[i] + a[j] * b[i-j]; } } @@ -278,6 +342,32 @@ static umodpoly operator*(const umodpoly& a, const umodpoly& b) return c; } +static upoly operator*(const upoly& a, const cl_I& x) +{ + if ( zerop(x) ) { + upoly r; + return r; + } + upoly r(a.size()); + for ( size_t i=0; i=ldeg; --deg ) { + up[deg] = the(ex_to(e.coeff(x, deg)).to_cl_N()); + } + for ( ; deg>=0; --deg ) { + up[deg] = 0; + } + canonicalize(up); +} + +static void umodpoly_from_upoly(umodpoly& ump, const upoly& e, const cl_modint_ring& R) +{ + int deg = degree(e); + ump.resize(deg+1); + for ( ; deg>=0; --deg ) { + ump[deg] = R->canonhom(e[deg]); + } + canonicalize(ump); +} + static void umodpoly_from_ex(umodpoly& ump, const ex& e, const ex& x, const cl_modint_ring& R) { // assert: e is in Z[x] @@ -304,12 +419,24 @@ 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) +{ + if ( a.empty() ) return 0; + ex e; + for ( int i=degree(a); i>=0; --i ) { + e += numeric(a[i]) * pow(x, i); + } + return e; +} -static ex umod_to_ex(const umodpoly& a, const ex& x) +static ex umodpoly_to_ex(const umodpoly& a, const ex& x) { if ( a.empty() ) return 0; cl_modint_ring R = a[0].ring(); @@ -327,6 +454,38 @@ static ex umod_to_ex(const umodpoly& a, const ex& x) return e; } +static upoly umodpoly_to_upoly(const umodpoly& a) +{ + upoly e(a.size()); + if ( a.empty() ) return e; + cl_modint_ring R = a[0].ring(); + cl_I mod = R->modulus; + cl_I halfmod = (mod-1) >> 1; + for ( int i=degree(a); i>=0; --i ) { + cl_I n = R->retract(a[i]); + if ( n > halfmod ) { + e[i] = n-mod; + } else { + e[i] = n; + } + } + 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 ( size_t 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. @@ -503,7 +662,7 @@ static bool squarefree(const umodpoly& a) umodpoly b; deriv(a, b); if ( b.empty() ) { - return true; + return false; } umodpoly c; gcd(a, b, c); @@ -516,6 +675,8 @@ static bool squarefree(const umodpoly& a) //////////////////////////////////////////////////////////////////////////////// // modular matrix +typedef vector mvec; + class modular_matrix { friend ostream& operator<<(ostream& o, const modular_matrix& m); @@ -530,90 +691,75 @@ public: cl_MI operator()(size_t row, size_t col) const { return m[row*c + col]; } void mul_col(size_t col, const cl_MI x) { - mvec::iterator i = m.begin() + col; for ( size_t rc=0; rc::iterator i = m.begin() + row*c; for ( size_t cc=0; cc::iterator i1 = m.begin() + row1*c; - vector::iterator i2 = m.begin() + row2*c; for ( size_t cc=0; cc::iterator i1 = m.begin() + row1*c; - vector::iterator i2 = m.begin() + row2*c; for ( size_t cc=0; cc& newrow) { - mvec::iterator i1 = m.begin() + row*c; - mvec::const_iterator i2 = newrow.begin(), end = newrow.end(); - for ( ; i2 != end; ++i1, ++i2 ) { - *i1 = *i2; + for (std::size_t i2 = 0; i2 < newrow.size(); ++i2) { + std::size_t i1 = row*c + i2; + m[i1] = newrow[i2]; } } mvec::const_iterator row_begin(size_t row) const { return m.begin()+row*c; } @@ -645,15 +791,19 @@ modular_matrix operator*(const modular_matrix& m1, const modular_matrix& m2) ostream& operator<<(ostream& o, const modular_matrix& m) { - vector::const_iterator i = m.m.begin(), end = m.m.end(); - size_t wrap = 1; - for ( ; i != end; ++i ) { - o << *i << " "; - if ( !(wrap++ % m.c) ) { - o << endl; + cl_modint_ring R = m(0,0).ring(); + o << "{"; + for ( size_t i=0; iretract(m(i,j)) << ","; + } + o << R->retract(m(i,m.colsize()-1)) << "}"; + if ( i != m.rowsize()-1 ) { + o << ","; } } - o << endl; + o << "}"; return o; } #endif // def DEBUGFACTOR @@ -661,41 +811,39 @@ ostream& operator<<(ostream& o, const modular_matrix& m) // END modular matrix //////////////////////////////////////////////////////////////////////////////// -static void q_matrix(const umodpoly& a, modular_matrix& Q) +/** Calculates the Q matrix for a polynomial. Used by Berlekamp's algorithm. + * + * @param[in] a_ modular polynomial + * @param[out] Q Q matrix + */ +static void q_matrix(const umodpoly& a_, modular_matrix& Q) { + umodpoly a = a_; + normalize_in_field(a); + int n = degree(a); unsigned int q = cl_I_to_uint(a[0].ring()->modulus); -// fast and buggy -// vector r(n, a.R->zero()); -// r[0] = a.R->one(); -// Q.set_row(0, r); -// unsigned int max = (n-1) * q; -// for ( size_t m=1; m<=max; ++m ) { -// cl_MI rn_1 = r.back(); -// for ( size_t i=n-1; i>0; --i ) { -// r[i] = r[i-1] - rn_1 * a[i]; -// } -// r[0] = -rn_1 * a[0]; -// if ( (m % q) == 0 ) { -// Q.set_row(m/q, r); -// } -// } -// slow and (hopefully) correct - cl_MI one = a[0].ring()->one(); - cl_MI zero = a[0].ring()->zero(); - for ( int i=0; izero()); + r[0] = a[0].ring()->one(); + Q.set_row(0, r); + unsigned int max = (n-1) * q; + for ( size_t m=1; m<=max; ++m ) { + cl_MI rn_1 = r.back(); + for ( size_t i=n-1; i>0; --i ) { + r[i] = r[i-1] - (rn_1 * a[i]); + } + r[0] = -rn_1 * a[0]; + if ( (m % q) == 0 ) { + Q.set_row(m/q, r); + } } } +/** Determine the nullspace of a matrix M-1. + * + * @param[in,out] M matrix, will be modified + * @param[out] basis calculated nullspace of M-1 + */ static void nullspace(modular_matrix& M, vector& basis) { const size_t n = M.rowsize(); @@ -740,17 +888,28 @@ static void nullspace(modular_matrix& M, vector& basis) } } +/** Berlekamp's modular factorization. + * + * The implementation follows the algorithm in chapter 8 of [GCL]. + * + * @param[in] a modular polynomial + * @param[out] upv vector containing modular factors. if upv was not empty the + * new elements are added at the end + */ static void berlekamp(const umodpoly& a, upvec& upv) { cl_modint_ring R = a[0].ring(); umodpoly one(1, R->one()); + // find nullspace of Q matrix modular_matrix Q(degree(a), degree(a), R->zero()); q_matrix(a, Q); vector nu; nullspace(Q, nu); + const unsigned int k = nu.size(); if ( k == 1 ) { + // irreducible return; } @@ -762,6 +921,7 @@ static void berlekamp(const umodpoly& a, upvec& upv) list::iterator u = factors.begin(); + // calculate all gcd's while ( true ) { for ( unsigned int s=0; s q ) { - c.resize(q+1, R->zero()); - c[q] = R->one(); - return; +/** Calculates a^(1/prime). + * + * @param[in] a polynomial + * @param[in] prime prime number -> exponent 1/prime + * @param[in] ap resulting polynomial + */ +static void expt_1_over_p(const umodpoly& a, unsigned int prime, umodpoly& ap) +{ + 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]; +/** Modular square free factorization. + * + * @param[in] a polynomial + * @param[out] factors modular factors + * @param[out] mult corresponding multiplicities (exponents) + */ +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; - 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; -} - -static void same_degree_factor(const umodpoly& a, upvec& result) -{ - 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; 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 degrees; + upvec ddfactors; + distinct_degree_factor(a, degrees, ddfactors); - upvec I(m); - umodpoly one(1, R->one()); - 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 { + berlekamp(ddfactors[i], upv); } } } -static void cantor_zassenhaus(const umodpoly& a, upvec& result) -{ -} +// Yes, we can (choose). +#define USE_SAME_DEGREE_FACTOR +/** Modular univariate factorization. + * + * In principle, we have two algorithms at our disposal: Berlekamp's algorithm + * and same degree factorization (SDF). SDF seems to be slightly faster in + * almost all cases so it is activated as default. + * + * @param[in] p modular polynomial + * @param[out] upv vector containing modular factors. if upv was not empty the + * new elements are added at the end + */ static void factor_modular(const umodpoly& p, upvec& upv) { - //same_degree_factor(p, upv); +#ifdef USE_SAME_DEGREE_FACTOR + same_degree_factor(p, upv); +#else berlekamp(p, upv); - return; +#endif } -static void exteuclid(const umodpoly& a, const umodpoly& b, umodpoly& g, umodpoly& s, umodpoly& t) +/** Calculates modular polynomials s and t such that a*s+b*t==1. + * Assertion: a and b are relatively prime and not zero. + * + * @param[in] a polynomial + * @param[in] b polynomial + * @param[out] s polynomial + * @param[out] t polynomial + */ +static void exteuclid(const umodpoly& a, const umodpoly& b, umodpoly& s, umodpoly& t) { if ( degree(a) < degree(b) ) { - exteuclid(b, a, g, t, s); + exteuclid(b, a, t, s); return; } + umodpoly one(1, a[0].ring()->one()); umodpoly c = a; normalize_in_field(c); umodpoly d = b; normalize_in_field(d); - umodpoly c1 = one; - umodpoly c2; + s = one; + t.clear(); umodpoly d1; umodpoly d2 = one; - while ( !d.empty() ) { - umodpoly q; + umodpoly q; + while ( true ) { div(c, d, q); umodpoly r = c - q * d; - umodpoly r1 = c1 - q * d1; - umodpoly r2 = c2 - q * d2; + umodpoly r1 = s - q * d1; + umodpoly r2 = t - q * d2; c = d; - c1 = d1; - c2 = d2; + s = d1; + t = d2; + if ( r.empty() ) break; d = r; d1 = r1; d2 = r2; } - g = c; normalize_in_field(g); - s = c1; - for ( int i=0; i<=degree(s); ++i ) { - s[i] = s[i] * recip(a[degree(a)] * c[degree(c)]); + cl_MI fac = recip(lcoeff(a) * lcoeff(c)); + umodpoly::iterator i = s.begin(), end = s.end(); + for ( ; i!=end; ++i ) { + *i = *i * fac; } canonicalize(s); - s = s * g; - t = c2; - for ( int i=0; i<=degree(t); ++i ) { - t[i] = t[i] * recip(b[degree(b)] * c[degree(c)]); + fac = recip(lcoeff(b) * lcoeff(c)); + i = t.begin(), end = t.end(); + for ( ; i!=end; ++i ) { + *i = *i * fac; } canonicalize(t); - t = t * g; } -static ex replace_lc(const ex& poly, const ex& x, const ex& lc) +/** Replaces the leading coefficient in a polynomial by a given number. + * + * @param[in] poly polynomial to change + * @param[in] lc new leading coefficient + * @return changed polynomial + */ +static upoly replace_lc(const upoly& poly, const cl_I& lc) { - ex r = expand(poly + (lc - poly.lcoeff(x)) * pow(x, poly.degree(x))); + if ( poly.empty() ) return poly; + upoly r = poly; + r.back() = 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) +/** Calculates the bound for the modulus. + * See [Mig]. + */ +static inline cl_I calc_bound(const ex& a, const ex& x, int maxdeg) +{ + cl_I maxcoeff = 0; + cl_R coeff = 0; + for ( int i=a.degree(x); i>=a.ldegree(x); --i ) { + cl_I aa = abs(the(ex_to(a.coeff(x, i)).to_cl_N())); + if ( aa > maxcoeff ) maxcoeff = aa; + coeff = coeff + square(aa); + } + cl_I coeffnorm = ceiling1(the(cln::sqrt(coeff))); + cl_I B = coeffnorm * expt_pos(cl_I(2), cl_I(maxdeg)); + return ( B > maxcoeff ) ? B : maxcoeff; +} + +/** Calculates the bound for the modulus. + * See [Mig]. + */ +static inline cl_I calc_bound(const upoly& a, int maxdeg) +{ + cl_I maxcoeff = 0; + cl_R coeff = 0; + for ( int i=degree(a); i>=0; --i ) { + cl_I aa = abs(a[i]); + if ( aa > maxcoeff ) maxcoeff = aa; + coeff = coeff + square(aa); + } + cl_I coeffnorm = ceiling1(the(cln::sqrt(coeff))); + cl_I B = coeffnorm * expt_pos(cl_I(2), cl_I(maxdeg)); + return ( B > maxcoeff ) ? B : maxcoeff; +} + +/** Hensel lifting as used by factor_univariate(). + * + * The implementation follows the algorithm in chapter 6 of [GCL]. + * + * @param[in] a_ primitive univariate polynomials + * @param[in] p prime number that does not divide lcoeff(a) + * @param[in] u1_ modular factor of a (mod p) + * @param[in] w1_ modular factor of a (mod p), relatively prime to u1_, + * fulfilling u1_*w1_ == a mod p + * @param[out] u lifted factor + * @param[out] w lifted factor, u*w = a + */ +static void hensel_univar(const upoly& a_, unsigned int p, const umodpoly& u1_, const umodpoly& w1_, upoly& u, upoly& w) { - ex a = a_; + upoly a = a_; const cl_modint_ring& R = u1_[0].ring(); // calc bound B - ex maxcoeff; - for ( int i=a.degree(x); i>=a.ldegree(x); --i ) { - maxcoeff += pow(abs(a.coeff(x, i)),2); - } - cl_I normmc = ceiling1(the(cln::sqrt(ex_to(maxcoeff).to_cl_N()))); - cl_I maxdegree = (degree(u1_) > degree(w1_)) ? degree(u1_) : degree(w1_); - cl_I B = normmc * expt_pos(cl_I(2), maxdegree); + int maxdeg = (degree(u1_) > degree(w1_)) ? degree(u1_) : degree(w1_); + cl_I maxmodulus = 2*calc_bound(a, maxdeg); // step 1 - ex alpha = a.lcoeff(x); - ex gamma = gamma_; - if ( gamma == 0 ) { - gamma = alpha; - } - numeric gamma_ui = ex_to(abs(gamma)); - a = a * gamma; + cl_I alpha = lcoeff(a); + a = a * alpha; umodpoly nu1 = u1_; normalize_in_field(nu1); umodpoly nw1 = w1_; normalize_in_field(nw1); - ex phi; - phi = gamma * umod_to_ex(nu1, x); + upoly phi; + phi = umodpoly_to_upoly(nu1) * alpha; umodpoly u1; - umodpoly_from_ex(u1, phi, x, R); - phi = alpha * umod_to_ex(nw1, x); + umodpoly_from_upoly(u1, phi, R); + phi = umodpoly_to_upoly(nw1) * alpha; umodpoly w1; - umodpoly_from_ex(w1, phi, x, R); + umodpoly_from_upoly(w1, phi, R); // step 2 - umodpoly g; umodpoly s; umodpoly t; - exteuclid(u1, w1, g, s, t); - if ( unequal_one(g) ) { - throw logic_error("gcd(u1,w1) != 1"); - } + exteuclid(u1, w1, s, t); // step 3 - ex u = replace_lc(umod_to_ex(u1, x), x, gamma); - ex w = replace_lc(umod_to_ex(w1, x), x, alpha); - ex e = expand(a - u * w); - numeric modulus = p; - const numeric maxmodulus = 2*numeric(B)*gamma_ui; + 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; // step 4 - while ( !e.is_zero() && modulus < maxmodulus ) { - ex c = e / modulus; - phi = expand(umod_to_ex(s, x) * c); + while ( !e.empty() && modulus < maxmodulus ) { + upoly c = e / modulus; + phi = umodpoly_to_upoly(s) * c; umodpoly sigmatilde; - umodpoly_from_ex(sigmatilde, phi, x, R); - phi = expand(umod_to_ex(t, x) * c); + umodpoly_from_upoly(sigmatilde, phi, R); + phi = umodpoly_to_upoly(t) * c; umodpoly tautilde; - umodpoly_from_ex(tautilde, phi, x, R); + umodpoly_from_upoly(tautilde, phi, R); umodpoly r, q; remdiv(sigmatilde, w1, r, q); umodpoly sigma = r; - phi = expand(umod_to_ex(tautilde, x) + umod_to_ex(q, x) * umod_to_ex(u1, x)); + phi = umodpoly_to_upoly(tautilde) + umodpoly_to_upoly(q) * umodpoly_to_upoly(u1); umodpoly tau; - umodpoly_from_ex(tau, phi, x, R); - u = expand(u + umod_to_ex(tau, x) * modulus); - w = expand(w + umod_to_ex(sigma, x) * modulus); - e = expand(a - u * w); + umodpoly_from_upoly(tau, phi, R); + u = u + umodpoly_to_upoly(tau) * modulus; + w = w + umodpoly_to_upoly(sigma) * modulus; + e = a - u * w; modulus = modulus * p; } // step 5 - if ( e.is_zero() ) { - ex delta = u.content(x); - u = u / delta; - w = w / gamma * delta; - return lst(u, w); + if ( e.empty() ) { + cl_I g = u[0]; + for ( size_t i=1; i primes; @@ -1127,93 +1356,183 @@ static unsigned int next_prime(unsigned int p) throw logic_error("next_prime: should not reach this point!"); } -class Partition +/** Manages the splitting a vector of of modular factors into two partitions. + */ +class factor_partition { public: - Partition(size_t n_) : n(n_) + /** Takes the vector of modular factors and initializes the first partition */ + factor_partition(const upvec& factors_) : factors(factors_) { - k.resize(n, 1); - k[0] = 0; - sum = n-1; + n = factors.size(); + k.resize(n, 0); + k[0] = 1; + cache.resize(n-1); + one.resize(1, factors.front()[0].ring()->one()); + len = 1; + last = 0; + 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=1; --i ) { - if ( k[i] ) { - --k[i]; - --sum; - return sum > 0; + if ( last == n-1 ) { + int rem = len - 1; + int p = last - 1; + while ( rem ) { + if ( k[p] ) { + --rem; + --p; + continue; + } + last = p - 1; + while ( k[last] == 0 ) { --last; } + if ( last == 0 && n == 2*len ) return false; + k[last++] = 0; + for ( size_t i=0; i<=len-rem; ++i ) { + k[last] = 1; + ++last; + } + fill(k.begin()+last, k.end(), 0); + --last; + split(); + return true; } - ++k[i]; - ++sum; + last = len; + ++len; + if ( len > n/2 ) return false; + fill(k.begin(), k.begin()+len, 1); + fill(k.begin()+len+1, k.end(), 0); } - return false; + else { + k[last++] = 0; + k[last] = 1; + } + split(); + return true; } + /** Get first partition */ + umodpoly& left() { return lr[0]; } + /** Get second partition */ + umodpoly& right() { return lr[1]; } private: - size_t n, sum; - vector 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= d ) { + lr[group] = lr[group] * cache[pos][d-1]; + } + else { + if ( cache[pos].size() == 0 ) { + cache[pos].push_back(factors[pos] * factors[pos+1]); + } + size_t j = pos + cache[pos].size() + 1; + d -= cache[pos].size(); + while ( d ) { + umodpoly buf = cache[pos].back() * factors[j]; + cache[pos].push_back(buf); + --d; + ++j; + } + lr[group] = lr[group] * cache[pos].back(); + } + } + else { + lr[group] = lr[group] * factors[pos]; + } + } while ( i < n ); + } + void split() + { + lr[0] = one; + lr[1] = one; + if ( n > 6 ) { + split_cached(); } else { - a = a * factors[i]; + for ( size_t i=0; i > cache; + upvec factors; + umodpoly one; + size_t n; + size_t len; + size_t last; + vector k; +}; +/** Contains a pair of univariate polynomial and its modular factors. + * Used by factor_univariate(). + */ struct ModFactors { - ex poly; + upoly poly; upvec factors; }; -static ex factor_univariate(const ex& poly, const ex& x) +/** Univariate polynomial factorization. + * + * Modular factorization is tried for several primes to minimize the number of + * modular factors. Then, Hensel lifting is performed. + * + * @param[in] poly expanded square free univariate polynomial + * @param[in] x symbol + * @param[in,out] prime prime number to start trying modular factorization with, + * output value is the prime number actually used + */ +static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime) { - ex unit, cont, prim; - poly.unitcontprim(x, unit, cont, prim); + ex unit, cont, prim_ex; + poly.unitcontprim(x, unit, cont, prim_ex); + upoly prim; + 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; - numeric lcoeff = ex_to(prim.lcoeff(x)); + + const numeric& cont_n = ex_to(cont); + cl_I i_cont; + if (cont_n.is_integer()) { + i_cont = the(cont_n.to_cl_N()); + } else { + // poly \in Q[x] => poly = q ipoly, ipoly \in Z[x], q \in Q + // factor(poly) \equiv q factor(ipoly) + i_cont = cl_I(1); + } + cl_I lc = lcoeff(prim)*i_cont; upvec factors; while ( trials < 2 ) { + umodpoly modpoly; while ( true ) { - p = next_prime(p); - if ( irem(lcoeff, p) != 0 ) { - R = find_modint_ring(p); - umodpoly modpoly; - umodpoly_from_ex(modpoly, prim, x, R); + prime = next_prime(prime); + if ( !zerop(rem(lc, prime)) ) { + R = find_modint_ring(prime); + 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 ) { @@ -1223,17 +1542,16 @@ static ex factor_univariate(const ex& poly, const ex& x) if ( minfactors == 0 || trialfactors.size() < minfactors ) { factors = trialfactors; - minfactors = factors.size(); - lastp = p; + minfactors = trialfactors.size(); + lastp = prime; trials = 1; } else { ++trials; } } - p = lastp; - R = find_modint_ring(p); - cl_univpoly_modint_ring UPR = find_univpoly_ring(R); + prime = lastp; + R = find_modint_ring(prime); // lift all factor combinations stack tocheck; @@ -1241,24 +1559,24 @@ 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() ) { - if ( part.size_first() == 1 ) { - if ( part.size_second() == 1 ) { - result *= answer.op(0) * answer.op(1); + // call Hensel lifting + hensel_univar(tocheck.top().poly, prime, part.left(), part.right(), f1, f2); + if ( !f1.empty() ) { + // successful, update the stack and the result + if ( part.size_left() == 1 ) { + if ( part.size_right() == 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==). + */ struct EvalPoint { ex x; int evalpoint; }; +#ifdef DEBUGFACTOR +ostream& operator<<(ostream& o, const vector& v) +{ + for ( size_t i=0; i 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) +/** Utility function for multivariate Hensel lifting. + * + * Solves the equation + * s_1*b_1 + ... + s_r*b_r == 1 mod p^k + * with deg(s_i) < deg(a_i) + * and with given b_1 = a_1 * ... * a_{i-1} * a_{i+1} * ... * a_r + * + * The implementation follows the algorithm in chapter 6 of [GCL]. + * + * @param[in] a vector of modular univariate polynomials + * @param[in] x symbol + * @param[in] p prime number + * @param[in] k p^k is modulus + * @return vector of polynomials (s_i) + */ +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)); @@ -1338,10 +1695,10 @@ upvec multiterm_eea_lift(const upvec& a, const ex& x, unsigned int p, unsigned i upvec s; for ( size_t j=1; j mdarg(2); - mdarg[0] = umod_to_ex(q[j-1], x); - mdarg[1] = umod_to_ex(a[j-1], x); + mdarg[0] = umodpoly_to_ex(q[j-1], x); + mdarg[1] = umodpoly_to_ex(a[j-1], x); vector empty; - vector exsigma = multivar_diophant(mdarg, x, umod_to_ex(beta, x), empty, 0, p, k); + vector exsigma = multivar_diophant(mdarg, x, umodpoly_to_ex(beta, x), empty, 0, p, k); umodpoly sigma1; umodpoly_from_ex(sigma1, exsigma[0], x, R); umodpoly sigma2; @@ -1353,10 +1710,12 @@ upvec multiterm_eea_lift(const upvec& a, const ex& x, unsigned int p, unsigned i return s; } -/** - * Assert: a not empty. +/** Changes the modulus of a modular polynomial. Used by eea_lift(). + * + * @param[in] R new modular ring + * @param[in,out] a polynomial to change (in situ) */ -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(); @@ -1367,7 +1726,21 @@ 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_) +/** Utility function for multivariate Hensel lifting. + * + * Solves s*a + t*b == 1 mod p^k given a,b. + * + * The implementation follows the algorithm in chapter 6 of [GCL]. + * + * @param[in] a polynomial + * @param[in] b polynomial + * @param[in] x symbol + * @param[in] p prime number + * @param[in] k p^k is modulus + * @param[out] s_ output polynomial + * @param[out] t_ output polynomial + */ +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; @@ -1375,13 +1748,9 @@ void eea_lift(const umodpoly& a, const umodpoly& b, const ex& x, unsigned int p, umodpoly bmod = b; change_modulus(R, bmod); - umodpoly g; umodpoly smod; umodpoly tmod; - exteuclid(amod, bmod, g, smod, tmod); - if ( unequal_one(g) ) { - throw logic_error("gcd(amod,bmod) != 1"); - } + exteuclid(amod, bmod, smod, tmod); cl_modint_ring Rpk = find_modint_ring(expt_pos(cl_I(p),k)); umodpoly s = smod; @@ -1414,7 +1783,22 @@ 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) +/** Utility function for multivariate Hensel lifting. + * + * Solves the equation + * s_1*b_1 + ... + s_r*b_r == x^m mod p^k + * with given b_1 = a_1 * ... * a_{i-1} * a_{i+1} * ... * a_r + * + * The implementation follows the algorithm in chapter 6 of [GCL]. + * + * @param a vector with univariate polynomials mod p^k + * @param x symbol + * @param m exponent of x^m in the equation to solve + * @param p prime number + * @param k p^k is modulus + * @return vector of polynomials (s_i) + */ +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)); @@ -1423,34 +1807,31 @@ upvec univar_diophant(const upvec& a, const ex& x, unsigned int m, unsigned int if ( r > 2 ) { upvec s = multiterm_eea_lift(a, x, p, k); for ( size_t j=0; j multivar_diophant(const vector& a_, const ex& x, const ex& c, const vector& I, unsigned int d, unsigned int p, unsigned int k) +/** Utility function for multivariate Hensel lifting. + * + * Returns the polynomials s_i that fulfill + * s_1*b_1 + ... + s_r*b_r == c mod + * with given b_1 = a_1 * ... * a_{i-1} * a_{i+1} * ... * a_r + * + * The implementation follows the algorithm in chapter 6 of [GCL]. + * + * @param a_ vector of multivariate factors mod p^k + * @param x symbol (equiv. x_1 in [GCL]) + * @param c polynomial mod p^k + * @param I vector of evaluation points + * @param d maximum total degree of result + * @param p prime number + * @param k p^k is modulus + * @return vector of polynomials (s_i) + */ +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_; @@ -1519,22 +1925,20 @@ vector multivar_diophant(const vector& a_, const ex& x, const ex& c, con ex e = make_modular(buf, R); ex monomial = 1; - for ( size_t m=1; m<=d; ++m ) { - while ( !e.is_zero() && e.has(xnu) ) { - monomial *= (xnu - alphanu); - monomial = expand(monomial); - ex cm = e.diff(ex_to(xnu), m).subs(xnu==alphanu) / factorial(m); - cm = make_modular(cm, R); - if ( !cm.is_zero() ) { - vector delta_s = multivar_diophant(anew, x, cm, Inew, d, p, k); - ex buf = e; - for ( size_t j=0; j(xnu), m).subs(xnu==alphanu) / factorial(m); + cm = make_modular(cm, R); + if ( !cm.is_zero() ) { + vector delta_s = multivar_diophant(anew, x, cm, Inew, d, p, k); + ex buf = e; + for ( size_t j=0; j multivar_diophant(const vector& a_, const ex& x, const ex& c, con modcm = cl_MI(R, poscm); for ( size_t j=0; j 1 ) { z = c.op(i+1); @@ -1584,18 +1988,24 @@ 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 upvec& u, const vector& lcU) +/** Multivariate Hensel lifting. + * The implementation follows the algorithm in chapter 6 of [GCL]. + * Since we don't have a data type for modular multivariate polynomials, the + * respective operations are done in a GiNaC::ex and the function + * make_modular() is then called to make the coefficient modular p^l. + * + * @param a multivariate polynomial primitive in x + * @param x symbol (equiv. x_1 in [GCL]) + * @param I vector of evaluation points (x_2==a_2,x_3==a_3,...) + * @param p prime number (should not divide lcoeff(a mod I)) + * @param l p^l is the modulus of the lifted univariate field + * @param u vector of modular (mod p^l) factors of a mod I + * @param lcU correct leading coefficient of the univariate factors of a mod I + * @return list GiNaC::lst with lifted factors (multivariate factors of a), + * empty if Hensel lifting did not succeed + */ +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)); @@ -1619,7 +2029,7 @@ ex hensel_multivar(const ex& a, const ex& x, const vector& I, unsigne const size_t n = u.size(); vector U(n); for ( size_t i=0; i& I, unsigne } } +/** Takes a factorized expression and puts the factors in a lst. The exponents + * of the factors are discarded, e.g. 7*x^2*(y+1)^4 --> {7,x,y+1}. The first + * element of the list is always the numeric coefficient. + */ static ex put_factors_into_lst(const ex& e) { lst result; - if ( is_a(e) ) { result.append(e); return result; @@ -1702,13 +2115,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) ) { @@ -1720,11 +2131,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); @@ -1733,25 +2142,19 @@ 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) +/** Checks a set of numbers for whether each number has a unique prime factor. + * + * @param[in] f list of numbers to check + * @return true: if number set is bad, false: if set is okay (has unique + * prime factors) + */ +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]; @@ -1768,13 +2171,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 [Wan] 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(); @@ -1783,286 +2203,213 @@ 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())); - if ( cont == 1 ) break; - } - ex pp = expand(normal(p / cont)); + // make polynomial primitive + ex unit, cont, pp; + poly.unitcontprim(x, unit, cont, pp); if ( !is_a(cont) ) { -#ifdef DEBUGFACTOR - return ::factor(cont) * ::factor(pp); -#else - return factor(cont) * factor(pp); -#endif + return factor_sqrfree(cont) * factor_sqrfree(pp); } - /* factor leading coefficient */ - pp = pp.collect(x); - ex vn = pp.lcoeff(x); - pp = pp.expand(); + // factor leading coefficient + ex vn = pp.collect(x).lcoeff(x); ex vnlst; if ( is_a(vn) ) { vnlst = lst(vn); } else { -#ifdef DEBUGFACTOR - ex vnfactors = ::factor(vn); -#else ex vnfactors = factor(vn); -#endif 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 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); - } -#ifdef DEBUGFACTOR - ufac = ::factor(u); -#else - ufac = factor(u); -#endif + // 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(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; - } - 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); - } - else { - D[i-1] = D[i-1] * pow(vnlst.op(2*(j-2)+1), count); - } - } - else { - ftilde[j-1] = ftilde[j-1] / prefac; - break; - } - ++j; - } - } - } - - bool some_factor_unused = false; - for ( size_t i=0; i C(factor_count); - if ( delta == 1 ) { - for ( size_t i=0; i(vn) ) { + // easy case + for ( size_t i=1; i ftilde(vnlst.nops()-1); + for ( size_t i=0; i(ft); + } + // calculate D and C + 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 { - ui = ufaclst.op(2*(i-1)+1); + } + else { + 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; + } + 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); + } + } + C[i] = D[i] * prefac; } - while ( true ) { - ex d = gcd(ui.lcoeff(x), Dtilde); - C[i] = D[i] * ( ui.lcoeff(x) / d ); - ui = ui * ( Dtilde[i] / d ); - delta = delta / ( Dtilde[i] / d ); - if ( delta == 1 ) break; - ui = delta * ui; - C[i] = delta * C[i]; - pp = pp * pow(delta, D.size()-1); + } + // check if something went wrong + bool some_factor_unused = false; + for ( size_t i=0; i epv; s = syms.begin(); @@ -2073,37 +2420,32 @@ static ex factor_multivariate(const ex& poly, const exset& syms) epv.push_back(ep); } - // calc bound B - ex maxcoeff; - for ( int i=u.degree(x); 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 (int)maxdegree ) { - maxdegree = ufaclst[2*i+1].degree(x); + // calc bound p^l + int maxdeg = 0; + for ( int i=1; i<=factor_count; ++i ) { + if ( ufaclst.op(i).degree(x) > maxdeg ) { + maxdeg = ufaclst[i].degree(x); } } - cl_I B = normmc * expt_pos(cl_I(2), maxdegree); + cl_I B = 2*calc_bound(u, x, maxdeg); cl_I l = 1; cl_I pl = prime; while ( pl < B ) { l = l + 1; pl = pl * prime; } - - upvec uvec; + + // set up modular factors (mod p^l) 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 0 ) { + // pull out direct factors int ld = poly.ldegree(x); ex res = factor_univariate(expand(poly/pow(x, ld)), x); return res * pow(x,ld); @@ -2153,17 +2501,16 @@ static ex factor_sqrfree(const ex& poly) return res; } +/** Map used by factor() when factor_options::all is given to access all + * subexpressions and to call factor() on them. + */ struct apply_factor_map : public map_function { unsigned options; apply_factor_map(unsigned options_) : options(options_) { } ex operator()(const ex& e) { if ( e.info(info_flags::polynomial) ) { -#ifdef DEBUGFACTOR - return ::factor(e, options); -#else return factor(e, options); -#endif } if ( is_a(e) ) { ex s1, s2; @@ -2177,11 +2524,7 @@ struct apply_factor_map : public map_function { } s1 = s1.eval(); s2 = s2.eval(); -#ifdef DEBUGFACTOR - return ::factor(s1, options) + s2.map(*this); -#else return factor(s1, options) + s2.map(*this); -#endif } return e.map(*this); } @@ -2189,11 +2532,11 @@ struct apply_factor_map : public map_function { } // anonymous namespace -#ifdef DEBUGFACTOR -ex factor(const ex& poly, unsigned options = 0) -#else +/** Interface function to the outside world. It checks the arguments, tries a + * square free factorization, and then calls factor_sqrfree to do the hard + * work. + */ ex factor(const ex& poly, unsigned options) -#endif { // check arguments if ( !poly.info(info_flags::polynomial) ) { @@ -2218,7 +2561,7 @@ ex factor(const ex& poly, unsigned options) } // make poly square free - ex sfpoly = sqrfree(poly, syms); + ex sfpoly = sqrfree(poly.expand(), syms); // factorize the square free components if ( is_a(sfpoly) ) { @@ -2265,3 +2608,7 @@ ex factor(const ex& poly, unsigned options) } } // namespace GiNaC + +#ifdef DEBUGFACTOR +#include "test.h" +#endif