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)
void compactEliminatedVariablesIfProfitable()
RawSquareFreeIdeal & getIdeal()
Definition EulerState.h:49
bool toColonSubState(const Word *pivot)
int getSign() const
Definition EulerState.h:44
static EulerState * construct(const Ideal &idealParam, Arena *arena)
void transpose()
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
bool _useManyDivSimplify
EulerState * processState(EulerState &state)
bool _useUniqueDivSimplify
auto_ptr< PivotStrategy > _pivotStrategy
const mpz_class & computeEulerCharacteristic(const Ideal &ideal)
mpz_class _euler
vector< size_t > _divCountsTmp
bool _initialAutoTranspose
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)
This header file includes common definitions and is included as the first line of code in every imple...
unsigned long Word
The native unsigned type for the CPU.
Definition stdinc.h:93
#define ASSERT(X)
Definition stdinc.h:86