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"
   7
   8 namespace GiNaC
   9 {
  10
  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);
  38
  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         }
  47
  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         }
  72
  73         pp = ex_collect_to_ex(ec, rest_vars).expand().smod(numeric(p));
  74         c = g;
  75 }
  76
  77 } // namespace GiNaC
  78