1 /** @file divide_in_z_p.cpp
3 * Implementation of polynomial division in Z/Zp. */
6 * GiNaC Copyright (C) 1999-2019 Johannes Gutenberg University Mainz, Germany
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 #include "operators.h"
26 #include "smod_helpers.h"
31 * Exact polynomial division of a, b \in Z_p[x_0, \ldots, x_n]
32 * It doesn't check whether the inputs are proper polynomials, so be careful
33 * of what you pass in.
35 * @param a first multivariate polynomial (dividend)
36 * @param b second multivariate polynomial (divisor)
37 * @param q quotient (returned)
38 * @param var variables X iterator to first element of vector of symbols
40 * @return "true" when exact division succeeds (the quotient is returned in
41 * q), "false" otherwise.
42 * @note @a p = 0 means the base ring is Z
44 bool divide_in_z_p(const ex &a, const ex &b, ex &q, const exvector& vars, const long p)
46 static const ex _ex1(1);
48 throw(std::overflow_error("divide_in_z: division by zero"));
49 if (b.is_equal(_ex1)) {
53 if (is_exactly_a<numeric>(a)) {
54 if (is_exactly_a<numeric>(b)) {
55 // p == 0 means division in Z
57 const numeric tmp = ex_to<numeric>(a/b);
58 if (tmp.is_integer()) {
64 q = (a*recip(ex_to<numeric>(b), p)).smod(numeric(p));
76 const ex &x = vars.back();
79 int adeg = a.degree(x), bdeg = b.degree(x);
83 // Polynomial long division (recursive)
89 ex blcoeff = eb.coeff(x, bdeg);
91 v.reserve(std::max(rdeg - bdeg + 1, 0));
92 exvector rest_vars(vars);
94 while (rdeg >= bdeg) {
95 ex term, rcoeff = r.coeff(x, rdeg);
96 if (!divide_in_z_p(rcoeff, blcoeff, term, rest_vars, p))
98 term = (term*power(x, rdeg - bdeg)).expand();
100 r = (r - term*eb).expand();
102 r = r.smod(numeric(p));
104 q = dynallocate<add>(v);