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 "debug.hpp"
  12 #include "compiler.h"
  13
  14 namespace GiNaC
  15 {
  16
  17 typedef std::vector<cln::cl_I> upoly;
  18 typedef std::vector<cln::cl_MI> umodpoly;
  19
  20 template<typename T> static std::size_t degree(const T& p)
  21 {
  22         return p.size() - 1;
  23 }
  24
  25 template<typename T> static typename T::value_type lcoeff(const T& p)
  26 {
  27         bug_on(p.empty(), "lcoeff of a zero polynomial is undefined");
  28         return p[p.size() - 1];
  29 }
  30
  31 template<typename T> static typename T::reference lcoeff(T& p)
  32 {
  33         bug_on(p.empty(), "lcoeff of a zero polynomial is undefined");
  34         return p[p.size() - 1];
  35 }
  36
  37 template<typename T> static typename T::value_type max_coeff(const T& p)
  38 {
  39         bug_on(p.empty(), "max_coeff of a zero polynomial is undefined");
  40         typename T::value_type curr = p[p.size() - 1];
  41         for (std::size_t i = p.size(); i-- != 0; ) {
  42                 if (p[i] > curr)
  43                         curr = p[i];
  44         }
  45         return curr;
  46 }
  47
  48
  49
  50 // Canonicalize the polynomial, i.e. remove leading zero coefficients
  51 template<typename T> static void
  52 canonicalize(T& p, const typename T::size_type hint =
  53                    std::numeric_limits<typename T::size_type>::max())
  54 {
  55         if (p.empty())
  56                 return;
  57
  58         std::size_t i = p.size() - 1;
  59         // Be fast if the polynomial is already canonicalized
  60         if (!zerop(p[i]))
  61                 return;
  62
  63         if (hint < p.size())
  64                 i = hint;
  65
  66         bool is_zero = false;
  67         do {
  68                 if (!zerop(p[i])) {
  69                         ++i;
  70                         break;
  71                 }
  72                 if (i == 0) {
  73                         is_zero = true;
  74                         break;
  75                 }
  76                 --i;
  77         } while (true);
  78
  79         if (is_zero) {
  80                 p.clear();
  81                 return;
  82         }
  83
  84         bug_on(!zerop(p.at(i)), "p[" << i << "] = " << p[i] << " != 0 would be erased.");
  85
  86         typename T::const_iterator it = p.begin() + i;
  87         for (std::size_t k = i; it != p.end(); ++it, ++k) {
  88                 bug_on(!zerop(*it), "p[" << k << "] = " << p[k] << " != 0 "
  89                                    "would be erased.");
  90         }
  91
  92         p.erase(p.begin() + i, p.end());
  93
  94         bug_on(!p.empty() && zerop(lcoeff(p)), "oops, lcoeff(p) = 0");
  95 }
  96
  97 // Multiply a univariate polynomial p \in R[x] with a constant c \in R
  98 template<typename T> static T&
  99 operator*=(T& p, const typename T::value_type& c)
 100 {
 101         if (p.empty())
 102                 return p;
 103         if (zerop(c)) {
 104                 p.clear();
 105                 return p;
 106         }
 107         if (c == the_one(p[0]))
 108                 return p;
 109
 110         for (std::size_t i = p.size(); i-- != 0; )
 111                 p[i] = p[i]*c;
 112         canonicalize(p);
 113         return p;
 114 }
 115
 116 // Convert Z[x] -> Z/p[x]
 117
 118 static void
 119 make_umodpoly(umodpoly& up, const upoly& p, const cln::cl_modint_ring& R)
 120 {
 121         for (std::size_t i = p.size(); i-- != 0; )
 122                 up[i] = R->canonhom(p[i]);
 123         canonicalize(up);
 124 }
 125
 126 } // namespace GiNaC
 127
 128 #endif // GINAC_NEW_UPOLY_HPP
 129