January 2009 - GiNaC-devel

git updated on server
by Richard B. Kreckel 27 Jan '09

27 Jan '09

Hi! I changed the git installation on www.ginac.de from self-built 1.5.4.3 to the version provided by openSuSE 11.1 RPM packages, which is git 1.6.0.2. As far as I can see, everything should be working at least as well as before. Please contact me, if something behaves in an unexpected way. Cheers -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>

1 0

[PATCH] Improve (fix?) smod: now it really converts into symmetric representation...
by Alexei Sheplyakov 20 Jan '09

20 Jan '09

... instead of clumsy convention inspirited by some proprietary CAS. --- doc/tutorial/ginac.texi | 2 +- ginac/numeric.cpp | 19 +++++++++++-------- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi index 3207a96..fe7a264 100644 --- a/doc/tutorial/ginac.texi +++ b/doc/tutorial/ginac.texi @@ -1479,7 +1479,7 @@ evaluated immediately: @tab modulus in positive representation (in the range @code{[0, abs(b)-1]} with the sign of b, or zero) @cindex @code{mod()} @item @code{smod(a, b)} -@tab modulus in symmetric representation (in the range @code{[-iquo(abs(b)-1, 2), iquo(abs(b), 2)]}) +@tab modulus in symmetric representation (in the range @code{[-iquo(abs(b), 2), iquo(abs(b), 2)]}) @cindex @code{smod()} @item @code{irem(a, b)} @tab integer remainder (has the sign of @math{a}, or is zero) diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index 91721e1..50d873a 100644 --- a/ginac/numeric.cpp +++ b/ginac/numeric.cpp @@ -2337,15 +2337,18 @@ const numeric mod(const numeric &a, const numeric &b) /** Modulus (in symmetric representation). - * Equivalent to Maple's mods. * - * @return a mod b in the range [-iquo(abs(b)-1,2), iquo(abs(b),2)]. */ -const numeric smod(const numeric &a, const numeric &b) -{ - if (a.is_integer() && b.is_integer()) { - const cln::cl_I b2 = cln::ceiling1(cln::the<cln::cl_I>(b.to_cl_N()) >> 1) - 1; - return numeric(cln::mod(cln::the<cln::cl_I>(a.to_cl_N()) + b2, - cln::the<cln::cl_I>(b.to_cl_N())) - b2); + * @return a mod b in the range [-iquo(abs(b),2), iquo(abs(b),2)]. */ +const numeric smod(const numeric &a_, const numeric &b_) +{ + if (a_.is_integer() && b_.is_integer()) { + const cln::cl_I a = cln::the<cln::cl_I>(a_.to_cl_N()); + const cln::cl_I b = cln::the<cln::cl_I>(b_.to_cl_N()); + const cln::cl_I b2 = b >> 1; + const cln::cl_I m = cln::mod(a, b); + const cln::cl_I m_b = m - b; + const cln::cl_I ret = m > b2 ? m_b : m; + return numeric(ret); } else return *_num0_p; } -- 1.5.6.5

2 1

symbolic matrix
by Luis Rivera 19 Jan '09

19 Jan '09

Hello, I have some problem accessing the elements of a symbolic matrix. For example the code below creates a 2x2 symbolic matrix 'A' but I can't access the individual elements. The output I get when I run the small program below is: [luis@urubamba ~/GiNaC]#87 ./AA A(0,0) = [[A00,A01],[A10,A11]] matrix::operator(): index out of range What am I making wrong? Thanks for any hint. Luis #############################33 #include #include #include using namespace std; using namespace GiNaC; int main() { int j,k,n=2; try{ matrix A; A = symbolic_matrix(n,n, "A"); cout << "A(0,0) = " << A(0,0) << endl; cout << "A(0,1) = " << A(0,1) << endl; cout << "A(1,0) = " << A(1,0) << endl; cout << "A(1,1) = " << A(1,1) << endl; } catch (exception &z) { cerr << z.what() << endl; return 1; } return 0; } ######################### Créez votre adresse électronique prenom.nom(a)laposte.net 1 Go d'espace de stockage, anti-spam et anti-virus intégrés.

1 0

[PATCH] Implement modular multivariate GCD (based on chinese remaindering algorithm).
by Alexei Sheplyakov 19 Jan '09

19 Jan '09

Both algorithm and implementation are a bit inefficient: * algorithm does not make use of sparseness of inputs (and most of multivariate polynomials are quite sparse) * GiNaC's expressions are essentially immutable. Thus some simple operations (i.e. multiplying a polynomial by an element of the base ring) are prohibitively expensive. * All numbers (i.e. GiNaC::numeric objects) are heap allocated. Still it's much faster (~5x on bivariate polynomials, ~100x on 3-variate ones) than (subresultant) PRS algorithm, so gcd() uses modular algorithm by default. --- ginac/Makefile.am | 18 ++++ ginac/normal.cpp | 10 ++- ginac/normal.h | 17 +++- ginac/polynomial/chinrem_gcd.cpp | 14 +++ ginac/polynomial/chinrem_gcd.h | 19 ++++ ginac/polynomial/collect_vargs.cpp | 164 ++++++++++++++++++++++++++++++ ginac/polynomial/collect_vargs.h | 34 ++++++ ginac/polynomial/debug.hpp | 1 + ginac/polynomial/divide_in_z_p.cpp | 88 ++++++++++++++++ ginac/polynomial/divide_in_z_p.h | 27 +++++ ginac/polynomial/euclid_gcd_wrap.h | 56 ++++++++++ ginac/polynomial/eval_point_finder.h | 51 +++++++++ ginac/polynomial/mgcd.cpp | 96 +++++++++++++++++ ginac/polynomial/newton_interpolate.h | 36 +++++++ ginac/polynomial/optimal_vars_finder.cpp | 134 ++++++++++++++++++++++++ ginac/polynomial/optimal_vars_finder.h | 20 ++++ ginac/polynomial/pgcd.cpp | 136 +++++++++++++++++++++++++ ginac/polynomial/pgcd.h | 27 +++++ ginac/polynomial/poly_cra.h | 38 +++++++ ginac/polynomial/primes_factory.h | 60 +++++++++++ ginac/polynomial/primpart_content.cpp | 78 ++++++++++++++ ginac/polynomial/smod_helpers.h | 72 +++++++++++++ 22 files changed, 1192 insertions(+), 4 deletions(-) create mode 100644 ginac/polynomial/chinrem_gcd.cpp create mode 100644 ginac/polynomial/chinrem_gcd.h create mode 100644 ginac/polynomial/collect_vargs.cpp create mode 100644 ginac/polynomial/collect_vargs.h create mode 100644 ginac/polynomial/divide_in_z_p.cpp create mode 100644 ginac/polynomial/divide_in_z_p.h create mode 100644 ginac/polynomial/euclid_gcd_wrap.h create mode 100644 ginac/polynomial/eval_point_finder.h create mode 100644 ginac/polynomial/mgcd.cpp create mode 100644 ginac/polynomial/newton_interpolate.h create mode 100644 ginac/polynomial/optimal_vars_finder.cpp create mode 100644 ginac/polynomial/optimal_vars_finder.h create mode 100644 ginac/polynomial/pgcd.cpp create mode 100644 ginac/polynomial/pgcd.h create mode 100644 ginac/polynomial/poly_cra.h create mode 100644 ginac/polynomial/primes_factory.h create mode 100644 ginac/polynomial/primpart_content.cpp create mode 100644 ginac/polynomial/smod_helpers.h diff --git a/ginac/Makefile.am b/ginac/Makefile.am index 7ab16d6..76e32b2 100644 --- a/ginac/Makefile.am +++ b/ginac/Makefile.am @@ -35,6 +35,24 @@ polynomial/interpolate_padic_uvar.h \ polynomial/sr_gcd_uvar.h \ polynomial/heur_gcd_uvar.h \ polynomial/gcd_uvar.cpp \ +polynomial/chinrem_gcd.cpp \ +polynomial/chinrem_gcd.h \ +polynomial/collect_vargs.cpp \ +polynomial/collect_vargs.h \ +polynomial/divide_in_z_p.cpp \ +polynomial/divide_in_z_p.h \ +polynomial/euclid_gcd_wrap.h \ +polynomial/eval_point_finder.h \ +polynomial/mgcd.cpp \ +polynomial/newton_interpolate.h \ +polynomial/optimal_vars_finder.cpp \ +polynomial/optimal_vars_finder.h \ +polynomial/pgcd.cpp \ +polynomial/pgcd.h \ +polynomial/poly_cra.h \ +polynomial/primes_factory.h \ +polynomial/primpart_content.cpp \ +polynomial/smod_helpers.h \ polynomial/debug.hpp libginac_la_LDFLAGS = -version-info $(LT_VERSION_INFO) -release $(LT_RELEASE) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 8f5ba73..54fff47 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -44,6 +44,7 @@ #include "pseries.h" #include "symbol.h" #include "utils.h" +#include "polynomial/chinrem_gcd.h" namespace GiNaC { @@ -1623,8 +1624,15 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio } #endif } + if (options & gcd_options::use_sr_gcd) { + g = sr_gcd(aex, bex, var); + } else { + exvector vars; + for (std::size_t n = sym_stats.size(); n-- != 0; ) + vars.push_back(sym_stats[n].sym); + g = chinrem_gcd(aex, bex, vars); + } - g = sr_gcd(aex, bex, var); if (g.is_equal(_ex1)) { // Keep cofactors factored if possible if (ca) diff --git a/ginac/normal.h b/ginac/normal.h index 5a735e3..a352be8 100644 --- a/ginac/normal.h +++ b/ginac/normal.h @@ -37,8 +37,12 @@ struct gcd_options { enum { /** - * Usually GiNaC tries heuristic GCD algorithm before PRS. - * Some people don't like this, so here's a flag to disable it. + * Usually GiNaC tries heuristic GCD first, because typically + * it's much faster than anything else. Even if heuristic + * algorithm fails, the overhead is negligible w.r.t. cost + * of computing the GCD by some other method. However, some + * people dislike it, so here's a flag which tells GiNaC + * to NOT use the heuristic algorithm. */ no_heur_gcd = 2, /** @@ -48,7 +52,14 @@ struct gcd_options * factored polynomials. DON'T SET THIS unless you *really* * know what are you doing! */ - no_part_factored = 4 + no_part_factored = 4, + /** + * By default GiNaC uses modular GCD algorithm. Typically + * it's much faster than PRS (pseudo remainder sequence) + * algorithm. This flag forces GiNaC to use PRS algorithm + */ + use_sr_gcd = 8 + }; }; diff --git a/ginac/polynomial/chinrem_gcd.cpp b/ginac/polynomial/chinrem_gcd.cpp new file mode 100644 index 0000000..4048ee0 --- /dev/null +++ b/ginac/polynomial/chinrem_gcd.cpp @@ -0,0 +1,14 @@ +#include "chinrem_gcd.h" +#include "optimal_vars_finder.h" + +namespace GiNaC +{ + +ex chinrem_gcd(const ex& A, const ex& B) +{ + const exvector vars = gcd_optimal_variables_order(A, B); + ex g = chinrem_gcd(A, B, vars); + return g; +} +} // namespace GiNaC + diff --git a/ginac/polynomial/chinrem_gcd.h b/ginac/polynomial/chinrem_gcd.h new file mode 100644 index 0000000..6c58a9f --- /dev/null +++ b/ginac/polynomial/chinrem_gcd.h @@ -0,0 +1,19 @@ +#ifndef GINAC_CHINREM_GCD_H +#define GINAC_CHINREM_GCD_H +#include "ex.h" + +namespace GiNaC +{ + +extern ex chinrem_gcd(const ex& A_, const ex& B_, const exvector& vars); +extern ex chinrem_gcd(const ex& A, const ex& B); + +struct chinrem_gcd_failed +{ + virtual ~chinrem_gcd_failed() { } +}; + +} + +#endif /* GINAC_CHINREM_GCD_H */ + diff --git a/ginac/polynomial/collect_vargs.cpp b/ginac/polynomial/collect_vargs.cpp new file mode 100644 index 0000000..5649e4a --- /dev/null +++ b/ginac/polynomial/collect_vargs.cpp @@ -0,0 +1,164 @@ +#include <iterator> +#include <map> +#include <algorithm> +#include <stdexcept> +#include <string> +#include <ginac/ginac.h> +#include "collect_vargs.h" +#include <cln/integer.h> +#include "smod_helpers.h" +#include "debug.hpp" + +namespace GiNaC +{ + +typedef std::map<exp_vector_t, ex> ex_collect_priv_t; + +static void +collect_vargs(ex_collect_priv_t& ec, ex e, const exvector& vars); +static void +collect_term(ex_collect_priv_t& ec, const ex& e, const exvector& vars); +static void wipe_out_zeros(ex_collect_priv_t& ec); + +template<typename T, typename CoeffCMP> +struct compare_terms +{ + const CoeffCMP& coeff_cmp; + explicit compare_terms(const CoeffCMP& coeff_cmp_) : coeff_cmp(coeff_cmp_) + { } + inline bool operator()(const T& t1, const T& t2) const + { + bool exponent_is_less = + std::lexicographical_compare(t1.first.rbegin(), + t1.first.rend(), + t2.first.rbegin(), + t2.first.rend()); + if (exponent_is_less) + return true; + + if ((t1.first == t2.first) && + coeff_cmp(t2.second, t2.second)) + return true; + return false; + } +}; + +template<typename T, typename CoeffCMP> +static struct compare_terms<T, CoeffCMP> +make_compare_terms(const T& dummy, const CoeffCMP& coeff_cmp) +{ + return compare_terms<T, CoeffCMP>(coeff_cmp); +} + +void collect_vargs(ex_collect_t& ec, const ex& e, const exvector& vars) +{ + ex_collect_priv_t ecp; + collect_vargs(ecp, e, vars); + ec.reserve(ecp.size()); + std::copy(ecp.begin(), ecp.end(), std::back_inserter(ec)); + std::sort(ec.begin(), ec.end(), + make_compare_terms(*ec.begin(), ex_is_less())); +} + +static void +collect_vargs(ex_collect_priv_t& ec, ex e, const exvector& vars) +{ + e = e.expand(); + if (e.is_zero()) { + ec.clear(); + return; + } + + if (!is_a<add>(e)) { + collect_term(ec, e, vars); + return; + } + + for (const_iterator i = e.begin(); i != e.end(); ++i) + collect_term(ec, *i, vars); + + wipe_out_zeros(ec); +} + +static void +collect_term(ex_collect_priv_t& ec, const ex& e, const exvector& vars) +{ + if (e.is_zero()) + return; + static const ex ex1(1); + exp_vector_t key(vars.size()); + ex pre_coeff = e; + for (std::size_t i = 0; i < vars.size(); ++i) { + const int var_i_pow = pre_coeff.degree(vars[i]); + key[i] = var_i_pow; + pre_coeff = pre_coeff.coeff(vars[i], var_i_pow); + } + ex_collect_priv_t::iterator i = ec.find(key); + if (i != ec.end()) + i->second += pre_coeff; + else + ec.insert(ex_collect_priv_t::value_type(key, pre_coeff)); +} + +static void wipe_out_zeros(ex_collect_priv_t& m) +{ + ex_collect_priv_t::iterator i = m.begin(); + while (i != m.end()) { + // be careful to not invalide iterator, use post-increment + // for that, see e.g. + // http://coding.derkeiler.com/Archive/C_CPP/comp.lang.cpp/2004-02/0502.html + if (i->second.is_zero()) + m.erase(i++); + else + ++i; + } +} + +GiNaC::ex +ex_collect_to_ex(const ex_collect_t& ec, const GiNaC::exvector& vars) +{ + exvector ev; + ev.reserve(ec.size()); + for (std::size_t i = 0; i < ec.size(); ++i) { + exvector tv; + tv.reserve(vars.size() + 1); + for (std::size_t j = 0; j < vars.size(); ++j) { + if (ec[i].first[j] != 0) + tv.push_back(power(vars[j], ec[i].first[j])); + } + tv.push_back(ec[i].second); + ex tmp = (new mul(tv))->setflag(status_flags::dynallocated); + ev.push_back(tmp); + } + ex ret = (new add(ev))->setflag(status_flags::dynallocated); + return ret; +} + +ex lcoeff_wrt(ex e, const exvector& x) +{ + static const ex ex0(0); + e = e.expand(); + if (e.is_zero()) + return ex0; + + ex_collect_t ec; + collect_vargs(ec, e, x); + return ec.rbegin()->second; +} + +cln::cl_I integer_lcoeff(const ex& e, const exvector& vars) +{ + ex_collect_t ec; + collect_vargs(ec, e, vars); + if (ec.size() == 0) + return cln::cl_I(0); + ex lc = ec.rbegin()->second; + bug_on(!is_a<numeric>(lc), "leading coefficient is not an integer"); + bug_on(!lc.info(info_flags::integer), + "leading coefficient is not an integer"); + + return to_cl_I(lc); +} + +} // namespace GiNaC + diff --git a/ginac/polynomial/collect_vargs.h b/ginac/polynomial/collect_vargs.h new file mode 100644 index 0000000..4dd16ba --- /dev/null +++ b/ginac/polynomial/collect_vargs.h @@ -0,0 +1,34 @@ +#ifndef GINAC_COLLECT_VARGS_HPP +#define GINAC_COLLECT_VARGS_HPP +#include "ex.h" +#include <cln/integer.h> +#include <vector> +#include <utility> // std::pair + +namespace GiNaC +{ + +typedef std::vector<int> exp_vector_t; +typedef std::vector<std::pair<exp_vector_t, ex> > ex_collect_t; + +extern void +collect_vargs(ex_collect_t& ec, const ex& e, const exvector& x); +extern ex +ex_collect_to_ex(const ex_collect_t& ec, const exvector& x); + +/** + * Leading coefficient of a multivariate polynomial e, considering it + * as a multivariate polynomial in x_0, \ldots x_{n-1} with coefficients + * being univariate polynomials in R[x_n] (where R is some ring) + */ +extern ex lcoeff_wrt(ex e, const exvector& x); + +/** + * Leading coefficient c \in R (where R = Z or Z_p) of a multivariate + * polynomial e \in R[x_0, \ldots, x_n] + */ +extern cln::cl_I integer_lcoeff(const ex& e, const exvector& vars); + +} // namespace GiNaC + +#endif /* GINAC_COLLECT_VARGS_HPP */ diff --git a/ginac/polynomial/debug.hpp b/ginac/polynomial/debug.hpp index 901ef87..48659fd 100644 --- a/ginac/polynomial/debug.hpp +++ b/ginac/polynomial/debug.hpp @@ -4,6 +4,7 @@ #include <string> #include <stdexcept> #include <sstream> +#include "compiler.h" #define DEBUG_PREFIX __func__ << ':' << __LINE__ << ": " #define EXCEPTION_PREFIX std::string(__func__) + std::string(": ") + diff --git a/ginac/polynomial/divide_in_z_p.cpp b/ginac/polynomial/divide_in_z_p.cpp new file mode 100644 index 0000000..dff30a9 --- /dev/null +++ b/ginac/polynomial/divide_in_z_p.cpp @@ -0,0 +1,88 @@ +#include <ginac/ginac.h> +#include "smod_helpers.h" + +namespace GiNaC +{ +/** + * Exact polynomial division of a, b \in Z_p[x_0, \ldots, x_n] + * It doesn't check whether the inputs are proper polynomials, so be careful + * of what you pass in. + * + * @param a first multivariate polynomial (dividend) + * @param b second multivariate polynomial (divisor) + * @param q quotient (returned) + * @param var variables X iterator to first element of vector of symbols + * + * @return "true" when exact division succeeds (the quotient is returned in + * q), "false" otherwise. + * @note @a p = 0 means the base ring is Z + */ +bool divide_in_z_p(const ex &a, const ex &b, ex &q, const exvector& vars, const long p) +{ + static const ex _ex1(1); + if (b.is_zero()) + throw(std::overflow_error("divide_in_z: division by zero")); + if (b.is_equal(_ex1)) { + q = a; + return true; + } + if (is_exactly_a<numeric>(a)) { + if (is_exactly_a<numeric>(b)) { + // p == 0 means division in Z + if (p == 0) { + const numeric tmp = ex_to<numeric>(a/b); + if (tmp.is_integer()) { + q = tmp; + return true; + } else + return false; + } else { + q = (a*recip(ex_to<numeric>(b), p)).smod(numeric(p)); + return true; + } + } else + return false; + } + if (a.is_equal(b)) { + q = _ex1; + return true; + } + + // Main symbol + const ex &x = vars.back(); + + // Compare degrees + int adeg = a.degree(x), bdeg = b.degree(x); + if (bdeg > adeg) + return false; + + // Polynomial long division (recursive) + ex r = a.expand(); + if (r.is_zero()) + return true; + int rdeg = adeg; + ex eb = b.expand(); + ex blcoeff = eb.coeff(x, bdeg); + exvector v; + v.reserve(std::max(rdeg - bdeg + 1, 0)); + exvector rest_vars(vars); + rest_vars.pop_back(); + while (rdeg >= bdeg) { + ex term, rcoeff = r.coeff(x, rdeg); + if (!divide_in_z_p(rcoeff, blcoeff, term, rest_vars, p)) + break; + term = (term*power(x, rdeg - bdeg)).expand(); + v.push_back(term); + r = (r - term*eb).expand(); + if (p != 0) + r = r.smod(numeric(p)); + if (r.is_zero()) { + q = (new add(v))->setflag(status_flags::dynallocated); + return true; + } + rdeg = r.degree(x); + } + return false; +} + +} diff --git a/ginac/polynomial/divide_in_z_p.h b/ginac/polynomial/divide_in_z_p.h new file mode 100644 index 0000000..f790fcd --- /dev/null +++ b/ginac/polynomial/divide_in_z_p.h @@ -0,0 +1,27 @@ +#ifndef GINAC_CHINREM_GCD_DIVIDE_IN_Z_P_H +#define GINAC_CHINREM_GCD_DIVIDE_IN_Z_P_H + +#include <ginac/ginac.h> +namespace GiNaC +{ + +/** + * Exact polynomial division of a, b \in Z_p[x_0, \ldots, x_n] + * It doesn't check whether the inputs are proper polynomials, so be careful + * of what you pass in. + * + * @param a first multivariate polynomial (dividend) + * @param b second multivariate polynomial (divisor) + * @param q quotient (returned) + * @param var variables X iterator to first element of vector of symbols + * + * @return "true" when exact division succeeds (the quotient is returned in + * q), "false" otherwise. + */ +extern bool +divide_in_z_p(const ex &a, const ex &b, ex &q, const exvector& vars, const long p); + +} // namespace GiNaC + +#endif /* GINAC_CHINREM_GCD_DIVIDE_IN_Z_P_H */ + diff --git a/ginac/polynomial/euclid_gcd_wrap.h b/ginac/polynomial/euclid_gcd_wrap.h new file mode 100644 index 0000000..d77f71c --- /dev/null +++ b/ginac/polynomial/euclid_gcd_wrap.h @@ -0,0 +1,56 @@ +#ifndef GINAC_PGCD_EUCLID_GCD_H +#define GINAC_PGCD_EUCLID_GCD_H +#include "upoly.hpp" +#include "gcd_euclid.tcc" +#include "smod_helpers.h" +#include <ginac/ginac.h> +#include "debug.hpp" + +namespace GiNaC +{ + +static void ex2upoly(umodpoly& u, ex e, const ex& var, const long p) +{ + e = e.expand(); + cln::cl_modint_ring R = cln::find_modint_ring(cln::cl_I(p)); + u.resize(e.degree(var) + 1); + for (int i = 0; i <= e.degree(var); ++i) { + ex ce = e.coeff(var, i); + bug_on(!is_a<numeric>(ce), "i = " << i << ", " << + "coefficient is not a number: " << ce); + const cln::cl_I c = to_cl_I(ce); + u[i] = R->canonhom(c); + } +} + +static ex umodpoly2ex(const umodpoly& a, const ex& var, const long p) +{ + cln::cl_modint_ring R = cln::find_modint_ring(cln::cl_I(p)); + const numeric pnum(p); + exvector ev(a.size()); + for (std::size_t i = a.size(); i-- != 0; ) { + const cln::cl_I c = smod(R->retract(a[i]), p); + const ex term = numeric(c)*power(var, i); + ev.push_back(term); + } + ex ret = (new add(ev))->setflag(status_flags::dynallocated); + return ret; +} + +static ex euclid_gcd(ex A, ex B, const ex& var, const long p) +{ + A = A.expand(); + B = B.expand(); + + umodpoly a, b; + ex2upoly(a, A, var, p); + ex2upoly(b, B, var, p); + umodpoly g; + gcd_euclid(g, a, b); + ex ge = umodpoly2ex(g, var, p); + return ge; +} + +} // namespace GiNaC + +#endif diff --git a/ginac/polynomial/eval_point_finder.h b/ginac/polynomial/eval_point_finder.h new file mode 100644 index 0000000..8ceadd1 --- /dev/null +++ b/ginac/polynomial/eval_point_finder.h @@ -0,0 +1,51 @@ +#ifndef GINAC_PGCD_EVAL_POINT_FINDER_H +#define GINAC_PGCD_EVAL_POINT_FINDER_H +#include <ginac/ginac.h> +#include <set> + +namespace GiNaC +{ + +/// Find a `good' evaluation point b \in Z_p for a pair of multivariate +/// polynomials A, B \in Z_p[x_n][x_0, \ldots, x_n]. Here `good' means +/// that b is not a root of GCD of contents of A and B. N.B. content +/// is univariate polynomial \in Z_p[x_n] +struct eval_point_finder +{ + typedef long value_type; + const value_type p; + std::set<value_type> points; + const random_modint modint_generator; + bool operator()(value_type& b, const ex& g, const ex& x); + eval_point_finder(const value_type& p_) : p(p_), modint_generator(p) + { } +}; + +bool eval_point_finder::operator()(value_type& b, const ex& lc, const ex& x) +{ + random_modint modint_generator(p); + // Search for a new element of field + while (points.size() < p - 1) { + value_type b_ = modint_generator(); + // check if this evaluation point was already used + if (points.find(b_) != points.end()) + continue; + + // mark found value as used, even if it's a root of lc + // (so we don't need to do the check once more) + points.insert(b_); + // Now make sure it's NOT the root of GCD's leading coeffient + if (lc.subs(x == numeric(b_)).smod(numeric(p)).is_zero()) + continue; + // Nice, it's our next evaluation point + b = b_; + return true; + } + // All possible evaluation points were used. + return false; +} + +} + +#endif /* GINAC_PGCD_EVAL_POINT_FINDER_H */ + diff --git a/ginac/polynomial/mgcd.cpp b/ginac/polynomial/mgcd.cpp new file mode 100644 index 0000000..b39475c --- /dev/null +++ b/ginac/polynomial/mgcd.cpp @@ -0,0 +1,96 @@ +#include "chinrem_gcd.h" +#include <cln/integer.h> +#include "pgcd.h" +#include "collect_vargs.h" +#include "primes_factory.h" +#include "divide_in_z_p.h" +#include "poly_cra.h" + +namespace GiNaC +{ + +static cln::cl_I extract_integer_content(ex& Apr, const ex& A) +{ + static const cln::cl_I n1(1); + const numeric icont_ = A.integer_content(); + const cln::cl_I icont = cln::the<cln::cl_I>(icont_.to_cl_N()); + if (icont != 1) { + Apr = (A/icont_).expand(); + return icont; + } else { + Apr = A; + return n1; + } +} + +ex chinrem_gcd(const ex& A_, const ex& B_, const exvector& vars) +{ + ex A, B; + const cln::cl_I a_icont = extract_integer_content(A, A_); + const cln::cl_I b_icont = extract_integer_content(B, B_); + const cln::cl_I c = cln::gcd(a_icont, b_icont); + + const cln::cl_I a_lc = integer_lcoeff(A, vars); + const cln::cl_I b_lc = integer_lcoeff(B, vars); + const cln::cl_I g_lc = cln::gcd(a_lc, b_lc); + + const ex& x(vars.back()); + int n = std::min(A.degree(x), B.degree(x)); + const cln::cl_I A_max_coeff = to_cl_I(A.max_coefficient()); + const cln::cl_I B_max_coeff = to_cl_I(B.max_coefficient()); + const cln::cl_I lcoeff_limit = (cln::cl_I(1) << n)*cln::abs(g_lc)* + std::min(A_max_coeff, B_max_coeff); + + cln::cl_I q = 0; + ex H = 0; + + long p; + primes_factory pfactory; + while (true) { + bool has_primes = pfactory(p, g_lc); + if (!has_primes) + throw chinrem_gcd_failed(); + + const numeric pnum(p); + ex Ap = A.smod(pnum); + ex Bp = B.smod(pnum); + ex Cp = pgcd(Ap, Bp, vars, p); + + const cln::cl_I g_lcp = smod(g_lc, p); + const cln::cl_I Cp_lc = integer_lcoeff(Cp, vars); + const cln::cl_I nlc = smod(recip(Cp_lc, p)*g_lcp, p); + Cp = (Cp*numeric(nlc)).expand().smod(pnum); + int cp_deg = Cp.degree(x); + if (cp_deg == 0) + return numeric(g_lc); + if (zerop(q)) { + H = Cp; + n = cp_deg; + q = p; + } else { + if (cp_deg == n) { + ex H_next = chinese_remainder(H, q, Cp, p); + q = q*cln::cl_I(p); + H = H_next; + } else if (cp_deg < n) { + // all previous homomorphisms are unlucky + q = p; + H = Cp; + n = cp_deg; + } else { + // dp_deg > d_deg: current prime is bad + } + } + if (q < lcoeff_limit) + continue; // don't bother to do division checks + ex C, dummy1, dummy2; + extract_integer_content(C, H); + if (divide_in_z_p(A, C, dummy1, vars, 0) && + divide_in_z_p(B, C, dummy2, vars, 0)) + return (numeric(c)*C).expand(); + // else: try more primes + } +} + +} // namespace GiNaC + diff --git a/ginac/polynomial/newton_interpolate.h b/ginac/polynomial/newton_interpolate.h new file mode 100644 index 0000000..e80687d --- /dev/null +++ b/ginac/polynomial/newton_interpolate.h @@ -0,0 +1,36 @@ +#ifndef GINAC_PGCD_NEWTON_INTERPOLATE_H +#define GINAC_PGCD_NEWTON_INTERPOLATE_H +#include "ex.h" +#include "numeric.h" +#include "smod_helpers.h" +namespace GiNaC +{ + +/** + * Newton interpolation -- incremental form. + * + * Given a polynomials e1 \in Z_p[Y], prev \in Z_p[x][Y], and evaluation + * points pt1, prevpts \in Z_p compute a polynomial P \in Z_[x][Y] such + * that P(x = pt1) = e1, P(x = pt \in prevpts) = prev(x = pt). + * + * @var{prevpts} enumerates previous evaluation points, should have a form + * (x - pt_2) (x - pt_3) \ldots (x - pt_n). + * @var{prev} encodes the result of previous interpolations. + */ +ex newton_interp(const ex& e1, const long pt1, + const ex& prev, const ex& prevpts, + const ex& x, const long p) +{ + const numeric pnum(p); + const numeric nc = ex_to<numeric>(prevpts.subs(x == numeric(pt1)).smod(pnum)); + const numeric nc_1 = recip(nc, p); + // result = prev + prevpts (e1 - prev|_{x = pt_1})/prevpts|_{x = pt_1) + ex tmp = prev.subs(x == numeric(pt1)).smod(pnum); + tmp = (((e1 - tmp).expand().smod(pnum))*nc_1).smod(pnum); + tmp = (prevpts*tmp).expand().smod(pnum); + tmp = (prev + tmp).expand().smod(pnum); + return tmp; +} + +} +#endif /* GINAC_PGCD_NEWTON_INTERPOLATE_H */ diff --git a/ginac/polynomial/optimal_vars_finder.cpp b/ginac/polynomial/optimal_vars_finder.cpp new file mode 100644 index 0000000..4903699 --- /dev/null +++ b/ginac/polynomial/optimal_vars_finder.cpp @@ -0,0 +1,134 @@ +#include <algorithm> +#include <cstddef> +#include "optimal_vars_finder.h" +#include "add.h" +#include "mul.h" +#include "power.h" +#include "symbol.h" +#include "numeric.h" + +namespace GiNaC +{ +namespace +{ +// XXX: copy-pasted from normal.cpp. +/* + * Statistical information about symbols in polynomials + */ + +/** This structure holds information about the highest and lowest degrees + * in which a symbol appears in two multivariate polynomials "a" and "b". + * A vector of these structures with information about all symbols in + * two polynomials can be created with the function get_symbol_stats(). + * + * @see get_symbol_stats */ +struct sym_desc +{ + /** Reference to symbol */ + ex sym; + + /** Highest degree of symbol in polynomial "a" */ + int deg_a; + + /** Highest degree of symbol in polynomial "b" */ + int deg_b; + + /** Lowest degree of symbol in polynomial "a" */ + int ldeg_a; + + /** Lowest degree of symbol in polynomial "b" */ + int ldeg_b; + + /** Maximum of deg_a and deg_b (Used for sorting) */ + int max_deg; + + /** Maximum number of terms of leading coefficient of symbol in both polynomials */ + std::size_t max_lcnops; + + /** Commparison operator for sorting */ + bool operator<(const sym_desc &x) const + { + if (max_deg == x.max_deg) + return max_lcnops < x.max_lcnops; + else + return max_deg < x.max_deg; + } +}; + +// Vector of sym_desc structures +typedef std::vector<sym_desc> sym_desc_vec; + +// Add symbol the sym_desc_vec (used internally by get_symbol_stats()) +static void add_symbol(const ex &s, sym_desc_vec &v) +{ + sym_desc_vec::const_iterator it = v.begin(), itend = v.end(); + while (it != itend) { + if (it->sym.is_equal(s)) // If it's already in there, don't add it a second time + return; + ++it; + } + sym_desc d; + d.sym = s; + v.push_back(d); +} + +// Collect all symbols of an expression (used internally by get_symbol_stats()) +static void collect_symbols(const ex &e, sym_desc_vec &v) +{ + if (is_a<symbol>(e)) { + add_symbol(e, v); + } else if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) { + for (size_t i=0; i<e.nops(); i++) + collect_symbols(e.op(i), v); + } else if (is_exactly_a<power>(e)) { + collect_symbols(e.op(0), v); + } +} + +/** + * @brief Find the order of variables which is optimal for GCD computation. + * + * Collects statistical information about the highest and lowest degrees + * of all variables that appear in two polynomials. Sorts the variables + * by minimum degree (lowest to highest). The information gathered by + * this function is used by GCD routines to find out the main variable + * for GCD computation. + * + * @param a first multivariate polynomial + * @param b second multivariate polynomial + * @param v vector of sym_desc structs (filled in) */ +static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v) +{ + collect_symbols(a, v); + collect_symbols(b, v); + sym_desc_vec::iterator it = v.begin(), itend = v.end(); + while (it != itend) { + int deg_a = a.degree(it->sym); + int deg_b = b.degree(it->sym); + it->deg_a = deg_a; + it->deg_b = deg_b; + it->max_deg = std::max(deg_a, deg_b); + it->max_lcnops = std::max(a.lcoeff(it->sym).nops(), b.lcoeff(it->sym).nops()); + it->ldeg_a = a.ldegree(it->sym); + it->ldeg_b = b.ldegree(it->sym); + ++it; + } + std::sort(v.begin(), v.end()); +} +// XXX: end copy-paste from normal.cpp + +} // anonymous namespace of helper functions + +exvector gcd_optimal_variables_order(const ex& a, const ex& b) +{ + sym_desc_vec v; + get_symbol_stats(a, b, v); + exvector vars; + vars.reserve(v.size()); + for (std::size_t i = v.size(); i-- != 0; ) + vars.push_back(v[i].sym); + return vars; +} + +} // namespace GiNaC + diff --git a/ginac/polynomial/optimal_vars_finder.h b/ginac/polynomial/optimal_vars_finder.h new file mode 100644 index 0000000..9715a47 --- /dev/null +++ b/ginac/polynomial/optimal_vars_finder.h @@ -0,0 +1,20 @@ +#ifndef GINAC_CHINREM_GCD_OPTIMAL_SYMBOL_FINDER_H +#define GINAC_CHINREM_GCD_OPTIMAL_SYMBOL_FINDER_H +#include <vector> +#include "ex.h" + +namespace GiNaC +{ +/** + * @brief Find the order of variables which is optimal for GCD computation. + * + * Collects statistical information about the highest and lowest degrees + * of all variables that appear in two polynomials. Sorts the variables + * by minimum degree (lowest to highest). The information gathered by + * this function is used by GCD routines to find out the main variable + * for GCD computation. + */ +extern exvector gcd_optimal_variables_order(const ex& A, const ex& B); +} +#endif /* GINAC_CHINREM_GCD_OPTIMAL_SYMBOL_FINDER_H */ + diff --git a/ginac/polynomial/pgcd.cpp b/ginac/polynomial/pgcd.cpp new file mode 100644 index 0000000..ddc3b8c --- /dev/null +++ b/ginac/polynomial/pgcd.cpp @@ -0,0 +1,136 @@ +#include "pgcd.h" +#include "collect_vargs.h" +#include "smod_helpers.h" +#include "euclid_gcd_wrap.h" +#include "eval_point_finder.h" +#include "newton_interpolate.h" +#include "divide_in_z_p.h" + +namespace GiNaC +{ + +extern void +primpart_content(ex& pp, ex& c, ex e, const exvector& vars, const long p); + +// Computes the GCD of two polynomials over a prime field. +// Based on Algorithm 7.2 from "Algorithms for Computer Algebra" +// A and B are considered as Z_p[x_n][x_0, \ldots, x_{n-1}], that is, +// as a polynomials in variables x_0, \ldots x_{n-1} having coefficients +// from the ring Z_p[x_n] +ex pgcd(const ex& A, const ex& B, const exvector& vars, const long p) +{ + static const ex ex1(1); + if (A.is_zero()) + return B; + + if (B.is_zero()) + return A; + + if (is_a<numeric>(A)) + return ex1; + + if (is_a<numeric>(B)) + return ex1; + + // Checks for univariate polynomial + if (vars.size() == 1) { + ex ret = euclid_gcd(A, B, vars[0], p); // Univariate GCD + return ret; + } + const ex& mainvar(vars.back()); + + // gcd of the contents + ex H = 0, Hprev = 0; // GCD candidate + ex newton_poly = 1; // for Newton Interpolation + + // Contents and primparts of A and B + ex Aprim, contA; + primpart_content(Aprim, contA, A, vars, p); + ex Bprim, contB; + primpart_content(Bprim, contB, B, vars, p); + // gcd of univariate polynomials + const ex cont_gcd = euclid_gcd(contA, contB, mainvar, p); + + exvector restvars = vars; + restvars.pop_back(); + const ex AL = lcoeff_wrt(Aprim, restvars); + const ex BL = lcoeff_wrt(Bprim, restvars); + // gcd of univariate polynomials + const ex lc_gcd = euclid_gcd(AL, BL, mainvar, p); + + // The estimate of degree of the gcd of Ab and Bb + int gcd_deg = std::min(degree(Aprim, mainvar), degree(Bprim, mainvar)); + eval_point_finder::value_type b; + + eval_point_finder find_eval_point(p); + const numeric pn(p); + do { + // Find a `good' evaluation point b. + bool has_more_pts = find_eval_point(b, lc_gcd, mainvar); + // If there are no more possible evaluation points, bail out + if (!has_more_pts) + throw pgcd_failed(); + + const numeric bn(b); + // Evaluate the polynomials in b + ex Ab = Aprim.subs(mainvar == bn).smod(pn); + ex Bb = Bprim.subs(mainvar == bn).smod(pn); + ex Cb = pgcd(Ab, Bb, restvars, p); + + // Set the correct the leading coefficient + const cln::cl_I lcb_gcd = + smod(to_cl_I(lc_gcd.subs(mainvar == bn)), p); + const cln::cl_I Cblc = integer_lcoeff(Cb, restvars); + const cln::cl_I correct_lc = smod(lcb_gcd*recip(Cblc, p), p); + Cb = (Cb*numeric(correct_lc)).smod(pn); + + // Test for unlucky homomorphisms + const int img_gcd_deg = Cb.degree(restvars.back()); + if (img_gcd_deg < gcd_deg) { + // The degree decreased, previous homomorphisms were + // bad, so we have to start it all over. + H = Cb; + newton_poly = mainvar - numeric(b); + Hprev = 0; + gcd_deg = img_gcd_deg; + continue; + } + if (img_gcd_deg > gcd_deg) { + // The degree of images GCD is too high, this + // evaluation point is bad. Skip it. + continue; + } + + // Image has the same degree as the previous one + // (or at least not higher than the limit) + Hprev = H; + H = newton_interp(Cb, b, H, newton_poly, mainvar, p); + newton_poly = newton_poly*(mainvar - b); + + // try to reduce the number of division tests. + const ex H_lcoeff = lcoeff_wrt(H, restvars); + + if (H_lcoeff.is_equal(lc_gcd)) { + if ((Hprev-H).expand().smod(pn).is_zero()) + continue; + ex C /* primitive part of H */, contH /* dummy */; + primpart_content(C, contH, H, vars, p); + // Normalize GCD so that leading coefficient is 1 + const cln::cl_I Clc = recip(integer_lcoeff(C, vars), p); + C = (C*numeric(Clc)).expand().smod(pn); + + ex dummy1, dummy2; + + if (divide_in_z_p(Aprim, C, dummy1, vars, p) && + divide_in_z_p(Bprim, C, dummy2, vars, p)) + return (cont_gcd*C).expand().smod(pn); + else if (img_gcd_deg == 0) + return cont_gcd; + // else continue building the candidate + } + } while(true); + throw pgcd_failed(); +} + +} // namespace GiNaC + diff --git a/ginac/polynomial/pgcd.h b/ginac/polynomial/pgcd.h new file mode 100644 index 0000000..5e72038 --- /dev/null +++ b/ginac/polynomial/pgcd.h @@ -0,0 +1,27 @@ +#ifndef GINAC_CHINREM_GCD_PGCD_H +#define GINAC_CHINREM_GCD_PGCD_H +#include "ex.h" + +namespace GiNaC +{ + +/// Exception to be thrown when modular GCD algorithm fails +struct pgcd_failed +{ + virtual ~pgcd_failed() { } +}; + +/** + * @brief Compute the GCD of two polynomials over a prime field Z_p + * + * @param vars variables + * @param p designates the coefficient field Z_p + * @param A polynomial \in Z_p[vars] + * @param B second polynomial \in Z_p[vars] + */ +extern ex +pgcd(const ex& A, const ex& B, const exvector& vars, const long p); + +} // namespace GiNaC + +#endif diff --git a/ginac/polynomial/poly_cra.h b/ginac/polynomial/poly_cra.h new file mode 100644 index 0000000..cbcff7b --- /dev/null +++ b/ginac/polynomial/poly_cra.h @@ -0,0 +1,38 @@ +#ifndef GINAC_POLY_CRA_H +#define GINAC_POLY_CRA_H +#include "ex.h" +#include <cln/integer.h> +#include "smod_helpers.h" + +namespace GiNaC +{ + +/** + * @brief Chinese reamainder algorithm for polynomials. + * + * Given two polynomials \f$e_1 \in Z_{q_1}[x_1, \ldots, x_n]\f$ and + * \f$e_2 \in Z_{q_2}[x_1, \ldots, x_n]\f$, compute the polynomial + * \f$r \in Z_{q_1 q_2}[x_1, \ldots, x_n]\f$ such that \f$ r mod q_1 = e_1\f$ + * and \f$ r mod q_2 = e_2 \f$ + */ +ex chinese_remainder(const ex& e1, const cln::cl_I& q1, + const ex& e2, const long q2) +{ + // res = v_1 + v_2 q_1 + // v_1 = e_1 mod q_1 + // v_2 = (e_2 - v_1)/q_1 mod q_2 + const numeric q2n(q2); + const numeric q1n(q1); + ex v1 = e1.smod(q1n); + ex u = v1.smod(q2n); + ex v2 = (e2.smod(q2n) - v1.smod(q2n)).expand().smod(q2n); + const numeric q1_1(recip(q1, q2)); // 1/q_1 mod q_2 + v2 = (v2*q1_1).smod(q2n); + ex ret = (v1 + v2*q1_1).expand(); + return ret; +} + +} // namespace GiNaC + +#endif /* GINAC_POLY_CRA_H */ + diff --git a/ginac/polynomial/primes_factory.h b/ginac/polynomial/primes_factory.h new file mode 100644 index 0000000..093c973 --- /dev/null +++ b/ginac/polynomial/primes_factory.h @@ -0,0 +1,60 @@ +#ifndef GINAC_CHINREM_GCD_PRIMES_FACTORY_H +#define GINAC_CHINREM_GCD_PRIMES_FACTORY_H +#include <cln/integer.h> +#include <cln/numtheory.h> +#include <limits> +#include "smod_helpers.h" +#include "debug.hpp" + +namespace GiNaC +{ + +/** + * Find a `big' prime p such that lc mod p != 0. Helper class used by modular + * GCD algorithm. + */ +class primes_factory +{ +private: + // These primes need to be large enough, so that the number of images + // we need to reconstruct the GCD (in Z) is reasonable. On the other + // hand, they should be as small as possible, so that operations on + // coefficients are efficient. Practically this means we coefficients + // should be native integers. (N.B.: as of now chinrem_gcd uses cl_I + // or even numeric. Eventually this will be fixed). + cln::cl_I last; + // This ensures coefficients are immediate. + static const int immediate_bits = 8*sizeof(void *) - __alignof__(void *); + static const long opt_hint = (1L << (immediate_bits >> 1)) - 1; +public: + primes_factory() + { + last = cln::nextprobprime(cln::cl_I(opt_hint)); + } + + bool operator()(long& p, const cln::cl_I& lc) + { + static const cln::cl_I maxval(std::numeric_limits<long>::max()); + while (last < maxval) { + long p_ = cln::cl_I_to_long(last); + last = cln::nextprobprime(last + 1); + + if (!zerop(smod(lc, p_))) { + p = p_; + return true; + } + } + return false; + } + + bool has_primes() const + { + static const cln::cl_I maxval(std::numeric_limits<long>::max()); + return last < maxval; + } +}; + +} // namespace GiNaC + +#endif /* GINAC_CHINREM_GCD_PRIMES_FACTORY_H */ + diff --git a/ginac/polynomial/primpart_content.cpp b/ginac/polynomial/primpart_content.cpp new file mode 100644 index 0000000..9d36a95 --- /dev/null +++ b/ginac/polynomial/primpart_content.cpp @@ -0,0 +1,78 @@ +#include "ex.h" +#include "numeric.h" +#include "collect_vargs.h" +#include "euclid_gcd_wrap.h" +#include "divide_in_z_p.h" +#include "debug.hpp" + +namespace GiNaC +{ + +/** + * Compute the primitive part and the content of a modular multivariate + * polynomial e \in Z_p[x_n][x_0, \ldots, x_{n-1}], i.e. e is considered + * a polynomial in variables x_0, \ldots, x_{n-1} with coefficients being + * modular polynomials Z_p[x_n] + * @param e polynomial to operate on + * @param pp primitive part of @a e, will be computed by this function + * @param c content (in the sense described above) of @a e, will be + * computed by this function + * @param vars variables x_0, \ldots, x_{n-1}, x_n + * @param p modulus + */ +void primpart_content(ex& pp, ex& c, ex e, const exvector& vars, + const long p) +{ + static const ex ex1(1); + static const ex ex0(0); + e = e.expand(); + if (e.is_zero()) { + pp = ex0; + c = ex1; + return; + } + exvector rest_vars = vars; + rest_vars.pop_back(); + ex_collect_t ec; + collect_vargs(ec, e, rest_vars); + + if (ec.size() == 1) { + // the input polynomial factorizes into + // p_1(x_n) p_2(x_0, \ldots, x_{n-1}) + c = ec.rbegin()->second; + ec.rbegin()->second = ex1; + pp = ex_collect_to_ex(ec, vars).expand().smod(numeric(p)); + return; + } + + // Start from the leading coefficient (which is stored as a last + // element of the terms array) + ex_collect_t::reverse_iterator i = ec.rbegin(); + ex g = i->second; + // there are at least two terms, so it's safe to... + ++i; + while (i != ec.rend() && !g.is_equal(ex1)) { + g = euclid_gcd(i->second, g, vars.back(), p); + ++i; + } + if (g.is_equal(ex1)) { + pp = e; + c = ex1; + return; + } + exvector mainvar; + mainvar.push_back(vars.back()); + for (i = ec.rbegin(); i != ec.rend(); ++i) { + ex tmp(0); + if (!divide_in_z_p(i->second, g, tmp, mainvar, p)) + throw std::logic_error(std::string(__func__) + + ": bogus division failure"); + i->second = tmp; + } + + pp = ex_collect_to_ex(ec, rest_vars).expand().smod(numeric(p)); + c = g; +} + +} // namespace GiNaC + diff --git a/ginac/polynomial/smod_helpers.h b/ginac/polynomial/smod_helpers.h new file mode 100644 index 0000000..81370ba --- /dev/null +++ b/ginac/polynomial/smod_helpers.h @@ -0,0 +1,72 @@ +#ifndef GINAC_POLYNOMIAL_SMOD_HELPERS_H +#define GINAC_POLYNOMIAL_SMOD_HELPERS_H +#include <cln/integer.h> +#include <cln/integer_io.h> +#include "ex.h" +#include "numeric.h" +#include "debug.hpp" + +namespace GiNaC +{ + +/// Z -> Z_p (in the symmetric representation) +static inline cln::cl_I smod(const cln::cl_I& a, long p) +{ + const cln::cl_I p2 = cln::cl_I(p >> 1); + const cln::cl_I m = cln::mod(a, p); + const cln::cl_I m_p = m - cln::cl_I(p); + const cln::cl_I ret = m > p2 ? m_p : m; + return ret; +} + +static inline cln::cl_I recip(const cln::cl_I& a, long p_) +{ + cln::cl_I p(p_); + cln::cl_I u, v; + const cln::cl_I g = xgcd(a, p, &u, &v); + cln::cl_I ret = smod(u, p_); + cln::cl_I chck = smod(a*ret, p_); + bug_on(chck != 1, "miscomputed recip(" << a << " (mod " << p_ << "))"); + return ret; + +} + +static inline numeric recip(const numeric& a_, long p) +{ + const cln::cl_I a = cln::the<cln::cl_I>(a_.to_cl_N()); + const cln::cl_I ret = recip(a, p); + return numeric(ret); +} + +static inline cln::cl_I to_cl_I(const ex& e) +{ + bug_on(!is_a<numeric>(e), "argument should be an integer"); + bug_on(!e.info(info_flags::integer), + "argument should be an integer"); + return cln::the<cln::cl_I>(ex_to<numeric>(e).to_cl_N()); +} + +struct random_modint +{ + typedef long value_type; + const value_type p; + const value_type p_2; + + random_modint(const value_type& p_) : p(p_), p_2((p >> 1)) + { } + value_type operator()() const + { + do { + cln::cl_I tmp_ = cln::random_I(p); + value_type tmp = cln::cl_I_to_long(tmp_); + if (tmp > p_2) + tmp -= p; + return tmp; + } while (true); + } + +}; + +} // namespace GiNaC + +#endif // GINAC_POLYNOMIAL_SMOD_HELPERS_H -- 1.5.6.5

1 0

[PATCH] Fix copy-paste error in GINACLIB_MICRO_VERSION.
by Alexei Sheplyakov 19 Jan '09

19 Jan '09

--- configure.ac | 2 +- 1 files changed, 1 insertions(+), 1 deletions(-) diff --git a/configure.ac b/configure.ac index 2c2c8e9..bf27156 100644 --- a/configure.ac +++ b/configure.ac @@ -38,7 +38,7 @@ AM_INIT_AUTOMAKE([gnu 1.7 dist-bzip2]) dnl Process GiNaC version information AC_DEFINE([GINACLIB_MAJOR_VERSION], ginac_major_version, [Major version of GiNaC]) AC_DEFINE([GINACLIB_MINOR_VERSION], ginac_minor_version, [Minor version of GiNaC]) -AC_DEFINE([GINACLIB_MICRO_VERSION], ginac_minor_version, [Micro version of GiNaC]) +AC_DEFINE([GINACLIB_MICRO_VERSION], ginac_micro_version, [Micro version of GiNaC]) dnl GiNaC archive file version information. dnl -- 1.5.6.5

1 0

[SCM] GiNaC -- a C++ library for symbolic computations branch, master, updated. release_1-4-0-170-g38f3143
by git＠ginac.de 08 Jan '09

08 Jan '09

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "GiNaC -- a C++ library for symbolic computations". The branch, master has been updated via 38f3143eee71d942f40ab47f5a561dc9736203fc (commit) from 4c47ecd9caa39ba1a31f5294e395fcbdf2006431 (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit 38f3143eee71d942f40ab47f5a561dc9736203fc Author: Jens Vollinga <jensv(a)balin.nikhef.nl> Date: Thu Jan 8 10:26:50 2009 +0100 Fixed bug in binomial(). binomial(n,0) erroneously returned 0. ----------------------------------------------------------------------- Summary of changes: check/exam_inifcns.cpp | 15 +++++++++++++++ ginac/inifcns.cpp | 2 +- 2 files changed, 16 insertions(+), 1 deletions(-) hooks/post-receive -- GiNaC -- a C++ library for symbolic computations

1 0