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();
64  Ideal::const_iterator it = state.getIdeal().begin();
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 
86 void 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 
95  if (_computeUnivariate) {
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);
102  _outputUnivariate.add(plus, _tmp);
103  } else
104  _outputMultivariate.add(plus, term);
105 }
106 
108 (CoefBigTermConsumer& consumer,
109  bool inCanonicalOrder) {
110  if (_computeUnivariate)
111  _outputUnivariate.feedTo(consumer, inCanonicalOrder);
112  else
113  _outputMultivariate.feedTo(_translator, consumer, inCanonicalOrder);
114 }
115 
117  _printDebug = value;
118 }
119 
121  _computeUnivariate = value;
122 }
123 
125  return _totalBaseCasesEver;
126 }
127 
129  return _totalTermsOutputEver;
130 }
131 
133  if (_computeUnivariate)
135  else
137 }
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
Cont::const_iterator const_iterator
Definition: Ideal.h:43
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
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,...
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