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-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  */
  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