Implement modular multivariate GCD (based on chinese remaindering algorithm).
authorAlexei Sheplyakov <varg@metalica.kh.ua>
Mon, 19 Jan 2009 06:45:38 +0000 (08:45 +0200)
committerJens Vollinga <jensv@balin.nikhef.nl>
Tue, 3 Feb 2009 12:16:23 +0000 (13:16 +0100)
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.

22 files changed:
ginac/Makefile.am
ginac/normal.cpp
ginac/normal.h
ginac/polynomial/chinrem_gcd.cpp [new file with mode: 0644]
ginac/polynomial/chinrem_gcd.h [new file with mode: 0644]
ginac/polynomial/collect_vargs.cpp [new file with mode: 0644]
ginac/polynomial/collect_vargs.h [new file with mode: 0644]
ginac/polynomial/debug.hpp
ginac/polynomial/divide_in_z_p.cpp [new file with mode: 0644]
ginac/polynomial/divide_in_z_p.h [new file with mode: 0644]
ginac/polynomial/euclid_gcd_wrap.h [new file with mode: 0644]
ginac/polynomial/eval_point_finder.h [new file with mode: 0644]
ginac/polynomial/mgcd.cpp [new file with mode: 0644]
ginac/polynomial/newton_interpolate.h [new file with mode: 0644]
ginac/polynomial/optimal_vars_finder.cpp [new file with mode: 0644]
ginac/polynomial/optimal_vars_finder.h [new file with mode: 0644]
ginac/polynomial/pgcd.cpp [new file with mode: 0644]
ginac/polynomial/pgcd.h [new file with mode: 0644]
ginac/polynomial/poly_cra.h [new file with mode: 0644]
ginac/polynomial/primes_factory.h [new file with mode: 0644]
ginac/polynomial/primpart_content.cpp [new file with mode: 0644]
ginac/polynomial/smod_helpers.h [new file with mode: 0644]

index 7ab16d6..76e32b2 100644 (file)
@@ -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)
index 8f5ba73..54fff47 100644 (file)
@@ -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)
index 5a735e3..a352be8 100644 (file)
@@ -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 (file)
index 0000000..4048ee0
--- /dev/null
@@ -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 (file)
index 0000000..6c58a9f
--- /dev/null
@@ -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 (file)
index 0000000..5649e4a
--- /dev/null
@@ -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 (file)
index 0000000..4dd16ba
--- /dev/null
@@ -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 */
index 901ef87..48659fd 100644 (file)
@@ -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 (file)
index 0000000..dff30a9
--- /dev/null
@@ -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 (file)
index 0000000..f790fcd
--- /dev/null
@@ -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 (file)
index 0000000..d77f71c
--- /dev/null
@@ -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 (file)
index 0000000..8ceadd1
--- /dev/null
@@ -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 (file)
index 0000000..b39475c
--- /dev/null
@@ -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 (file)
index 0000000..e80687d
--- /dev/null
@@ -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 (file)
index 0000000..4903699
--- /dev/null
@@ -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 (file)
index 0000000..9715a47
--- /dev/null
@@ -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 (file)
index 0000000..ddc3b8c
--- /dev/null
@@ -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 (file)
index 0000000..5e72038
--- /dev/null
@@ -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 (file)
index 0000000..cbcff7b
--- /dev/null
@@ -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 (file)
index 0000000..093c973
--- /dev/null
@@ -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 (file)
index 0000000..9d36a95
--- /dev/null
@@ -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 (file)
index 0000000..81370ba
--- /dev/null
@@ -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