Frobby  0.9.1
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 
35 namespace Ops = SquareFreeTermOps;
36 
37 //typedef vector<size_t> DivCounts;
38 typedef size_t* DivCounts;
39 
40 namespace {
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 
159 searchForMore:
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 
212  if (_useUniqueDivSimplify &&
213  optimizeOneDivCounts(state, divCountsTmp, _termTmp))
214  continue;
215  if (_useManyDivSimplify &&
216  optimizeSimpleFromDivCounts(_euler, state, divCountsTmp, _termTmp))
217  continue;
218  if (_useAllPairsSimplify) {
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 
239 void 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 
253 const 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();
265  EulerState* state = EulerState::construct(ideal, &(Arena::getArena()));
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)
277  _pivotStrategy = newDefaultPivotStrategy();
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();
288  EulerState* state = EulerState::construct(ideal, &(Arena::getArena()));
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 }
PivotEulerAlg::_useAllPairsSimplify
bool _useAllPairsSimplify
Definition: PivotEulerAlg.h:70
RawSquareFreeIdeal::getGeneratorCount
size_t getGeneratorCount() const
Definition: RawSquareFreeIdeal.h:124
stdinc.h
EulerState::getSign
int getSign() const
Definition: EulerState.h:44
PivotEulerAlg::PivotEulerAlg
PivotEulerAlg()
Definition: PivotEulerAlg.cpp:243
RawSquareFreeIdeal::getVarCount
size_t getVarCount() const
Definition: RawSquareFreeIdeal.h:125
SquareFreeTermOps::getWordCount
size_t getWordCount(size_t varCount)
Definition: RawSquareFreeTerm.cpp:91
Task.h
EulerState::construct
static EulerState * construct(const Ideal &idealParam, Arena *arena)
Definition: EulerState.cpp:28
DivCounts
size_t * DivCounts
Definition: PivotEulerAlg.cpp:38
PivotEulerAlg::_autoTranspose
bool _autoTranspose
Definition: PivotEulerAlg.h:71
EulerState::toColonSubStateNoReminimizeNecessary
void toColonSubStateNoReminimizeNecessary(size_t pivotVar)
Definition: EulerState.cpp:140
newDefaultPivotStrategy
auto_ptr< PivotStrategy > newDefaultPivotStrategy()
Definition: PivotStrategy.cpp:959
SquareFreeTermOps::assign
void assign(Word *a, const Word *b, size_t varCount)
Definition: RawSquareFreeTerm.cpp:186
LocalArray::begin
T * begin() const
Definition: LocalArray.h:54
EulerState::getIdeal
RawSquareFreeIdeal & getIdeal()
Definition: EulerState.h:49
RawSquareFreeIdeal::getLcmOfNonMultiples
void getLcmOfNonMultiples(Word *lcm, size_t var) const
Sets lcm to be the least common multple of those generators that var does not divide.
Definition: RawSquareFreeIdeal.cpp:304
RawSquareFreeIdeal.h
RawSquareFreeTerm.h
PivotEulerAlg::autoTranspose
bool autoTranspose(EulerState &state)
Definition: PivotEulerAlg.cpp:311
SquareFreeTermOps
Definition: RawSquareFreeTerm.cpp:24
EulerState::getParent
EulerState * getParent()
Definition: EulerState.h:48
EulerState::transpose
void transpose()
Definition: EulerState.cpp:190
PivotEulerAlg::getPivot
void getPivot(const EulerState &state, Word *pivot)
Definition: PivotEulerAlg.cpp:239
EulerState::toColonSubState
bool toColonSubState(const Word *pivot)
Definition: EulerState.cpp:118
LocalArray
Emulates stack allocation of an array using an Arena.
Definition: LocalArray.h:36
Ideal::getGeneratorCount
size_t getGeneratorCount() const
Definition: Ideal.h:57
Arena::getArena
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
SquareFreeTermOps::setExponent
void setExponent(Word *a, size_t var, bool value)
Definition: RawSquareFreeTerm.h:170
PivotEulerAlg::processState
EulerState * processState(EulerState &state)
Definition: PivotEulerAlg.cpp:199
Ideal.h
RawSquareFreeIdeal::hasFullSupport
bool hasFullSupport(const Word *ignore) const
Returns true if for every variable it either divides ignore or it divides some (not necessarily minim...
Definition: RawSquareFreeIdeal.cpp:545
EulerState::flipSign
void flipSign()
Definition: EulerState.h:43
SquareFreeTermOps::lcmInPlace
void lcmInPlace(Word *res, const Word *resEnd, const Word *a)
Definition: RawSquareFreeTerm.cpp:264
EulerState::compactEliminatedVariablesIfProfitable
void compactEliminatedVariablesIfProfitable()
Definition: EulerState.cpp:196
TaskEngine.h
SquareFreeTermOps::toZeroAtSupport
void toZeroAtSupport(const Word *a, size_t *inc, size_t varCount)
For every variable var that divides a, set inc[var] to zero.
Definition: RawSquareFreeTerm.cpp:376
PivotEulerAlg::_euler
mpz_class _euler
Definition: PivotEulerAlg.h:64
PivotEulerAlg.h
Ideal::getVarCount
size_t getVarCount() const
Definition: Ideal.h:56
RawSquareFreeIdeal
A bit packed square free ideal placed in a pre-allocated buffer.
Definition: RawSquareFreeIdeal.h:37
PivotEulerAlg::computeEulerCharacteristic
const mpz_class & computeEulerCharacteristic(const Ideal &ideal)
Definition: PivotEulerAlg.cpp:253
EulerState
Definition: EulerState.h:25
PivotEulerAlg::_initialAutoTranspose
bool _initialAutoTranspose
Definition: PivotEulerAlg.h:72
PivotEulerAlg::_useUniqueDivSimplify
bool _useUniqueDivSimplify
Definition: PivotEulerAlg.h:68
RawSquareFreeIdeal::getMultiple
size_t getMultiple(size_t var) const
Returns the index of the first generator that var divides or getGeneratorCount() if no such generator...
Definition: RawSquareFreeIdeal.cpp:393
RawSquareFreeIdeal::getGenerator
Word * getGenerator(size_t index)
Returns the generator at index.
Definition: RawSquareFreeIdeal.h:350
Word
unsigned long Word
The native unsigned type for the CPU.
Definition: stdinc.h:93
RawSquareFreeIdeal::getVarDividesCounts
void getVarDividesCounts(vector< size_t > &counts) const
Sets counts[var] to the number of generators that var divides.
Definition: RawSquareFreeIdeal.cpp:366
SquareFreeTermOps::lcm
void lcm(Word *res, const Word *resEnd, const Word *a, const Word *b)
Definition: RawSquareFreeTerm.cpp:251
SquareFreeTermOps::getExponent
bool getExponent(const Word *a, size_t var)
returns true if var divides a and false otherwise.
Definition: RawSquareFreeTerm.h:151
PivotEulerAlg::_useManyDivSimplify
bool _useManyDivSimplify
Definition: PivotEulerAlg.h:69
PivotEulerAlg::_pivotStrategy
auto_ptr< PivotStrategy > _pivotStrategy
Definition: PivotEulerAlg.h:73
EulerState.h
PivotStrategy.h
Arena.h
EulerState::getVarCount
size_t getVarCount() const
Definition: EulerState.h:52
LocalArray.h
Ideal
Represents a monomial ideal with int exponents.
Definition: Ideal.h:27
PivotEulerAlg::computeEuler
void computeEuler(EulerState *state)
Definition: PivotEulerAlg.cpp:297
PivotEulerAlg::_termTmp
Word * _termTmp
Definition: PivotEulerAlg.h:65
PivotEulerAlg::_divCountsTmp
vector< size_t > _divCountsTmp
Definition: PivotEulerAlg.h:66
SquareFreeTermOps::hasFullSupport
bool hasFullSupport(const Word *a, size_t varCount)
Definition: RawSquareFreeTerm.h:183
EulerState::removeGenerator
void removeGenerator(size_t index)
Definition: EulerState.h:55
ASSERT
#define ASSERT(X)
Definition: stdinc.h:86
RawSquareFreeIdeal::getNonMultiple
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...
Definition: RawSquareFreeIdeal.cpp:404
Arena::freeAndAllAfter
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
EulerState::getEliminatedVars
const Word * getEliminatedVars() const
Definition: EulerState.h:51
ElementDeleter.h