Frobby 0.9.5
TermGrader.cpp
Go to the documentation of this file.
1/* Frobby: Software for monomial ideal computations.
2 Copyright (C) 2007 Bjarke Hammersholt Roune (www.broune.com)
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2 of the License, or
7 (at your option) any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see http://www.gnu.org/licenses/.
16*/
17#include "stdinc.h"
18#include "TermGrader.h"
19
20#include "Projection.h"
21#include "TermTranslator.h"
22#include "Term.h"
23
24TermGrader::TermGrader(const vector<mpz_class>& varDegrees,
25 const TermTranslator& translator):
26 _grades(varDegrees.size()) {
27
28 // Set up _signs.
29 _signs.resize(varDegrees.size());
30 for (size_t var = 0; var < varDegrees.size(); ++var) {
31 if (varDegrees[var] > 0)
32 _signs[var] = 1;
33 else if (varDegrees[var] < 0)
34 _signs[var] = -1;
35 }
36
37 // Set up _grades.
38 for (size_t var = 0; var < varDegrees.size(); ++var) {
39 size_t maxId = translator.getMaxId(var);
40 _grades[var].resize(maxId + 1);
41
42 for (Exponent e = 0; e <= maxId; ++e)
43 _grades[var][e] = varDegrees[var] * translator.getExponent(var, e);
44 }
45}
46
47mpz_class TermGrader::getDegree(const Term& term) const {
48 mpz_class degree;
49 getDegree(term, degree);
50 return degree;
51}
52
53void TermGrader::getDegree(const Term& term, mpz_class& degree) const {
54 ASSERT(term.getVarCount() == _grades.size());
55 degree = 0;
56 for (size_t var = 0; var < term.getVarCount(); ++var)
57 degree += getGrade(var, term[var]);
58}
59
60void TermGrader::getDegree(const Term& term,
61 const Projection& projection,
62 mpz_class& degree) const {
63 ASSERT(term.getVarCount() == projection.getRangeVarCount());
64 degree = 0;
65 for (size_t var = 0; var < term.getVarCount(); ++var)
66 degree += getGrade(projection.inverseProjectVar(var), term[var]);
67}
68
69void TermGrader::getUpperBound(const Term& divisor,
70 const Term& dominator,
71 mpz_class& bound) const {
72 ASSERT(divisor.getVarCount() == getVarCount());
73 ASSERT(dominator.getVarCount() == getVarCount());
74 ASSERT(divisor.divides(dominator));
75
76 bound = 0;
77 size_t varCount = getVarCount();
78 for (size_t var = 0; var < varCount; ++var) {
79 int sign = getGradeSign(var);
80 if (sign == 0)
81 continue;
82
83 Exponent div = divisor[var];
84 Exponent dom = dominator[var];
85
86 Exponent optimalExponent;
87 if (div == dom)
88 optimalExponent = div; // Nothing to decide in this case.
89 else if (sign > 0) {
90 // In this case we normally prefer a high exponent.
91 //
92 // When computing irreducible decomposition or Alexander dual,
93 // we add pure powers of maximal degree that map to zero, in
94 // which case we want to avoid using that degree. This happens
95 // for dom == getMaxExponent(var).
96 if (dom == getMaxExponent(var)) {
97 ASSERT(getGrade(var, dom - 1) > getGrade(var, dom));
98 optimalExponent = dom - 1; // OK as div < dom.
99 } else
100 optimalExponent = dom;
101 } else {
102 ASSERT(sign < 0);
103
104 // In this case we normally prefer a low exponent. However, as
105 // above, we need to consider that the highest exponent could
106 // map to zero, which may be better.
107 if (dom == getMaxExponent(var)) {
108 ASSERT(getGrade(var, dom) > getGrade(var, div));
109 optimalExponent = dom;
110 } else
111 optimalExponent = div;
112 }
113
114 bound += getGrade(var, optimalExponent);
115 }
116}
117
118mpz_class TermGrader::getUpperBound(const Term& divisor,
119 const Term& dominator) const {
120 mpz_class bound;
121 getUpperBound(divisor, dominator, bound);
122 return bound;
123}
124
126(size_t var,
127 Exponent from,
128 Exponent to,
129 Exponent& index,
130 const mpz_class& maxDegree) const {
131 ASSERT(var < getVarCount());
132 ASSERT(from < _grades[var].size());
133 ASSERT(to < _grades[var].size());
134
135 if (from > to)
136 return false;
137
138 Exponent e = from;
139 while (true) {
140 const mpz_class& exp = _grades[var][e];
141 if (exp <= maxDegree) {
142 index = e;
143 return true;
144 }
145
146 if (e == to)
147 return false;
148 ++e;
149 }
150}
151
153(size_t var,
154 Exponent from,
155 Exponent to,
156 Exponent& index,
157 const mpz_class& maxDegree) const {
158 ASSERT(var < getVarCount());
159 ASSERT(from < _grades[var].size());
160 ASSERT(to < _grades[var].size());
161
162 if (from > to)
163 return false;
164
165 Exponent e = to;
166 while (true) {
167 const mpz_class& exp = _grades[var][e];
168 if (exp <= maxDegree) {
169 index = e;
170 return true;
171 }
172
173 if (e == from)
174 return false;
175 --e;
176 }
177}
178
180(size_t var, const mpz_class& value, bool strict) const {
181 ASSERT(var < getVarCount());
182 ASSERT(!_grades[var].empty());
183
184 bool first = true;
185 size_t best = 0;
186
187 for (size_t e = 1; e < _grades[var].size(); ++e) {
188 const mpz_class& exp = _grades[var][e];
189
190 if (exp <= value && (first || exp > _grades[var][best])) {
191 best = e;
192 first = false;
193 }
194 }
195
196 return best;
197}
198
200 const mpz_class& value, bool strict) const {
201 ASSERT(from <= to);
202
203 // If sign is negative, reverse the roles of < and > below.
204 int sign = getGradeSign(var);
205 if (sign == 0)
206 return 0;
207 bool positive = sign > 0;
208
209 // We are expecting that the correct value will usually be close to
210 // from, so we start with an exponential search starting at from and
211 // then move to a binary search when the endpoints become close.
212
213 Exponent low = from;
214 Exponent high = to;
215
216 // We carry on as though strict is true, and adjust the value
217 // below. The invariant is that degree(low) <= value < degree(high +
218 // 1), if that is true to begin with. You can check that both the
219 // cases value < degree(from) and degree(high + 1) <= value work out
220 // also.
221 while (true) {
222 ASSERT(low <= high);
223 size_t gap = high - low;
224 if (gap == 0)
225 break;
226
227 Exponent lowDelta = low - from;
228
229 // pivot is the point we do binary or exponential search on.
230 Exponent pivot;
231 if (lowDelta < gap) {
232 // In this case we have not moved much from the lower endpoint,
233 // so we double the distance, and add one in case lowDelta is
234 // zero.
235 pivot = low + lowDelta + 1;
236 } else {
237 // We use binary search. This formula sets pivot to be the
238 // average of low and high rounded up, while avoiding the
239 // possible overflow inherent in adding low and high.
240 pivot = low + (gap + 1) / 2;
241 }
242 ASSERT(low < pivot);
243 ASSERT(pivot <= high);
244
245 if (positive ? getGrade(var, pivot) <= value : getGrade(var, pivot) >= value) {
246 low = pivot;
247 }
248 else {
249 high = pivot - 1;
250 }
251 }
252 ASSERT(low == high);
253
254#ifdef DEBUG
255 Exponent reference = getLargestLessThan2(var, value, strict);
256 if (reference < from)
257 reference = from;
258 if (reference > to)
259 reference = to;
260 ASSERT(low == reference);
261#endif
262
263 return low;
264}
265
267 const Projection& projection,
268 mpz_class& degree) const {
269 ASSERT(term.getVarCount() == projection.getRangeVarCount());
270 degree = 0;
271 for (size_t var = 0; var < term.getVarCount(); ++var)
272 degree += getGrade(projection.inverseProjectVar(var), term[var] + 1);
273}
274
275const mpz_class& TermGrader::getGrade(size_t var, Exponent exponent) const {
276 ASSERT(var < _grades.size());
277 ASSERT(exponent < _grades[var].size());
278
279 return _grades[var][exponent];
280}
281
283 ASSERT(!_grades[var].empty());
284 return _grades[var].size() - 1;
285}
286
288 return _grades.size();
289}
290
291void TermGrader::print(ostream& out) const {
292 out << "TermGrader (\n";
293 for (size_t var = 0; var < _grades.size(); ++var) {
294 out << " var " << var << ':';
295 for (size_t e = 0; e < _grades[var].size(); ++e)
296 out << ' ' << _grades[var][e];
297 out << '\n';
298 }
299 out << ")\n";
300}
301
302int TermGrader::getGradeSign(size_t var) const {
303 ASSERT(var < _grades.size());
304 return _signs[var];
305}
306
307ostream& operator<<(ostream& out, const TermGrader& grader) {
308 grader.print(out);
309 return out;
310}
ostream & operator<<(ostream &out, const TermGrader &grader)
size_t inverseProjectVar(size_t rangeVar) const
size_t getRangeVarCount() const
A TermGrader assigns a value, the degree, to each monomial.
Definition TermGrader.h:27
void getIncrementedDegree(const Term &term, const Projection &projection, mpz_class &degree) const
void print(ostream &out) const
vector< int > _signs
Definition TermGrader.h:120
Exponent getLargestLessThan2(size_t var, const mpz_class &value, bool strict=true) const
Returns the index of the largest stored exponent of var that is less than value.
mpz_class getDegree(const Term &term) const
Returns the degree of term.
size_t getVarCount() const
Exponent getMaxExponent(size_t var) const
TermGrader(const vector< mpz_class > &varDegrees, const TermTranslator &translator)
bool getMaxIndexLessThan(size_t var, Exponent from, Exponent to, Exponent &index, const mpz_class &maxDegree) const
Finds maximal index in [from, to] to such that degree(t) <= maxDegree.
void getUpperBound(const Term &divisor, const Term &dominator, mpz_class &bound) const
Assigns to bound the degree of the largest term v such that divisor divides v and v divides dominator...
bool getMinIndexLessThan(size_t var, Exponent from, Exponent to, Exponent &index, const mpz_class &maxDegree) const
Finds minimal index in [from, to] to such that degree(t) <= maxDegree.
int getGradeSign(size_t var) const
Returns 1 if the grade strictly increases with the exponent of var, returns -1 if it strictly decreas...
const mpz_class & getGrade(size_t var, Exponent exponent) const
vector< vector< mpz_class > > _grades
Definition TermGrader.h:119
TermTranslator handles translation between terms whose exponents are infinite precision integers and ...
const mpz_class & getExponent(size_t variable, Exponent exponent) const
This method translates from IDs to arbitrary precision integers.
Exponent getMaxId(size_t variable) const
The assigned IDs are those in the range [0, getMaxId(var)].
Term represents a product of variables which does not include a coefficient.
Definition Term.h:49
size_t getVarCount() const
Definition Term.h:85
static bool divides(const Exponent *a, const Exponent *b, size_t varCount)
Returns whether a divides b.
Definition Term.h:152
This header file includes common definitions and is included as the first line of code in every imple...
unsigned int Exponent
Definition stdinc.h:89
#define ASSERT(X)
Definition stdinc.h:86