Frobby 0.9.5
HilbertBasecase.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 "HilbertBasecase.h"
19
20#include "Ideal.h"
21#include "Term.h"
22#include "error.h"
23
25 _idealCacheDeleter(_idealCache),
26 _stepsPerformed(0) {
27}
28
32
33// Performs one or more steps of computation. Steps are performed on
34// entry until entry becomes a base case, or until it needs to be
35// split into two.
36//
37// If entry becomes a base case, then false is returned, entry can be
38// discarded and newEntry is unchanged.
39//
40// If entry needs to be split into two, then true is returned, entry
41// becomes one of the subcomputations, and newEntry becomes the other.
44
45 size_t varCount = entry.ideal->getVarCount();
46
47 // This loops keeps iterating as long as entry is not a base case
48 // and there is some computation that can be done on entry that does
49 // not require splitting it into two.
50 while (true) {
51 // Here _term is used to contain support counts to choose best pivot and
52 // to detect base case.
53
54 // We start off checking different ways that entry can be a base
55 // case.
57 if (_term.getSizeOfSupport() + entry.extraSupport != varCount)
58 return false;
59
60 if (_term.isSquareFree()) {
61 if ((entry.ideal->getGeneratorCount() % 2) == 1)
62 entry.negate = !entry.negate;
63 if (entry.negate)
64 --_sum;
65 else
66 ++_sum;
67 return false;
68 }
69
70 if (entry.ideal->getGeneratorCount() == 2) {
71 if (entry.negate)
72 --_sum;
73 else
74 ++_sum;
75 return false;
76 }
77
78 // This is a simplification step, and if we can perform it, we
79 // just start over with the new entry we get that way. This is
80 // necessary because the base case checks below assume that this
81 // simplification has been performed.
82 size_t ridden = eliminate1Counts(*entry.ideal, _term, entry.negate);
83 if (ridden != 0) {
84 entry.extraSupport += ridden;
85 continue;
86 }
87
88 if (entry.ideal->getGeneratorCount() == 3) {
89 if (entry.negate)
90 _sum -= 2;
91 else
92 _sum += 2;
93 return false;
94 }
95
96 if (entry.ideal->getGeneratorCount() == 4 &&
98 _term.getSizeOfSupport() == 4) {
99 if (entry.negate)
100 ++_sum;
101 else
102 --_sum;
103 return false;
104 }
105
106 // At this point entry is not a base case, and it cannot be
107 // simplified, so we have to split it into two.
108
109 size_t bestPivotVar = _term.getFirstMaxExponent();
110
111 // Handle outer slice.
112 auto_ptr<Ideal> outer = getNewIdeal();
113 outer->clearAndSetVarCount(varCount);
114 outer->insertNonMultiples(bestPivotVar, 1, *entry.ideal);
115
116 // outer is subtracted instead of added due to having added the
117 // pivot to the ideal.
118 newEntry.negate = !entry.negate;
119 newEntry.extraSupport = entry.extraSupport + 1;
120 newEntry.ideal = outer.get();
121
122 // Handle inner slice in-place on entry.
123 entry.ideal->colonReminimize(bestPivotVar, 1);
124 ++entry.extraSupport;
125
126 outer.release();
127 return true;
128 }
129}
130
132 ASSERT(_todo.empty());
133
134 try { // Here to clear _todo in case of an exception
135 // _sum is updated as a side-effect of calling stepComputation.
136 _sum = 0;
137
138 // _term is reused for several different purposes in order to avoid
139 // having to allocate and deallocate the underlying data structure.
140 _term.reset(originalIdeal.getVarCount());
141
142 // entry is the Entry currently being processed. Additional entries
143 // are added to _todo, though this only happens if there are two,
144 // since otherwise entry can just be updated to the next value
145 // directly and so we avoid the overhead of using _todo when we can.
146 Entry entry;
147 entry.negate = false;
148 entry.extraSupport = 0;
149 entry.ideal = &originalIdeal;
150
151 // This should normally point to entry.ideal, but since we do not
152 // have ownership of originalIdeal, it starts out pointing nowhere.
153 auto_ptr<Ideal> entryIdealDeleter;
154
155 while (true) {
156 // Do an inner loop since there is no reason to add entry to _todo
157 // and then immediately take it off again.
158 Entry newEntry;
159 while (stepComputation(entry, newEntry)) {
160 auto_ptr<Ideal> newEntryIdealDeleter(newEntry.ideal);
161 _todo.push_back(newEntry);
162 newEntryIdealDeleter.release();
163 }
164
165 if (_todo.empty())
166 break;
167
168 if (entryIdealDeleter.get() != 0)
169 freeIdeal(entryIdealDeleter);
170 entry = _todo.back();
171 _todo.pop_back();
172
173 ASSERT(entryIdealDeleter.get() == 0);
174 entryIdealDeleter.reset(entry.ideal);
175 }
176 ASSERT(_todo.empty());
177
178 // originalIdeal is in some state that depends on the particular
179 // steps the algorithm took. This information should not escape
180 // HilbertBasecase, and we ensure this by clearing originalIdeal.
181 originalIdeal.clear();
182 } catch (...) {
183 for (vector<Entry>::iterator it = _todo.begin(); it != _todo.end(); ++it)
184 delete it->ideal;
185 _todo.clear();
186 throw;
187 }
188}
189
191 return _sum;
192}
193
195 const Ideal& ideal,
196 const Term& counts) {
197 if (counts[var] == 0)
198 return false;
199 Ideal::const_iterator stop = ideal.end();
200
201 size_t varCount = counts.getVarCount();
202 for (size_t other = 0; other < varCount; ++other) {
203 if (other == var || counts[other] == 0)
204 continue;
205
206 // todo: the answer is always no, I think, if var appears in less
207 // generators than other does, since then there must be some
208 // generator that other appears in that var does not.
209
210 bool can = true;
211 for (Ideal::const_iterator it = ideal.begin(); it != stop; ++it) {
212 if ((*it)[var] == 0 && (*it)[other] > 0) {
213 can = false;
214 break;
215 }
216 }
217 if (can)
218 return true;
219 }
220 return false;
221}
222
224 Term& counts,
225 bool& negate) {
226 size_t varCount = ideal.getVarCount();
227 size_t adj = 0;
228 for (size_t var = 0; var < varCount; ++var) {
229 if (counts[var] != 1)
230 continue;
231
232 Ideal::const_iterator it = ideal.getMultiple(var);
233 ASSERT(it != ideal.end());
234
235 for (size_t other = 0; other < varCount; ++other) {
236 if ((*it)[other] > 0) {
237 ++adj;
238 if (counts[other] == 1)
239 counts[other] = 0;
240 } else
241 counts[other] = 0;
242 }
243
244 for (size_t other = 0; other < varCount; ++other) {
245 if (counts[other] > 0) {
246 if (!ideal.colonReminimize(other, 1)) {
247 ideal.clear();
248 return 1;
249 }
250 }
251 }
252
253 it = ideal.getMultiple(var);
254 if (it == ideal.end()) {
255 ideal.clear();
256 return 1;
257 }
258 ideal.remove(it);
259 negate = !negate;
260
261 return adj;
262 }
263
264 for (size_t var = 0; var < varCount; ++var) {
265 if (canSimplify(var, ideal, counts)) {
266 if (!ideal.colonReminimize(var, 1))
267 ideal.clear();
268 return adj + 1;
269 }
270 }
271
272 return adj;
273}
274
276 if (_idealCache.empty())
277 return auto_ptr<Ideal>(new Ideal());
278
279 auto_ptr<Ideal> ideal(_idealCache.back());
280 _idealCache.pop_back();
281
282 return ideal;
283}
284
285void HilbertBasecase::freeIdeal(auto_ptr<Ideal> ideal) {
286 ASSERT(ideal.get() != 0);
287
288 ideal->clear();
290}
void exceptionSafePushBack(Container &container, auto_ptr< Element > pointer)
bool canSimplify(size_t var, const Ideal &ideal, const Term &counts)
vector< Ideal * > _idealCache
void computeCoefficient(Ideal &ideal)
void freeIdeal(auto_ptr< Ideal > ideal)
bool stepComputation(Entry &entry, Entry &newEntry)
size_t eliminate1Counts(Ideal &ideal, Term &counts, bool &negate)
vector< Entry > _todo
auto_ptr< Ideal > getNewIdeal()
const mpz_class & getLastCoefficient()
Represents a monomial ideal with int exponents.
Definition Ideal.h:27
void remove(const_iterator it)
Definition Ideal.cpp:576
size_t getGeneratorCount() const
Definition Ideal.h:57
const_iterator getMultiple(size_t var) const
Definition Ideal.cpp:668
void clear()
Definition Ideal.cpp:641
const_iterator end() const
Definition Ideal.h:49
const_iterator begin() const
Definition Ideal.h:48
void getSupportCounts(Exponent *counts) const
counts[var] will be the number of generators divisible by var.
Definition Ideal.cpp:227
bool colonReminimize(const Exponent *colon)
Definition Ideal.cpp:550
Cont::const_iterator const_iterator
Definition Ideal.h:43
size_t getVarCount() const
Definition Ideal.h:56
Term represents a product of variables which does not include a coefficient.
Definition Term.h:49
void reset(size_t newVarCount)
Definition Term.h:551
static size_t getSizeOfSupport(const Exponent *a, size_t varCount)
Returns the number of variables such that divides .
Definition Term.h:402
size_t getVarCount() const
Definition Term.h:85
static bool isSquareFree(const Exponent *a, size_t varCount)
Returns whether a is square free, i.e. for each where .
Definition Term.h:331
static size_t getFirstMaxExponent(const Exponent *a, size_t varCount)
Returns a var such that a[var] >= a[i] for all i.
Definition Term.h:381
This header file includes common definitions and is included as the first line of code in every imple...
#define ASSERT(X)
Definition stdinc.h:86