Frobby 0.9.5
Minimizer.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 "Minimizer.h"
19
20#include "TermPredicate.h"
21#include "Term.h"
22#include <algorithm>
23
28namespace {
29 typedef vector<Exponent*>::iterator TermIterator;
30}
31
32TermIterator simpleMinimize(TermIterator begin, TermIterator end, size_t varCount) {
33 if (begin == end)
34 return end;
35
36 std::sort(begin, end, LexComparator(varCount));
37
38 TermIterator newEnd = begin;
39 ++newEnd; // The first one is always kept
40 TermIterator dominator = newEnd;
41 for (; dominator != end; ++dominator) {
42 bool remove = false;
43 for (TermIterator divisor = begin; divisor != newEnd; ++divisor) {
44 if (Term::divides(*divisor, *dominator, varCount)) {
45 remove = true;
46 break;
47 }
48 }
49
50 if (!remove) {
51 *newEnd = *dominator;
52 ++newEnd;
53 }
54 }
55 return newEnd;
56}
57
58TermIterator twoVarMinimize(TermIterator begin, TermIterator end) {
59 if (begin == end)
60 return end;
61
62 std::sort(begin, end, LexComparator(2));
63
64 TermIterator last = begin;
65 TermIterator it = begin;
66 ++it;
67 for (; it != end; ++it) {
68 if ((*it)[1] < (*last)[1]) {
69 ++last;
70 *last = *it;
71 }
72 }
73
74 ++last;
75 return last;
76}
77
78class TreeNode {
79 typedef vector<Exponent*>::iterator iterator;
80
81public:
82 TreeNode(iterator begin, iterator end, size_t varCount):
83 _lessOrEqual(0),
84 _greater(0),
85 _var(0),
86 _pivot(0),
87 _varCount(varCount),
88 _begin(begin),
89 _end(end) {
90 }
91
93 delete _lessOrEqual;
94 delete _greater;
95 }
96
97 void makeTree() {
98 ASSERT(_greater == 0);
99 ASSERT(_lessOrEqual == 0);
100
101 if (distance(_begin, _end) > 20) {
102
103 Term lcm(_varCount);
104 for (iterator it = _begin; it != _end; ++it)
105 lcm.lcm(lcm, *it);
106
107 while (true) {
108 size_t maxVar = 0;
109
110 for (size_t var = 0; var < _varCount; ++var)
111 if (lcm[var] > lcm[maxVar])
112 maxVar = var;
113 if (lcm[maxVar] == 0) {
114 break; // we are not making any progress anyway
115 }
116
117 ASSERT(lcm[maxVar] >= 1);
118 _var = maxVar;
119 _pivot = lcm[maxVar] / 4;
120 lcm[maxVar] = 0; // So we do not process this same var again.
121
122 iterator left = _begin;
123 iterator right = _end - 1;
124 while (left != right) {
125 while ((*left)[_var] <= _pivot && left != right)
126 ++left;
127 while ((*right)[_var] > _pivot && left != right)
128 --right;
129 swap(*left, *right);
130 }
131 ASSERT((*right)[_var] > _pivot);
132 if ((*_begin)[_var] > _pivot)
133 continue;
134
135 iterator middle = right;
136 while ((*middle)[maxVar] > _pivot)
137 --middle;
138 ++middle;
139
140 ASSERT(middle != _begin);
141
142 _lessOrEqual = new TreeNode(_begin, middle, _varCount);
143 _greater = new TreeNode(middle, _end, _varCount);
144 _end = _begin;
145
148 return;
149 }
150 }
151
153
154 }
155
156 bool isRedundant(const Exponent* term) {
157 if (_begin != _end) {
158 ASSERT(_lessOrEqual == 0);
159 ASSERT(_greater == 0);
160
161 for (iterator it = _begin; it != _end; ++it)
162 if (Term::dominates(term, *it, _varCount))
163 return true;
164 return false;
165 } else {
166 ASSERT(_lessOrEqual != 0);
167 ASSERT(_greater != 0);
168
169 if (term[_var] > _pivot && _greater->isRedundant(term))
170 return true;
171 if (_lessOrEqual->isRedundant(term))
172 return true;
173 return false;
174 }
175 }
176
177 void collect(vector<Exponent*>& terms) {
178 if (_begin != _end) {
179 ASSERT(_lessOrEqual == 0);
180 ASSERT(_greater == 0);
181 terms.insert(terms.end(), _begin, _end);
182 return;
183 }
184 ASSERT(_lessOrEqual != 0);
185 ASSERT(_greater != 0);
186
187 size_t oldSize = terms.size();
188 _greater->collect(terms);
189 for (size_t i = oldSize; i < terms.size();) {
190 if (_lessOrEqual->isRedundant(terms[i])) {
191 swap(terms[i], terms.back());
192 terms.pop_back();
193 } else
194 ++i;
195 }
196
197 _lessOrEqual->collect(terms);
198 }
199
200 void print(FILE* out) {
201 if (_begin == _end) {
202 ASSERT(_lessOrEqual != 0);
203 ASSERT(_greater != 0);
204
205 fprintf(out, "NODE (pivot=%lu^%lu, varCount = %lu\n",
206 (unsigned long)_var,
207 (unsigned long)_pivot,
208 (unsigned long)_varCount);
209 fputs("-lessOrEqual: ", out);
210 _lessOrEqual->print(out);
211 fputs("-greater: ", out);
212 _greater->print(out);
213 fputs(")\n", out);
214 } else {
215 ASSERT(_lessOrEqual == 0);
216 ASSERT(_greater == 0);
217
218 fprintf(out, "NODE (_varCount = %lu terms:\n", (unsigned long)_varCount);
219 for (iterator it = _begin; it != _end; ++it) {
220 fputc(' ', out);
221 Term::print(out, *it, _varCount);
222 fprintf(out, " %p\n", (void*)*it);
223 }
224 fputs(")\n", out);
225 }
226 }
227
228private:
231 size_t _var;
233 size_t _varCount;
234
237};
238
240 if (_varCount == 2)
241 return twoVarMinimize(begin, end);
242 if (distance(begin, end) < 1000 || _varCount == 0)
243 return simpleMinimize(begin, end, _varCount);
244
245 vector<Exponent*> terms;
246 terms.clear();
247 terms.reserve(distance(begin, end));
248
249 TreeNode node(begin, end, _varCount);
250 node.makeTree();
251 node.collect(terms);
252
253 return copy(terms.begin(), terms.end(), begin);
254}
255
256pair<Minimizer::iterator, bool> Minimizer::colonReminimize
257(iterator begin, iterator end, const Exponent* colon) {
258 ASSERT(isMinimallyGenerated(begin, end));
259
260 if (Term::getSizeOfSupport(colon, _varCount) == 1) {
261 size_t var = Term::getFirstNonZeroExponent(colon, _varCount);
262 return colonReminimize(begin, end, var, colon[var]);
263 }
264
265 iterator blockBegin = end;
266 for (iterator it = begin; it != blockBegin;) {
267 bool block = true;
268 bool strictDivision = true;
269 for (size_t var = 0; var < _varCount; ++var) {
270 if (colon[var] >= (*it)[var]) {
271 if ((*it)[var] > 0)
272 block = false;
273 if (colon[var] > 0)
274 strictDivision = false;
275 (*it)[var] = 0;
276 } else
277 (*it)[var] -= colon[var];
278 }
279
280 if (strictDivision) {
281 swap(*begin, *it);
282 ++begin;
283 ++it;
284 } else if (block) {
285 --blockBegin;
286 swap(*it, *blockBegin);
287 } else
288 ++it;
289 }
290
291 if (begin == blockBegin)
292 return make_pair(end, false);
293
294 iterator newEnd = minimize(begin, blockBegin);
295
296 for (iterator it = blockBegin; it != end; ++it) {
297 if (!dominatesAny(begin, blockBegin, *it)) {
298 *newEnd = *it;
299 ++newEnd;
300 }
301 }
302
303 ASSERT(isMinimallyGenerated(begin, newEnd));
304 return make_pair(newEnd, true);
305}
306
308(const_iterator begin, const_iterator end) {
309 if (distance(begin, end) < 1000 || _varCount == 0) {
310 for (const_iterator divisor = begin; divisor != end; ++divisor)
311 for (const_iterator dominator = begin; dominator != end; ++dominator)
312 if (Term::divides(*divisor, *dominator, _varCount) &&
313 divisor != dominator)
314 return false;
315 return true;
316 }
317
318 vector<Exponent*> terms(begin, end);
319 TreeNode node(terms.begin(), terms.end(), _varCount);
320 node.makeTree();
321
322 vector<Exponent*> terms2;
323 node.collect(terms2);
324
325 return terms.size() == terms2.size();
326}
327
329(iterator begin, iterator end, const Exponent* term) {
330 for (; begin != end; ++begin)
331 if (Term::dominates(term, *begin, _varCount))
332 return true;
333 return false;
334}
335
337(iterator begin, iterator end, const Exponent* term) {
338 for (; begin != end; ++begin)
339 if (Term::divides(term, *begin, _varCount))
340 return true;
341 return false;
342}
343
344pair<Minimizer::iterator, bool> Minimizer::colonReminimize
345(iterator begin, iterator end, size_t var, Exponent exponent) {
346
347 // Sort in descending order according to exponent of var while
348 // ignoring everything that is strictly divisible by
349 // var^exponent. We put the zero entries at the right end
350 // immediately, before calling sort, because there are likely to be
351 // many of them, and we can do so while we are anyway looking for
352 // the strictly divisible monomials. The combination of these
353 // significantly reduce the number of monomials that need to be
354 // sorted.
355 iterator zeroBegin = end;
356 for (iterator it = begin; it != zeroBegin;) {
357 if ((*it)[var] > exponent) {
358 (*it)[var] -= exponent; // apply colon
359 swap(*it, *begin);
360 ++begin;
361 ++it;
362 } else if ((*it)[var] == 0) {
363 // no need to apply colon in this case
364 --zeroBegin;
365 swap(*it, *zeroBegin);
366 } else
367 ++it;
368 }
369
370 if (begin == zeroBegin)
371 return make_pair(end, false);
372
373 // Sort the part of the array that we have not handled yet.
374 std::sort(begin, zeroBegin,
376
377 // We group terms into blocks according to term[var].
378 iterator previousBlockEnd = begin;
379 iterator newEnd = begin;
380
381 Exponent block = (*begin)[var];
382
383 for (iterator it = begin; it != end; ++it) {
384 // Detect if we are moving on to next block.
385 if ((*it)[var] != block) {
386 block = (*it)[var];
387 previousBlockEnd = newEnd;
388 }
389
390 ASSERT((*it)[var] <= exponent);
391 (*it)[var] = 0;
392
393 bool remove = false;
394
395 for (iterator divisor = begin; divisor != previousBlockEnd; ++divisor) {
396 if (Term::divides(*divisor, *it, _varCount)) {
397 remove = true;
398 break;
399 }
400 }
401
402 if (!remove) {
403 *newEnd = *it;
404 ++newEnd;
405 }
406 }
407
408 return make_pair(newEnd, true);
409}
void swap(HilbertSlice &a, HilbertSlice &b)
TermIterator simpleMinimize(TermIterator begin, TermIterator end, size_t varCount)
Definition Minimizer.cpp:32
TermIterator twoVarMinimize(TermIterator begin, TermIterator end)
Definition Minimizer.cpp:58
A predicate that sorts terms according to lexicographic order.
bool isMinimallyGenerated(const_iterator begin, const_iterator end)
pair< iterator, bool > colonReminimize(iterator begin, iterator end, const Exponent *colon)
iterator minimize(iterator begin, iterator end) const
bool dominatesAny(iterator begin, iterator end, const Exponent *term)
size_t _varCount
Definition Minimizer.h:44
vector< Exponent * >::iterator iterator
Definition Minimizer.h:24
bool dividesAny(iterator begin, iterator end, const Exponent *term)
vector< Exponent * >::const_iterator const_iterator
Definition Minimizer.h:25
A predicate that sorts terms in weakly descending order according to degree of the specified variable...
Term represents a product of variables which does not include a coefficient.
Definition Term.h:49
static bool dominates(const Exponent *a, const Exponent *b, size_t varCount)
Returns whether a dominates b, i.e. whether b divides a.
Definition Term.h:172
static void print(FILE *file, const Exponent *e, size_t varCount)
Writes e to file in a format suitable for debug output.
Definition Term.cpp:110
size_t getFirstNonZeroExponent() const
Definition Term.h:354
size_t getSizeOfSupport() const
Definition Term.h:411
static bool divides(const Exponent *a, const Exponent *b, size_t varCount)
Returns whether a divides b.
Definition Term.h:152
iterator _begin
size_t _varCount
TreeNode(iterator begin, iterator end, size_t varCount)
Definition Minimizer.cpp:82
vector< Exponent * >::iterator iterator
Definition Minimizer.cpp:79
void print(FILE *out)
TreeNode * _greater
Exponent _pivot
void makeTree()
Definition Minimizer.cpp:97
iterator _end
TreeNode * _lessOrEqual
bool isRedundant(const Exponent *term)
void collect(vector< Exponent * > &terms)
size_t _var
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