X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?a=blobdiff_plain;ds=sidebyside;f=ginac%2Ffactor.cpp;h=2b77d1b0110f3e327216a11a20510b901a6b93ee;hb=HEAD;hp=33b8f50125ed853f36ff9a0ceb3480504419e61f;hpb=dfaa3383381a37f871cdd8ae7b19db579dddd83d;p=ginac.git diff --git a/ginac/factor.cpp b/ginac/factor.cpp index 33b8f501..efe438ea 100644 --- a/ginac/factor.cpp +++ b/ginac/factor.cpp @@ -33,7 +33,7 @@ */ /* - * GiNaC Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2024 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 @@ -65,11 +65,12 @@ #include "normal.h" #include "add.h" +#include #include -#include #include #include #include +#include #ifdef DEBUGFACTOR #include #endif @@ -80,6 +81,9 @@ using namespace cln; namespace GiNaC { +// anonymous namespace to hide all utility functions +namespace { + #ifdef DEBUGFACTOR #define DCOUT(str) cout << #str << endl #define DCOUTVAR(var) cout << #var << ": " << var << endl @@ -133,9 +137,6 @@ ostream& operator<<(ostream& o, const vector>& v) #define DCOUT2(str,var) #endif // def DEBUGFACTOR -// anonymous namespace to hide all utility functions -namespace { - //////////////////////////////////////////////////////////////////////////////// // modular univariate polynomial code @@ -143,7 +144,8 @@ typedef std::vector umodpoly; typedef std::vector upoly; typedef vector upvec; -// COPY FROM UPOLY.HPP + +// COPY FROM UPOLY.H // CHANGED size_t -> int !!! template static int degree(const T& p) @@ -156,6 +158,11 @@ template static typename T::value_type lcoeff(const T& p) return p[p.size() - 1]; } +/** Make the polynomial unit normal (having unit normal leading coefficient). + * + * @param[in, out] a polynomial to make unit normal + * @return true if polynomial a was already unit normal, false otherwise + */ static bool normalize_in_field(umodpoly& a) { if (a.size() == 0) @@ -170,61 +177,22 @@ static bool normalize_in_field(umodpoly& a) return false; } +/** Remove leading zero coefficients from polynomial. + * + * @param[in, out] p polynomial from which the zero leading coefficients will be removed + * @param[in] hint assume all coefficients of order ≥ hint are zero + */ template static void canonicalize(T& p, const typename T::size_type hint = std::numeric_limits::max()) { - if (p.empty()) - return; + std::size_t i = min(p.size(), hint); - std::size_t i = p.size() - 1; - // Be fast if the polynomial is already canonicalized - if (!zerop(p[i])) - return; - - if (hint < p.size()) - i = hint; - - bool is_zero = false; - do { - if (!zerop(p[i])) { - ++i; - break; - } - if (i == 0) { - is_zero = true; - break; - } - --i; - } while (true); - - if (is_zero) { - p.clear(); - return; - } + while ( i-- && zerop(p[i]) ) { } - p.erase(p.begin() + i, p.end()); + p.erase(p.begin() + i + 1, p.end()); } -// END COPY FROM UPOLY.HPP - -static void expt_pos(umodpoly& a, unsigned int q) -{ - if ( a.empty() ) return; - cl_MI zero = a[0].ring()->zero(); - 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 */ }; +// END COPY FROM UPOLY.H template struct uvar_poly_p { @@ -487,12 +455,12 @@ static umodpoly umodpoly_to_umodpoly(const umodpoly& a, const cl_modint_ring& R, return e; } -/** Divides all coefficients of the polynomial a by the integer x. +/** Divides all coefficients of the polynomial a by the positive integer x. * All coefficients are supposed to be divisible by x. If they are not, the - * the cast will raise an exception. + * division will raise an exception. * * @param[in,out] a polynomial of which the coefficients will be reduced by x - * @param[in] x integer that divides the coefficients + * @param[in] x positive integer that divides the coefficients */ static void reduce_coeff(umodpoly& a, const cl_I& x) { @@ -502,7 +470,7 @@ static void reduce_coeff(umodpoly& a, const cl_I& x) for (auto & i : a) { // cln cannot perform this division in the modular field cl_I c = R->retract(i); - i = cl_MI(R, the(c / x)); + i = cl_MI(R, exquopos(c, x)); } } @@ -532,7 +500,7 @@ static void rem(const umodpoly& a, const umodpoly& b, umodpoly& r) } while ( k-- ); fill(r.begin()+n, r.end(), a[0].ring()->zero()); - canonicalize(r); + canonicalize(r, n); } /** Calculates quotient of a/b. @@ -643,7 +611,6 @@ static void deriv(const umodpoly& a, umodpoly& d) static bool unequal_one(const umodpoly& a) { - if ( a.empty() ) return true; return ( a.size() != 1 || a[0] != a[0].ring()->one() ); } @@ -669,6 +636,26 @@ static bool squarefree(const umodpoly& a) return equal_one(c); } +/** Computes w^q mod a. + * Uses theorem 2.1 from A.K.Lenstra's PhD thesis; see exercise 8.13 in [GCL]. + * + * @param[in] w polynomial + * @param[in] a modulus polynomial + * @param[in] q common modulus of w and a + * @param[out] r result + */ +static void expt_pos_Q(const umodpoly& w, const umodpoly& a, unsigned int q, umodpoly& r) +{ + if ( w.empty() ) return; + cl_MI zero = w[0].ring()->zero(); + int deg = degree(w); + umodpoly buf(deg*q+1, zero); + for ( size_t i=0; i<=deg; ++i ) { + buf[i*q] = w[i]; + } + rem(buf, a, r); +} + // END modular univariate polynomial code //////////////////////////////////////////////////////////////////////////////// @@ -679,7 +666,9 @@ typedef vector mvec; class modular_matrix { +#ifdef DEBUGFACTOR friend ostream& operator<<(ostream& o, const modular_matrix& m); +#endif public: modular_matrix(size_t r_, size_t c_, const cl_MI& init) : r(r_), c(c_) { @@ -812,6 +801,8 @@ ostream& operator<<(ostream& o, const modular_matrix& m) //////////////////////////////////////////////////////////////////////////////// /** Calculates the Q matrix for a polynomial. Used by Berlekamp's algorithm. + * + * The implementation follows algorithm 8.5 of [GCL]. * * @param[in] a_ modular polynomial * @param[out] Q Q matrix @@ -890,7 +881,7 @@ static void nullspace(modular_matrix& M, vector& basis) /** Berlekamp's modular factorization. * - * The implementation follows the algorithm in chapter 8 of [GCL]. + * The implementation follows algorithm 8.4 of [GCL]. * * @param[in] a modular polynomial * @param[out] upv vector containing modular factors. if upv was not empty the @@ -913,8 +904,7 @@ static void berlekamp(const umodpoly& a, upvec& upv) return; } - list factors; - factors.push_back(a); + list factors = {a}; unsigned int size = 1; unsigned int r = 1; unsigned int q = cl_I_to_uint(R->modulus); @@ -934,8 +924,7 @@ static void berlekamp(const umodpoly& a, upvec& upv) div(*u, g, uo); if ( equal_one(uo) ) { throw logic_error("berlekamp: unexpected divisor."); - } - else { + } else { *u = uo; } factors.push_back(g); @@ -965,9 +954,9 @@ static void berlekamp(const umodpoly& a, upvec& upv) /** Calculates a^(1/prime). * - * @param[in] a polynomial - * @param[in] prime prime number -> exponent 1/prime - * @param[in] ap resulting polynomial + * @param[in] a polynomial + * @param[in] prime prime number -> exponent 1/prime + * @param[out] ap resulting polynomial */ static void expt_1_over_p(const umodpoly& a, unsigned int prime, umodpoly& ap) { @@ -1018,8 +1007,7 @@ static void modsqrfree(const umodpoly& a, upvec& factors, vector& mult) mult[i] *= prime; } } - } - else { + } else { umodpoly ap; expt_1_over_p(a, prime, ap); size_t previ = mult.size(); @@ -1034,7 +1022,7 @@ static void modsqrfree(const umodpoly& a, upvec& factors, vector& mult) /** Distinct degree factorization (DDF). * - * The implementation follows the algorithm in chapter 8 of [GCL]. + * The implementation follows algorithm 8.8 of [GCL]. * * @param[in] a_ modular polynomial * @param[out] degrees vector containing the degrees of the factors of the @@ -1051,23 +1039,17 @@ static void distinct_degree_factor(const umodpoly& a_, vector& degrees, upv int nhalf = degree(a)/2; int i = 1; - umodpoly w(2); - w[0] = R->zero(); - w[1] = R->one(); + umodpoly w = {R->zero(), R->one()}; umodpoly x = w; while ( i <= nhalf ) { - expt_pos(w, q); umodpoly buf; - rem(w, a, buf); + expt_pos_Q(w, a, q, buf); w = buf; - umodpoly wx = w - x; - gcd(a, wx, buf); + gcd(a, w - x, 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; @@ -1104,8 +1086,7 @@ static void same_degree_factor(const umodpoly& a, upvec& upv) for ( size_t i=0; 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); + radicand = radicand + 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; + return ceiling1(the(cln::sqrt(radicand))); } -/** Calculates the bound for the modulus. - * See [Mig]. +/** Calculates bound for the product of absolute values (modulus) of the roots. + * Uses Landau's inequality, see [Mig]. */ -static inline cl_I calc_bound(const upoly& a, int maxdeg) +static inline cl_I calc_bound(const upoly& a) { - cl_I maxcoeff = 0; - cl_R coeff = 0; + cl_R radicand = 0; for ( int i=degree(a); i>=0; --i ) { cl_I aa = abs(a[i]); - if ( aa > maxcoeff ) maxcoeff = aa; - coeff = coeff + square(aa); + radicand = radicand + 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; + return ceiling1(the(cln::sqrt(radicand))); } /** Hensel lifting as used by factor_univariate(). * - * The implementation follows the algorithm in chapter 6 of [GCL]. + * The implementation follows algorithm 6.1 of [GCL]. * * @param[in] a_ primitive univariate polynomials * @param[in] p prime number that does not divide lcoeff(a) @@ -1248,7 +1221,7 @@ static void hensel_univar(const upoly& a_, unsigned int p, const umodpoly& u1_, // calc bound B int maxdeg = (degree(u1_) > degree(w1_)) ? degree(u1_) : degree(w1_); - cl_I maxmodulus = 2*calc_bound(a, maxdeg); + cl_I maxmodulus = ash(calc_bound(a), maxdeg+1); // = 2 * calc_bound(a) * 2^maxdeg // step 1 cl_I alpha = lcoeff(a); @@ -1311,46 +1284,41 @@ static void hensel_univar(const upoly& a_, unsigned int p, const umodpoly& u1_, if ( alpha != 1 ) { w = w / alpha; } - } - else { + } else { u.clear(); } } -/** Returns a new prime number. +/** Returns a new small prime number. * - * @param[in] p prime number - * @return next prime number after p + * @param[in] n an integer + * @return smallest prime greater than n */ -static unsigned int next_prime(unsigned int p) -{ - static vector primes; - if ( primes.size() == 0 ) { - primes.push_back(3); primes.push_back(5); primes.push_back(7); - } - if ( p >= primes.back() ) { - unsigned int candidate = primes.back() + 2; - while ( true ) { - size_t n = primes.size()/2; - for ( size_t i=0; i primes = {2, 3, 5, 7}; + unsigned int candidate = primes.back(); + while (primes.back() <= n) { + candidate += 2; + bool is_prime = true; + for (size_t i=1; primes[i]*primes[i]<=candidate; ++i) { + if (candidate % primes[i] == 0) { + is_prime = false; + break; } - primes.push_back(candidate); - if ( candidate > p ) break; } - return candidate; + if (is_prime) + primes.push_back(candidate); } for (auto & it : primes) { - if ( it > p ) { + if ( it > n ) { return it; } } throw logic_error("next_prime: should not reach this point!"); } -/** Manages the splitting a vector of of modular factors into two partitions. +/** Manages the splitting of a vector of modular factors into two partitions. */ class factor_partition { @@ -1402,8 +1370,7 @@ public: if ( len > n/2 ) return false; fill(k.begin(), k.begin()+len, 1); fill(k.begin()+len+1, k.end(), 0); - } - else { + } else { k[last++] = 0; k[last] = 1; } @@ -1426,8 +1393,7 @@ private: if ( d ) { if ( cache[pos].size() >= d ) { lr[group] = lr[group] * cache[pos][d-1]; - } - else { + } else { if ( cache[pos].size() == 0 ) { cache[pos].push_back(factors[pos] * factors[pos+1]); } @@ -1441,8 +1407,7 @@ private: } lr[group] = lr[group] * cache[pos].back(); } - } - else { + } else { lr[group] = lr[group] * factors[pos]; } } while ( i < n ); @@ -1453,8 +1418,7 @@ private: lr[1] = one; if ( n > 6 ) { split_cached(); - } - else { + } else { for ( size_t i=0; i multivar_diophant(const vector& a_, const ex& x, const ex& * 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]. + * The implementation follows algorithm 6.3 of [GCL]. * * @param[in] a vector of modular univariate polynomials * @param[in] x symbol @@ -1723,7 +1686,7 @@ static void change_modulus(const cl_modint_ring& R, umodpoly& a) * * Solves s*a + t*b == 1 mod p^k given a,b. * - * The implementation follows the algorithm in chapter 6 of [GCL]. + * The implementation follows algorithm 6.3 of [GCL]. * * @param[in] a polynomial * @param[in] b polynomial @@ -1782,7 +1745,7 @@ static void eea_lift(const umodpoly& a, const umodpoly& b, const ex& x, unsigned * 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]. + * The implementation follows algorithm 6.3 of [GCL]. * * @param a vector with univariate polynomials mod p^k * @param x symbol @@ -1805,8 +1768,7 @@ static upvec univar_diophant(const upvec& a, const ex& x, unsigned int m, unsign rem(bmod, a[j], buf); result.push_back(buf); } - } - else { + } else { umodpoly s, t; eea_lift(a[1], a[0], x, p, k, s, t); umodpoly bmod = umodpoly_to_umodpoly(s, R, m); @@ -1828,7 +1790,7 @@ static upvec univar_diophant(const upvec& a, const ex& x, unsigned int m, unsign struct make_modular_map : public map_function { cl_modint_ring R; make_modular_map(const cl_modint_ring& R_) : R(R_) { } - ex operator()(const ex& e) + ex operator()(const ex& e) override { if ( is_a(e) || is_a(e) ) { return e.map(*this); @@ -1840,8 +1802,7 @@ struct make_modular_map : public map_function { numeric n(R->retract(emod)); if ( n > halfmod ) { return n-mod; - } - else { + } else { return n; } } @@ -1868,7 +1829,7 @@ static ex make_modular(const ex& e, const cl_modint_ring& R) * 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]. + * The implementation follows algorithm 6.2 of [GCL]. * * @param a_ vector of multivariate factors mod p^k * @param x symbol (equiv. x_1 in [GCL]) @@ -1884,7 +1845,8 @@ static vector multivar_diophant(const vector& a_, const ex& x, const ex& { vector a = a_; - const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k)); + const cl_I modulus = expt_pos(cl_I(p),k); + const cl_modint_ring R = find_modint_ring(modulus); const size_t r = a.size(); const size_t nu = I.size() + 1; @@ -1934,8 +1896,7 @@ static vector multivar_diophant(const vector& a_, const ex& x, const ex& e = make_modular(buf, R); } } - } - else { + } else { upvec amod; for ( size_t i=0; i multivar_diophant(const vector& a_, const ex& x, const ex& if ( is_a(c) ) { nterms = c.nops(); z = c.op(0); - } - else { + } else { nterms = 1; z = c; } @@ -1959,16 +1919,13 @@ static vector multivar_diophant(const vector& a_, const ex& x, const ex& cl_I cm = the(ex_to(z.lcoeff(x)).to_cl_N()); upvec delta_s = univar_diophant(amod, x, m, p, k); cl_MI modcm; - cl_I poscm = cm; - while ( poscm < 0 ) { - poscm = poscm + expt_pos(cl_I(p),k); - } + cl_I poscm = plusp(cm) ? cm : mod(cm, modulus); modcm = cl_MI(R, poscm); for ( size_t j=0; j 1 ) { + if ( nterms > 1 && i+1 != nterms ) { z = c.op(i+1); } } @@ -1982,7 +1939,7 @@ static vector multivar_diophant(const vector& a_, const ex& x, const ex& } /** Multivariate Hensel lifting. - * The implementation follows the algorithm in chapter 6 of [GCL]. + * The implementation follows algorithm 6.4 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. @@ -2082,74 +2039,69 @@ static ex hensel_multivar(const ex& a, const ex& x, const vector& I, acand *= U[i]; } if ( expand(a-acand).is_zero() ) { - lst res; - for ( size_t i=0; i {7,x,y+1}. The first - * element of the list is always the numeric coefficient. + * element of the result is always the numeric coefficient. */ -static ex put_factors_into_lst(const ex& e) +static exvector put_factors_into_vec(const ex& e) { - lst result; + exvector result; if ( is_a(e) ) { - result.append(e); + result.push_back(e); return result; } if ( is_a(e) ) { - result.append(1); - result.append(e.op(0)); + result.push_back(1); + result.push_back(e.op(0)); return result; } if ( is_a(e) || is_a(e) ) { ex icont(e.integer_content()); - result.append(icont); - result.append(e/icont); + result.push_back(icont); + result.push_back(e/icont); return result; } if ( is_a(e) ) { ex nfac = 1; + result.push_back(nfac); for ( size_t i=0; i(op) ) { nfac = op; } if ( is_a(op) ) { - result.append(op.op(0)); + result.push_back(op.op(0)); } if ( is_a(op) || is_a(op) ) { - result.append(op); + result.push_back(op); } } - result.prepend(nfac); + result[0] = nfac; return result; } - throw runtime_error("put_factors_into_lst: bad term."); + throw runtime_error("put_factors_into_vec: bad term."); } /** Checks a set of numbers for whether each number has a unique prime factor. * - * @param[in] f list of numbers to check + * @param[in] f 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) +static bool checkdivisors(const exvector& f) { - const int k = f.nops(); + const int k = f.size(); numeric q, r; vector d(k); - d[0] = ex_to(abs(f.op(0))); + d[0] = ex_to(abs(f[0])); for ( int i=1; i(abs(f.op(i))); + q = ex_to(abs(f[i])); for ( int j=i-1; j>=0; --j ) { r = d[j]; do { @@ -2174,25 +2126,24 @@ static bool checkdivisors(const lst& f) * * @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] x first symbol that appears in u + * @param[in] syms_wox remaining symbols that appear in u + * @param[in] f vector 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, +static void generate_set(const ex& u, const ex& vn, const ex& x, const exset& syms_wox, const exvector& f, numeric& modulus, ex& u0, vector& a) { - const ex& x = *syms.begin(); while ( true ) { ++modulus; // generate a set of integers ... u0 = u; ex vna = vn; ex vnatry; - exset::const_iterator s = syms.begin(); - ++s; + auto s = syms_wox.begin(); for ( size_t i=0; i(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; + exvector fnum = f; + fnum[0] = fnum[0] * u0.content(x); + for ( size_t i=1; i(fnum[i]) ) { + s = syms_wox.begin(); for ( size_t j=0; j(cont) ) { - return factor_sqrfree(cont) * factor_sqrfree(pp); - } - - // factor leading coefficient - ex vn = pp.collect(x).lcoeff(x); - ex vnlst; - if ( is_a(vn) ) { - vnlst = lst(vn); - } - else { - ex vnfactors = factor(vn); - vnlst = put_factors_into_lst(vnfactors); - } - - const unsigned int maxtrials = 3; - numeric modulus = (vnlst.nops() > 3) ? vnlst.nops() : 3; - vector a(syms.size()-1, 0); - - // try now to factorize until we are successful - while ( true ) { +struct factorization_ctx { + const ex poly, x; // polynomial, first symbol x... + const exset syms_wox; // ...remaining symbols w/o x + ex unit, cont, pp; // unit * cont * pp == poly + ex vn; exvector vnlst; // leading coeff, factors of leading coeff + numeric modulus; // incremented each time we try + /** returns factors or empty if it did not succeed */ + ex try_next_evaluation_homomorphism() + { + constexpr unsigned maxtrials = 3; + vector a(syms_wox.size(), 0); unsigned int trialcount = 0; unsigned int prime; int factor_count = 0; int min_factor_count = -1; ex u, delta; - ex ufac, ufaclst; + ex ufac; + exvector ufaclst; // try several evaluation points to reduce the number of factors while ( trialcount < maxtrials ) { // generate a set of valid evaluation points - generate_set(pp, vn, syms, ex_to(vnlst), modulus, u, a); + generate_set(pp, vn, x, syms_wox, vnlst, modulus, u, a); ufac = factor_univariate(u, x, prime); - ufaclst = put_factors_into_lst(ufac); - factor_count = ufaclst.nops()-1; - delta = ufaclst.op(0); + ufaclst = put_factors_into_vec(ufac); + factor_count = ufaclst.size()-1; + delta = ufaclst[0]; if ( factor_count <= 1 ) { // irreducible - return poly; + return lst{pp}; } if ( min_factor_count < 0 ) { // first time here @@ -2319,20 +2239,18 @@ static ex factor_multivariate(const ex& poly, const exset& syms) vector C(factor_count); if ( is_a(vn) ) { // easy case - for ( size_t i=1; i ftilde(vnlst.nops()-1); + vector ftilde(vnlst.size()-1); for ( size_t i=0; i D(factor_count, 1); if ( delta == 1 ) { for ( int i=0; i(ufaclst.op(i+1).lcoeff(x)); + numeric prefac = ex_to(ufaclst[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]); + numeric q; + while ( irem(prefac, ftilde[j], q) == 0 ) { + prefac = q; ++count; } if ( count ) { used_flag[j] = true; - D[i] = D[i] * pow(vnlst.op(j+1), count); + D[i] = D[i] * pow(vnlst[j+1], count); } } C[i] = D[i] * prefac; } - } - else { + } else { for ( int i=0; i(ufaclst.op(i+1).lcoeff(x)); + numeric prefac = ex_to(ufaclst[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]); + numeric q; + while ( irem(prefac, ftilde[j], q) == 0 ) { + prefac = q; ++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); + ufaclst[i+1] = ufaclst[i+1] * (ftilde[j]/g); ++count; } if ( count ) { used_flag[j] = true; - D[i] = D[i] * pow(vnlst.op(j+1), count); + D[i] = D[i] * pow(vnlst[j+1], count); } } C[i] = D[i] * prefac; @@ -2392,60 +2311,114 @@ static ex factor_multivariate(const ex& poly, const exset& syms) } } if ( some_factor_unused ) { - continue; + return lst{}; // next try } } - + // multiply the remaining content of the univariate polynomial into the // first factor if ( delta != 1 ) { C[0] = C[0] * delta; - ufaclst.let_op(1) = ufaclst.op(1) * delta; + ufaclst[1] = ufaclst[1] * delta; } // set up evaluation points - EvalPoint ep; vector epv; - s = syms.begin(); - ++s; + auto s = syms_wox.begin(); for ( size_t i=0; i maxdeg ) { + if ( ufaclst[i].degree(x) > maxdeg ) { maxdeg = ufaclst[i].degree(x); } } - cl_I B = 2*calc_bound(u, x, maxdeg); + cl_I B = ash(calc_bound(u, x), maxdeg+1); // = 2 * calc_bound(u,x) * 2^maxdeg cl_I l = 1; cl_I pl = prime; while ( pl < B ) { 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 ctx_in_x; + for (auto x : syms) { + exset syms_wox; // remaining syms w/o x + copy_if(syms.begin(), syms.end(), + inserter(syms_wox, syms_wox.end()), [x](const ex& y){ return y != x; }); + + factorization_ctx ctx{poly, x, syms_wox}; + + // make polynomial primitive + poly.unitcontprim(x, ctx.unit, ctx.cont, ctx.pp); + if ( !is_a(ctx.cont) ) { + // content is a polynomial in one or more of remaining syms, let's start over + return ctx.unit * factor_sqrfree(ctx.cont) * factor_sqrfree(ctx.pp); + } + + // find factors of leading coefficient + ctx.vn = ctx.pp.collect(x).lcoeff(x); + ctx.vnlst = put_factors_into_vec(factor(ctx.vn)); + + ctx.modulus = (ctx.vnlst.size() > 3) ? ctx.vnlst.size() : numeric(3); + + ctx_in_x.push_back(ctx); + } + + // try an evaluation homomorphism for each context in round-robin + auto ctx = ctx_in_x.begin(); + while ( true ) { + + ex res = ctx->try_next_evaluation_homomorphism(); + + if ( res != lst{} ) { + // found the factors + ex result = ctx->cont * ctx->unit; for ( size_t i=0; ix, unit, cont, pp); + result *= unit * cont * pp; } return result; } + + // switch context for next symbol + if (++ctx == ctx_in_x.end()) { + ctx = ctx_in_x.begin(); + } } } @@ -2453,7 +2426,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms) */ struct find_symbols_map : public map_function { exset syms; - ex operator()(const ex& e) + ex operator()(const ex& e) override { if ( is_a(e) ) { syms.insert(e); @@ -2478,13 +2451,12 @@ static ex factor_sqrfree(const ex& poly) if ( findsymbols.syms.size() == 1 ) { // univariate case const ex& x = *(findsymbols.syms.begin()); - if ( poly.ldegree(x) > 0 ) { + int ld = poly.ldegree(x); + if ( ld > 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); - } - else { + } else { ex res = factor_univariate(poly, x); return res; } @@ -2501,7 +2473,7 @@ static ex factor_sqrfree(const ex& poly) struct apply_factor_map : public map_function { unsigned options; apply_factor_map(unsigned options_) : options(options_) { } - ex operator()(const ex& e) + ex operator()(const ex& e) override { if ( e.info(info_flags::polynomial) ) { return factor(e, options); @@ -2511,26 +2483,51 @@ struct apply_factor_map : public map_function { for ( size_t i=0; i void +factor_iter(const ex &e, F yield) +{ + if (is_a(e)) { + for (const auto &f : e) { + if (is_a(f)) { + yield(f.op(0), f.op(1)); + } else { + yield(f, ex(1)); + } + } + } else { + if (is_a(e)) { + yield(e.op(0), e.op(1)); + } else { + yield(e, ex(1)); + } + } +} -/** 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. +/** This function factorizes a polynomial. It checks the arguments, + * tries a square free factorization, and then calls factor_sqrfree + * to do the hard work. + * + * This function expands its argument, so for polynomials with + * explicit factors it's better to call it on each one separately + * (or use factor() which does just that). */ -ex factor(const ex& poly, unsigned options) +static ex factor1(const ex& poly, unsigned options) { // check arguments if ( !poly.info(info_flags::polynomial) ) { @@ -2557,51 +2554,35 @@ ex factor(const ex& poly, unsigned options) ex sfpoly = sqrfree(poly.expand(), syms); // factorize the square free components - if ( is_a(sfpoly) ) { - // case: (polynomial)^exponent - const ex& base = sfpoly.op(0); - if ( !is_a(base) ) { - // simple case: (monomial)^exponent - return sfpoly; - } - ex f = factor_sqrfree(base); - return pow(f, sfpoly.op(1)); - } - if ( is_a(sfpoly) ) { - // case: multiple factors - ex res = 1; - for ( size_t i=0; i(t) ) { - const ex& base = t.op(0); - if ( !is_a(base) ) { - res *= t; - } - else { - ex f = factor_sqrfree(base); - res *= pow(f, t.op(1)); - } - } - else if ( is_a(t) ) { - ex f = factor_sqrfree(t); - res *= f; - } - else { - res *= t; + ex res = 1; + factor_iter(sfpoly, + [&](const ex &f, const ex &k) { + if ( is_a(f) ) { + res *= pow(factor_sqrfree(f), k); + } else { + // simple case: (monomial)^exponent + res *= pow(f, k); } - } - return res; - } - if ( is_a(sfpoly) ) { - return poly; - } - // case: (polynomial) - ex f = factor_sqrfree(sfpoly); - return f; + }); + return res; } -} // namespace GiNaC +} // anonymous namespace -#ifdef DEBUGFACTOR -#include "test.h" -#endif +/** Interface function to the outside world. It uses factor1() + * on each of the explicitly present factors of poly. + */ +ex factor(const ex& poly, unsigned options) +{ + ex result = 1; + factor_iter(poly, + [&](const ex &f1, const ex &k1) { + factor_iter(factor1(f1, options), + [&](const ex &f2, const ex &k2) { + result *= pow(f2, k1*k2); + }); + }); + return result; +} + +} // namespace GiNaC