]> www.ginac.de Git - ginac.git/blob - ginac/polynomial/primpart_content.cpp
1 #include "ex.h"
2 #include "numeric.h"
3 #include "collect_vargs.h"
4 #include "euclid_gcd_wrap.h"
5 #include "divide_in_z_p.h"
6 #include "debug.h"
8 namespace GiNaC
9 {
11 /**
12  * Compute the primitive part and the content of a modular multivariate
13  * polynomial e \in Z_p[x_n][x_0, \ldots, x_{n-1}], i.e. e is considered
14  * a polynomial in variables x_0, \ldots, x_{n-1} with coefficients being
15  * modular polynomials Z_p[x_n]
16  * @param e polynomial to operate on
17  * @param pp primitive part of @a e, will be computed by this function
18  * @param c content (in the sense described above) of @a e, will be
19  *        computed by this function
20  * @param vars variables x_0, \ldots, x_{n-1}, x_n
21  * @param p modulus
22  */
23 void primpart_content(ex& pp, ex& c, ex e, const exvector& vars,
24                       const long p)
25 {
26         static const ex ex1(1);
27         static const ex ex0(0);
28         e = e.expand();
29         if (e.is_zero()) {
30                 pp = ex0;
31                 c = ex1;
32                 return;
33         }
34         exvector rest_vars = vars;
35         rest_vars.pop_back();
36         ex_collect_t ec;
37         collect_vargs(ec, e, rest_vars);
39         if (ec.size() == 1) {
40                 // the input polynomial factorizes into
41                 // p_1(x_n) p_2(x_0, \ldots, x_{n-1})
42                 c = ec.rbegin()->second;
43                 ec.rbegin()->second = ex1;
44                 pp = ex_collect_to_ex(ec, vars).expand().smod(numeric(p));
45                 return;
46         }
48         // Start from the leading coefficient (which is stored as a last
49         // element of the terms array)
50         ex_collect_t::reverse_iterator i = ec.rbegin();
51         ex g = i->second;
52         // there are at least two terms, so it's safe to...
53         ++i;
54         while (i != ec.rend() && !g.is_equal(ex1)) {
55                 g = euclid_gcd(i->second, g, vars.back(), p);
56                 ++i;
57         }
58         if (g.is_equal(ex1)) {
59                 pp = e;
60                 c = ex1;
61                 return;
62         }
63         exvector mainvar;
64         mainvar.push_back(vars.back());
65         for (i = ec.rbegin(); i != ec.rend(); ++i) {
66                 ex tmp(0);
67                 if (!divide_in_z_p(i->second, g, tmp, mainvar, p))
68                         throw std::logic_error(std::string(__func__) +
69                                         ": bogus division failure");
70                 i->second = tmp;
71         }
73         pp = ex_collect_to_ex(ec, rest_vars).expand().smod(numeric(p));
74         c = g;
75 }
77 } // namespace GiNaC