Frobby  0.9.0
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 
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 }
Cont::const_iterator const_iterator
Definition: Ideal.h:43
bool genericBaseCase(const BigattiState &state)
Returns ture if state is a base case slice while also considering genericity.
Used in enumerateScarfComplex and necessary to have here to define _states.
size_t getTermCount() const
void add(const mpz_class &coef, const Term &term)
Add coef*term to the polynomial.
size_t getVarCount() const
Definition: Term.h:85
bool strictlyContains(const Exponent *term) const
Definition: Ideal.cpp:73
const Term & getMultiply() const
#define ASSERT(X)
Definition: stdinc.h:85
const_iterator begin() const
Definition: Ideal.h:48
void feedTo(const TermTranslator &translator, CoefBigTermConsumer &consumer, bool inCanonicalOrder) const
size_t getTermCount() const
Represents a monomial ideal with int exponents.
Definition: Ideal.h:27
const_iterator end() const
Definition: Ideal.h:49
void enumerateScarfComplex(const BigattiState &state, bool allFaces)
The ideal in state must be weakly generic.
bool isZeroIdeal() const
Definition: Ideal.cpp:86
BigattiBaseCase(const TermTranslator &translator)
Initialize this object to handle the computation of Hilbert-Poincare series numerator polynomials in ...
const mpz_class & getExponent(size_t variable, Exponent exponent) const
This method translates from IDs to arbitrary precision integers.
size_t getTotalTermsOutputEver() const
Returns the total number of terms this object has output.
size_t getVarCount() const
bool disjointSupport() const
Returns true if all pairs of generators have disjoint support.
Definition: Ideal.cpp:142
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:272
void feedOutputTo(CoefBigTermConsumer &consumer, bool inCanonicalOrder)
Feed the output Hilbert-Poincare numerator polynomial computed so far to the consumer.
size_t _totalTermsOutputEver
For statistics.
const Ideal & getIdeal() const
size_t getGeneratorCount() const
Definition: Ideal.h:57
void output(bool plus, const Term &term)
Add +term or -term to the output polynomial when plus is true or false respectively.
static size_t getSizeOfSupport(const Exponent *a, size_t varCount)
Returns the number of variables such that divides .
Definition: Term.h:394
size_t getVarCount() const
Definition: Ideal.h:56
void feedTo(CoefBigTermConsumer &consumer, bool inCanonicalOrder=false) const
UniHashPolynomial _outputUnivariate
The part of the coarsely graded Hilbert-Poincare numerator polynomial computed so far...
void setComputeUnivariate(bool value)
Use the fine grading if value is false, otherwise grade each variable by the same variable t...
bool simpleBaseCase(const BigattiState &state)
Computes the Hilbert-Poincare series of state and returns true if state is a particularly simple and ...
TermTranslator handles translation between terms whose exponents are infinite precision integers and ...
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...
const TermTranslator & _translator
Used to translate the output from ints.
void setPrintDebug(bool value)
Starts to print debug output on what happens if value is true.
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
vector< size_t > _maxCount
HashPolynomial _outputMultivariate
The part of the finely graded Hilbert-Poincare numerator polynomial computed so far.
bool univariateAllFaces(const BigattiState &state)
size_t _totalBaseCasesEver
For statistics.
vector< State > _states
Used in enumerateScarfCompex.
size_t getTotalBaseCasesEver() const
Returns the total number of base cases this object has seen.
bool isWeaklyGeneric() const
Definition: Ideal.cpp:121
bool baseCase(const BigattiState &state)
Returns true if state is a base case slice without considering genericity.
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:213
size_t getTotalTermsInOutput() const
Returns the number of terms in the output polynomial right now.
unsigned int Exponent
Definition: stdinc.h:88
Term represents a product of variables which does not include a coefficient.
Definition: Term.h:49
bool _computeUnivariate
Use the fine grading if false, otherwise grade each variable by the same variable t...
void getLcm(Exponent *lcm) const
Sets lcm to the least common multiple of all generators.
Definition: Ideal.cpp:157