1 /** @file primpart_content.cpp
2 *
3 * Function to find primitive part of a multivariate polynomial. */
4
5 /*
6 * GiNaC Copyright (C) 1999-2022 Johannes Gutenberg University Mainz, Germany
7 *
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.
12 *
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.
17 *
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
21 */
22
23 #include "ex.h"
24 #include "numeric.h"
25 #include "collect_vargs.h"
26 #include "euclid_gcd_wrap.h"
27 #include "divide_in_z_p.h"
28 #include "debug.h"
29
30 namespace GiNaC {
31
32 /**
33 * Compute the primitive part and the content of a modular multivariate
34 * polynomial e \in Z_p[x_n][x_0, \ldots, x_{n-1}], i.e. e is considered
35 * a polynomial in variables x_0, \ldots, x_{n-1} with coefficients being
36 * modular polynomials Z_p[x_n]
37 * @param e polynomial to operate on
38 * @param pp primitive part of @a e, will be computed by this function
39 * @param c content (in the sense described above) of @a e, will be
40 * computed by this function
41 * @param vars variables x_0, \ldots, x_{n-1}, x_n
42 * @param p modulus
43 */
primpart_content(ex & pp,ex & c,ex e,const exvector & vars,const long p)44 void primpart_content(ex& pp, ex& c, ex e, const exvector& vars,
45 const long p)
46 {
47 static const ex ex1(1);
48 static const ex ex0(0);
49 e = e.expand();
50 if (e.is_zero()) {
51 pp = ex0;
52 c = ex1;
53 return;
54 }
55 exvector rest_vars = vars;
56 rest_vars.pop_back();
57 ex_collect_t ec;
58 collect_vargs(ec, e, rest_vars);
59
60 if (ec.size() == 1) {
61 // the input polynomial factorizes into
62 // p_1(x_n) p_2(x_0, \ldots, x_{n-1})
63 c = ec.rbegin()->second;
64 ec.rbegin()->second = ex1;
65 pp = ex_collect_to_ex(ec, rest_vars).expand().smod(numeric(p));
66 return;
67 }
68
69 // Start from the leading coefficient (which is stored as a last
70 // element of the terms array)
71 auto i = ec.rbegin();
72 ex g = i->second;
73 // there are at least two terms, so it's safe to...
74 ++i;
75 while (i != ec.rend() && !g.is_equal(ex1)) {
76 g = euclid_gcd(i->second, g, vars.back(), p);
77 ++i;
78 }
79 if (g.is_equal(ex1)) {
80 pp = e;
81 c = ex1;
82 return;
83 }
84 exvector mainvar;
85 mainvar.push_back(vars.back());
86 for (i = ec.rbegin(); i != ec.rend(); ++i) {
87 ex tmp(0);
88 if (!divide_in_z_p(i->second, g, tmp, mainvar, p))
89 throw std::logic_error(std::string(__func__) +
90 ": bogus division failure");
91 i->second = tmp;
92 }
93
94 pp = ex_collect_to_ex(ec, rest_vars).expand().smod(numeric(p));
95 c = g;
96 }
97
98 } // namespace GiNaC
99