X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Ffactor.cpp;h=b71291356722f81a921b22796a75c854f2a753bb;hp=fb73897218db06ed31b40f32bf52aabca7e33acf;hb=dcce72e85f936117b39774a839ec1d76add66aec;hpb=2b0ad5c381dc081cc4066b0c5f939f5169ad9ab3 diff --git a/ginac/factor.cpp b/ginac/factor.cpp index fb738972..b7129135 100644 --- a/ginac/factor.cpp +++ b/ginac/factor.cpp @@ -1,12 +1,35 @@ /** @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. */ /* @@ -61,17 +84,6 @@ namespace GiNaC { #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 - -// anonymous namespace to hide all utility functions -namespace { - -typedef vector mvec; -#ifdef DEBUGFACTOR ostream& operator<<(ostream& o, const vector& v) { vector::const_iterator i = v.begin(), end = v.end(); @@ -98,6 +110,13 @@ ostream& operator<<(ostream& o, const vector& v) } 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(); @@ -107,7 +126,14 @@ ostream& operator<<(ostream& o, const vector< vector >& v) } return o; } -#endif +#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 @@ -621,6 +647,8 @@ static bool squarefree(const umodpoly& a) //////////////////////////////////////////////////////////////////////////////// // modular matrix +typedef vector mvec; + class modular_matrix { friend ostream& operator<<(ostream& o, const modular_matrix& m); @@ -770,6 +798,11 @@ ostream& operator<<(ostream& o, const modular_matrix& m) // END modular matrix //////////////////////////////////////////////////////////////////////////////// +/** 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_; @@ -793,6 +826,11 @@ static void q_matrix(const umodpoly& a_, modular_matrix& Q) } } +/** 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(); @@ -837,11 +875,20 @@ 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; @@ -849,6 +896,7 @@ static void berlekamp(const umodpoly& a, upvec& upv) const unsigned int k = nu.size(); if ( k == 1 ) { + // irreducible return; } @@ -860,6 +908,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 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; @@ -909,6 +968,12 @@ static void expt_1_over_p(const umodpoly& a, unsigned int prime, umodpoly& ap) } } +/** 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); @@ -954,6 +1019,18 @@ static void modsqrfree(const umodpoly& a, upvec& factors, vector& mult) } } +#endif // deactivation of square free factorization + +/** Distinct degree factorization (DDF). + * + * The implementation follows the algorithm in chapter 8 of [GCL]. + * + * @param[in] a_ modular polynomial + * @param[out] degrees vector containing the degrees of the factors of the + * corresponding polynomials in ddfactors. + * @param[out] ddfactors vector containing polynomials which factors have the + * degree given in degrees. + */ static void distinct_degree_factor(const umodpoly& a_, vector& degrees, upvec& ddfactors) { umodpoly a = a_; @@ -995,6 +1072,16 @@ static void distinct_degree_factor(const umodpoly& a_, vector& degrees, upv } } +/** Modular same degree factorization. + * Same degree factorization is a kind of misnomer. It performs distinct degree + * factorization, but instead of using the Cantor-Zassenhaus algorithm it + * (sub-optimally) uses Berlekamp's algorithm for the factors of the same + * degree. + * + * @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 same_degree_factor(const umodpoly& a, upvec& upv) { cl_modint_ring R = a[0].ring(); @@ -1013,8 +1100,19 @@ static void same_degree_factor(const umodpoly& a, upvec& upv) } } +// 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) { #ifdef USE_SAME_DEGREE_FACTOR @@ -1024,7 +1122,7 @@ static void factor_modular(const umodpoly& p, upvec& upv) #endif } -/** Calculates polynomials s and t such that a*s+b*t==1. +/** 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 @@ -1074,6 +1172,12 @@ static void exteuclid(const umodpoly& a, const umodpoly& b, umodpoly& s, umodpol canonicalize(t); } +/** 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) { if ( poly.empty() ) return poly; @@ -1082,6 +1186,9 @@ static upoly replace_lc(const upoly& poly, const cl_I& lc) return r; } +/** 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; @@ -1096,6 +1203,9 @@ static inline cl_I calc_bound(const ex& a, const ex& x, int 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; @@ -1110,6 +1220,18 @@ static inline cl_I calc_bound(const upoly& a, int 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) { upoly a = a_; @@ -1186,6 +1308,11 @@ static void hensel_univar(const upoly& a_, unsigned int p, const umodpoly& u1_, } } +/** Returns a new prime number. + * + * @param[in] p prime number + * @return next prime number after p + */ static unsigned int next_prime(unsigned int p) { static vector primes; @@ -1216,9 +1343,12 @@ static unsigned int next_prime(unsigned int p) throw logic_error("next_prime: should not reach this point!"); } +/** Manages the splitting a vector of of modular factors into two partitions. + */ class factor_partition { public: + /** Takes the vector of modular factors and initializes the first partition */ factor_partition(const upvec& factors_) : factors(factors_) { n = factors.size(); @@ -1234,9 +1364,8 @@ public: size_t size() const { return n; } size_t size_left() const { return n-len; } size_t size_right() const { return len; } -#ifdef DEBUGFACTOR - void get() const { DCOUTVAR(k); } -#endif + /** Initializes the next partition. + Returns true, if there is one, false otherwise. */ bool next() { if ( last == n-1 ) { @@ -1274,7 +1403,9 @@ public: split(); return true; } + /** Get first partition */ umodpoly& left() { return lr[0]; } + /** Get second partition */ umodpoly& right() { return lr[1]; } private: void split_cached() @@ -1333,12 +1464,25 @@ private: vector k; }; +/** Contains a pair of univariate polynomial and its modular factors. + * Used by factor_univariate(). + */ struct ModFactors { upoly poly; upvec factors; }; +/** 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_ex; @@ -1398,8 +1542,10 @@ static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime) const size_t n = tocheck.top().factors.size(); factor_partition part(tocheck.top().factors); while ( true ) { + // 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); @@ -1453,7 +1599,10 @@ static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime) } } else { + // not successful if ( !part.next() ) { + // if no more combinations left, return polynomial as + // irreducible result *= upoly_to_ex(tocheck.top().poly, x); tocheck.pop(); break; @@ -1465,21 +1614,51 @@ static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime) return unit * cont * result; } +/** Second interface to factor_univariate() to be used if the information about + * the prime is not needed. + */ static inline ex factor_univariate(const ex& poly, const ex& x) { unsigned int prime; return factor_univariate(poly, x, prime); } +/** Represents an evaluation point (==). + */ 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); +/** 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(); @@ -1508,8 +1687,10 @@ static upvec multiterm_eea_lift(const upvec& a, const ex& x, unsigned int p, uns 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) */ static void change_modulus(const cl_modint_ring& R, umodpoly& a) { @@ -1522,6 +1703,20 @@ static void change_modulus(const cl_modint_ring& R, umodpoly& a) canonicalize(a); } +/** 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); @@ -1565,6 +1760,21 @@ static void eea_lift(const umodpoly& a, const umodpoly& b, const ex& x, unsigned s_ = s; t_ = t; } +/** 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)); @@ -1595,6 +1805,10 @@ static upvec univar_diophant(const upvec& a, const ex& x, unsigned int m, unsign return result; } +/** Map used by function make_modular(). + * Finds every coefficient in a polynomial and replaces it by is value in the + * given modular ring R (symmetric representation). + */ struct make_modular_map : public map_function { cl_modint_ring R; make_modular_map(const cl_modint_ring& R_) : R(R_) { } @@ -1619,12 +1833,36 @@ struct make_modular_map : public map_function { } }; +/** Helps mimicking modular multivariate polynomial arithmetic. + * + * @param e expression of which to make the coefficients equal to their value + * in the modular ring R (symmetric representation) + * @param R modular ring + * @return resulting expression + */ static ex make_modular(const ex& e, const cl_modint_ring& R) { make_modular_map map(R); return map(e.expand()); } +/** 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) { @@ -1727,17 +1965,24 @@ static vector multivar_diophant(const vector& a_, const ex& x, const ex& 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)); @@ -1833,10 +2078,13 @@ static ex hensel_multivar(const ex& a, const ex& x, const vector& I, } } +/** 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; @@ -1871,20 +2119,11 @@ 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(x))); if ( !is_a(g) ) { continue; } if ( !is_a(vn) ) { - /* ... and for which the evaluated factors have each an unique prime factor */ + // ... 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=p.ldegree(x); --i ) { @@ -1997,7 +2248,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms) return factor_sqrfree(cont) * factor_sqrfree(pp); } - /* factor leading coefficient */ + // factor leading coefficient ex vn = pp.collect(x).lcoeff(x); ex vnlst; if ( is_a(vn) ) { @@ -2012,7 +2263,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms) numeric modulus = (vnlst.nops() > 3) ? vnlst.nops() : 3; vector a(syms.size()-1, 0); - /* try now to factorize until we are successful */ + // try now to factorize until we are successful while ( true ) { unsigned int trialcount = 0; @@ -2022,10 +2273,10 @@ static ex factor_multivariate(const ex& poly, const exset& syms) ex u, delta; ex ufac, ufaclst; - /* try several evaluation points to reduce the number of modular factors */ + // try several evaluation points to reduce the number of factors while ( trialcount < maxtrials ) { - /* generate a set of valid evaluation points */ + // generate a set of valid evaluation points generate_set(pp, vn, syms, ex_to(vnlst), modulus, u, a); ufac = factor_univariate(u, x, prime); @@ -2034,19 +2285,19 @@ static ex factor_multivariate(const ex& poly, const exset& syms) delta = ufaclst.op(0); if ( factor_count <= 1 ) { - /* irreducible */ + // irreducible return poly; } if ( min_factor_count < 0 ) { - /* first time here */ + // first time here min_factor_count = factor_count; } else if ( min_factor_count == factor_count ) { - /* one less to try */ + // one less to try ++trialcount; } else if ( min_factor_count > factor_count ) { - /* new minimum, reset trial counter */ + // new minimum, reset trial counter min_factor_count = factor_count; trialcount = 0; } @@ -2055,11 +2306,16 @@ static ex factor_multivariate(const ex& poly, const exset& syms) // determine true leading coefficients for the Hensel lifting vector C(factor_count); if ( is_a(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 ) { @@ -2115,7 +2371,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms) C[i] = D[i] * prefac; } } - + // check if something went wrong bool some_factor_unused = false; for ( size_t i=0; i epv; s = syms.begin(); @@ -2157,12 +2416,15 @@ static ex factor_multivariate(const ex& poly, const exset& syms) l = l + 1; pl = pl * prime; } + + // set up modular factors (mod p^l) cl_modint_ring R = find_modint_ring(expt_pos(cl_I(prime),l)); 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); @@ -2215,6 +2483,9 @@ 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_) { } @@ -2243,6 +2514,10 @@ struct apply_factor_map : public map_function { } // anonymous namespace +/** 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) { // check arguments