1 #ifndef GINAC_NEW_UPOLY_HPP
2 #define GINAC_NEW_UPOLY_HPP
5 #include <cln/integer.h>
6 #include <cln/modinteger.h>
7 #include <cln/integer_io.h>
11 #include "ring_traits.hpp"
18 typedef std::vector<cln::cl_I> upoly;
19 typedef std::vector<cln::cl_MI> umodpoly;
21 template<typename T> static std::size_t degree(const T& p)
26 template<typename T> static typename T::value_type lcoeff(const T& p)
28 bug_on(p.empty(), "lcoeff of a zero polynomial is undefined");
29 return p[p.size() - 1];
32 template<typename T> static typename T::reference lcoeff(T& p)
34 bug_on(p.empty(), "lcoeff of a zero polynomial is undefined");
35 return p[p.size() - 1];
38 template<typename T> static typename T::value_type max_coeff(const T& p)
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; ) {
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())
59 std::size_t i = p.size() - 1;
60 // Be fast if the polynomial is already canonicalized
85 bug_on(!zerop(p.at(i)), "p[" << i << "] = " << p[i] << " != 0 would be erased.");
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 "
93 p.erase(p.begin() + i, p.end());
95 bug_on(!p.empty() && zerop(lcoeff(p)), "oops, lcoeff(p) = 0");
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)
108 if (c == the_one(p[0]))
111 for (std::size_t i = p.size(); i-- != 0; )
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)
125 if (r.size() != p.size())
127 bool divisible = true;
128 for (std::size_t i = p.size(); i-- != 0; ) {
129 divisible = div(r[i], p[i], c);
136 template<typename T> bool divide(T& p, const typename T::value_type& c)
141 bool divisible = divide(r, p, c);
147 // Convert Z[x] -> Z/p[x]
150 make_umodpoly(umodpoly& up, const upoly& p, const cln::cl_modint_ring& R)
152 for (std::size_t i = p.size(); i-- != 0; )
153 up[i] = R->canonhom(p[i]);
159 #endif // GINAC_NEW_UPOLY_HPP