Finalize 1.7.8 release.
[ginac.git] / ginac / polynomial / upoly.h
1 /** @file upoly.h
2  *
3  *  Interface to polynomials with integer and modular coefficients. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2019 Johannes Gutenberg University Mainz, Germany
7  *
8  *  This program is free software; you can redistribute it and/or modify
9  *  it under the terms of the GNU General Public License as published by
10  *  the Free Software Foundation; either version 2 of the License, or
11  *  (at your option) any later version.
12  *
13  *  This program is distributed in the hope that it will be useful,
14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  *  GNU General Public License for more details.
17  *
18  *  You should have received a copy of the GNU General Public License
19  *  along with this program; if not, write to the Free Software
20  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
21  */
22
23 #ifndef GINAC_UPOLY_H
24 #define GINAC_UPOLY_H
25
26 #include "ring_traits.h"
27 #include "debug.h"
28 #include "compiler.h"
29
30 #include <cln/integer.h>
31 #include <cln/integer_io.h>
32 #include <cln/modinteger.h>
33 #include <cstddef>
34 #include <iterator>
35 #include <limits>
36 #include <stdexcept>
37 #include <vector>
38
39 namespace GiNaC {
40
41 typedef std::vector<cln::cl_I> upoly;
42 typedef std::vector<cln::cl_MI> umodpoly;
43
44 template<typename T> static std::size_t degree(const T& p)
45 {
46         return p.size() - 1;
47 }
48
49 template<typename T> static typename T::value_type lcoeff(const T& p)
50 {
51         bug_on(p.empty(), "lcoeff of a zero polynomial is undefined");
52         return p[p.size() - 1];
53 }
54
55 template<typename T> static typename T::reference lcoeff(T& p)
56 {
57         bug_on(p.empty(), "lcoeff of a zero polynomial is undefined");
58         return p[p.size() - 1];
59 }
60
61 template<typename T> static typename T::value_type max_coeff(const T& p)
62 {
63         bug_on(p.empty(), "max_coeff of a zero polynomial is undefined");
64         typename T::value_type curr = p[p.size() - 1];
65         for (std::size_t i = p.size(); i-- != 0; ) {
66                 if (p[i] > curr)
67                         curr = p[i];
68         }
69         return curr;
70 }
71
72
73
74 // Canonicalize the polynomial, i.e. remove leading zero coefficients
75 template<typename T> static void
76 canonicalize(T& p, const typename T::size_type hint =
77                    std::numeric_limits<typename T::size_type>::max())
78 {
79         if (p.empty())
80                 return;
81
82         std::size_t i = p.size() - 1;
83         // Be fast if the polynomial is already canonicalized
84         if (!zerop(p[i]))
85                 return;
86
87         if (hint < p.size())
88                 i = hint;
89
90         bool is_zero = false;
91         do {
92                 if (!zerop(p[i])) {
93                         ++i;
94                         break;
95                 }
96                 if (i == 0) {
97                         is_zero = true;
98                         break;
99                 }
100                 --i;
101         } while (true);
102
103         if (is_zero) {
104                 p.clear();
105                 return;
106         }
107
108         bug_on(!zerop(p.at(i)), "p[" << i << "] = " << p[i] << " != 0 would be erased.");
109
110         typename T::const_iterator it = p.begin() + i;
111         for (std::size_t k = i; it != p.end(); ++it, ++k) {
112                 bug_on(!zerop(*it), "p[" << k << "] = " << p[k] << " != 0 "
113                                    "would be erased.");
114         }
115
116         p.erase(p.begin() + i, p.end());
117
118         bug_on(!p.empty() && zerop(lcoeff(p)), "oops, lcoeff(p) = 0");
119 }
120
121 // Multiply a univariate polynomial p \in R[x] with a constant c \in R
122 template<typename T> static T&
123 operator*=(T& p, const typename T::value_type& c)
124 {
125         if (p.empty())
126                 return p;
127         if (zerop(c)) {
128                 p.clear();
129                 return p;
130         }
131         if (c == the_one(p[0]))
132                 return p;
133
134         for (std::size_t i = p.size(); i-- != 0; )
135                 p[i] = p[i]*c;
136         canonicalize(p);
137         return p;
138 }
139
140 /// Divide the polynomial @a p by the ring element @a c, put the result
141 /// into @a r. If @a p is not divisible by @a c @a r is undefined.
142 template<typename T> bool divide(T& r, const T& p, const typename T::value_type& c)
143 {
144         if (p.empty()) {
145                 r.clear();
146                 return true;
147         }
148         if (r.size() != p.size())
149                 r.resize(p.size());
150         bool divisible = true;
151         for (std::size_t i = p.size(); i-- != 0; ) {
152                 divisible = div(r[i], p[i], c);
153                 if (!divisible)
154                         break;
155         }
156         return divisible;
157 }
158
159 template<typename T> bool divide(T& p, const typename T::value_type& c)
160 {
161         if (p.empty())
162                 return true;
163         T r(p.size());
164         bool divisible = divide(r, p, c);
165         if (divisible)
166                 swap(p, r);
167         return divisible;
168 }
169
170 // Convert Z[x] -> Z/p[x]
171
172 static inline void
173 make_umodpoly(umodpoly& up, const upoly& p, const cln::cl_modint_ring& R)
174 {
175         for (std::size_t i = p.size(); i-- != 0; )
176                 up[i] = R->canonhom(p[i]);
177         canonicalize(up);
178 }
179
180 } // namespace GiNaC
181
182 #endif // GINAC_UPOLY_H