Implemented modular GCD algorithm for univariate polynomials.
[ginac.git] / 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