From 546bababce8ef2b3c5ed3011cb7e49bd78016844 Mon Sep 17 00:00:00 2001 From: Alexei Sheplyakov Date: Sat, 20 Sep 2008 23:18:46 +0400 Subject: [PATCH 1/1] Implemented modular GCD algorithm for univariate polynomials. --- check/Makefile.am | 9 +- check/exam_cra.cpp | 120 ++++++++++++++++++++++ check/exam_mod_gcd.cpp | 85 ++++++++++++++++ ginac/Makefile.am | 16 ++- ginac/polynomial/cra_garner.cpp | 88 +++++++++++++++++ ginac/polynomial/cra_garner.hpp | 12 +++ ginac/polynomial/debug.hpp | 29 ++++++ ginac/polynomial/gcd_euclid.tcc | 45 +++++++++ ginac/polynomial/mod_gcd.cpp | 165 +++++++++++++++++++++++++++++++ ginac/polynomial/mod_gcd.hpp | 11 +++ ginac/polynomial/normalize.tcc | 93 +++++++++++++++++ ginac/polynomial/remainder.tcc | 116 ++++++++++++++++++++++ ginac/polynomial/ring_traits.hpp | 32 ++++++ ginac/polynomial/upoly.hpp | 129 ++++++++++++++++++++++++ ginac/polynomial/upoly_io.cpp | 53 ++++++++++ ginac/polynomial/upoly_io.hpp | 12 +++ 16 files changed, 1012 insertions(+), 3 deletions(-) create mode 100644 check/exam_cra.cpp create mode 100644 check/exam_mod_gcd.cpp create mode 100644 ginac/polynomial/cra_garner.cpp create mode 100644 ginac/polynomial/cra_garner.hpp create mode 100644 ginac/polynomial/debug.hpp create mode 100644 ginac/polynomial/gcd_euclid.tcc create mode 100644 ginac/polynomial/mod_gcd.cpp create mode 100644 ginac/polynomial/mod_gcd.hpp create mode 100644 ginac/polynomial/normalize.tcc create mode 100644 ginac/polynomial/remainder.tcc create mode 100644 ginac/polynomial/ring_traits.hpp create mode 100644 ginac/polynomial/upoly.hpp create mode 100644 ginac/polynomial/upoly_io.cpp create mode 100644 ginac/polynomial/upoly_io.hpp diff --git a/check/Makefile.am b/check/Makefile.am index 87b489d1..5c24aee6 100644 --- a/check/Makefile.am +++ b/check/Makefile.am @@ -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 index 00000000..11e1ab70 --- /dev/null +++ b/check/exam_cra.cpp @@ -0,0 +1,120 @@ +#include +#include +#include +#include +#include +#include "polynomial/cra_garner.hpp" +#include +#include +#include +using namespace cln; +using namespace std; + +/// Generate a sequences of primes p_i such that \prod_i p_i < limit +static std::vector +make_random_moduli(const cln::cl_I& limit); + +static std::vector +calc_residues(const cln::cl_I& x, const std::vector& moduli); + +static void dump(const std::vector& 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 moduli = make_random_moduli(lim); + cln::cl_I x = random_I(lim) + 1; + + if (x > (lim >> 1)) + x = x - lim; + + std::vector 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 map_t; + map_t the_map; + // Run 1024 tests with native 32-bit numbers + the_map[cln::cl_I(std::numeric_limits::max())] = 1024; + + // Run 512 tests with native 64-bit integers + if (sizeof(long) > sizeof(int)) + the_map[cln::cl_I(std::numeric_limits::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 +calc_residues(const cln::cl_I& x, const std::vector& moduli) +{ + std::vector residues(moduli.size()); + for (std::size_t i = moduli.size(); i-- != 0; ) + residues[i] = mod(x, moduli[i]); + return residues; +} + +static std::vector +make_random_moduli(const cln::cl_I& limit) +{ + std::vector 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& 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 index 00000000..4407d90a --- /dev/null +++ b/check/exam_mod_gcd.cpp @@ -0,0 +1,85 @@ +#include "polynomial/upoly.hpp" +#include "polynomial/upoly_io.hpp" +#include "polynomial/mod_gcd.hpp" +#include "ginac.h" +#include +#include +#include +#include +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 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::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(ex_to(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; +} + diff --git a/ginac/Makefile.am b/ginac/Makefile.am index db09d388..b683d9b4 100644 --- a/ginac/Makefile.am +++ b/ginac/Makefile.am @@ -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 index 00000000..b400adbe --- /dev/null +++ b/ginac/polynomial/cra_garner.cpp @@ -0,0 +1,88 @@ +#include +#include +#include +#include +#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& dst, + const vector& 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& dst, + const vector& residues, + const vector& moduli, + const vector& 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& mixed_radix_coeffs, + const vector& 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& residues, + const vector& moduli) +{ + + vector recips(moduli.size() - 1); + compute_recips(recips, moduli); + + vector 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 index 00000000..eff8e9ba --- /dev/null +++ b/ginac/polynomial/cra_garner.hpp @@ -0,0 +1,12 @@ +#ifndef CL_INTEGER_CRA +#define CL_INTEGER_CRA +#include +#include + +namespace cln +{ +extern cl_I integer_cra(const std::vector& residues, + const std::vector& moduli); +} +#endif // CL_INTEGER_CRA + diff --git a/ginac/polynomial/debug.hpp b/ginac/polynomial/debug.hpp new file mode 100644 index 00000000..901ef873 --- /dev/null +++ b/ginac/polynomial/debug.hpp @@ -0,0 +1,29 @@ +#ifndef GINAC_MOD_GCD_DEBUG_HPP +#define GINAC_MOD_GCD_DEBUG_HPP +#include +#include +#include +#include + +#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 index 00000000..deae2522 --- /dev/null +++ b/ginac/polynomial/gcd_euclid.tcc @@ -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 index 00000000..8aa84808 --- /dev/null +++ b/ginac/polynomial/mod_gcd.cpp @@ -0,0 +1,165 @@ +#include "upoly.hpp" +#include "gcd_euclid.tcc" +#include "cra_garner.hpp" +#include +#include +#include +#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 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 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 index 00000000..43e1e929 --- /dev/null +++ b/ginac/polynomial/mod_gcd.hpp @@ -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 index 00000000..e2463075 --- /dev/null +++ b/ginac/polynomial/normalize.tcc @@ -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 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 index 00000000..4db2843b --- /dev/null +++ b/ginac/polynomial/remainder.tcc @@ -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 +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 index 00000000..10fb172a --- /dev/null +++ b/ginac/polynomial/ring_traits.hpp @@ -0,0 +1,32 @@ +#ifndef GINAC_RING_TRAITS_HPP +#define GINAC_RING_TRAITS_HPP +#include +#include + +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 +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 index 00000000..3564ae02 --- /dev/null +++ b/ginac/polynomial/upoly.hpp @@ -0,0 +1,129 @@ +#ifndef GINAC_NEW_UPOLY_HPP +#define GINAC_NEW_UPOLY_HPP +#include +#include +#include +#include +#include +#include +#include +#include +#include "debug.hpp" +#include "compiler.h" + +namespace GiNaC +{ + +typedef std::vector upoly; +typedef std::vector umodpoly; + +template static std::size_t degree(const T& p) +{ + return p.size() - 1; +} + +template 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 static typename T::reference lcoeff(T& p) +{ + bug_on(p.empty(), "lcoeff of a zero polynomial is undefined"); + return p[p.size() - 1]; +} + +template 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 static void +canonicalize(T& p, const typename T::size_type hint = + std::numeric_limits::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 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 index 00000000..35b623a2 --- /dev/null +++ b/ginac/polynomial/upoly_io.cpp @@ -0,0 +1,53 @@ +#include +#include +#include "upoly.hpp" +#include "upoly_io.hpp" +#include + +namespace GiNaC +{ +using std::ostream; +using std::string; + +template 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 index 00000000..a9fb2591 --- /dev/null +++ b/ginac/polynomial/upoly_io.hpp @@ -0,0 +1,12 @@ +#ifndef GINAC_UPOLY_IO_HPP +#define GINAC_UPOLY_IO_HPP +#include +#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 -- 2.44.0