ginac/polynomial/upoly.hpp

   1 #ifndef GINAC_NEW_UPOLY_HPP
   2 #define GINAC_NEW_UPOLY_HPP
   3 #include <vector>
   4 #include <cstddef>
   5 #include <cln/integer.h>
   6 #include <cln/modinteger.h>
   7 #include <cln/integer_io.h>
   8 #include <stdexcept>
   9 #include <iterator>
  10 #include <limits>
  11 #include "ring_traits.hpp"
  12 #include "debug.hpp"
  13 #include "compiler.h"
  14
  15 namespace GiNaC
  16 {
  17
  18 typedef std::vector<cln::cl_I> upoly;
  19 typedef std::vector<cln::cl_MI> umodpoly;
  20
  21 template<typename T> static std::size_t degree(const T& p)
  22 {
  23         return p.size() - 1;
  24 }
  25
  26 template<typename T> static typename T::value_type lcoeff(const T& p)
  27 {
  28         bug_on(p.empty(), "lcoeff of a zero polynomial is undefined");
  29         return p[p.size() - 1];
  30 }
  31
  32 template<typename T> static typename T::reference lcoeff(T& p)
  33 {
  34         bug_on(p.empty(), "lcoeff of a zero polynomial is undefined");
  35         return p[p.size() - 1];
  36 }
  37
  38 template<typename T> static typename T::value_type max_coeff(const T& p)
  39 {
  40         bug_on(p.empty(), "max_coeff of a zero polynomial is undefined");
  41         typename T::value_type curr = p[p.size() - 1];
  42         for (std::size_t i = p.size(); i-- != 0; ) {
  43                 if (p[i] > curr)
  44                         curr = p[i];
  45         }
  46         return curr;
  47 }
  48
  49
  50
  51 // Canonicalize the polynomial, i.e. remove leading zero coefficients
  52 template<typename T> static void
  53 canonicalize(T& p, const typename T::size_type hint =
  54                    std::numeric_limits<typename T::size_type>::max())
  55 {
  56         if (p.empty())
  57                 return;
  58
  59         std::size_t i = p.size() - 1;
  60         // Be fast if the polynomial is already canonicalized
  61         if (!zerop(p[i]))
  62                 return;
  63
  64         if (hint < p.size())
  65                 i = hint;
  66
  67         bool is_zero = false;
  68         do {
  69                 if (!zerop(p[i])) {
  70                         ++i;
  71                         break;
  72                 }
  73                 if (i == 0) {
  74                         is_zero = true;
  75                         break;
  76                 }
  77                 --i;
  78         } while (true);
  79
  80         if (is_zero) {
  81                 p.clear();
  82                 return;
  83         }
  84
  85         bug_on(!zerop(p.at(i)), "p[" << i << "] = " << p[i] << " != 0 would be erased.");
  86
  87         typename T::const_iterator it = p.begin() + i;
  88         for (std::size_t k = i; it != p.end(); ++it, ++k) {
  89                 bug_on(!zerop(*it), "p[" << k << "] = " << p[k] << " != 0 "
  90                                    "would be erased.");
  91         }
  92
  93         p.erase(p.begin() + i, p.end());
  94
  95         bug_on(!p.empty() && zerop(lcoeff(p)), "oops, lcoeff(p) = 0");
  96 }
  97
  98 // Multiply a univariate polynomial p \in R[x] with a constant c \in R
  99 template<typename T> static T&
 100 operator*=(T& p, const typename T::value_type& c)
 101 {
 102         if (p.empty())
 103                 return p;
 104         if (zerop(c)) {
 105                 p.clear();
 106                 return p;
 107         }
 108         if (c == the_one(p[0]))
 109                 return p;
 110
 111         for (std::size_t i = p.size(); i-- != 0; )
 112                 p[i] = p[i]*c;
 113         canonicalize(p);
 114         return p;
 115 }
 116
 117 /// Divide the polynomial @a p by the ring element @a c, put the result
 118 /// into @a r. If @a p is not divisible by @a c @a r is undefined.
 119 template<typename T> bool divide(T& r, const T& p, const typename T::value_type& c)
 120 {
 121         if (p.empty()) {
 122                 r.clear();
 123                 return true;
 124         }
 125         if (r.size() != p.size())
 126                 r.resize(p.size());
 127         bool divisible = true;
 128         for (std::size_t i = p.size(); i-- != 0; ) {
 129                 divisible = div(r[i], p[i], c);
 130                 if (!divisible)
 131                         break;
 132         }
 133         return divisible;
 134 }
 135
 136 template<typename T> bool divide(T& p, const typename T::value_type& c)
 137 {
 138         if (p.empty())
 139                 return true;
 140         T r(p.size());
 141         bool divisible = divide(r, p, c);
 142         if (divisible)
 143                 swap(p, r);
 144         return divisible;
 145 }
 146
 147 // Convert Z[x] -> Z/p[x]
 148
 149 static void
 150 make_umodpoly(umodpoly& up, const upoly& p, const cln::cl_modint_ring& R)
 151 {
 152         for (std::size_t i = p.size(); i-- != 0; )
 153                 up[i] = R->canonhom(p[i]);
 154         canonicalize(up);
 155 }
 156
 157 } // namespace GiNaC
 158
 159 #endif // GINAC_NEW_UPOLY_HPP
 160