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
73 virtual void doneConsuming() {
74 if (_univar)
76 else
78 }
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.
Definition: CommonParams.h:66
bool getPrintDebug() const
Returns whether to print information about what the algorithm is doing to standard error as it runs.
Definition: CommonParams.h:61
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
Definition: IdealTree.cpp:161
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
virtual void doneConsuming()
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
unsigned int Exponent
Definition: stdinc.h:89
#define ASSERT(X)
Definition: stdinc.h:86