ginac/polynomial/divide_in_z_p.cpp

   1 /** @file divide_in_z_p.cpp
   2  *
   3  *  Implementation of polynomial division in Z/Zp. */
   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 "add.h"
  24 #include "operators.h"
  25 #include "power.h"
  26 #include "smod_helpers.h"
  27
  28 namespace GiNaC {
  29
  30 /**
  31  * Exact polynomial division of a, b \in Z_p[x_0, \ldots, x_n]
  32  * It doesn't check whether the inputs are proper polynomials, so be careful
  33  * of what you pass in.
  34  *
  35  * @param a  first multivariate polynomial (dividend)
  36  * @param b  second multivariate polynomial (divisor)
  37  * @param q  quotient (returned)
  38  * @param var variables X iterator to first element of vector of symbols
  39  *
  40  * @return "true" when exact division succeeds (the quotient is returned in
  41  *          q), "false" otherwise.
  42  * @note @a p = 0 means the base ring is Z
  43  */
  44 bool divide_in_z_p(const ex &a, const ex &b, ex &q, const exvector& vars, const long p)
  45 {
  46         static const ex _ex1(1);
  47         if (b.is_zero())
  48                 throw(std::overflow_error("divide_in_z: division by zero"));
  49         if (b.is_equal(_ex1)) {
  50                 q = a;
  51                 return true;
  52         }
  53         if (is_exactly_a<numeric>(a)) {
  54                 if (is_exactly_a<numeric>(b)) {
  55                         // p == 0 means division in Z
  56                         if (p == 0) {
  57                                 const numeric tmp = ex_to<numeric>(a/b);
  58                                 if (tmp.is_integer()) {
  59                                         q = tmp;
  60                                         return true;
  61                                 } else
  62                                         return false;
  63                         } else {
  64                                 q = (a*recip(ex_to<numeric>(b), p)).smod(numeric(p));
  65                                 return true;
  66                         }
  67                 } else
  68                         return false;
  69         }
  70         if (a.is_equal(b)) {
  71                 q = _ex1;
  72                 return true;
  73         }
  74
  75         // Main symbol
  76         const ex &x = vars.back();
  77
  78         // Compare degrees
  79         int adeg = a.degree(x), bdeg = b.degree(x);
  80         if (bdeg > adeg)
  81                 return false;
  82
  83         // Polynomial long division (recursive)
  84         ex r = a.expand();
  85         if (r.is_zero())
  86                 return true;
  87         int rdeg = adeg;
  88         ex eb = b.expand();
  89         ex blcoeff = eb.coeff(x, bdeg);
  90         exvector v;
  91         v.reserve(std::max(rdeg - bdeg + 1, 0));
  92         exvector rest_vars(vars);
  93         rest_vars.pop_back();
  94         while (rdeg >= bdeg) {
  95                 ex term, rcoeff = r.coeff(x, rdeg);
  96                 if (!divide_in_z_p(rcoeff, blcoeff, term, rest_vars, p))
  97                         break;
  98                 term = (term*power(x, rdeg - bdeg)).expand();
  99                 v.push_back(term);
 100                 r = (r - term*eb).expand();
 101                 if (p != 0)
 102                         r = r.smod(numeric(p));
 103                 if (r.is_zero()) {
 104                         q = dynallocate<add>(v);
 105                         return true;
 106                 }
 107                 rdeg = r.degree(x);
 108         }
 109         return false;
 110 }
 111
 112 } // namespace GiNaC