2db82588f3a744315a6393f3428cb7ccfbb27678
[ginac.git] / ginac / polynomial / primpart_content.cpp
1 /** @file primpart_content.cpp
2  *
3  *  Function to find primitive part of a multivariate polynomial. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2019 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  */
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