Implemented modular GCD algorithm for univariate polynomials.
authorAlexei Sheplyakov <varg@theor.jinr.ru>
Sat, 20 Sep 2008 19:18:46 +0000 (23:18 +0400)
committerAlexei Sheplyakov <varg@theor.jinr.ru>
Sat, 20 Sep 2008 19:18:46 +0000 (23:18 +0400)
16 files changed:
check/Makefile.am
check/exam_cra.cpp [new file with mode: 0644]
check/exam_mod_gcd.cpp [new file with mode: 0644]
ginac/Makefile.am
ginac/polynomial/cra_garner.cpp [new file with mode: 0644]
ginac/polynomial/cra_garner.hpp [new file with mode: 0644]
ginac/polynomial/debug.hpp [new file with mode: 0644]
ginac/polynomial/gcd_euclid.tcc [new file with mode: 0644]
ginac/polynomial/mod_gcd.cpp [new file with mode: 0644]
ginac/polynomial/mod_gcd.hpp [new file with mode: 0644]
ginac/polynomial/normalize.tcc [new file with mode: 0644]
ginac/polynomial/remainder.tcc [new file with mode: 0644]
ginac/polynomial/ring_traits.hpp [new file with mode: 0644]
ginac/polynomial/upoly.hpp [new file with mode: 0644]
ginac/polynomial/upoly_io.cpp [new file with mode: 0644]
ginac/polynomial/upoly_io.hpp [new file with mode: 0644]

index 87b489d..5c24aee 100644 (file)
@@ -27,7 +27,9 @@ EXAMS = exam_paranoia \
        exam_archive  \
        exam_structure  \
        exam_hashmap  \
-       exam_misc
+       exam_misc \
+       exam_mod_gcd \
+       exam_cra
 
 TIMES = time_dennyfliegner \
        time_gammaseries \
@@ -140,6 +142,11 @@ exam_hashmap_LDADD = ../ginac/libginac.la
 exam_misc_SOURCES = exam_misc.cpp
 exam_misc_LDADD = ../ginac/libginac.la
 
+exam_mod_gcd_SOURCES = exam_mod_gcd.cpp
+exam_mod_gcd_LDADD = ../ginac/libginac.la
+
+exam_cra_SOURCES = exam_cra.cpp
+exam_cra_LDADD = ../ginac/libginac.la
 
 time_dennyfliegner_SOURCES = time_dennyfliegner.cpp \
                             randomize_serials.cpp timer.cpp timer.h
