Frobby 0.9.5
BigattiBaseCase.cpp
Go to the documentation of this file.
1/* Frobby: Software for monomial ideal computations.
2 Copyright (C) 2009 University of Aarhus
3 Contact Bjarke Hammersholt Roune for license information (www.broune.com)
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program. If not, see http://www.gnu.org/licenses/.
17*/
18#include "stdinc.h"
19#include "BigattiBaseCase.h"
20
21#include "BigattiState.h"
22#include "TermTranslator.h"
23#include <algorithm>
24
26 _maxCount(translator.getVarCount()),
27 _lcm(translator.getVarCount()),
28 _outputMultivariate(translator.getVarCount()),
29 _computeUnivariate(false),
30 _translator(translator),
31 _totalBaseCasesEver(0),
32 _totalTermsOutputEver(0),
33 _printDebug(false) {
34}
35
37 if (baseCase(state))
38 return true;
39
40 if (!state.getIdeal().isWeaklyGeneric())
41 return false;
42
43 enumerateScarfComplex(state, false);
44
46 return true;
47}
48
50 ASSERT(_maxCount.size() == state.getVarCount());
51
52 if (simpleBaseCase(state))
53 return true;
54
55 if (state.getIdeal().getGeneratorCount() > state.getVarCount())
56 return false;
57
58 state.getIdeal().getLcm(_lcm);
60 return false;
61
62 fill(_maxCount.begin(), _maxCount.end(), 0);
63 Ideal::const_iterator end = state.getIdeal().end();
65 for (; it != end; ++it) {
66 bool hasMax = false;
67 for (size_t var = 0; var < state.getVarCount(); ++var) {
68 ASSERT((*it)[var] <= _lcm[var]);
69 if ((*it)[var] == _lcm[var] && _lcm[var] > 0) {
70 hasMax = true;
71 _maxCount[var] += 1;
72 if (_maxCount[var] > 1)
73 return false;
74 }
75 }
76 if (!hasMax)
77 return false;
78 }
79
80 enumerateScarfComplex(state, true);
81
83 return true;
84}
85
86void BigattiBaseCase::output(bool plus, const Term& term) {
87 if (_printDebug) {
88 fputs("Debug: Outputting term ", stderr);
89 fputc(plus ? '+' : '-', stderr);
90 term.print(stderr);
91 fputs(".\n", stderr);
92 }
93
96 if (term.getVarCount() == 0)
97 _tmp = 0;
98 else
99 _tmp = _translator.getExponent(0, term);
100 for (size_t var = 1; var < term.getVarCount(); ++var)
101 _tmp += _translator.getExponent(var, term);
103 } else
104 _outputMultivariate.add(plus, term);
105}
106
108(CoefBigTermConsumer& consumer,
109 bool inCanonicalOrder) {
111 _outputUnivariate.feedTo(consumer, inCanonicalOrder);
112 else
113 _outputMultivariate.feedTo(_translator, consumer, inCanonicalOrder);
114}
115
117 _printDebug = value;
118}
119
123
127
131
138
140 const Ideal& ideal = state.getIdeal();
141 size_t genCount = ideal.getGeneratorCount();
142 const Term& multiply = state.getMultiply();
143
144 if (genCount > 2)
145 return false;
146
147 output(true, multiply);
148 if (genCount == 0)
149 return true;
150
151 _lcm.product(multiply, ideal[0]);
152 output(false, _lcm);
153 if (genCount == 1)
154 return true;
155
156 ASSERT(genCount == 2);
157 _lcm.product(multiply, ideal[1]);
158 output(false, _lcm);
159
160 _lcm.lcm(ideal[0], ideal[1]);
161 _lcm.product(_lcm, multiply);
162 output(true, _lcm);
163
165 return true;
166}
167
170 const Ideal& ideal = state.getIdeal();
171 const Term& multiply = state.getMultiply();
172
173 if (!ideal.disjointSupport())
174 return false;
175
176 if (ideal.getGeneratorCount() > 30)
177 return false; // Coefficients may not fit in 32 bits
178
179 Term max(ideal.getVarCount());
180 ideal.getLcm(max);
181 max.product(max, multiply);
182
183 _tmp = 0;
184 for (size_t var = 0; var < max.getVarCount(); ++var)
185 _tmp += _translator.getExponent(var, max);
186 if (_tmp > 1024*1024)
187 return false; // Too high memory requirement.
188 ASSERT(_tmp.fits_uint_p());
189 size_t maxDegree = _tmp.get_ui();
190
191 size_t approxWorkForScarfComplex = (1 << ideal.getGeneratorCount());
192 if (approxWorkForScarfComplex < maxDegree)
193 return false; // Scarf complex is on par or faster.
194
195 // At this point we have made sure that this instance makes sense to
196 // solve using this method. We have also determined that we can do
197 // everything in machine ints.
198
199 vector<int> poly;
200 poly.reserve(maxDegree);
201 poly.push_back(1);
202
203 // TODO: sort with smallest first to decrease work in early
204 // iterations.
205 for (size_t i = 0; i < ideal.getGeneratorCount(); ++i) {
206 ASSERT(poly.back() != 0);
207 const Exponent* gen = ideal[i];
208
209 // calculate degree
210 int degree = 0;
211 for (size_t var = 0; var < max.getVarCount(); ++var)
212 degree += _translator.getExponent(var, multiply[var] + gen[var]).get_ui()
213 - _translator.getExponent(var, multiply[var]).get_ui();
214
215 // replace poly P by (1-t^degree)*P = P-P*t^degree.
216 size_t oldSize = poly.size();
217 poly.resize(oldSize + degree);
218 for (size_t e = oldSize; e > 0;) {
219 --e;
220 poly[e+degree] -= poly[e];
221 }
222 }
223
224 int degree = 0;
225 for (size_t var = 0; var < max.getVarCount(); ++var)
226 degree += _translator.getExponent(var, multiply).get_ui();
227
228 for (size_t e = 0; e < poly.size(); ++e) {
229 if (_printDebug) {
230 fprintf(stderr, "Debug: Outputting term %i*t^%u.\n",
231 poly[e], (unsigned int)(e + degree));
232 }
233
235 _outputUnivariate.add(poly[e], e+degree);
236 }
237
238 return true;
239}
240
242 bool allFaces) {
243 if (allFaces &&
245 univariateAllFaces(state)) {
246 return;
247 }
248
249 const Ideal& ideal = state.getIdeal();
250
251 // Set up _states with enough entries of the right size.
252 size_t needed = ideal.getGeneratorCount() + 1;
253 if (_states.size() < needed)
254 _states.resize(needed);
255 for (size_t i = 0; i < _states.size(); ++i)
256 _states[i].term.reset(state.getVarCount());
257
258 // Set up the initial state
259 ASSERT(!ideal.isZeroIdeal());
260 _states[0].plus = true;
261 _states[0].pos = ideal.begin();
262 ASSERT(_states[0].term.isIdentity());
263
264 // Cache this to avoid repeated calls to end().
265 Ideal::const_iterator stop = ideal.end();
266
267 // Iterate until all states are done. The active entries of _states
268 // are those from index 0 up to and including _current.
269 size_t current = 0;
270 while (true) {
271 ASSERT(current < _states.size());
272 State& currentState = _states[current];
273
274 if (currentState.pos == stop) {
275 // This is a base case since we have considered all minimal
276 // generators.
277 _lcm.product(currentState.term, state.getMultiply());
278 output(currentState.plus, _lcm);
279
280 // We are done with this entry, so go back to the previous
281 // active entry.
282 if (current == 0)
283 break; // Nothing remains to be done.
284 --current;
285 } else {
286 // Split into two cases according to whether we put the minimal
287 // generator at pos into the face or not.
288 ASSERT(current + 1 < _states.size());
289 State& next = _states[current + 1];
290
291 next.term.lcm(currentState.term, *currentState.pos);
292 ++currentState.pos;
293
294 if (allFaces || !ideal.strictlyContains(next.term)) {
295 // If allFaces is true we do not need to check the condition
296 // since we know it should always hold.
297 ASSERT(!ideal.strictlyContains(next.term));
298
299 next.plus = !currentState.plus;
300 next.pos = currentState.pos;
301 ++current;
302 }
303 }
304 }
305}
BigattiBaseCase(const TermTranslator &translator)
Initialize this object to handle the computation of Hilbert-Poincare series numerator polynomials in ...
size_t _totalBaseCasesEver
For statistics.
size_t getTotalTermsOutputEver() const
Returns the total number of terms this object has output.
void setPrintDebug(bool value)
Starts to print debug output on what happens if value is true.
bool genericBaseCase(const BigattiState &state)
Returns ture if state is a base case slice while also considering genericity.
bool univariateAllFaces(const BigattiState &state)
size_t getTotalBaseCasesEver() const
Returns the total number of base cases this object has seen.
void setComputeUnivariate(bool value)
Use the fine grading if value is false, otherwise grade each variable by the same variable t.
vector< size_t > _maxCount
bool simpleBaseCase(const BigattiState &state)
Computes the Hilbert-Poincare series of state and returns true if state is a particularly simple and ...
bool _computeUnivariate
Use the fine grading if false, otherwise grade each variable by the same variable t.
bool baseCase(const BigattiState &state)
Returns true if state is a base case slice without considering genericity.
void enumerateScarfComplex(const BigattiState &state, bool allFaces)
The ideal in state must be weakly generic.
HashPolynomial _outputMultivariate
The part of the finely graded Hilbert-Poincare numerator polynomial computed so far.
const TermTranslator & _translator
Used to translate the output from ints.
size_t _totalTermsOutputEver
For statistics.
vector< State > _states
Used in enumerateScarfCompex.
void output(bool plus, const Term &term)
Add +term or -term to the output polynomial when plus is true or false respectively.
UniHashPolynomial _outputUnivariate
The part of the coarsely graded Hilbert-Poincare numerator polynomial computed so far.
void feedOutputTo(CoefBigTermConsumer &consumer, bool inCanonicalOrder)
Feed the output Hilbert-Poincare numerator polynomial computed so far to the consumer.
size_t getTotalTermsInOutput() const
Returns the number of terms in the output polynomial right now.
const Term & getMultiply() const
size_t getVarCount() const
const Ideal & getIdeal() const
void add(const mpz_class &coef, const Term &term)
Add coef*term to the polynomial.
size_t getTermCount() const
void feedTo(const TermTranslator &translator, CoefBigTermConsumer &consumer, bool inCanonicalOrder) const
Represents a monomial ideal with int exponents.
Definition Ideal.h:27
bool isZeroIdeal() const
Definition Ideal.cpp:86
size_t getGeneratorCount() const
Definition Ideal.h:57
void getLcm(Exponent *lcm) const
Sets lcm to the least common multiple of all generators.
Definition Ideal.cpp:157
bool strictlyContains(const Exponent *term) const
Definition Ideal.cpp:73
bool disjointSupport() const
Returns true if all pairs of generators have disjoint support.
Definition Ideal.cpp:142
const_iterator end() const
Definition Ideal.h:49
const_iterator begin() const
Definition Ideal.h:48
bool isWeaklyGeneric() const
Definition Ideal.cpp:121
Cont::const_iterator const_iterator
Definition Ideal.h:43
size_t getVarCount() const
Definition Ideal.h:56
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.
Term represents a product of variables which does not include a coefficient.
Definition Term.h:49
static size_t getSizeOfSupport(const Exponent *a, size_t varCount)
Returns the number of variables such that divides .
Definition Term.h:402
static void product(Exponent *res, const Exponent *a, const Exponent *b, size_t varCount)
Sets res equal to the product of a and b.
Definition Term.h:280
size_t getVarCount() const
Definition Term.h:85
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
static void lcm(Exponent *res, const Exponent *a, const Exponent *b, size_t varCount)
Sets res equal to the least commom multiple of a and b.
Definition Term.h:221
size_t getTermCount() const
void feedTo(CoefBigTermConsumer &consumer, bool inCanonicalOrder=false) const
void add(bool plus, const mpz_class &exponent)
Add +t^exponent or -t^exponent to the polynomial depending on whether plus is true or false,...
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
Used in enumerateScarfComplex and necessary to have here to define _states.
Ideal::const_iterator pos