Both algorithm and implementation are a bit inefficient: * algorithm does not make use of sparseness of inputs (and most of multivariate polynomials are quite sparse) * GiNaC's expressions are essentially immutable. Thus some simple operations (i.e. multiplying a polynomial by an element of the base ring) are prohibitively expensive. * All numbers (i.e. GiNaC::numeric objects) are heap allocated. Still it's much faster (~5x on bivariate polynomials, ~100x on 3-variate ones) than (subresultant) PRS algorithm, so gcd() uses modular algorithm by default. --- ginac/Makefile.am | 18 ++++ ginac/normal.cpp | 10 ++- ginac/normal.h | 17 +++- ginac/polynomial/chinrem_gcd.cpp | 14 +++ ginac/polynomial/chinrem_gcd.h | 19 ++++ ginac/polynomial/collect_vargs.cpp | 164 ++++++++++++++++++++++++++++++ ginac/polynomial/collect_vargs.h | 34 ++++++ ginac/polynomial/debug.hpp | 1 + ginac/polynomial/divide_in_z_p.cpp | 88 ++++++++++++++++ ginac/polynomial/divide_in_z_p.h | 27 +++++ ginac/polynomial/euclid_gcd_wrap.h | 56 ++++++++++ ginac/polynomial/eval_point_finder.h | 51 +++++++++ ginac/polynomial/mgcd.cpp | 96 +++++++++++++++++ ginac/polynomial/newton_interpolate.h | 36 +++++++ ginac/polynomial/optimal_vars_finder.cpp | 134 ++++++++++++++++++++++++ ginac/polynomial/optimal_vars_finder.h | 20 ++++ ginac/polynomial/pgcd.cpp | 136 +++++++++++++++++++++++++ ginac/polynomial/pgcd.h | 27 +++++ ginac/polynomial/poly_cra.h | 38 +++++++ ginac/polynomial/primes_factory.h | 60 +++++++++++ ginac/polynomial/primpart_content.cpp | 78 ++++++++++++++ ginac/polynomial/smod_helpers.h | 72 +++++++++++++ 22 files changed, 1192 insertions(+), 4 deletions(-) create mode 100644 ginac/polynomial/chinrem_gcd.cpp create mode 100644 ginac/polynomial/chinrem_gcd.h create mode 100644 ginac/polynomial/collect_vargs.cpp create mode 100644 ginac/polynomial/collect_vargs.h create mode 100644 ginac/polynomial/divide_in_z_p.cpp create mode 100644 ginac/polynomial/divide_in_z_p.h create mode 100644 ginac/polynomial/euclid_gcd_wrap.h create mode 100644 ginac/polynomial/eval_point_finder.h create mode 100644 ginac/polynomial/mgcd.cpp create mode 100644 ginac/polynomial/newton_interpolate.h create mode 100644 ginac/polynomial/optimal_vars_finder.cpp create mode 100644 ginac/polynomial/optimal_vars_finder.h create mode 100644 ginac/polynomial/pgcd.cpp create mode 100644 ginac/polynomial/pgcd.h create mode 100644 ginac/polynomial/poly_cra.h create mode 100644 ginac/polynomial/primes_factory.h create mode 100644 ginac/polynomial/primpart_content.cpp create mode 100644 ginac/polynomial/smod_helpers.h diff --git a/ginac/Makefile.am b/ginac/Makefile.am index 7ab16d6..76e32b2 100644 --- a/ginac/Makefile.am +++ b/ginac/Makefile.am @@ -35,6 +35,24 @@ polynomial/interpolate_padic_uvar.h \ polynomial/sr_gcd_uvar.h \ polynomial/heur_gcd_uvar.h \ polynomial/gcd_uvar.cpp \ +polynomial/chinrem_gcd.cpp \ +polynomial/chinrem_gcd.h \ +polynomial/collect_vargs.cpp \ +polynomial/collect_vargs.h \ +polynomial/divide_in_z_p.cpp \ +polynomial/divide_in_z_p.h \ +polynomial/euclid_gcd_wrap.h \ +polynomial/eval_point_finder.h \ +polynomial/mgcd.cpp \ +polynomial/newton_interpolate.h \ +polynomial/optimal_vars_finder.cpp \ +polynomial/optimal_vars_finder.h \ +polynomial/pgcd.cpp \ +polynomial/pgcd.h \ +polynomial/poly_cra.h \ +polynomial/primes_factory.h \ +polynomial/primpart_content.cpp \ +polynomial/smod_helpers.h \ polynomial/debug.hpp libginac_la_LDFLAGS = -version-info $(LT_VERSION_INFO) -release $(LT_RELEASE) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 8f5ba73..54fff47 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -44,6 +44,7 @@ #include "pseries.h" #include "symbol.h" #include "utils.h" +#include "polynomial/chinrem_gcd.h" namespace GiNaC { @@ -1623,8 +1624,15 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio } #endif } + if (options & gcd_options::use_sr_gcd) { + g = sr_gcd(aex, bex, var); + } else { + exvector vars; + for (std::size_t n = sym_stats.size(); n-- != 0; ) + vars.push_back(sym_stats[n].sym); + g = chinrem_gcd(aex, bex, vars); + } - g = sr_gcd(aex, bex, var); if (g.is_equal(_ex1)) { // Keep cofactors factored if possible if (ca) diff --git a/ginac/normal.h b/ginac/normal.h index 5a735e3..a352be8 100644 --- a/ginac/normal.h +++ b/ginac/normal.h @@ -37,8 +37,12 @@ struct gcd_options { enum { /** - * Usually GiNaC tries heuristic GCD algorithm before PRS. - * Some people don't like this, so here's a flag to disable it. + * Usually GiNaC tries heuristic GCD first, because typically + * it's much faster than anything else. Even if heuristic + * algorithm fails, the overhead is negligible w.r.t. cost + * of computing the GCD by some other method. However, some + * people dislike it, so here's a flag which tells GiNaC + * to NOT use the heuristic algorithm. */ no_heur_gcd = 2, /** @@ -48,7 +52,14 @@ struct gcd_options * factored polynomials. DON'T SET THIS unless you *really* * know what are you doing! */ - no_part_factored = 4 + no_part_factored = 4, + /** + * By default GiNaC uses modular GCD algorithm. Typically + * it's much faster than PRS (pseudo remainder sequence) + * algorithm. This flag forces GiNaC to use PRS algorithm + */ + use_sr_gcd = 8 + }; }; diff --git a/ginac/polynomial/chinrem_gcd.cpp b/ginac/polynomial/chinrem_gcd.cpp new file mode 100644 index 0000000..4048ee0 --- /dev/null +++ b/ginac/polynomial/chinrem_gcd.cpp @@ -0,0 +1,14 @@ +#include "chinrem_gcd.h" +#include "optimal_vars_finder.h" + +namespace GiNaC +{ + +ex chinrem_gcd(const ex& A, const ex& B) +{ + const exvector vars = gcd_optimal_variables_order(A, B); + ex g = chinrem_gcd(A, B, vars); + return g; +} +} // namespace GiNaC + diff --git a/ginac/polynomial/chinrem_gcd.h b/ginac/polynomial/chinrem_gcd.h new file mode 100644 index 0000000..6c58a9f --- /dev/null +++ b/ginac/polynomial/chinrem_gcd.h @@ -0,0 +1,19 @@ +#ifndef GINAC_CHINREM_GCD_H +#define GINAC_CHINREM_GCD_H +#include "ex.h" + +namespace GiNaC +{ + +extern ex chinrem_gcd(const ex& A_, const ex& B_, const exvector& vars); +extern ex chinrem_gcd(const ex& A, const ex& B); + +struct chinrem_gcd_failed +{ + virtual ~chinrem_gcd_failed() { } +}; + +} + +#endif /* GINAC_CHINREM_GCD_H */ + diff --git a/ginac/polynomial/collect_vargs.cpp b/ginac/polynomial/collect_vargs.cpp new file mode 100644 index 0000000..5649e4a --- /dev/null +++ b/ginac/polynomial/collect_vargs.cpp @@ -0,0 +1,164 @@ +#include <iterator> +#include <map> +#include <algorithm> +#include <stdexcept> +#include <string> +#include <ginac/ginac.h> +#include "collect_vargs.h" +#include <cln/integer.h> +#include "smod_helpers.h" +#include "debug.hpp" + +namespace GiNaC +{ + +typedef std::map<exp_vector_t, ex> ex_collect_priv_t; + +static void +collect_vargs(ex_collect_priv_t& ec, ex e, const exvector& vars); +static void +collect_term(ex_collect_priv_t& ec, const ex& e, const exvector& vars); +static void wipe_out_zeros(ex_collect_priv_t& ec); + +template<typename T, typename CoeffCMP> +struct compare_terms +{ + const CoeffCMP& coeff_cmp; + explicit compare_terms(const CoeffCMP& coeff_cmp_) : coeff_cmp(coeff_cmp_) + { } + inline bool operator()(const T& t1, const T& t2) const + { + bool exponent_is_less = + std::lexicographical_compare(t1.first.rbegin(), + t1.first.rend(), + t2.first.rbegin(), + t2.first.rend()); + if (exponent_is_less) + return true; + + if ((t1.first == t2.first) && + coeff_cmp(t2.second, t2.second)) + return true; + return false; + } +}; + +template<typename T, typename CoeffCMP> +static struct compare_terms<T, CoeffCMP> +make_compare_terms(const T& dummy, const CoeffCMP& coeff_cmp) +{ + return compare_terms<T, CoeffCMP>(coeff_cmp); +} + +void collect_vargs(ex_collect_t& ec, const ex& e, const exvector& vars) +{ + ex_collect_priv_t ecp; + collect_vargs(ecp, e, vars); + ec.reserve(ecp.size()); + std::copy(ecp.begin(), ecp.end(), std::back_inserter(ec)); + std::sort(ec.begin(), ec.end(), + make_compare_terms(*ec.begin(), ex_is_less())); +} + +static void +collect_vargs(ex_collect_priv_t& ec, ex e, const exvector& vars) +{ + e = e.expand(); + if (e.is_zero()) { + ec.clear(); + return; + } + + if (!is_a<add>(e)) { + collect_term(ec, e, vars); + return; + } + + for (const_iterator i = e.begin(); i != e.end(); ++i) + collect_term(ec, *i, vars); + + wipe_out_zeros(ec); +} + +static void +collect_term(ex_collect_priv_t& ec, const ex& e, const exvector& vars) +{ + if (e.is_zero()) + return; + static const ex ex1(1); + exp_vector_t key(vars.size()); + ex pre_coeff = e; + for (std::size_t i = 0; i < vars.size(); ++i) { + const int var_i_pow = pre_coeff.degree(vars[i]); + key[i] = var_i_pow; + pre_coeff = pre_coeff.coeff(vars[i], var_i_pow); + } + ex_collect_priv_t::iterator i = ec.find(key); + if (i != ec.end()) + i->second += pre_coeff; + else + ec.insert(ex_collect_priv_t::value_type(key, pre_coeff)); +} + +static void wipe_out_zeros(ex_collect_priv_t& m) +{ + ex_collect_priv_t::iterator i = m.begin(); + while (i != m.end()) { + // be careful to not invalide iterator, use post-increment + // for that, see e.g. + // http://coding.derkeiler.com/Archive/C_CPP/comp.lang.cpp/2004-02/0502.html + if (i->second.is_zero()) + m.erase(i++); + else + ++i; + } +} + +GiNaC::ex +ex_collect_to_ex(const ex_collect_t& ec, const GiNaC::exvector& vars) +{ + exvector ev; + ev.reserve(ec.size()); + for (std::size_t i = 0; i < ec.size(); ++i) { + exvector tv; + tv.reserve(vars.size() + 1); + for (std::size_t j = 0; j < vars.size(); ++j) { + if (ec[i].first[j] != 0) + tv.push_back(power(vars[j], ec[i].first[j])); + } + tv.push_back(ec[i].second); + ex tmp = (new mul(tv))->setflag(status_flags::dynallocated); + ev.push_back(tmp); + } + ex ret = (new add(ev))->setflag(status_flags::dynallocated); + return ret; +} + +ex lcoeff_wrt(ex e, const exvector& x) +{ + static const ex ex0(0); + e = e.expand(); + if (e.is_zero()) + return ex0; + + ex_collect_t ec; + collect_vargs(ec, e, x); + return ec.rbegin()->second; +} + +cln::cl_I integer_lcoeff(const ex& e, const exvector& vars) +{ + ex_collect_t ec; + collect_vargs(ec, e, vars); + if (ec.size() == 0) + return cln::cl_I(0); + ex lc = ec.rbegin()->second; + bug_on(!is_a<numeric>(lc), "leading coefficient is not an integer"); + bug_on(!lc.info(info_flags::integer), + "leading coefficient is not an integer"); + + return to_cl_I(lc); +} + +} // namespace GiNaC + diff --git a/ginac/polynomial/collect_vargs.h b/ginac/polynomial/collect_vargs.h new file mode 100644 index 0000000..4dd16ba --- /dev/null +++ b/ginac/polynomial/collect_vargs.h @@ -0,0 +1,34 @@ +#ifndef GINAC_COLLECT_VARGS_HPP +#define GINAC_COLLECT_VARGS_HPP +#include "ex.h" +#include <cln/integer.h> +#include <vector> +#include <utility> // std::pair + +namespace GiNaC +{ + +typedef std::vector<int> exp_vector_t; +typedef std::vector<std::pair<exp_vector_t, ex> > ex_collect_t; + +extern void +collect_vargs(ex_collect_t& ec, const ex& e, const exvector& x); +extern ex +ex_collect_to_ex(const ex_collect_t& ec, const exvector& x); + +/** + * Leading coefficient of a multivariate polynomial e, considering it + * as a multivariate polynomial in x_0, \ldots x_{n-1} with coefficients + * being univariate polynomials in R[x_n] (where R is some ring) + */ +extern ex lcoeff_wrt(ex e, const exvector& x); + +/** + * Leading coefficient c \in R (where R = Z or Z_p) of a multivariate + * polynomial e \in R[x_0, \ldots, x_n] + */ +extern cln::cl_I integer_lcoeff(const ex& e, const exvector& vars); + +} // namespace GiNaC + +#endif /* GINAC_COLLECT_VARGS_HPP */ diff --git a/ginac/polynomial/debug.hpp b/ginac/polynomial/debug.hpp index 901ef87..48659fd 100644 --- a/ginac/polynomial/debug.hpp +++ b/ginac/polynomial/debug.hpp @@ -4,6 +4,7 @@ #include <string> #include <stdexcept> #include <sstream> +#include "compiler.h" #define DEBUG_PREFIX __func__ << ':' << __LINE__ << ": " #define EXCEPTION_PREFIX std::string(__func__) + std::string(": ") + diff --git a/ginac/polynomial/divide_in_z_p.cpp b/ginac/polynomial/divide_in_z_p.cpp new file mode 100644 index 0000000..dff30a9 --- /dev/null +++ b/ginac/polynomial/divide_in_z_p.cpp @@ -0,0 +1,88 @@ +#include <ginac/ginac.h> +#include "smod_helpers.h" + +namespace GiNaC +{ +/** + * Exact polynomial division of a, b \in Z_p[x_0, \ldots, x_n] + * It doesn't check whether the inputs are proper polynomials, so be careful + * of what you pass in. + * + * @param a first multivariate polynomial (dividend) + * @param b second multivariate polynomial (divisor) + * @param q quotient (returned) + * @param var variables X iterator to first element of vector of symbols + * + * @return "true" when exact division succeeds (the quotient is returned in + * q), "false" otherwise. + * @note @a p = 0 means the base ring is Z + */ +bool divide_in_z_p(const ex &a, const ex &b, ex &q, const exvector& vars, const long p) +{ + static const ex _ex1(1); + if (b.is_zero()) + throw(std::overflow_error("divide_in_z: division by zero")); + if (b.is_equal(_ex1)) { + q = a; + return true; + } + if (is_exactly_a<numeric>(a)) { + if (is_exactly_a<numeric>(b)) { + // p == 0 means division in Z + if (p == 0) { + const numeric tmp = ex_to<numeric>(a/b); + if (tmp.is_integer()) { + q = tmp; + return true; + } else + return false; + } else { + q = (a*recip(ex_to<numeric>(b), p)).smod(numeric(p)); + return true; + } + } else + return false; + } + if (a.is_equal(b)) { + q = _ex1; + return true; + } + + // Main symbol + const ex &x = vars.back(); + + // Compare degrees + int adeg = a.degree(x), bdeg = b.degree(x); + if (bdeg > adeg) + return false; + + // Polynomial long division (recursive) + ex r = a.expand(); + if (r.is_zero()) + return true; + int rdeg = adeg; + ex eb = b.expand(); + ex blcoeff = eb.coeff(x, bdeg); + exvector v; + v.reserve(std::max(rdeg - bdeg + 1, 0)); + exvector rest_vars(vars); + rest_vars.pop_back(); + while (rdeg >= bdeg) { + ex term, rcoeff = r.coeff(x, rdeg); + if (!divide_in_z_p(rcoeff, blcoeff, term, rest_vars, p)) + break; + term = (term*power(x, rdeg - bdeg)).expand(); + v.push_back(term); + r = (r - term*eb).expand(); + if (p != 0) + r = r.smod(numeric(p)); + if (r.is_zero()) { + q = (new add(v))->setflag(status_flags::dynallocated); + return true; + } + rdeg = r.degree(x); + } + return false; +} + +} diff --git a/ginac/polynomial/divide_in_z_p.h b/ginac/polynomial/divide_in_z_p.h new file mode 100644 index 0000000..f790fcd --- /dev/null +++ b/ginac/polynomial/divide_in_z_p.h @@ -0,0 +1,27 @@ +#ifndef GINAC_CHINREM_GCD_DIVIDE_IN_Z_P_H +#define GINAC_CHINREM_GCD_DIVIDE_IN_Z_P_H + +#include <ginac/ginac.h> +namespace GiNaC +{ + +/** + * Exact polynomial division of a, b \in Z_p[x_0, \ldots, x_n] + * It doesn't check whether the inputs are proper polynomials, so be careful + * of what you pass in. + * + * @param a first multivariate polynomial (dividend) + * @param b second multivariate polynomial (divisor) + * @param q quotient (returned) + * @param var variables X iterator to first element of vector of symbols + * + * @return "true" when exact division succeeds (the quotient is returned in + * q), "false" otherwise. + */ +extern bool +divide_in_z_p(const ex &a, const ex &b, ex &q, const exvector& vars, const long p); + +} // namespace GiNaC + +#endif /* GINAC_CHINREM_GCD_DIVIDE_IN_Z_P_H */ + diff --git a/ginac/polynomial/euclid_gcd_wrap.h b/ginac/polynomial/euclid_gcd_wrap.h new file mode 100644 index 0000000..d77f71c --- /dev/null +++ b/ginac/polynomial/euclid_gcd_wrap.h @@ -0,0 +1,56 @@ +#ifndef GINAC_PGCD_EUCLID_GCD_H +#define GINAC_PGCD_EUCLID_GCD_H +#include "upoly.hpp" +#include "gcd_euclid.tcc" +#include "smod_helpers.h" +#include <ginac/ginac.h> +#include "debug.hpp" + +namespace GiNaC +{ + +static void ex2upoly(umodpoly& u, ex e, const ex& var, const long p) +{ + e = e.expand(); + cln::cl_modint_ring R = cln::find_modint_ring(cln::cl_I(p)); + u.resize(e.degree(var) + 1); + for (int i = 0; i <= e.degree(var); ++i) { + ex ce = e.coeff(var, i); + bug_on(!is_a<numeric>(ce), "i = " << i << ", " << + "coefficient is not a number: " << ce); + const cln::cl_I c = to_cl_I(ce); + u[i] = R->canonhom(c); + } +} + +static ex umodpoly2ex(const umodpoly& a, const ex& var, const long p) +{ + cln::cl_modint_ring R = cln::find_modint_ring(cln::cl_I(p)); + const numeric pnum(p); + exvector ev(a.size()); + for (std::size_t i = a.size(); i-- != 0; ) { + const cln::cl_I c = smod(R->retract(a[i]), p); + const ex term = numeric(c)*power(var, i); + ev.push_back(term); + } + ex ret = (new add(ev))->setflag(status_flags::dynallocated); + return ret; +} + +static ex euclid_gcd(ex A, ex B, const ex& var, const long p) +{ + A = A.expand(); + B = B.expand(); + + umodpoly a, b; + ex2upoly(a, A, var, p); + ex2upoly(b, B, var, p); + umodpoly g; + gcd_euclid(g, a, b); + ex ge = umodpoly2ex(g, var, p); + return ge; +} + +} // namespace GiNaC + +#endif diff --git a/ginac/polynomial/eval_point_finder.h b/ginac/polynomial/eval_point_finder.h new file mode 100644 index 0000000..8ceadd1 --- /dev/null +++ b/ginac/polynomial/eval_point_finder.h @@ -0,0 +1,51 @@ +#ifndef GINAC_PGCD_EVAL_POINT_FINDER_H +#define GINAC_PGCD_EVAL_POINT_FINDER_H +#include <ginac/ginac.h> +#include <set> + +namespace GiNaC +{ + +/// Find a `good' evaluation point b \in Z_p for a pair of multivariate +/// polynomials A, B \in Z_p[x_n][x_0, \ldots, x_n]. Here `good' means +/// that b is not a root of GCD of contents of A and B. N.B. content +/// is univariate polynomial \in Z_p[x_n] +struct eval_point_finder +{ + typedef long value_type; + const value_type p; + std::set<value_type> points; + const random_modint modint_generator; + bool operator()(value_type& b, const ex& g, const ex& x); + eval_point_finder(const value_type& p_) : p(p_), modint_generator(p) + { } +}; + +bool eval_point_finder::operator()(value_type& b, const ex& lc, const ex& x) +{ + random_modint modint_generator(p); + // Search for a new element of field + while (points.size() < p - 1) { + value_type b_ = modint_generator(); + // check if this evaluation point was already used + if (points.find(b_) != points.end()) + continue; + + // mark found value as used, even if it's a root of lc + // (so we don't need to do the check once more) + points.insert(b_); + // Now make sure it's NOT the root of GCD's leading coeffient + if (lc.subs(x == numeric(b_)).smod(numeric(p)).is_zero()) + continue; + // Nice, it's our next evaluation point + b = b_; + return true; + } + // All possible evaluation points were used. + return false; +} + +} + +#endif /* GINAC_PGCD_EVAL_POINT_FINDER_H */ + diff --git a/ginac/polynomial/mgcd.cpp b/ginac/polynomial/mgcd.cpp new file mode 100644 index 0000000..b39475c --- /dev/null +++ b/ginac/polynomial/mgcd.cpp @@ -0,0 +1,96 @@ +#include "chinrem_gcd.h" +#include <cln/integer.h> +#include "pgcd.h" +#include "collect_vargs.h" +#include "primes_factory.h" +#include "divide_in_z_p.h" +#include "poly_cra.h" + +namespace GiNaC +{ + +static cln::cl_I extract_integer_content(ex& Apr, const ex& A) +{ + static const cln::cl_I n1(1); + const numeric icont_ = A.integer_content(); + const cln::cl_I icont = cln::the<cln::cl_I>(icont_.to_cl_N()); + if (icont != 1) { + Apr = (A/icont_).expand(); + return icont; + } else { + Apr = A; + return n1; + } +} + +ex chinrem_gcd(const ex& A_, const ex& B_, const exvector& vars) +{ + ex A, B; + const cln::cl_I a_icont = extract_integer_content(A, A_); + const cln::cl_I b_icont = extract_integer_content(B, B_); + const cln::cl_I c = cln::gcd(a_icont, b_icont); + + const cln::cl_I a_lc = integer_lcoeff(A, vars); + const cln::cl_I b_lc = integer_lcoeff(B, vars); + const cln::cl_I g_lc = cln::gcd(a_lc, b_lc); + + const ex& x(vars.back()); + int n = std::min(A.degree(x), B.degree(x)); + const cln::cl_I A_max_coeff = to_cl_I(A.max_coefficient()); + const cln::cl_I B_max_coeff = to_cl_I(B.max_coefficient()); + const cln::cl_I lcoeff_limit = (cln::cl_I(1) << n)*cln::abs(g_lc)* + std::min(A_max_coeff, B_max_coeff); + + cln::cl_I q = 0; + ex H = 0; + + long p; + primes_factory pfactory; + while (true) { + bool has_primes = pfactory(p, g_lc); + if (!has_primes) + throw chinrem_gcd_failed(); + + const numeric pnum(p); + ex Ap = A.smod(pnum); + ex Bp = B.smod(pnum); + ex Cp = pgcd(Ap, Bp, vars, p); + + const cln::cl_I g_lcp = smod(g_lc, p); + const cln::cl_I Cp_lc = integer_lcoeff(Cp, vars); + const cln::cl_I nlc = smod(recip(Cp_lc, p)*g_lcp, p); + Cp = (Cp*numeric(nlc)).expand().smod(pnum); + int cp_deg = Cp.degree(x); + if (cp_deg == 0) + return numeric(g_lc); + if (zerop(q)) { + H = Cp; + n = cp_deg; + q = p; + } else { + if (cp_deg == n) { + ex H_next = chinese_remainder(H, q, Cp, p); + q = q*cln::cl_I(p); + H = H_next; + } else if (cp_deg < n) { + // all previous homomorphisms are unlucky + q = p; + H = Cp; + n = cp_deg; + } else { + // dp_deg > d_deg: current prime is bad + } + } + if (q < lcoeff_limit) + continue; // don't bother to do division checks + ex C, dummy1, dummy2; + extract_integer_content(C, H); + if (divide_in_z_p(A, C, dummy1, vars, 0) && + divide_in_z_p(B, C, dummy2, vars, 0)) + return (numeric(c)*C).expand(); + // else: try more primes + } +} + +} // namespace GiNaC + diff --git a/ginac/polynomial/newton_interpolate.h b/ginac/polynomial/newton_interpolate.h new file mode 100644 index 0000000..e80687d --- /dev/null +++ b/ginac/polynomial/newton_interpolate.h @@ -0,0 +1,36 @@ +#ifndef GINAC_PGCD_NEWTON_INTERPOLATE_H +#define GINAC_PGCD_NEWTON_INTERPOLATE_H +#include "ex.h" +#include "numeric.h" +#include "smod_helpers.h" +namespace GiNaC +{ + +/** + * Newton interpolation -- incremental form. + * + * Given a polynomials e1 \in Z_p[Y], prev \in Z_p[x][Y], and evaluation + * points pt1, prevpts \in Z_p compute a polynomial P \in Z_[x][Y] such + * that P(x = pt1) = e1, P(x = pt \in prevpts) = prev(x = pt). + * + * @var{prevpts} enumerates previous evaluation points, should have a form + * (x - pt_2) (x - pt_3) \ldots (x - pt_n). + * @var{prev} encodes the result of previous interpolations. + */ +ex newton_interp(const ex& e1, const long pt1, + const ex& prev, const ex& prevpts, + const ex& x, const long p) +{ + const numeric pnum(p); + const numeric nc = ex_to<numeric>(prevpts.subs(x == numeric(pt1)).smod(pnum)); + const numeric nc_1 = recip(nc, p); + // result = prev + prevpts (e1 - prev|_{x = pt_1})/prevpts|_{x = pt_1) + ex tmp = prev.subs(x == numeric(pt1)).smod(pnum); + tmp = (((e1 - tmp).expand().smod(pnum))*nc_1).smod(pnum); + tmp = (prevpts*tmp).expand().smod(pnum); + tmp = (prev + tmp).expand().smod(pnum); + return tmp; +} + +} +#endif /* GINAC_PGCD_NEWTON_INTERPOLATE_H */ diff --git a/ginac/polynomial/optimal_vars_finder.cpp b/ginac/polynomial/optimal_vars_finder.cpp new file mode 100644 index 0000000..4903699 --- /dev/null +++ b/ginac/polynomial/optimal_vars_finder.cpp @@ -0,0 +1,134 @@ +#include <algorithm> +#include <cstddef> +#include "optimal_vars_finder.h" +#include "add.h" +#include "mul.h" +#include "power.h" +#include "symbol.h" +#include "numeric.h" + +namespace GiNaC +{ +namespace +{ +// XXX: copy-pasted from normal.cpp. +/* + * Statistical information about symbols in polynomials + */ + +/** This structure holds information about the highest and lowest degrees + * in which a symbol appears in two multivariate polynomials "a" and "b". + * A vector of these structures with information about all symbols in + * two polynomials can be created with the function get_symbol_stats(). + * + * @see get_symbol_stats */ +struct sym_desc +{ + /** Reference to symbol */ + ex sym; + + /** Highest degree of symbol in polynomial "a" */ + int deg_a; + + /** Highest degree of symbol in polynomial "b" */ + int deg_b; + + /** Lowest degree of symbol in polynomial "a" */ + int ldeg_a; + + /** Lowest degree of symbol in polynomial "b" */ + int ldeg_b; + + /** Maximum of deg_a and deg_b (Used for sorting) */ + int max_deg; + + /** Maximum number of terms of leading coefficient of symbol in both polynomials */ + std::size_t max_lcnops; + + /** Commparison operator for sorting */ + bool operator<(const sym_desc &x) const + { + if (max_deg == x.max_deg) + return max_lcnops < x.max_lcnops; + else + return max_deg < x.max_deg; + } +}; + +// Vector of sym_desc structures +typedef std::vector<sym_desc> sym_desc_vec; + +// Add symbol the sym_desc_vec (used internally by get_symbol_stats()) +static void add_symbol(const ex &s, sym_desc_vec &v) +{ + sym_desc_vec::const_iterator it = v.begin(), itend = v.end(); + while (it != itend) { + if (it->sym.is_equal(s)) // If it's already in there, don't add it a second time + return; + ++it; + } + sym_desc d; + d.sym = s; + v.push_back(d); +} + +// Collect all symbols of an expression (used internally by get_symbol_stats()) +static void collect_symbols(const ex &e, sym_desc_vec &v) +{ + if (is_a<symbol>(e)) { + add_symbol(e, v); + } else if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) { + for (size_t i=0; i<e.nops(); i++) + collect_symbols(e.op(i), v); + } else if (is_exactly_a<power>(e)) { + collect_symbols(e.op(0), v); + } +} + +/** + * @brief Find the order of variables which is optimal for GCD computation. + * + * Collects statistical information about the highest and lowest degrees + * of all variables that appear in two polynomials. Sorts the variables + * by minimum degree (lowest to highest). The information gathered by + * this function is used by GCD routines to find out the main variable + * for GCD computation. + * + * @param a first multivariate polynomial + * @param b second multivariate polynomial + * @param v vector of sym_desc structs (filled in) */ +static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v) +{ + collect_symbols(a, v); + collect_symbols(b, v); + sym_desc_vec::iterator it = v.begin(), itend = v.end(); + while (it != itend) { + int deg_a = a.degree(it->sym); + int deg_b = b.degree(it->sym); + it->deg_a = deg_a; + it->deg_b = deg_b; + it->max_deg = std::max(deg_a, deg_b); + it->max_lcnops = std::max(a.lcoeff(it->sym).nops(), b.lcoeff(it->sym).nops()); + it->ldeg_a = a.ldegree(it->sym); + it->ldeg_b = b.ldegree(it->sym); + ++it; + } + std::sort(v.begin(), v.end()); +} +// XXX: end copy-paste from normal.cpp + +} // anonymous namespace of helper functions + +exvector gcd_optimal_variables_order(const ex& a, const ex& b) +{ + sym_desc_vec v; + get_symbol_stats(a, b, v); + exvector vars; + vars.reserve(v.size()); + for (std::size_t i = v.size(); i-- != 0; ) + vars.push_back(v[i].sym); + return vars; +} + +} // namespace GiNaC + diff --git a/ginac/polynomial/optimal_vars_finder.h b/ginac/polynomial/optimal_vars_finder.h new file mode 100644 index 0000000..9715a47 --- /dev/null +++ b/ginac/polynomial/optimal_vars_finder.h @@ -0,0 +1,20 @@ +#ifndef GINAC_CHINREM_GCD_OPTIMAL_SYMBOL_FINDER_H +#define GINAC_CHINREM_GCD_OPTIMAL_SYMBOL_FINDER_H +#include <vector> +#include "ex.h" + +namespace GiNaC +{ +/** + * @brief Find the order of variables which is optimal for GCD computation. + * + * Collects statistical information about the highest and lowest degrees + * of all variables that appear in two polynomials. Sorts the variables + * by minimum degree (lowest to highest). The information gathered by + * this function is used by GCD routines to find out the main variable + * for GCD computation. + */ +extern exvector gcd_optimal_variables_order(const ex& A, const ex& B); +} +#endif /* GINAC_CHINREM_GCD_OPTIMAL_SYMBOL_FINDER_H */ + diff --git a/ginac/polynomial/pgcd.cpp b/ginac/polynomial/pgcd.cpp new file mode 100644 index 0000000..ddc3b8c --- /dev/null +++ b/ginac/polynomial/pgcd.cpp @@ -0,0 +1,136 @@ +#include "pgcd.h" +#include "collect_vargs.h" +#include "smod_helpers.h" +#include "euclid_gcd_wrap.h" +#include "eval_point_finder.h" +#include "newton_interpolate.h" +#include "divide_in_z_p.h" + +namespace GiNaC +{ + +extern void +primpart_content(ex& pp, ex& c, ex e, const exvector& vars, const long p); + +// Computes the GCD of two polynomials over a prime field. +// Based on Algorithm 7.2 from "Algorithms for Computer Algebra" +// A and B are considered as Z_p[x_n][x_0, \ldots, x_{n-1}], that is, +// as a polynomials in variables x_0, \ldots x_{n-1} having coefficients +// from the ring Z_p[x_n] +ex pgcd(const ex& A, const ex& B, const exvector& vars, const long p) +{ + static const ex ex1(1); + if (A.is_zero()) + return B; + + if (B.is_zero()) + return A; + + if (is_a<numeric>(A)) + return ex1; + + if (is_a<numeric>(B)) + return ex1; + + // Checks for univariate polynomial + if (vars.size() == 1) { + ex ret = euclid_gcd(A, B, vars[0], p); // Univariate GCD + return ret; + } + const ex& mainvar(vars.back()); + + // gcd of the contents + ex H = 0, Hprev = 0; // GCD candidate + ex newton_poly = 1; // for Newton Interpolation + + // Contents and primparts of A and B + ex Aprim, contA; + primpart_content(Aprim, contA, A, vars, p); + ex Bprim, contB; + primpart_content(Bprim, contB, B, vars, p); + // gcd of univariate polynomials + const ex cont_gcd = euclid_gcd(contA, contB, mainvar, p); + + exvector restvars = vars; + restvars.pop_back(); + const ex AL = lcoeff_wrt(Aprim, restvars); + const ex BL = lcoeff_wrt(Bprim, restvars); + // gcd of univariate polynomials + const ex lc_gcd = euclid_gcd(AL, BL, mainvar, p); + + // The estimate of degree of the gcd of Ab and Bb + int gcd_deg = std::min(degree(Aprim, mainvar), degree(Bprim, mainvar)); + eval_point_finder::value_type b; + + eval_point_finder find_eval_point(p); + const numeric pn(p); + do { + // Find a `good' evaluation point b. + bool has_more_pts = find_eval_point(b, lc_gcd, mainvar); + // If there are no more possible evaluation points, bail out + if (!has_more_pts) + throw pgcd_failed(); + + const numeric bn(b); + // Evaluate the polynomials in b + ex Ab = Aprim.subs(mainvar == bn).smod(pn); + ex Bb = Bprim.subs(mainvar == bn).smod(pn); + ex Cb = pgcd(Ab, Bb, restvars, p); + + // Set the correct the leading coefficient + const cln::cl_I lcb_gcd = + smod(to_cl_I(lc_gcd.subs(mainvar == bn)), p); + const cln::cl_I Cblc = integer_lcoeff(Cb, restvars); + const cln::cl_I correct_lc = smod(lcb_gcd*recip(Cblc, p), p); + Cb = (Cb*numeric(correct_lc)).smod(pn); + + // Test for unlucky homomorphisms + const int img_gcd_deg = Cb.degree(restvars.back()); + if (img_gcd_deg < gcd_deg) { + // The degree decreased, previous homomorphisms were + // bad, so we have to start it all over. + H = Cb; + newton_poly = mainvar - numeric(b); + Hprev = 0; + gcd_deg = img_gcd_deg; + continue; + } + if (img_gcd_deg > gcd_deg) { + // The degree of images GCD is too high, this + // evaluation point is bad. Skip it. + continue; + } + + // Image has the same degree as the previous one + // (or at least not higher than the limit) + Hprev = H; + H = newton_interp(Cb, b, H, newton_poly, mainvar, p); + newton_poly = newton_poly*(mainvar - b); + + // try to reduce the number of division tests. + const ex H_lcoeff = lcoeff_wrt(H, restvars); + + if (H_lcoeff.is_equal(lc_gcd)) { + if ((Hprev-H).expand().smod(pn).is_zero()) + continue; + ex C /* primitive part of H */, contH /* dummy */; + primpart_content(C, contH, H, vars, p); + // Normalize GCD so that leading coefficient is 1 + const cln::cl_I Clc = recip(integer_lcoeff(C, vars), p); + C = (C*numeric(Clc)).expand().smod(pn); + + ex dummy1, dummy2; + + if (divide_in_z_p(Aprim, C, dummy1, vars, p) && + divide_in_z_p(Bprim, C, dummy2, vars, p)) + return (cont_gcd*C).expand().smod(pn); + else if (img_gcd_deg == 0) + return cont_gcd; + // else continue building the candidate + } + } while(true); + throw pgcd_failed(); +} + +} // namespace GiNaC + diff --git a/ginac/polynomial/pgcd.h b/ginac/polynomial/pgcd.h new file mode 100644 index 0000000..5e72038 --- /dev/null +++ b/ginac/polynomial/pgcd.h @@ -0,0 +1,27 @@ +#ifndef GINAC_CHINREM_GCD_PGCD_H +#define GINAC_CHINREM_GCD_PGCD_H +#include "ex.h" + +namespace GiNaC +{ + +/// Exception to be thrown when modular GCD algorithm fails +struct pgcd_failed +{ + virtual ~pgcd_failed() { } +}; + +/** + * @brief Compute the GCD of two polynomials over a prime field Z_p + * + * @param vars variables + * @param p designates the coefficient field Z_p + * @param A polynomial \in Z_p[vars] + * @param B second polynomial \in Z_p[vars] + */ +extern ex +pgcd(const ex& A, const ex& B, const exvector& vars, const long p); + +} // namespace GiNaC + +#endif diff --git a/ginac/polynomial/poly_cra.h b/ginac/polynomial/poly_cra.h new file mode 100644 index 0000000..cbcff7b --- /dev/null +++ b/ginac/polynomial/poly_cra.h @@ -0,0 +1,38 @@ +#ifndef GINAC_POLY_CRA_H +#define GINAC_POLY_CRA_H +#include "ex.h" +#include <cln/integer.h> +#include "smod_helpers.h" + +namespace GiNaC +{ + +/** + * @brief Chinese reamainder algorithm for polynomials. + * + * Given two polynomials \f$e_1 \in Z_{q_1}[x_1, \ldots, x_n]\f$ and + * \f$e_2 \in Z_{q_2}[x_1, \ldots, x_n]\f$, compute the polynomial + * \f$r \in Z_{q_1 q_2}[x_1, \ldots, x_n]\f$ such that \f$ r mod q_1 = e_1\f$ + * and \f$ r mod q_2 = e_2 \f$ + */ +ex chinese_remainder(const ex& e1, const cln::cl_I& q1, + const ex& e2, const long q2) +{ + // res = v_1 + v_2 q_1 + // v_1 = e_1 mod q_1 + // v_2 = (e_2 - v_1)/q_1 mod q_2 + const numeric q2n(q2); + const numeric q1n(q1); + ex v1 = e1.smod(q1n); + ex u = v1.smod(q2n); + ex v2 = (e2.smod(q2n) - v1.smod(q2n)).expand().smod(q2n); + const numeric q1_1(recip(q1, q2)); // 1/q_1 mod q_2 + v2 = (v2*q1_1).smod(q2n); + ex ret = (v1 + v2*q1_1).expand(); + return ret; +} + +} // namespace GiNaC + +#endif /* GINAC_POLY_CRA_H */ + diff --git a/ginac/polynomial/primes_factory.h b/ginac/polynomial/primes_factory.h new file mode 100644 index 0000000..093c973 --- /dev/null +++ b/ginac/polynomial/primes_factory.h @@ -0,0 +1,60 @@ +#ifndef GINAC_CHINREM_GCD_PRIMES_FACTORY_H +#define GINAC_CHINREM_GCD_PRIMES_FACTORY_H +#include <cln/integer.h> +#include <cln/numtheory.h> +#include <limits> +#include "smod_helpers.h" +#include "debug.hpp" + +namespace GiNaC +{ + +/** + * Find a `big' prime p such that lc mod p != 0. Helper class used by modular + * GCD algorithm. + */ +class primes_factory +{ +private: + // These primes need to be large enough, so that the number of images + // we need to reconstruct the GCD (in Z) is reasonable. On the other + // hand, they should be as small as possible, so that operations on + // coefficients are efficient. Practically this means we coefficients + // should be native integers. (N.B.: as of now chinrem_gcd uses cl_I + // or even numeric. Eventually this will be fixed). + cln::cl_I last; + // This ensures coefficients are immediate. + static const int immediate_bits = 8*sizeof(void *) - __alignof__(void *); + static const long opt_hint = (1L << (immediate_bits >> 1)) - 1; +public: + primes_factory() + { + last = cln::nextprobprime(cln::cl_I(opt_hint)); + } + + bool operator()(long& p, const cln::cl_I& lc) + { + static const cln::cl_I maxval(std::numeric_limits<long>::max()); + while (last < maxval) { + long p_ = cln::cl_I_to_long(last); + last = cln::nextprobprime(last + 1); + + if (!zerop(smod(lc, p_))) { + p = p_; + return true; + } + } + return false; + } + + bool has_primes() const + { + static const cln::cl_I maxval(std::numeric_limits<long>::max()); + return last < maxval; + } +}; + +} // namespace GiNaC + +#endif /* GINAC_CHINREM_GCD_PRIMES_FACTORY_H */ + diff --git a/ginac/polynomial/primpart_content.cpp b/ginac/polynomial/primpart_content.cpp new file mode 100644 index 0000000..9d36a95 --- /dev/null +++ b/ginac/polynomial/primpart_content.cpp @@ -0,0 +1,78 @@ +#include "ex.h" +#include "numeric.h" +#include "collect_vargs.h" +#include "euclid_gcd_wrap.h" +#include "divide_in_z_p.h" +#include "debug.hpp" + +namespace GiNaC +{ + +/** + * Compute the primitive part and the content of a modular multivariate + * polynomial e \in Z_p[x_n][x_0, \ldots, x_{n-1}], i.e. e is considered + * a polynomial in variables x_0, \ldots, x_{n-1} with coefficients being + * modular polynomials Z_p[x_n] + * @param e polynomial to operate on + * @param pp primitive part of @a e, will be computed by this function + * @param c content (in the sense described above) of @a e, will be + * computed by this function + * @param vars variables x_0, \ldots, x_{n-1}, x_n + * @param p modulus + */ +void primpart_content(ex& pp, ex& c, ex e, const exvector& vars, + const long p) +{ + static const ex ex1(1); + static const ex ex0(0); + e = e.expand(); + if (e.is_zero()) { + pp = ex0; + c = ex1; + return; + } + exvector rest_vars = vars; + rest_vars.pop_back(); + ex_collect_t ec; + collect_vargs(ec, e, rest_vars); + + if (ec.size() == 1) { + // the input polynomial factorizes into + // p_1(x_n) p_2(x_0, \ldots, x_{n-1}) + c = ec.rbegin()->second; + ec.rbegin()->second = ex1; + pp = ex_collect_to_ex(ec, vars).expand().smod(numeric(p)); + return; + } + + // Start from the leading coefficient (which is stored as a last + // element of the terms array) + ex_collect_t::reverse_iterator i = ec.rbegin(); + ex g = i->second; + // there are at least two terms, so it's safe to... + ++i; + while (i != ec.rend() && !g.is_equal(ex1)) { + g = euclid_gcd(i->second, g, vars.back(), p); + ++i; + } + if (g.is_equal(ex1)) { + pp = e; + c = ex1; + return; + } + exvector mainvar; + mainvar.push_back(vars.back()); + for (i = ec.rbegin(); i != ec.rend(); ++i) { + ex tmp(0); + if (!divide_in_z_p(i->second, g, tmp, mainvar, p)) + throw std::logic_error(std::string(__func__) + + ": bogus division failure"); + i->second = tmp; + } + + pp = ex_collect_to_ex(ec, rest_vars).expand().smod(numeric(p)); + c = g; +} + +} // namespace GiNaC + diff --git a/ginac/polynomial/smod_helpers.h b/ginac/polynomial/smod_helpers.h new file mode 100644 index 0000000..81370ba --- /dev/null +++ b/ginac/polynomial/smod_helpers.h @@ -0,0 +1,72 @@ +#ifndef GINAC_POLYNOMIAL_SMOD_HELPERS_H +#define GINAC_POLYNOMIAL_SMOD_HELPERS_H +#include <cln/integer.h> +#include <cln/integer_io.h> +#include "ex.h" +#include "numeric.h" +#include "debug.hpp" + +namespace GiNaC +{ + +/// Z -> Z_p (in the symmetric representation) +static inline cln::cl_I smod(const cln::cl_I& a, long p) +{ + const cln::cl_I p2 = cln::cl_I(p >> 1); + const cln::cl_I m = cln::mod(a, p); + const cln::cl_I m_p = m - cln::cl_I(p); + const cln::cl_I ret = m > p2 ? m_p : m; + return ret; +} + +static inline cln::cl_I recip(const cln::cl_I& a, long p_) +{ + cln::cl_I p(p_); + cln::cl_I u, v; + const cln::cl_I g = xgcd(a, p, &u, &v); + cln::cl_I ret = smod(u, p_); + cln::cl_I chck = smod(a*ret, p_); + bug_on(chck != 1, "miscomputed recip(" << a << " (mod " << p_ << "))"); + return ret; + +} + +static inline numeric recip(const numeric& a_, long p) +{ + const cln::cl_I a = cln::the<cln::cl_I>(a_.to_cl_N()); + const cln::cl_I ret = recip(a, p); + return numeric(ret); +} + +static inline cln::cl_I to_cl_I(const ex& e) +{ + bug_on(!is_a<numeric>(e), "argument should be an integer"); + bug_on(!e.info(info_flags::integer), + "argument should be an integer"); + return cln::the<cln::cl_I>(ex_to<numeric>(e).to_cl_N()); +} + +struct random_modint +{ + typedef long value_type; + const value_type p; + const value_type p_2; + + random_modint(const value_type& p_) : p(p_), p_2((p >> 1)) + { } + value_type operator()() const + { + do { + cln::cl_I tmp_ = cln::random_I(p); + value_type tmp = cln::cl_I_to_long(tmp_); + if (tmp > p_2) + tmp -= p; + return tmp; + } while (true); + } + +}; + +} // namespace GiNaC + +#endif // GINAC_POLYNOMIAL_SMOD_HELPERS_H -- 1.5.6.5