diff --git a/check/exam_cra.cpp b/check/exam_cra.cpp
new file mode 100644 (file)
index 0000000..11e1ab7
--- /dev/null
@@ -0,0 +1,120 @@
+#include <iostream>
+#include <cln/integer.h>
+#include <cln/integer_io.h>
+#include <cln/random.h>
+#include <cln/numtheory.h>
+#include "polynomial/cra_garner.hpp"
+#include <map>
+#include <vector>
+#include <stdexcept>
+using namespace cln;
+using namespace std;
+
+/// Generate a sequences of primes p_i such that \prod_i p_i < limit
+static std::vector<cln::cl_I>
+make_random_moduli(const cln::cl_I& limit);
+
+static std::vector<cln::cl_I>
+calc_residues(const cln::cl_I& x, const std::vector<cln::cl_I>& moduli);
+
+static void dump(const std::vector<cln::cl_I>& v);
+
+/// Make @a n random relatively prime moduli, each < limit, make a 
+/// random number x < \prod_{i=0}{n-1}, calculate residues, and 
+/// compute x' by chinese remainder algorithm. Check if the result
+/// of computation matches the original value x.
+static void run_test_once(const cln::cl_I& lim)
+{
+       std::vector<cln::cl_I> moduli = make_random_moduli(lim);
+       cln::cl_I x = random_I(lim) + 1;
+
+       if (x > (lim >> 1))
+               x = x - lim;
+
+       std::vector<cln::cl_I> residues = calc_residues(x, moduli);
+       cln::cl_I x_test;
+
+       bool error = false;
+       try {
+               x_test = integer_cra(residues, moduli);
+       } catch (std::exception& oops) {
+               std::cerr << "Oops: " << oops.what() << std::endl;
+               error = true;
+       }
+
+       if (x != x_test)
+               error = true;
+
+       if (error) {
+               std::cerr << "Expected x = " << x << ", got " <<
+                       x_test << " instead" << std::endl;
+               std::cerr << "moduli = ";
+               dump(moduli);
+               std::cerr << std::endl;
+               std::cerr << "residues = ";
+               dump(residues);
+               std::cerr << std::endl;
+               throw std::logic_error("bug in integer_cra?");
+       }
+}
+
+static void run_test(const cln::cl_I& limit, const std::size_t ntimes)
+{
+       for (std::size_t i = 0; i < ntimes; ++i)
+               run_test_once(limit);
+}
+
+int main(int argc, char** argv)
+{
+       typedef std::map<cln::cl_I, std::size_t> map_t;
+       map_t the_map;
+       // Run 1024 tests with native 32-bit numbers
+       the_map[cln::cl_I(std::numeric_limits<int>::max())] = 1024;
+
+       // Run 512 tests with native 64-bit integers
+       if (sizeof(long) > sizeof(int)) 
+               the_map[cln::cl_I(std::numeric_limits<long>::max())] = 512;
+
+       // Run 32 tests with a bit bigger numbers
+       the_map[cln::cl_I("987654321098765432109876543210")] = 32;
+
+       std::cout << "examining Garner's integer chinese remainder algorithm " << std::flush;
+
+       for (map_t::const_iterator i = the_map.begin(); i != the_map.end(); ++i)
+               run_test(i->first, i->second);
+
+       return 0;
+}
+
+static std::vector<cln::cl_I>
+calc_residues(const cln::cl_I& x, const std::vector<cln::cl_I>& moduli)
+{
+       std::vector<cln::cl_I> residues(moduli.size());
+       for (std::size_t i = moduli.size(); i-- != 0; )
+               residues[i] = mod(x, moduli[i]);
+       return residues;
+}
+
+static std::vector<cln::cl_I>
+make_random_moduli(const cln::cl_I& limit)
+{
+       std::vector<cln::cl_I> moduli;
+       cln::cl_I prod(1);
+       cln::cl_I next = random_I(std::min(limit >> 1, cln::cl_I(128)));
+       do {
+               cln::cl_I tmp = nextprobprime(next);
+               next = tmp + random_I(cln::cl_I(10)) + 1;
+               prod = prod*tmp;
+               moduli.push_back(tmp);
+       } while (prod < limit);
+       return moduli;
+}
+
+static void dump(const std::vector<cln::cl_I>& v)
+{
+       std::cerr << "[ ";
+       for (std::size_t i = 0; i < v.size(); ++i)
+               std::cerr << v[i] << " ";
+       std::cerr << "]";
+}
+
diff --git a/check/exam_mod_gcd.cpp b/check/exam_mod_gcd.cpp
new file mode 100644 (file)
index 0000000..4407d90
--- /dev/null
@@ -0,0 +1,85 @@
+#include "polynomial/upoly.hpp"
+#include "polynomial/upoly_io.hpp"
+#include "polynomial/mod_gcd.hpp"
+#include "ginac.h"
+#include <cln/random.h>
+#include <string>
+#include <iostream>
+#include <map>
+using namespace GiNaC;
+
+static upoly ex_to_upoly(const ex& e, const symbol& x);
+static ex upoly_to_ex(const upoly& p, const symbol& x);
+
+// make a univariate polynomial \in Z[x] of degree deg
+static upoly make_random_upoly(const std::size_t deg);
+
+static void run_test_once(const std::size_t deg)
+{
+       static const symbol xsym("x");
+
+       const upoly a = make_random_upoly(deg);
+       const upoly b = make_random_upoly(deg);
+
+       upoly g;
+       mod_gcd(g, a, b);
+
+       ex ea = upoly_to_ex(a, xsym);
+       ex eb = upoly_to_ex(b, xsym);
+       ex eg = gcd(ea, eb);
+
+       const upoly g_check = ex_to_upoly(eg, xsym);
+       if (g != g_check) {
+               std::cerr << "a = " << a << std::endl;
+               std::cerr << "b = " << b << std::endl;
+               std::cerr << "mod_gcd(a, b) = " << g << std::endl;
+               std::cerr << "sr_gcd(a, b) = " << g_check << std::endl;
+               throw std::logic_error("bug in mod_gcd");
+       }
+}
+
+int main(int argc, char** argv)
+{
+       std::cout << "examining modular gcd. ";
+       std::map<std::size_t, std::size_t> n_map;
+       // run 256 tests with polynomials of degree 10
+       n_map[10] = 256;
+       // run 32 tests with polynomials of degree 100
+       n_map[100] = 32;
+       std::map<std::size_t, std::size_t>::const_iterator i = n_map.begin();
+       for (; i != n_map.end(); ++i) {
+               for (std::size_t k = 0; k < i->second; ++k)
+                       run_test_once(i->first);
+       }
+       return 0;
+}
+
+static upoly ex_to_upoly(const ex& e, const symbol& x)
+{
+       upoly p(e.degree(x) + 1);
+       for (int i = 0; i <= e.degree(x); ++i)
+               p[i] = cln::the<cln::cl_I>(ex_to<numeric>(e.coeff(x, i)).to_cl_N());
+       return p;
+}
+
+static ex upoly_to_ex(const upoly& p, const symbol& x)
+{
+       exvector tv(p.size());
+       for (std::size_t i = 0; i < p.size(); ++i)
+               tv[i] = pow(x, i)*numeric(p[i]);
+       return (new add(tv))->setflag(status_flags::dynallocated);
+}
+
+static upoly make_random_upoly(const std::size_t deg)
+{
+       static const cln::cl_I biggish("98765432109876543210");
+       upoly p(deg + 1);
+       for (std::size_t i = 0; i <= deg; ++i)
+               p[i] = cln::random_I(biggish);
+
+       // Make sure the leading coefficient is non-zero
+       while (zerop(p[deg])) 
+               p[deg] = cln::random_I(biggish);
+       return p;
+}
+
index db09d38..b683d9b 100644 (file)
@@ -17,8 +17,20 @@ libginac_la_SOURCES = add.cpp archive.cpp basic.cpp clifford.cpp color.cpp \
   parser/lexer.cpp \
   parser/lexer.hpp \
   parser/parser_compat.cpp \
