Frobby 0.9.5
ScarfHilbertAlgorithm.cpp
Go to the documentation of this file.
1/* Frobby: Software for monomial ideal computations.
2 Copyright (C) 2010 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"
20
21#include "CoefTermConsumer.h"
22#include "TermTranslator.h"
23#include "Deformer.h"
24#include "CoefBigTermConsumer.h"
25#include "HashPolynomial.h"
26#include "UniHashPolynomial.h"
27#include "ScarfParams.h"
28#include "IdealTree.h"
29#include "IdealOrderer.h"
30
32public:
34 const TermTranslator& translator,
35 CoefBigTermConsumer& consumer,
36 const IdealOrderer& order,
37 bool univar,
38 bool canonical,
39 bool doStrongDeformation):
40 _univar(univar),
41 _tmp(toDeform.getVarCount()),
42 _deformer(toDeform, order, doStrongDeformation),
43 _translator(translator),
44 _canonical(canonical),
45 _consumer(consumer),
46 _poly(toDeform.getVarCount()) {
47 }
48
49 virtual void consumeRing(const VarNames& names) {
50 ASSERT(names == _translator.getNames());
51 }
52
53 virtual void beginConsuming() {
54 }
55
56 virtual void consume(const mpz_class& coef, const Term& term) {
58 _tmp = term;
60
61 if (_univar) {
62 if (_tmp.getVarCount() == 0)
63 _tdeg = 0;
64 else
66 for (size_t var = 1; var < _tmp.getVarCount(); ++var)
68 _uniPoly.add(coef, _tdeg);
69 } else
70 _poly.add(coef, _tmp);
71 }
72
79
80private:
81 bool _univar;
89 mpz_class _tdeg;
90};
91
93(const TermTranslator& translator,
94 const ScarfParams& params,
95 auto_ptr<IdealOrderer> enumerationOrder,
96 auto_ptr<IdealOrderer> deformationOrder):
97 _translator(translator),
98 _params(params),
99 _enumerationOrder(enumerationOrder),
100 _deformationOrder(deformationOrder),
101 _totalStates(0),
102 _totalFaces(0) {
103 ASSERT(_enumerationOrder.get() != 0);
104 ASSERT(_deformationOrder.get() != 0);
105}
106
108 // Destructor defined so auto_ptr<T> in the header does not need
109 // definition of T.
110}
111
113 CoefBigTermConsumer& consumer,
114 bool univariate,
115 bool canonical) {
116 Ideal deformed(ideal);
117 UndeformConsumer undeformer(deformed,
119 consumer,
121 univariate,
122 canonical,
124
125 undeformer.consumeRing(_translator.getNames());
126 undeformer.beginConsuming();
127 ASSERT(_enumerationOrder.get() != 0);
128 _enumerationOrder->order(deformed);
129 enumerateScarfComplex(deformed, undeformer);
130 undeformer.doneConsuming();
131
133 fputs("*** Statistics ***\n", stderr);
134 fprintf(stderr, "Total states considered: %u\n",
135 static_cast<unsigned int>(_totalStates));
136 fprintf(stderr, "Total faces accepted: %u\n",
137 static_cast<unsigned int>(_totalFaces));
138 }
139}
140
142 size_t& activeStateCount) {
144
145 if (_params.getPrintDebug()) {
146 fputs("Enumerating faces of Scarf complex of:\n", stderr);
147 ideal.print(stderr);
148 }
149
150 // Set up _states with enough entries. The maximal number of active
151 // entries at any time is one for each generator plus one for the
152 // empty face. We need one more than this because we take a
153 // reference to the next state even when there is no next state.
154 size_t statesNeeded = ideal.getGeneratorCount() + 2;
155 if (_states.size() < statesNeeded) {
156 _states.resize(statesNeeded);
157 for (size_t i = 0; i < _states.size(); ++i) {
158 _states[i].term.reset(ideal.getVarCount());
159 _states[i].face.reserve(ideal.getVarCount());
160 }
161 }
162
163 // Set up the initial state
164 activeStateCount = 0;
165 if (ideal.containsIdentity())
166 return;
167
168 ++activeStateCount;
169 _states[0].plus = true;
170 _states[0].pos = ideal.begin();
171 ASSERT(_states[0].term.isIdentity());
172}
173
175 const IdealTree& tree,
176 State& state,
177 State& nextState) {
178 if (_params.getPrintDebug()) {
179 fputs("DEBUG:*Looking at element ", stderr);
180 if (state.pos == ideal.end())
181 fputs("end", stderr);
182 else
183 Term::print(stderr, *state.pos, ideal.getVarCount());
184 fputs(" with lcm(face)=", stderr);
185 state.term.print(stderr);
186 fputs(" and face=", stderr);
187 if (state.face.empty())
188 fputs("empty", stderr);
189 for (size_t i = 0; i < state.face.size(); ++i) {
190 fputs("\nDEBUG: ", stderr);
191 Term::print(stderr, state.face[i], ideal.getVarCount());
192 }
193 fputc('\n', stderr);
194 fflush(stderr);
195 }
196
197 Exponent* termToAdd;
198 while (true) {
199 ++_totalStates;
200 if (state.face.size() == ideal.getVarCount() || state.pos == ideal.end())
201 return false; // A base case
202
203 termToAdd = *state.pos;
204
205 // This accounts for the possibility of not adding termToAdd to the
206 // face. We do that in-place on state.
207 ++state.pos;
208
209 // The other possibility is to add termToAdd to the face. We record
210 // this only if that is still a face of the complex, i.e. if no
211 // generator strictly divides the lcm of the set.
212 nextState.term.lcm(state.term, termToAdd);
213 // The elements of face are more likely to become strict divisors
214 // than a random generator, so check those first.
215 for (size_t i = 0; i < state.face.size(); ++i) {
216 if (Term::strictlyDivides(state.face[i],
217 nextState.term,
218 ideal.getVarCount())) {
219 goto doNext;
220 }
221 }
222 if (tree.strictlyContains(nextState.term))
223 goto doNext;
224 ASSERT(!ideal.strictlyContains(nextState.term));
225 break;
226 doNext:;
227 ASSERT(ideal.strictlyContains(nextState.term));
228 }
229
230 nextState.plus = !state.plus;
231 nextState.pos = state.pos;
232 nextState.face = state.face;
233 nextState.face.push_back(termToAdd);
234
235 return true;
236}
237
239 CoefTermConsumer& consumer) {
240 if (_params.getPrintDebug()) {
241 fputs("DEBUG: Found base case with lcm(face)=", stderr);
242 state.term.print(stderr);
243 fputc('\n', stderr);
244 fflush(stderr);
245 }
246
247 consumer.consume(state.plus ? 1 : -1, state.term);
248
249 // Every face ends up as a base case exactly once, so this is a
250 // convenient place to count them.
251 ++_totalFaces;
252}
253
255 CoefTermConsumer& consumer) {
256 ASSERT(Ideal(ideal).isWeaklyGeneric());
257
258 IdealTree tree(ideal);
259
260 size_t activeStateCount = 0;
261 initializeEnumeration(ideal, activeStateCount);
262 while (activeStateCount > 0) {
263 ASSERT(activeStateCount < _states.size());
264 State& currentState = _states[activeStateCount - 1];
265 State& nextState = _states[activeStateCount];
266 if (doEnumerationStep(ideal, tree, currentState, nextState))
267 ++activeStateCount;
268 else {
269 doEnumerationBaseCase(currentState, consumer);
270 --activeStateCount;
271 }
272 }
273}
virtual void consume(const Polynomial &poly)
bool getPrintStatistics() const
Returns whether to print statistics on what the algorithm did to standard error after it has run.
bool getPrintDebug() const
Returns whether to print information about what the algorithm is doing to standard error as it runs.
Objects of this class encapsulate the process of applying a generic deformation to a monomial ideal.
Definition Deformer.h:31
void undeform(Term &term) const
Apply the reverse transformation on term than that applied to the Ideal passed to the constructor.
Definition Deformer.cpp:99
A sparse multivariate polynomial represented by a hash table mapping terms to coefficients.
void add(const mpz_class &coef, const Term &term)
Add coef*term to the polynomial.
void feedTo(const TermTranslator &translator, CoefBigTermConsumer &consumer, bool inCanonicalOrder) const
Objects of this class represents a monomial ideal.
Definition IdealTree.h:29
bool strictlyContains(const Exponent *term) const
Represents a monomial ideal with int exponents.
Definition Ideal.h:27
size_t getGeneratorCount() const
Definition Ideal.h:57
bool containsIdentity() const
Definition Ideal.cpp:65
bool strictlyContains(const Exponent *term) const
Definition Ideal.cpp:73
const_iterator end() const
Definition Ideal.h:49
void print(FILE *file) const
Definition Ideal.cpp:440
const_iterator begin() const
Definition Ideal.h:48
size_t getVarCount() const
Definition Ideal.h:56
void enumerateScarfComplex(const Ideal &ideal, CoefTermConsumer &consumer)
void doEnumerationBaseCase(const State &state, CoefTermConsumer &consumer)
const ScarfParams & _params
bool doEnumerationStep(const Ideal &ideal, const IdealTree &tree, State &state, State &nextState)
void runGeneric(const Ideal &ideal, CoefBigTermConsumer &consumer, bool univariate, bool canonical)
const TermTranslator & _translator
const auto_ptr< IdealOrderer > _deformationOrder
void initializeEnumeration(const Ideal &ideal, size_t &activeStateCount)
ScarfHilbertAlgorithm(const TermTranslator &translator, const ScarfParams &params, auto_ptr< IdealOrderer > enumerationOrder, auto_ptr< IdealOrderer > deformationOrder)
const auto_ptr< IdealOrderer > _enumerationOrder
bool getDeformToStronglyGeneric() const
Returns true if deforming to a strongly generic ideal.
Definition ScarfParams.h:33
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.
size_t getVarCount() const
const VarNames & getNames() const
Term represents a product of variables which does not include a coefficient.
Definition Term.h:49
static bool strictlyDivides(const Exponent *a, const Exponent *b, size_t varCount)
Returns whether a strictly divides b.
Definition Term.h:196
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
UniHashPolynomial _uniPoly
CoefBigTermConsumer & _consumer
virtual void consumeRing(const VarNames &names)
virtual void consume(const mpz_class &coef, const Term &term)
UndeformConsumer(Ideal &toDeform, const TermTranslator &translator, CoefBigTermConsumer &consumer, const IdealOrderer &order, bool univar, bool canonical, bool doStrongDeformation)
virtual void beginConsuming()
const TermTranslator & _translator
A sparse univariate polynomial represented by a hash table mapping terms to coefficients.
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,...
Defines the variables of a polynomial ring and facilities IO involving them.
Definition VarNames.h:40
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