[GiNaC-devel] [PATCH] Implement modular multivariate GCD (based on chinese remaindering algorithm).
Alexei Sheplyakov
varg at metalica.kh.ua
Mon Jan 19 07:45:38 CET 2009
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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 835 bytes
Desc: Digital signature
URL: <http://www.ginac.de/pipermail/ginac-devel/attachments/20090119/28b09c12/attachment.sig>
More information about the GiNaC-devel
mailing list