-  parser/debug.hpp
-  
+  parser/debug.hpp \
+polynomial/mod_gcd.cpp \
+polynomial/cra_garner.cpp \
+polynomial/gcd_euclid.tcc \
+polynomial/remainder.tcc \
+polynomial/normalize.tcc \
+polynomial/upoly.hpp \
+polynomial/ring_traits.hpp \
+polynomial/mod_gcd.hpp \
+polynomial/cra_garner.hpp \
+polynomial/upoly_io.hpp \
+polynomial/upoly_io.cpp \
+polynomial/debug.hpp
+
 libginac_la_LDFLAGS = -version-info $(LT_VERSION_INFO) -release $(LT_RELEASE)
 libginac_la_LIBADD = $(DL_LIBS)
 ginacincludedir = $(includedir)/ginac
diff --git a/ginac/polynomial/cra_garner.cpp b/ginac/polynomial/cra_garner.cpp
new file mode 100644 (file)
index 0000000..b400adb
--- /dev/null
@@ -0,0 +1,88 @@
+#include <cln/integer.h>
+#include <cln/modinteger.h>
+#include <vector>
+#include <cstddef>
+#include "cra_garner.hpp"
+
+namespace cln
+{
+using std::vector;
+using std::size_t;
+
+static cl_I 
+retract_symm(const cl_MI& x, const cl_modint_ring& R,
+            const cl_I& modulus)
+{
+       cl_I result = R->retract(x);
+       if (result > (modulus >> 1))
+               result = result - modulus;
+       return result;
+}
+
+static void
+compute_recips(vector<cl_MI>& dst, 
+              const vector<cl_I>& moduli)
+{
+       for (size_t k = 1; k < moduli.size(); ++k) {
+               cl_modint_ring R = find_modint_ring(moduli[k]);
+               cl_MI product = R->canonhom(moduli[0]);
+               for (size_t i = 1; i < k; ++i)
+                       product = product*moduli[i];
+               dst[k-1] = recip(product);
+       }
+}
+
+static void
+compute_mix_radix_coeffs(vector<cl_I>& dst,
+                        const vector<cl_I>& residues,
+                        const vector<cl_I>& moduli,
+                        const vector<cl_MI>& recips)
+{
+       dst[0] = residues[0];
+
+       do {
+               cl_modint_ring R = find_modint_ring(moduli[1]);
+               cl_MI tmp = R->canonhom(residues[0]);
+               cl_MI next = (R->canonhom(residues[1]) - tmp)*recips[0];
+               dst[1] = retract_symm(next, R, moduli[1]);
+       } while (0);
+
+       for (size_t k = 2; k < residues.size(); ++k) {
+               cl_modint_ring R = find_modint_ring(moduli[k]);
+               cl_MI tmp = R->canonhom(dst[k-1]);
+
+               for (size_t j = k - 1 /* NOT k - 2 */; j-- != 0; )
+                       tmp = tmp*moduli[j] + R->canonhom(dst[j]);
+
+               cl_MI next = (R->canonhom(residues[k]) - tmp)*recips[k-1];
+               dst[k] = retract_symm(next, R, moduli[k]);
+       }
+}
+
+static cl_I
+mixed_radix_2_ordinary(const vector<cl_I>& mixed_radix_coeffs,
+                      const vector<cl_I>& moduli)
+{
+       size_t k = mixed_radix_coeffs.size() - 1;
+       cl_I u = mixed_radix_coeffs[k];
+       for (; k-- != 0; )
+               u = u*moduli[k] + mixed_radix_coeffs[k];
+       return u;
+}
+
+cl_I integer_cra(const vector<cl_I>& residues,
+                const vector<cl_I>& moduli)
+{
+
+       vector<cl_MI> recips(moduli.size() - 1);
+       compute_recips(recips, moduli);
+
+       vector<cl_I> coeffs(moduli.size());
+       compute_mix_radix_coeffs(coeffs, residues, moduli, recips);
+       cl_I result = mixed_radix_2_ordinary(coeffs, moduli);
+
+       return result;
+}
+
+} // namespace cln
+
diff --git a/ginac/polynomial/cra_garner.hpp b/ginac/polynomial/cra_garner.hpp
new file mode 100644 (file)
index 0000000..eff8e9b
--- /dev/null
@@ -0,0 +1,12 @@
+#ifndef CL_INTEGER_CRA
+#define CL_INTEGER_CRA
+#include <vector>
+#include <cln/integer.h>
+
+namespace cln
+{
+extern cl_I integer_cra(const std::vector<cl_I>& residues,
+                       const std::vector<cl_I>& moduli);
+}
+#endif // CL_INTEGER_CRA
+
diff --git a/ginac/polynomial/debug.hpp b/ginac/polynomial/debug.hpp
new file mode 100644 (file)
index 0000000..901ef87
--- /dev/null
@@ -0,0 +1,29 @@
+#ifndef GINAC_MOD_GCD_DEBUG_HPP
+#define GINAC_MOD_GCD_DEBUG_HPP
+#include <iostream>
+#include <string>
+#include <stdexcept>
+#include <sstream>
+
+#define DEBUG_PREFIX __func__ << ':' << __LINE__ << ": "
+#define EXCEPTION_PREFIX std::string(__func__) + std::string(": ") +
+
+#define Dout2(stream, msg) \
+do {                                                                   \
+       stream << DEBUG_PREFIX << msg << std::endl << std::flush;       \
+} while (0)
+#define Dout(msg) Dout2(std::cout, msg)
+
+#define bug3_on(condition, the_exception, msg)                         \
+do {                                                                   \
+       if (unlikely(condition)) {                                      \
+               std::ostringstream err_stream;                          \
+               Dout2(err_stream, "BUG: " << msg);                      \
+               throw the_exception(err_stream.str());                  \
+       }                                                               \
+} while (0)
+
+#define bug_on(condition, msg) bug3_on(condition, std::logic_error, msg)
+
+#endif // GINAC_MOD_GCD_DEBUG_HPP
+
diff --git a/ginac/polynomial/gcd_euclid.tcc b/ginac/polynomial/gcd_euclid.tcc
new file mode 100644 (file)
index 0000000..deae252
--- /dev/null
@@ -0,0 +1,45 @@
+#ifndef GINAC_GCD_EUCLID
+#define GINAC_GCD_EUCLID
+#include "upoly.hpp"
+#include "remainder.tcc"
+#include "normalize.tcc"
+#include "debug.hpp"
+#include "upoly_io.hpp"
+
+namespace GiNaC
+{
+
+static bool
+gcd_euclid(umodpoly& c, umodpoly /* passed by value */ a, umodpoly b)
+{
+       if (a.size() == 0) {
+               c.clear();
+               return true;
+       }
+       if (b.size() == 0) {
+               c.clear();
+               return true;
+       }
+       bug_on(a[0].ring()->modulus != b[0].ring()->modulus,
+               "different moduli");
+
+       normalize_in_field(a);
+       normalize_in_field(b);
+       if (degree(a) < degree(b))
+               std::swap(a, b);
+
+       umodpoly r;
+       while (b.size() != 0) {
+               remainder_in_field(r, a, b); 
+               a = b;
+               b = r;
+       }
+       normalize_in_field(a);
+       c = a;
+       return false;
+}
+
+} // namespace GiNaC
+
+#endif // GINAC_GCD_EUCLID
+
diff --git a/ginac/polynomial/mod_gcd.cpp b/ginac/polynomial/mod_gcd.cpp
new file mode 100644 (file)
index 0000000..8aa8480
--- /dev/null
@@ -0,0 +1,165 @@
+#include "upoly.hpp"
+#include "gcd_euclid.tcc"
+#include "cra_garner.hpp"
+#include <cln/random.h>
+#include <cln/numtheory.h>
+#include <stdexcept>
+#include "debug.hpp"
+
+namespace GiNaC
+{
+/**
+ * @brief Remove the integer content from univariate polynomials A and B.
+ *
+ * As a byproduct compute the GCD of contents.
+ */
+static void remove_content(upoly& A, upoly& B, upoly::value_type& c)
+{
+       // remove the integer
+       upoly::value_type acont, bcont;
+       normalize_in_ring(A, &acont);
+       normalize_in_ring(B, &bcont);
+       c = gcd(acont, bcont);
+};
+
+/// Check if @a candidate divides both @a A and @a B
+static bool
+do_division_check(const upoly& A, const upoly& B, const upoly& candidate)
+{
+       upoly r1;
+       remainder_in_ring(r1, A, candidate);
+       if (r1.size() != 0)
+               return false;
+
+       upoly r2;
+       remainder_in_ring(r2, B, candidate);
+       if (r2.size() != 0)
+               return false;
+
+       return true;
+}
+
+/**
+ * Given two GCD candidates H \in Z/q[x], and C \in Z/p[x] (where p is a prime)
+ * compute the candidate in Z/(q*p) with CRA (chinise remainder algorithm)
+ *
+ * @param H \in Z/q[x] GCD candidate, will be updated by this function
+ * @param q modulus of H, will NOT be updated by this function
+ * @param C \in Z/p[x] GCD candidate, will be updated by this function
+ * @param p modulus of C
+ */
+static void
+update_the_candidate(upoly& H, const upoly::value_type& q,
+                    const umodpoly& C,
+                    const upoly::value_type& p,
+                    const cln::cl_modint_ring& R)
+{
+       typedef upoly::value_type ring_t;
+       std::vector<ring_t> moduli(2);
+       moduli[0] = q;
+       moduli[1] = p;
+       if (H.size() < C.size())
+               H.resize(C.size());
+
+       for (std::size_t  i = C.size(); i-- != 0; ) {
+               std::vector<ring_t> coeffs(2);
+               coeffs[0] = H[i];
+               coeffs[1] = R->retract(C[i]); 
+               H[i] = integer_cra(coeffs, moduli);
+       }
+}
+
+/// Convert Z/p[x] -> Z[x]
+static void retract(upoly& p, const umodpoly& cp,
+                   const cln::cl_modint_ring& Rp)
+{
+       p.resize(cp.size());
+       for (std::size_t i = cp.size(); i-- != 0; )
+               p[i] = Rp->retract(cp[i]);
+}
+
+
+/// Find the prime which is > p, and does NOT divide g
+static void find_next_prime(cln::cl_I& p, const cln::cl_I& g)
+{
+       do {
+               ++p;
+               p = nextprobprime(p);
+       } while (zerop(mod(g, p)));
+}
+
+/// Compute the GCD of univariate polynomials A, B \in Z[x]
+void mod_gcd(upoly& result, upoly A, upoly B)
+{
+       typedef upoly::value_type ring_t;
+       ring_t content_gcd;
+       // remove the integer content
+       remove_content(A, B, content_gcd);
+
+       // compute the coefficient bound for GCD(a, b)
+       ring_t g = gcd(lcoeff(A), lcoeff(B));
+       std::size_t max_gcd_degree = std::min(degree(A), degree(B));
+       ring_t limit = (ring_t(1) << max_gcd_degree)*g*
+                      std::min(max_coeff(A), max_coeff(B));
+       ring_t q(0);
+       upoly H;
+
+       ring_t p(2);
+       while (true) {
+               find_next_prime(p, g);
+
+               // Map the polynomials onto Z/p[x]
+               cln::cl_modint_ring Rp = cln::find_modint_ring(p);
+               cln::cl_MI gp = Rp->canonhom(g);
+               umodpoly ap(A.size()), bp(B.size());
+               make_umodpoly(ap, A, Rp);
+               make_umodpoly(bp, B, Rp);
+
+               // Compute the GCD in Z/p[x]
+               umodpoly cp;
+               gcd_euclid(cp, ap, bp);
+               bug_on(cp.size() == 0, "gcd(ap, bp) = 0, with ap = " <<
+                                       ap << ", and bp = " << bp);
+
+
+               // Normalize the candidate so that its leading coefficient
+               // is g mod p
+               umodpoly::value_type norm_factor = gp*recip(lcoeff(cp));
+               bug_on(zerop(norm_factor), "division in a field give 0");
+
+               lcoeff(cp) = gp;
+               for (std::size_t k = cp.size() - 1; k-- != 0; )
+                       cp[k] = cp[k]*norm_factor;
+
+
+               // check for unlucky homomorphisms
+               if (degree(cp) < max_gcd_degree) {
+                       q = p;
+                       max_gcd_degree = degree(cp);
+                       retract(H, cp, Rp);
+               } else {
+                       update_the_candidate(H, q, cp, p, Rp);
+                       q = q*p;
+               }
+
+               if (q > limit) {
+                       result = H;
+                       normalize_in_ring(result);
+                       // if H divides both A and B it's a GCD
+                       if (do_division_check(A, B, result)) {
+                               result *= content_gcd;
+                               break; 
+                       }
+                       // H does not divide A and/or B, look for
+                       // another one
+               } else if (degree(cp) == 0) {
+                       // Polynomials are relatively prime
+                       result.resize(1);
+                       result[0] = content_gcd;
+                       break;
+               }
+       }
+}
+
+} // namespace GiNaC
+
diff --git a/ginac/polynomial/mod_gcd.hpp b/ginac/polynomial/mod_gcd.hpp
new file mode 100644 (file)
index 0000000..43e1e92
--- /dev/null
@@ -0,0 +1,11 @@
+#ifndef GINAC_MOD_GCD_HPP
+#define GINAC_MOD_GCD_HPP
+#include "upoly.hpp"
+
+namespace GiNaC
+{
+extern void mod_gcd(upoly& result, upoly A, upoly B);
+}
+
+#endif // GINAC_MOD_GCD_HPP
+
diff --git a/ginac/polynomial/normalize.tcc b/ginac/polynomial/normalize.tcc
new file mode 100644 (file)
index 0000000..e246307
--- /dev/null
@@ -0,0 +1,93 @@
+#ifndef GINAC_UPOLY_NORMALIZE_TCC
+#define GINAC_UPOLY_NORMALIZE_TCC
+#include "upoly.hpp"
+#include "ring_traits.hpp"
+#include "debug.hpp"
+
+namespace GiNaC
+{
+
+/// Make the univariate polynomial @a a \in F[x] unit normal.
+/// F should be a field.
+/// Returns true if the polynomial @x is already unit normal, and false
+/// otherwise.
+static bool normalize_in_field(umodpoly& a, cln::cl_MI* content_ = 0)
+{
+       if (a.size() == 0)
+               return true;
+       if (lcoeff(a) == the_one(a[0])) {
+               if (content_)
+                       *content_ = the_one(a[0]);
+               return true;
+       }
+
+       const cln::cl_MI lc_1 = recip(lcoeff(a));
+       for (std::size_t k = a.size(); k-- != 0; )
+               a[k] = a[k]*lc_1;
+       if (content_)
+               *content_ = lc_1;
+       return false;
+}
+
+/// Make the univariate polynomial @a x unit normal. This version is used
+/// for rings which are not fields. 
+/// Returns true if the polynomial @x is already unit normal, and false
+/// otherwise.
+template<typename T> bool
+normalize_in_ring(T& x, typename T::value_type* content_ = 0, int* unit_ = 0)
+{
+       typedef typename T::value_type ring_t;
+       static const ring_t one(1);
+       if (x.empty())
+               return true;
+
+       bool something_changed = false;
+       if (minusp(lcoeff(x))) {
+               something_changed = true;
+               if (unit_)
+                       *unit_ = -1;
+               for (std::size_t i = x.size(); i-- != 0; )
+                       x[i] = - x[i];
+       }
+
+       if (degree(x) == 0) {
+               if (content_)
+                       *content_ = x[0];
+               if (x[0] == one)
+                       return something_changed;
+               x[0] = one;
+               return false; // initial polynomial was unit normal
+       }
+               
+       // Compute the gcd of coefficients
+       ring_t content = lcoeff(x);
+       // We want this function to be fast when applied to unit normal
+       // polynomials. Hence we start from the leading coefficient.
+       for (std::size_t i = x.size() - 1; i-- != 0; ) {
+               if (content == one) {
+                       if (content_)
+                               *content_ = one;
+                       return something_changed;
+               }
+               content = gcd(x[i], content);
+       }
+
+       if (content == one) {
+               if (content_)
+                       *content_ = one;
+               return something_changed;
+       }
+       
+       lcoeff(x) = one;
+       for (std::size_t i = x.size() - 1; i-- != 0; )
+               x[i] = exquo(x[i], content);
+
+       if (content_)
+               *content_ = content;
+       return false; // initial polynomial was not unit normal
+}
+
+} // namespace GiNaC
+       
+#endif // GINAC_UPOLY_NORMALIZE_TCC
+
diff --git a/ginac/polynomial/remainder.tcc b/ginac/polynomial/remainder.tcc
new file mode 100644 (file)
index 0000000..4db2843
--- /dev/null
@@ -0,0 +1,116 @@
+#ifndef GINAC_UPOLY_REMAINDER_TCC
+#define GINAC_UPOLY_REMAINDER_TCC
+#include "upoly.hpp"
+#include "ring_traits.hpp"
+#include "upoly_io.hpp"
+#include "debug.hpp"
+
+namespace GiNaC
+{
+/**
+ * @brief Polynomial remainder for univariate polynomials over fields
+ *
+ * Given two univariate polynomials \f$a, b \in F[x]\f$, where F is
+ * a finite field (presumably Z/p) computes the remainder @a r, which is
+ * defined as \f$a = b q + r\f$. Returns true if the remainder is zero
+ * and false otherwise.
+ */
+static bool
+remainder_in_field(umodpoly& r, const umodpoly& a, const umodpoly& b)
+{
+       typedef cln::cl_MI field_t;
+
+       if (degree(a) < degree(b)) {
+               r = a;
+               return false;
+       }
+       // The coefficient ring is a field, so any 0 degree polynomial
+       // divides any other polynomial.
+       if (degree(b) == 0) {
+               r.clear();
+               return true;
+       }
+
+       r = a;
+       const field_t b_lcoeff = lcoeff(b);
+       for (std::size_t k = a.size(); k-- >= b.size(); ) {
+
+               // r -= r_k/b_n x^{k - n} b(x)
+               if (zerop(r[k]))
+                       continue;
+
+               field_t qk = div(r[k], b_lcoeff);
+               bug_on(zerop(qk), "division in a field yield zero: "
+                                  << r[k] << '/' << b_lcoeff);
+
+               // Why C++ is so off-by-one prone?
+               for (std::size_t j = k, i = b.size(); i-- != 0; --j) {
+                       if (zerop(b[i]))
+                               continue;
+                       r[j] = r[j] - qk*b[i];
+               }
+               bug_on(!zerop(r[k]), "polynomial division in field failed: " <<
+                                     "r[" << k << "] = " << r[k] << ", " <<
+                                     "r = " << r << ", b = " << b);
+
+       }
+
+       // Canonicalize the remainder: remove leading zeros. Give a hint
+       // to canonicalize(): we know degree(remainder) < degree(b) 
+       // (because the coefficient ring is a field), so 
+       // c_{degree(b)} \ldots c_{degree(a)} are definitely zero.
+       std::size_t from = degree(b) - 1;
+       canonicalize(r, from);
+       return r.empty();
+}
+
+/**
+ * @brief Polynomial remainder for univariate polynomials over a ring. 
+ *
+ * Given two univariate polynomials \f$a, b \in R[x]\f$, where R is
+ * a ring (presumably Z) computes the remainder @a r, which is
+ * defined as \f$a = b q + r\f$. Returns true if the remainder is zero
+ * and false otherwise.
+ */
+template<typename T>
+bool remainder_in_ring(T& r, const T& a, const T& b)
+{
+       typedef typename T::value_type ring_t;
+       if (degree(a) < degree(b)) {
+               r = a;
+               return false;
+       }
+       // N.B: don't bother to optimize division by constant
+
+       r = a;
+       const ring_t b_lcoeff = lcoeff(b);
+       for (std::size_t k = a.size(); k-- >= b.size(); ) {
+
+               // r -= r_k/b_n x^{k - n} b(x)
+               if (zerop(r[k]))
+                       continue;
+
+               const ring_t qk = div(r[k], b_lcoeff);
+
+               // Why C++ is so off-by-one prone?
+               for (std::size_t j = k, i = b.size(); i-- != 0; --j) {
+                       if (zerop(b[i]))
+                               continue;
+                       r[j] = r[j] - qk*b[i];
+               }
+
+               if (!zerop(r[k])) {
+                       // division failed, don't bother to continue
+                       break;
+               }
+       }
+
+       // Canonicalize the remainder: remove leading zeros. We can't say
+       // anything about the degree of the remainder here.
+       canonicalize(r);
+       return r.empty();
+}
+} // namespace GiNaC
+
+#endif // GINAC_UPOLY_REMAINDER_TCC
+
diff --git a/ginac/polynomial/ring_traits.hpp b/ginac/polynomial/ring_traits.hpp
new file mode 100644 (file)
index 0000000..10fb172
--- /dev/null
@@ -0,0 +1,32 @@
+#ifndef GINAC_RING_TRAITS_HPP
+#define GINAC_RING_TRAITS_HPP
+#include <cln/integer.h>
+#include <cln/modinteger.h>
+
+namespace cln
+{
+static inline cln::cl_I div(const cln::cl_I& x, const cln::cl_I& y)
+{
+       return cln::exquo(x, y);
+}
+
+static inline cln::cl_I get_ring_elt(const cln::cl_I& sample, const int val)
+{
+       return cln::cl_I(val);
+}
+
+static inline cln::cl_MI get_ring_elt(const cln::cl_MI& sample, const int val)
+{
+       return sample.ring()->canonhom(val);
+}
+
+template<typename T>
+static inline T the_one(const T& sample)
+{
+       return get_ring_elt(sample, 1);
+}
+
+} // namespace cln
+
+#endif // GINAC_RING_TRAITS_HPP
+
diff --git a/ginac/polynomial/upoly.hpp b/ginac/polynomial/upoly.hpp
new file mode 100644 (file)
index 0000000..3564ae0
--- /dev/null
@@ -0,0 +1,129 @@
+#ifndef GINAC_NEW_UPOLY_HPP
+#define GINAC_NEW_UPOLY_HPP
+#include <vector>
+#include <cstddef>
+#include <cln/integer.h>
+#include <cln/modinteger.h>
+#include <cln/integer_io.h>
+#include <stdexcept>
+#include <iterator>
+#include <limits>
+#include "debug.hpp"
+#include "compiler.h"
+
+namespace GiNaC
+{
+
+typedef std::vector<cln::cl_I> upoly;
+typedef std::vector<cln::cl_MI> umodpoly;
+
+template<typename T> static std::size_t degree(const T& p)
+{
+       return p.size() - 1;
+}
+
+template<typename T> static typename T::value_type lcoeff(const T& p)
+{
+       bug_on(p.empty(), "lcoeff of a zero polynomial is undefined");
+       return p[p.size() - 1];
+}
+
+template<typename T> static typename T::reference lcoeff(T& p)
+{
+       bug_on(p.empty(), "lcoeff of a zero polynomial is undefined");
+       return p[p.size() - 1];
+}
+
+template<typename T> static typename T::value_type max_coeff(const T& p)
+{
+       bug_on(p.empty(), "max_coeff of a zero polynomial is undefined");
+       typename T::value_type curr = p[p.size() - 1];
+       for (std::size_t i = p.size(); i-- != 0; ) {
+               if (p[i] > curr)
+                       curr = p[i];
+       }
+       return curr;
+}
+
+
+
+// Canonicalize the polynomial, i.e. remove leading zero coefficients
+template<typename T> static void
+canonicalize(T& p, const typename T::size_type hint =
+                  std::numeric_limits<typename T::size_type>::max())
+{
+       if (p.empty())
+               return;
+
+       std::size_t i = p.size() - 1;
+       // Be fast if the polynomial is already canonicalized
+       if (!zerop(p[i]))
+               return;
+
+       if (hint < p.size())
+               i = hint;
+
+       bool is_zero = false;
+       do {
+               if (!zerop(p[i])) {
+                       ++i;
+                       break;
+               }
+               if (i == 0) {
+                       is_zero = true;
+                       break;
+               }
+               --i;
+       } while (true);
+
+       if (is_zero) {
+               p.clear();
+               return;
+       }
+
+       bug_on(!zerop(p.at(i)), "p[" << i << "] = " << p[i] << " != 0 would be erased.");
+
+       typename T::const_iterator it = p.begin() + i;
+       for (std::size_t k = i; it != p.end(); ++it, ++k) {
+               bug_on(!zerop(*it), "p[" << k << "] = " << p[k] << " != 0 "
+                                  "would be erased.");
+       }
+
+       p.erase(p.begin() + i, p.end());
+
+       bug_on(!p.empty() && zerop(lcoeff(p)), "oops, lcoeff(p) = 0");
+}
+
+// Multiply a univariate polynomial p \in R[x] with a constant c \in R
+template<typename T> static T&
+operator*=(T& p, const typename T::value_type& c)
+{
+       if (p.empty())
+               return p;
+       if (zerop(c)) {
+               p.clear();
+               return p;
+       }
+       if (c == the_one(p[0]))
+               return p;
+
+       for (std::size_t i = p.size(); i-- != 0; )
+               p[i] = p[i]*c;
+       canonicalize(p);
+       return p;
+}
+
+// Convert Z[x] -> Z/p[x]
+
+static void
+make_umodpoly(umodpoly& up, const upoly& p, const cln::cl_modint_ring& R)
+{
+       for (std::size_t i = p.size(); i-- != 0; )
+               up[i] = R->canonhom(p[i]);
+       canonicalize(up);
+}
+
+} // namespace GiNaC
+
+#endif // GINAC_NEW_UPOLY_HPP
+
diff --git a/ginac/polynomial/upoly_io.cpp b/ginac/polynomial/upoly_io.cpp
new file mode 100644 (file)
index 0000000..35b623a
--- /dev/null
@@ -0,0 +1,53 @@
+#include <iostream>
+#include <string>
+#include "upoly.hpp"
+#include "upoly_io.hpp"
+#include <cln/integer_io.h>
+
+namespace GiNaC
+{
+using std::ostream;
+using std::string;
+
+template<typename T> static void 
+print(const T& p, ostream& os, const string& varname = string("x"))
+{
+       if (p.size() == 0)
+               os << '0';
+
+       bool seen_nonzero = false;
+
+       for (std::size_t i = p.size(); i-- != 0;  ) {
+               if (zerop(p[i])) {
+                       if (seen_nonzero)
+                               continue;
+                       os << "+ [WARNING: 0]*" << varname << "^" << i << "]";
+                       continue;
+               }
+               seen_nonzero = true;
+               os << "+ (" << p[i] << ")";
+               if (i != 0)
+                       os << "*" << varname;
+               if (i > 1)
+                       os << '^' << i;
+               os << " ";
+       }
+}
+
+#define DEFINE_OPERATOR_OUT(type)                                      \
+std::ostream& operator<<(std::ostream& os, const type& p)              \
+{                                                                      \
+       print(p, os);                                                   \
+       return os;                                                      \
+}                                                                      \
+void dbgprint(const type& p)                                           \
+{                                                                      \
+       print(p, std::cerr);                                            \
+}
+
+DEFINE_OPERATOR_OUT(upoly);
+DEFINE_OPERATOR_OUT(umodpoly);
+#undef DEFINE_OPERATOR_OUT
+
+} // namespace GiNaC
+
diff --git a/ginac/polynomial/upoly_io.hpp b/ginac/polynomial/upoly_io.hpp
new file mode 100644 (file)
index 0000000..a9fb259
--- /dev/null
@@ -0,0 +1,12 @@
+#ifndef GINAC_UPOLY_IO_HPP
+#define GINAC_UPOLY_IO_HPP
+#include <iostream>
+#include "upoly.hpp"
+
+namespace GiNaC
+{
+extern std::ostream& operator<<(std::ostream&, const upoly& );
+extern std::ostream& operator<<(std::ostream&, const umodpoly& );
+}
+
+#endif // GINAC_UPOLY_IO_HPP