ginac/polynomial/divide_in_z_p.cpp

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