1 #include <ginac/ginac.h>
2 #include "smod_helpers.h"
7 * Exact polynomial division of a, b \in Z_p[x_0, \ldots, x_n]
8 * It doesn't check whether the inputs are proper polynomials, so be careful
11 * @param a first multivariate polynomial (dividend)
12 * @param b second multivariate polynomial (divisor)
13 * @param q quotient (returned)
14 * @param var variables X iterator to first element of vector of symbols
16 * @return "true" when exact division succeeds (the quotient is returned in
17 * q), "false" otherwise.
18 * @note @a p = 0 means the base ring is Z
20 bool divide_in_z_p(const ex &a, const ex &b, ex &q, const exvector& vars, const long p)
22 static const ex _ex1(1);
24 throw(std::overflow_error("divide_in_z: division by zero"));
25 if (b.is_equal(_ex1)) {
29 if (is_exactly_a<numeric>(a)) {
30 if (is_exactly_a<numeric>(b)) {
31 // p == 0 means division in Z
33 const numeric tmp = ex_to<numeric>(a/b);
34 if (tmp.is_integer()) {
40 q = (a*recip(ex_to<numeric>(b), p)).smod(numeric(p));
52 const ex &x = vars.back();
55 int adeg = a.degree(x), bdeg = b.degree(x);
59 // Polynomial long division (recursive)
65 ex blcoeff = eb.coeff(x, bdeg);
67 v.reserve(std::max(rdeg - bdeg + 1, 0));
68 exvector rest_vars(vars);
70 while (rdeg >= bdeg) {
71 ex term, rcoeff = r.coeff(x, rdeg);
72 if (!divide_in_z_p(rcoeff, blcoeff, term, rest_vars, p))
74 term = (term*power(x, rdeg - bdeg)).expand();
76 r = (r - term*eb).expand();
78 r = r.smod(numeric(p));
80 q = (new add(v))->setflag(status_flags::dynallocated);