PivotEulerAlg.cpp
Go to the documentation of this file.
/* Frobby: Software for monomial ideal computations.
Copyright (C) 2010 University of Aarhus
Contact Bjarke Hammersholt Roune for license information (www.broune.com)
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see http://www.gnu.org/licenses/.
*/
#include "stdinc.h"
#include "PivotEulerAlg.h"
#include "Ideal.h"
#include "RawSquareFreeTerm.h"
#include "RawSquareFreeIdeal.h"
#include "Task.h"
#include "TaskEngine.h"
#include "ElementDeleter.h"
#include "EulerState.h"
#include "PivotStrategy.h"
#include "Arena.h"
#include "LocalArray.h"
#include <sstream>
#include <vector>
namespace Ops = SquareFreeTermOps;
//typedef vector<size_t> DivCounts;
typedef size_t* DivCounts;
namespace {
bool baseCaseSimple1(mpz_class& accumulator,
const EulerState& state) {
const size_t varCount = state.getVarCount();
const RawSquareFreeIdeal& ideal = state.getIdeal();
const Word* eliminated = state.getEliminatedVars();
const size_t genCount = ideal.getGeneratorCount();
if (!ideal.hasFullSupport(eliminated))
return true;
if (genCount > 2)
return false;
if (genCount == 0)
accumulator += state.getSign();
else if (genCount == 2)
accumulator += state.getSign();
else {
ASSERT(genCount == 1);
if (!Ops::hasFullSupport(eliminated, varCount))
accumulator -= state.getSign();
}
return true;
}
bool optimizeSimpleFromDivCounts(mpz_class& accumulator,
EulerState& state,
DivCounts& divCounts,
Word* termTmp) {
const size_t varCount = state.getVarCount();
const size_t genCount = state.getIdeal().getGeneratorCount();
ASSERT(genCount > 2);
for (size_t var = 0; var < varCount; ++var) {
ASSERT(genCount == state.getIdeal().getGeneratorCount());
ASSERT(genCount > 2);
if (divCounts[var] < genCount - 2)
continue;
if (divCounts[var] == genCount - 1) {
const Word* nonMultiple =
state.getIdeal().getGenerator(state.getIdeal().getNonMultiple(var));
Ops::lcm(termTmp, nonMultiple, state.getEliminatedVars(), varCount);
Ops::setExponent(termTmp, var, 1);
if (Ops::hasFullSupport(termTmp, varCount))
accumulator += state.getSign();
if (state.toColonSubState(var))
return true;
divCounts[var] = 0;
} else if (divCounts[var] == genCount - 2) {
state.getIdeal().getLcmOfNonMultiples(termTmp, var);
Ops::lcmInPlace(termTmp, state.getEliminatedVars(), varCount);
Ops::setExponent(termTmp, var, 1);
if (Ops::hasFullSupport(termTmp, varCount))
accumulator -= state.getSign();
if (state.toColonSubState(var))
return true;
divCounts[var] = 0;
} else {
ASSERT(divCounts[var] == genCount);
state.toColonSubStateNoReminimizeNecessary(var);
divCounts[var] = 0;
}
}
return false;
}
bool baseCaseSimple2(mpz_class& accumulator,
const EulerState& state,
const DivCounts& divCounts) {
const size_t varCount = state.getVarCount();
const RawSquareFreeIdeal& ideal = state.getIdeal();
const size_t genCount = state.getIdeal().getGeneratorCount();
for (size_t var = 0; var < varCount; ++var)
if (divCounts[var] != 1 && divCounts[var] != genCount)
return false;
if ((ideal.getGeneratorCount() & 1) == 0)
accumulator += state.getSign();
else
accumulator -= state.getSign();
return true;
}
bool baseCasePreconditionSimplified(mpz_class& accumulator,
const EulerState& state) {
const RawSquareFreeIdeal& ideal = state.getIdeal();
if (ideal.getGeneratorCount() == 3) {
accumulator += state.getSign() + state.getSign();
return true;
}
return false;
}
bool optimizeOneDivCounts(EulerState& state,
DivCounts& divCounts,
Word* tmp) {
const size_t varCount = state.getVarCount();
const RawSquareFreeIdeal& ideal = state.getIdeal();
size_t var = 0;
for (; var < varCount; ++var) {
if (divCounts[var] != 1)
continue;
size_t index = ideal.getMultiple(var);
ASSERT(ideal.getGeneratorCount() > index);
Ops::assign(tmp, ideal.getGenerator(index), varCount);
state.removeGenerator(index);
state.flipSign();
goto searchForMore;
}
return false;
searchForMore:
for (++var; var < varCount; ++var) {
if (divCounts[var] != 1 || Ops::getExponent(tmp, var) == 1)
continue;
size_t index = ideal.getMultiple(var);
ASSERT(ideal.getGeneratorCount() > index);
Ops::lcmInPlace(tmp, ideal.getGenerator(index), varCount);
state.removeGenerator(index);
state.flipSign();
}
if (state.toColonSubState(tmp) || ideal.getGeneratorCount() <= 2)
return true;
Ops::toZeroAtSupport(tmp, &(divCounts[0]), varCount);
return false;
}
bool optimizeVarPairs(EulerState& state, Word* tmp, DivCounts& divCounts) {
const size_t varCount = state.getVarCount();
const RawSquareFreeIdeal& ideal = state.getIdeal();
const Word* eliminated = state.getEliminatedVars();
for (size_t var = 0; var < varCount; ++var) {
if (Ops::getExponent(eliminated, var) == 1)
continue;
ideal.getLcmOfNonMultiples(tmp, var);
Ops::lcmInPlace(tmp, state.getEliminatedVars(), varCount);
Ops::setExponent(tmp, var, true);
if (!Ops::hasFullSupport(tmp, varCount)) {
if (state.toColonSubState(var))
return true;
divCounts[var] = 0;
}
}
return false;
}
}
EulerState* PivotEulerAlg::processState(EulerState& state) {
state.compactEliminatedVariablesIfProfitable();
// ** First optimize state and return false if a base case is detected.
while (true) {
ASSERT(state.debugIsValid());
if (baseCaseSimple1(_euler, state))
return 0;
state.getIdeal().getVarDividesCounts(_divCountsTmp);
size_t* divCountsTmp = &(_divCountsTmp[0]);
if (_useUniqueDivSimplify &&
optimizeOneDivCounts(state, divCountsTmp, _termTmp))
continue;
if (_useManyDivSimplify &&
optimizeSimpleFromDivCounts(_euler, state, divCountsTmp, _termTmp))
continue;
if (_useAllPairsSimplify) {
if (optimizeVarPairs(state, _termTmp, divCountsTmp))
continue;
if (baseCasePreconditionSimplified(_euler, state))
return 0;
}
if (_autoTranspose && autoTranspose(state))
continue;
break;
}
// ** State is not a base case so perform a split while putting the
// two sub-states into state and newState.
size_t* divCountsTmp = &(_divCountsTmp[0]);
ASSERT(_pivotStrategy.get() != 0);
EulerState* next = _pivotStrategy->doPivot(state, divCountsTmp);
return next;
}
void PivotEulerAlg::getPivot(const EulerState& state, Word* pivot) {
ASSERT(false);
}
PivotEulerAlg::PivotEulerAlg():
_euler(0),
_termTmp(0),
_useUniqueDivSimplify(true),
_useManyDivSimplify(true),
_useAllPairsSimplify(false),
_autoTranspose(true),
_initialAutoTranspose(true) {
}
const mpz_class& PivotEulerAlg::computeEulerCharacteristic(const Ideal& ideal) {
if (_pivotStrategy.get() == 0)
_pivotStrategy = newDefaultPivotStrategy();
if (ideal.getGeneratorCount() == 0)
_euler = 0;
else if (ideal.getVarCount() == 0)
_euler = -1;
else {
const size_t maxDim = std::max(ideal.getVarCount(), ideal.getGeneratorCount());
LocalArray<Word> termTmp(Ops::getWordCount(maxDim));
_termTmp = termTmp.begin();
EulerState* state = EulerState::construct(ideal, &(Arena::getArena()));
computeEuler(state);
_termTmp = 0;
}
_pivotStrategy->computationCompleted(*this);
return _euler;
}
const mpz_class& PivotEulerAlg::computeEulerCharacteristic
(const RawSquareFreeIdeal& ideal) {
if (_pivotStrategy.get() == 0)
_pivotStrategy = newDefaultPivotStrategy();
if (ideal.getGeneratorCount() == 0)
_euler = 0;
else if (ideal.getVarCount() == 0)
_euler = -1;
else {
const size_t maxDim = std::max(ideal.getVarCount(), ideal.getGeneratorCount());
LocalArray<Word> termTmp(Ops::getWordCount(maxDim));
_termTmp = termTmp.begin();
EulerState* state = EulerState::construct(ideal, &(Arena::getArena()));
computeEuler(state);
_termTmp = 0;
}
_pivotStrategy->computationCompleted(*this);
return _euler;
}
void PivotEulerAlg::computeEuler(EulerState* state) {
_euler = 0;
if (_initialAutoTranspose)
autoTranspose(*state);
while (state != 0) {
EulerState* nextState = processState(*state);
if (nextState == 0) {
nextState = state->getParent();
Arena::getArena().freeAndAllAfter(state);
}
state = nextState;
}
}
bool PivotEulerAlg::autoTranspose(EulerState& state) {
if (!_pivotStrategy->shouldTranspose(state))
return false;
state.transpose();
return true;
}
Arena.h
ElementDeleter.h
EulerState.h
Ideal.h
LocalArray.h
DivCounts
size_t * DivCounts
Definition PivotEulerAlg.cpp:38
PivotEulerAlg.h
newDefaultPivotStrategy
unique_ptr< PivotStrategy > newDefaultPivotStrategy()
Definition PivotStrategy.cpp:959
PivotStrategy.h
RawSquareFreeIdeal.h
RawSquareFreeTerm.h
TaskEngine.h
Task.h
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
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
Definition EulerState.h:25
EulerState::removeGenerator
void removeGenerator(size_t index)
Definition EulerState.h:55
EulerState::flipSign
void flipSign()
Definition EulerState.h:43
EulerState::toColonSubStateNoReminimizeNecessary
void toColonSubStateNoReminimizeNecessary(size_t pivotVar)
Definition EulerState.cpp:140
EulerState::compactEliminatedVariablesIfProfitable
void compactEliminatedVariablesIfProfitable()
Definition EulerState.cpp:196
EulerState::getIdeal
RawSquareFreeIdeal & getIdeal()
Definition EulerState.h:49
EulerState::toColonSubState
bool toColonSubState(const Word *pivot)
Definition EulerState.cpp:118
EulerState::getSign
int getSign() const
Definition EulerState.h:44
EulerState::construct
static EulerState * construct(const Ideal &idealParam, Arena *arena)
Definition EulerState.cpp:28
EulerState::transpose
void transpose()
Definition EulerState.cpp:190
EulerState::getVarCount
size_t getVarCount() const
Definition EulerState.h:52
EulerState::getEliminatedVars
const Word * getEliminatedVars() const
Definition EulerState.h:51
EulerState::getParent
EulerState * getParent()
Definition EulerState.h:48
Ideal
Represents a monomial ideal with int exponents.
Definition Ideal.h:27
Ideal::getGeneratorCount
size_t getGeneratorCount() const
Definition Ideal.h:57
Ideal::getVarCount
size_t getVarCount() const
Definition Ideal.h:56
LocalArray
Emulates stack allocation of an array using an Arena.
Definition LocalArray.h:36
LocalArray::begin
T * begin() const
Definition LocalArray.h:54
PivotEulerAlg::_useAllPairsSimplify
bool _useAllPairsSimplify
Definition PivotEulerAlg.h:70
PivotEulerAlg::_useManyDivSimplify
bool _useManyDivSimplify
Definition PivotEulerAlg.h:69
PivotEulerAlg::PivotEulerAlg
PivotEulerAlg()
Definition PivotEulerAlg.cpp:243
PivotEulerAlg::processState
EulerState * processState(EulerState &state)
Definition PivotEulerAlg.cpp:199
PivotEulerAlg::_useUniqueDivSimplify
bool _useUniqueDivSimplify
Definition PivotEulerAlg.h:68
PivotEulerAlg::_autoTranspose
bool _autoTranspose
Definition PivotEulerAlg.h:71
PivotEulerAlg::computeEulerCharacteristic
const mpz_class & computeEulerCharacteristic(const Ideal &ideal)
Definition PivotEulerAlg.cpp:253
PivotEulerAlg::_euler
mpz_class _euler
Definition PivotEulerAlg.h:64
PivotEulerAlg::_divCountsTmp
vector< size_t > _divCountsTmp
Definition PivotEulerAlg.h:66
PivotEulerAlg::_termTmp
Word * _termTmp
Definition PivotEulerAlg.h:65
PivotEulerAlg::_initialAutoTranspose
bool _initialAutoTranspose
Definition PivotEulerAlg.h:72
PivotEulerAlg::autoTranspose
bool autoTranspose(EulerState &state)
Definition PivotEulerAlg.cpp:311
PivotEulerAlg::_pivotStrategy
unique_ptr< PivotStrategy > _pivotStrategy
Definition PivotEulerAlg.h:73
PivotEulerAlg::computeEuler
void computeEuler(EulerState *state)
Definition PivotEulerAlg.cpp:297
PivotEulerAlg::getPivot
void getPivot(const EulerState &state, Word *pivot)
Definition PivotEulerAlg.cpp:239
RawSquareFreeIdeal
A bit packed square free ideal placed in a pre-allocated buffer.
Definition RawSquareFreeIdeal.h:37
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
RawSquareFreeIdeal::getVarDividesCounts
void getVarDividesCounts(vector< size_t > &counts) const
Sets counts[var] to the number of generators that var divides.
Definition RawSquareFreeIdeal.cpp:366
RawSquareFreeIdeal::getVarCount
size_t getVarCount() const
Definition RawSquareFreeIdeal.h:125
RawSquareFreeIdeal::getGenerator
Word * getGenerator(size_t index)
Returns the generator at index.
Definition RawSquareFreeIdeal.h:350
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
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::getGeneratorCount
size_t getGeneratorCount() const
Definition RawSquareFreeIdeal.h:124
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
SquareFreeTermOps
Definition RawSquareFreeTerm.cpp:24
SquareFreeTermOps::setExponent
void setExponent(Word *a, size_t var, bool value)
Definition RawSquareFreeTerm.h:170
SquareFreeTermOps::getWordCount
size_t getWordCount(size_t varCount)
Definition RawSquareFreeTerm.cpp:91
SquareFreeTermOps::hasFullSupport
bool hasFullSupport(const Word *a, size_t varCount)
Definition RawSquareFreeTerm.h:183
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
SquareFreeTermOps::lcmInPlace
void lcmInPlace(Word *res, const Word *resEnd, const Word *a)
Definition RawSquareFreeTerm.cpp:264
SquareFreeTermOps::getExponent
bool getExponent(const Word *a, size_t var)
returns true if var divides a and false otherwise.
Definition RawSquareFreeTerm.h:151
SquareFreeTermOps::lcm
void lcm(Word *res, const Word *resEnd, const Word *a, const Word *b)
Definition RawSquareFreeTerm.cpp:251
SquareFreeTermOps::assign
void assign(Word *a, const Word *b, size_t varCount)
Definition RawSquareFreeTerm.cpp:186
stdinc.h
This header file includes common definitions and is included as the first line of code in every imple...
Word
unsigned long Word
The native unsigned type for the CPU.
Definition stdinc.h:93
ASSERT
#define ASSERT(X)
Definition stdinc.h:86