Frobby 0.9.5
PivotEulerAlg.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"
19#include "PivotEulerAlg.h"
20
21#include "Ideal.h"
22#include "RawSquareFreeTerm.h"
23#include "RawSquareFreeIdeal.h"
24#include "Task.h"
25#include "TaskEngine.h"
26#include "ElementDeleter.h"
27#include "EulerState.h"
28#include "PivotStrategy.h"
29#include "Arena.h"
30#include "LocalArray.h"
31
32#include <sstream>
33#include <vector>
34
35namespace Ops = SquareFreeTermOps;
36
37//typedef vector<size_t> DivCounts;
38typedef size_t* DivCounts;
39
40namespace {
41 bool baseCaseSimple1(mpz_class& accumulator,
42 const EulerState& state) {
43 const size_t varCount = state.getVarCount();
44 const RawSquareFreeIdeal& ideal = state.getIdeal();
45 const Word* eliminated = state.getEliminatedVars();
46 const size_t genCount = ideal.getGeneratorCount();
47
48 if (!ideal.hasFullSupport(eliminated))
49 return true;
50 if (genCount > 2)
51 return false;
52
53 if (genCount == 0)
54 accumulator += state.getSign();
55 else if (genCount == 2)
56 accumulator += state.getSign();
57 else {
58 ASSERT(genCount == 1);
59 if (!Ops::hasFullSupport(eliminated, varCount))
60 accumulator -= state.getSign();
61 }
62 return true;
63 }
64
65 bool optimizeSimpleFromDivCounts(mpz_class& accumulator,
66 EulerState& state,
67 DivCounts& divCounts,
68 Word* termTmp) {
69 const size_t varCount = state.getVarCount();
70 const size_t genCount = state.getIdeal().getGeneratorCount();
71 ASSERT(genCount > 2);
72
73 for (size_t var = 0; var < varCount; ++var) {
74 ASSERT(genCount == state.getIdeal().getGeneratorCount());
75 ASSERT(genCount > 2);
76
77 if (divCounts[var] < genCount - 2)
78 continue;
79
80 if (divCounts[var] == genCount - 1) {
81 const Word* nonMultiple =
82 state.getIdeal().getGenerator(state.getIdeal().getNonMultiple(var));
83 Ops::lcm(termTmp, nonMultiple, state.getEliminatedVars(), varCount);
84 Ops::setExponent(termTmp, var, 1);
85 if (Ops::hasFullSupport(termTmp, varCount))
86 accumulator += state.getSign();
87
88 if (state.toColonSubState(var))
89 return true;
90 divCounts[var] = 0;
91 } else if (divCounts[var] == genCount - 2) {
92 state.getIdeal().getLcmOfNonMultiples(termTmp, var);
93 Ops::lcmInPlace(termTmp, state.getEliminatedVars(), varCount);
94 Ops::setExponent(termTmp, var, 1);
95 if (Ops::hasFullSupport(termTmp, varCount))
96 accumulator -= state.getSign();
97
98 if (state.toColonSubState(var))
99 return true;
100 divCounts[var] = 0;
101 } else {
102 ASSERT(divCounts[var] == genCount);
103
105 divCounts[var] = 0;
106 }
107 }
108 return false;
109 }
110
111 bool baseCaseSimple2(mpz_class& accumulator,
112 const EulerState& state,
113 const DivCounts& divCounts) {
114 const size_t varCount = state.getVarCount();
115 const RawSquareFreeIdeal& ideal = state.getIdeal();
116 const size_t genCount = state.getIdeal().getGeneratorCount();
117
118 for (size_t var = 0; var < varCount; ++var)
119 if (divCounts[var] != 1 && divCounts[var] != genCount)
120 return false;
121
122 if ((ideal.getGeneratorCount() & 1) == 0)
123 accumulator += state.getSign();
124 else
125 accumulator -= state.getSign();
126 return true;
127 }
128
129 bool baseCasePreconditionSimplified(mpz_class& accumulator,
130 const EulerState& state) {
131 const RawSquareFreeIdeal& ideal = state.getIdeal();
132
133 if (ideal.getGeneratorCount() == 3) {
134 accumulator += state.getSign() + state.getSign();
135 return true;
136 }
137 return false;
138 }
139
140 bool optimizeOneDivCounts(EulerState& state,
141 DivCounts& divCounts,
142 Word* tmp) {
143 const size_t varCount = state.getVarCount();
144 const RawSquareFreeIdeal& ideal = state.getIdeal();
145
146 size_t var = 0;
147 for (; var < varCount; ++var) {
148 if (divCounts[var] != 1)
149 continue;
150 size_t index = ideal.getMultiple(var);
151 ASSERT(ideal.getGeneratorCount() > index);
152 Ops::assign(tmp, ideal.getGenerator(index), varCount);
153 state.removeGenerator(index);
154 state.flipSign();
155 goto searchForMore;
156 }
157 return false;
158
159searchForMore:
160 for (++var; var < varCount; ++var) {
161 if (divCounts[var] != 1 || Ops::getExponent(tmp, var) == 1)
162 continue;
163 size_t index = ideal.getMultiple(var);
164 ASSERT(ideal.getGeneratorCount() > index);
165 Ops::lcmInPlace(tmp, ideal.getGenerator(index), varCount);
166
167 state.removeGenerator(index);
168 state.flipSign();
169 }
170
171 if (state.toColonSubState(tmp) || ideal.getGeneratorCount() <= 2)
172 return true;
173
174 Ops::toZeroAtSupport(tmp, &(divCounts[0]), varCount);
175 return false;
176 }
177
178 bool optimizeVarPairs(EulerState& state, Word* tmp, DivCounts& divCounts) {
179 const size_t varCount = state.getVarCount();
180 const RawSquareFreeIdeal& ideal = state.getIdeal();
181 const Word* eliminated = state.getEliminatedVars();
182
183 for (size_t var = 0; var < varCount; ++var) {
184 if (Ops::getExponent(eliminated, var) == 1)
185 continue;
186 ideal.getLcmOfNonMultiples(tmp, var);
187 Ops::lcmInPlace(tmp, state.getEliminatedVars(), varCount);
188 Ops::setExponent(tmp, var, true);
189 if (!Ops::hasFullSupport(tmp, varCount)) {
190 if (state.toColonSubState(var))
191 return true;
192 divCounts[var] = 0;
193 }
194 }
195 return false;
196 }
197}
198
201
202 // ** First optimize state and return false if a base case is detected.
203 while (true) {
204 ASSERT(state.debugIsValid());
205
206 if (baseCaseSimple1(_euler, state))
207 return 0;
208
210 size_t* divCountsTmp = &(_divCountsTmp[0]);
211
213 optimizeOneDivCounts(state, divCountsTmp, _termTmp))
214 continue;
216 optimizeSimpleFromDivCounts(_euler, state, divCountsTmp, _termTmp))
217 continue;
219 if (optimizeVarPairs(state, _termTmp, divCountsTmp))
220 continue;
221 if (baseCasePreconditionSimplified(_euler, state))
222 return 0;
223 }
224 if (_autoTranspose && autoTranspose(state))
225 continue;
226 break;
227 }
228
229 // ** State is not a base case so perform a split while putting the
230 // two sub-states into state and newState.
231
232 size_t* divCountsTmp = &(_divCountsTmp[0]);
233 ASSERT(_pivotStrategy.get() != 0);
234 EulerState* next = _pivotStrategy->doPivot(state, divCountsTmp);
235
236 return next;
237}
238
239void PivotEulerAlg::getPivot(const EulerState& state, Word* pivot) {
240 ASSERT(false);
241}
242
244 _euler(0),
245 _termTmp(0),
246 _useUniqueDivSimplify(true),
247 _useManyDivSimplify(true),
248 _useAllPairsSimplify(false),
249 _autoTranspose(true),
250 _initialAutoTranspose(true) {
251}
252
253const mpz_class& PivotEulerAlg::computeEulerCharacteristic(const Ideal& ideal) {
254 if (_pivotStrategy.get() == 0)
256
257 if (ideal.getGeneratorCount() == 0)
258 _euler = 0;
259 else if (ideal.getVarCount() == 0)
260 _euler = -1;
261 else {
262 const size_t maxDim = std::max(ideal.getVarCount(), ideal.getGeneratorCount());
263 LocalArray<Word> termTmp(Ops::getWordCount(maxDim));
264 _termTmp = termTmp.begin();
266 computeEuler(state);
267 _termTmp = 0;
268 }
269 _pivotStrategy->computationCompleted(*this);
270
271 return _euler;
272}
273
275(const RawSquareFreeIdeal& ideal) {
276 if (_pivotStrategy.get() == 0)
278
279
280 if (ideal.getGeneratorCount() == 0)
281 _euler = 0;
282 else if (ideal.getVarCount() == 0)
283 _euler = -1;
284 else {
285 const size_t maxDim = std::max(ideal.getVarCount(), ideal.getGeneratorCount());
286 LocalArray<Word> termTmp(Ops::getWordCount(maxDim));
287 _termTmp = termTmp.begin();
289 computeEuler(state);
290 _termTmp = 0;
291 }
292 _pivotStrategy->computationCompleted(*this);
293
294 return _euler;
295}
296
298 _euler = 0;
300 autoTranspose(*state);
301 while (state != 0) {
302 EulerState* nextState = processState(*state);
303 if (nextState == 0) {
304 nextState = state->getParent();
306 }
307 state = nextState;
308 }
309}
310
312 if (!_pivotStrategy->shouldTranspose(state))
313 return false;
314 state.transpose();
315 return true;
316}
size_t * DivCounts
auto_ptr< PivotStrategy > newDefaultPivotStrategy()
static Arena & getArena()
Returns an arena object that can be used for non-thread safe scratch memory after static objects have...
Definition: Arena.h:126
void freeAndAllAfter(void *ptr)
Frees the buffer pointed to by ptr and all not yet freed allocations that have happened since that bu...
Definition: Arena.h:224
void removeGenerator(size_t index)
Definition: EulerState.h:55
void flipSign()
Definition: EulerState.h:43
void toColonSubStateNoReminimizeNecessary(size_t pivotVar)
Definition: EulerState.cpp:140
void compactEliminatedVariablesIfProfitable()
Definition: EulerState.cpp:196
RawSquareFreeIdeal & getIdeal()
Definition: EulerState.h:49
bool toColonSubState(const Word *pivot)
Definition: EulerState.cpp:118
int getSign() const
Definition: EulerState.h:44
static EulerState * construct(const Ideal &idealParam, Arena *arena)
Definition: EulerState.cpp:28
void transpose()
Definition: EulerState.cpp:190
size_t getVarCount() const
Definition: EulerState.h:52
const Word * getEliminatedVars() const
Definition: EulerState.h:51
EulerState * getParent()
Definition: EulerState.h:48
Represents a monomial ideal with int exponents.
Definition: Ideal.h:27
size_t getGeneratorCount() const
Definition: Ideal.h:57
size_t getVarCount() const
Definition: Ideal.h:56
Emulates stack allocation of an array using an Arena.
Definition: LocalArray.h:36
T * begin() const
Definition: LocalArray.h:54
bool _useAllPairsSimplify
Definition: PivotEulerAlg.h:70
bool _useManyDivSimplify
Definition: PivotEulerAlg.h:69
EulerState * processState(EulerState &state)
bool _useUniqueDivSimplify
Definition: PivotEulerAlg.h:68
auto_ptr< PivotStrategy > _pivotStrategy
Definition: PivotEulerAlg.h:73
const mpz_class & computeEulerCharacteristic(const Ideal &ideal)
mpz_class _euler
Definition: PivotEulerAlg.h:64
vector< size_t > _divCountsTmp
Definition: PivotEulerAlg.h:66
bool _initialAutoTranspose
Definition: PivotEulerAlg.h:72
bool autoTranspose(EulerState &state)
void computeEuler(EulerState *state)
void getPivot(const EulerState &state, Word *pivot)
A bit packed square free ideal placed in a pre-allocated buffer.
bool hasFullSupport(const Word *ignore) const
Returns true if for every variable it either divides ignore or it divides some (not necessarily minim...
void getVarDividesCounts(vector< size_t > &counts) const
Sets counts[var] to the number of generators that var divides.
size_t getVarCount() const
Word * getGenerator(size_t index)
Returns the generator at index.
size_t getNonMultiple(size_t var) const
Returns the index of the first generator that var does not divide or getGeneratorCount() if no such g...
size_t getMultiple(size_t var) const
Returns the index of the first generator that var divides or getGeneratorCount() if no such generator...
size_t getGeneratorCount() const
void getLcmOfNonMultiples(Word *lcm, size_t var) const
Sets lcm to be the least common multple of those generators that var does not divide.
void setExponent(Word *a, size_t var, bool value)
size_t getWordCount(size_t varCount)
bool hasFullSupport(const Word *a, size_t varCount)
void toZeroAtSupport(const Word *a, size_t *inc, size_t varCount)
For every variable var that divides a, set inc[var] to zero.
void lcmInPlace(Word *res, const Word *resEnd, const Word *a)
bool getExponent(const Word *a, size_t var)
returns true if var divides a and false otherwise.
void lcm(Word *res, const Word *resEnd, const Word *a, const Word *b)
void assign(Word *a, const Word *b, size_t varCount)
unsigned long Word
The native unsigned type for the CPU.
Definition: stdinc.h:93
#define ASSERT(X)
Definition: stdinc.h:86