ginac/polynomial/divide_in_z_p.cpp

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