exam_archive \
exam_structure \
exam_hashmap \
- exam_misc
+ exam_misc \
+ exam_mod_gcd \
+ exam_cra
TIMES = time_dennyfliegner \
time_gammaseries \
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
--- /dev/null
+#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 << "]";
+}
+
--- /dev/null
+#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;
+}
+
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
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
git://www.ginac.de / ginac.git / commitdiff