[mathicgb] 365/393: Clean-up on F4MatrixBuilder2, FixedSizeMonomialMap and MonomialMap, migrating them to use a monoid.
Doug Torrance
dtorrance-guest at moszumanska.debian.org
Fri Apr 3 15:59:35 UTC 2015
This is an automated email from the git hooks/post-receive script.
dtorrance-guest pushed a commit to branch upstream
in repository mathicgb.
commit 229b90237b41a0b7ce719985a278bbea2f1f5de0
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date: Mon Sep 16 13:52:40 2013 +0200
Clean-up on F4MatrixBuilder2, FixedSizeMonomialMap and MonomialMap, migrating them to use a monoid.
---
src/mathicgb/F4MatrixBuilder.cpp | 30 +-
src/mathicgb/F4MatrixBuilder.hpp | 14 +-
src/mathicgb/F4MatrixBuilder2.cpp | 910 ++++++++++++++++++------------------
src/mathicgb/F4MatrixBuilder2.hpp | 18 +-
src/mathicgb/FixedSizeMonomialMap.h | 130 ++++--
src/mathicgb/MonoMonoid.hpp | 74 ++-
src/mathicgb/MonomialMap.hpp | 37 +-
src/mathicgb/Poly.hpp | 6 +
src/mathicgb/QuadMatrixBuilder.cpp | 6 +-
9 files changed, 648 insertions(+), 577 deletions(-)
diff --git a/src/mathicgb/F4MatrixBuilder.cpp b/src/mathicgb/F4MatrixBuilder.cpp
index 46ebfe0..9052631 100755
--- a/src/mathicgb/F4MatrixBuilder.cpp
+++ b/src/mathicgb/F4MatrixBuilder.cpp
@@ -13,34 +13,32 @@ MATHICGB_DEFINE_LOG_DOMAIN(
MATHICGB_NAMESPACE_BEGIN
MATHICGB_NO_INLINE
-std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonomial>
-F4MatrixBuilder::findOrCreateColumn(
+auto F4MatrixBuilder::findOrCreateColumn(
const const_monomial monoA,
const const_monomial monoB,
TaskFeeder& feeder
-) {
+) -> std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonoRef> {
MATHICGB_ASSERT(!monoA.isNull());
MATHICGB_ASSERT(!monoB.isNull());
const auto col = ColReader(mMap).findProduct(monoA, monoB);
if (col.first != 0)
- return std::make_pair(*col.first, col.second);
+ return std::make_pair(*col.first, *col.second);
return createColumn(monoA, monoB, feeder);
}
MATHICGB_INLINE
-std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonomial>
-F4MatrixBuilder::findOrCreateColumn(
+auto F4MatrixBuilder::findOrCreateColumn(
const const_monomial monoA,
const const_monomial monoB,
const ColReader& colMap,
TaskFeeder& feeder
-) {
+) -> std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonoRef> {
MATHICGB_ASSERT(!monoA.isNull());
MATHICGB_ASSERT(!monoB.isNull());
const auto col = colMap.findProduct(monoA, monoB);
if (col.first == 0)
return findOrCreateColumn(monoA, monoB, feeder);
- return std::make_pair(*col.first, col.second);
+ return std::make_pair(*col.first, *col.second);
}
MATHICGB_NO_INLINE
@@ -206,7 +204,7 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
for (auto it = reader.begin(); it != end; ++it) {
const auto p = *it;
monomial copy = ring().allocMonomial();
- ring().monomialCopy(p.second, copy);
+ ring().monoid().copy(p.second, copy);
auto& monos = p.first.left() ?
matrix.leftColumnMonomials : matrix.rightColumnMonomials;
const auto index = p.first.index();
@@ -229,12 +227,11 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
mMap.clearNonConcurrent();
}
-std::pair<F4MatrixBuilder::LeftRightColIndex, ConstMonomial>
-F4MatrixBuilder::createColumn(
+auto F4MatrixBuilder::createColumn(
const const_monomial monoA,
const const_monomial monoB,
TaskFeeder& feeder
-) {
+) -> std::pair<F4MatrixBuilder::LeftRightColIndex, ConstMonoRef> {
MATHICGB_ASSERT(!monoA.isNull());
MATHICGB_ASSERT(!monoB.isNull());
@@ -243,7 +240,7 @@ F4MatrixBuilder::createColumn(
{
const auto found(ColReader(mMap).findProduct(monoA, monoB));
if (found.first != 0)
- return std::make_pair(*found.first, found.second);
+ return std::make_pair(*found.first, *found.second);
}
// The column really does not exist, so we need to create it
@@ -264,17 +261,18 @@ F4MatrixBuilder::createColumn(
++colCount;
MATHICGB_ASSERT(inserted.second);
MATHICGB_ASSERT(inserted.first.first != 0);
+ auto ptr = const_cast<exponent*>(Monoid::toOld(*inserted.first.second));
// schedule new task if we found a reducer
if (insertLeft) {
RowTask task = {};
task.addToTop = true;
task.poly = &mBasis.poly(reducerIndex);
- task.desiredLead = inserted.first.second.castAwayConst();
+ task.desiredLead = ptr;
feeder.add(task);
}
- return std::make_pair(*inserted.first.first, inserted.first.second);
+ return std::make_pair(*inserted.first.first, Monoid::toRef(ptr));
}
void F4MatrixBuilder::appendRowBottom(
@@ -410,7 +408,7 @@ void F4MatrixBuilder::appendRowBottom(
(itA.getMonomial(), mulA, colMap, feeder);
const auto colB = findOrCreateColumn
(itB.getMonomial(), mulB, colMap, feeder);
- const auto cmp = ring().monomialCompare(colA.second, colB.second);
+ const auto cmp = ring().monoid().compare(colA.second, colB.second);
if (cmp != LT) {
coeff = itA.getCoefficient();
col = colA.first;
diff --git a/src/mathicgb/F4MatrixBuilder.hpp b/src/mathicgb/F4MatrixBuilder.hpp
index 63a3750..6e6ac90 100755
--- a/src/mathicgb/F4MatrixBuilder.hpp
+++ b/src/mathicgb/F4MatrixBuilder.hpp
@@ -29,6 +29,14 @@ private:
typedef QuadMatrixBuilder::MonomialsType Monomials;
public:
+ typedef PolyRing::Monoid Monoid;
+ typedef Monoid::Mono Mono;
+ typedef Monoid::MonoRef MonoRef;
+ typedef Monoid::ConstMonoRef ConstMonoRef;
+ typedef Monoid::MonoPtr MonoPtr;
+ typedef Monoid::ConstMonoPtr ConstMonoPtr;
+
+
/// memoryQuantum is how much to increase the memory size by each time the
/// current amount of memory is exhausted. A value of 0 indicates to start
/// small and double the quantum at each exhaustion.
@@ -91,7 +99,7 @@ private:
/// reduce that column if possible. Here x is monoA if monoB is
/// null and otherwise x is the product of monoA and monoB.
MATHICGB_NO_INLINE
- std::pair<LeftRightColIndex, ConstMonomial>
+ std::pair<LeftRightColIndex, ConstMonoRef>
createColumn(
const_monomial monoA,
const_monomial monoB,
@@ -122,7 +130,7 @@ private:
);
MATHICGB_NO_INLINE
- std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonomial>
+ std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonoRef>
findOrCreateColumn(
const_monomial monoA,
const_monomial monoB,
@@ -130,7 +138,7 @@ private:
);
MATHICGB_INLINE
- std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonomial>
+ std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonoRef>
findOrCreateColumn(
const_monomial monoA,
const_monomial monoB,
diff --git a/src/mathicgb/F4MatrixBuilder2.cpp b/src/mathicgb/F4MatrixBuilder2.cpp
index 2758b45..55e3a41 100755
--- a/src/mathicgb/F4MatrixBuilder2.cpp
+++ b/src/mathicgb/F4MatrixBuilder2.cpp
@@ -20,10 +20,259 @@ MATHICGB_NAMESPACE_BEGIN
class F4MatrixBuilder2::Builder {
public:
- void buildMatrixAndClear( std::vector<RowTask>& mTodo, QuadMatrix& quadMatrix);
+ typedef SparseMatrix::ColIndex ColIndex;
+ typedef SparseMatrix::Scalar Scalar;
+ typedef MonomialMap<ColIndex> Map;
+ typedef SparseMatrix::RowIndex RowIndex;
+
+ /// Initializes the set of add-this-row tasks.
+ void initializeRowsToReduce(std::vector<RowTask>& tasks) {
+ // If aF-bG is an S-pair that is added as a bottom row in the matrix, and
+ // we do not already have a reducer for the common leading term of aF and
+ // bG, then we can instead set bG as a top/reducer/pivot row and put aF
+ // as a bottom row. This way, we do not need to add a reducer for that
+ // column if that monomial turns out to be used. We can also use the
+ // faster code path for adding a single multiple like aF to the matrix
+ // instead of the more complicated path that adds the difference between
+ // two multiples aF-bG.
+ //
+ // If bG is already the top/reducer/pivot row for its leading monomial,
+ // then we can simply replace a bottom row aF-bG with aF directly, since
+ // the first computation that will happen will be to calculate aF-bG.
+ // This way we get to use the much faster code in the matrix reducer for
+ // making that reduction/Gaussian-elimination step. More importantly, we
+ // have 1 row less to construct.
+ //
+ // The above two cases turn rows specified like aF-bG into rows specified
+ // more simply as just aF. Note that we can interchange aF and bG in these
+ // arguments with no problems.
+ //
+ // These two cases will apply to all
+ // S-pairs with a given leading monomial if and only if there is a single
+ // bG that appears in all of those S-pairs. The absolute best case is
+ // if that single common bG is at the same time the most preferable
+ // (sparsest or oldest) reducer for that monomial. It is not trivial to
+ // realize this (in fact I should write a paper it), but the particular
+ // way that S-pair elimination works in MathicGB ensures that this best
+ // case is always the case. So the code below will get rid of all the
+ // S-pairs aF-bG and replace them with single bottom rows aF.
+ // However, if the S-pair elimination code is ever changed, that may
+ // no longer be the case, so this code is still supposed to work even
+ // if some S-pairs are not gotten rid of. Still, there is an assert
+ // to flag it if the S-pair elimination is accidentally broken so that it
+ // no longer gives this guarantee.
+
+
+ // Bring S-pairs with the same leading monomial together by ordering
+ // them according in increasing order of monomial. Put the non-S-pairs
+ // together at the front.
+ auto cmp = [&](const RowTask& a, const RowTask& b) {
+ if (a.sPairPoly == nullptr)
+ return b.sPairPoly != nullptr;
+ else
+ return monoid().lessThan(*a.desiredLead, *b.desiredLead);
+ };
+ mgb::mtbb::parallel_sort(tasks.begin(), tasks.end(), cmp);
+
+ const auto taskCount = tasks.size();
+ for (size_t i = 0; i < taskCount;) {
+ if (tasks[i].sPairPoly == 0) {
+ ++i;
+ continue;
+ }
+ MATHICGB_ASSERT(tasks[i].desiredLead != nullptr);
+ MATHICGB_ASSERT(ColReader(mMap).find(*tasks[i].desiredLead).first == 0);
+
+ // Create column for the lead term that cancels in the S-pair
+ if (mIsColumnToLeft.size() >= std::numeric_limits<ColIndex>::max())
+ throw std::overflow_error("Too many columns in matrix.");
+ const auto newIndex = static_cast<ColIndex>(mIsColumnToLeft.size());
+ const auto inserted =
+ mMap.insert(std::make_pair(tasks[i].desiredLead, newIndex));
+ mIsColumnToLeft.push_back(true);
+ const auto& mono = inserted.first.second;
+
+ // Schedule the two parts of the S-pair as separate rows. This adds a row
+ // while creating the column in the hash table without adding a reducer
+ // removes a row, effectively making this operation equivalent to taking
+ // half of the S-pair and using it as a reducer.
+ auto desiredLead = monoid().alloc();
+ monoid().copy(*mono, *desiredLead);
+ RowTask newTask = {desiredLead.release(), tasks[i].sPairPoly, nullptr};
+
+ // Now we can strip off any part of an S-pair with the same cancelling lead
+ // term that equals a or b since those rows are in the matrix.
+ const auto* const a = tasks[i].poly;
+ const auto* const b = tasks[i].sPairPoly;
+
+ tasks[i].sPairPoly = 0;
+
+ tasks.push_back(newTask);
+ for (++i; i < taskCount; ++i) {
+ auto& task = tasks[i];
+ if (tasks[i].sPairPoly == 0)
+ continue;
+ MATHICGB_ASSERT(tasks[i].desiredLead != nullptr);
+ if (!monoid().equal(*tasks[i].desiredLead, *mono))
+ break;
+ if (tasks[i].poly == a || tasks[i].poly == b) {
+ tasks[i].poly = tasks[i].sPairPoly;
+ tasks[i].sPairPoly = 0;
+ } else if (tasks[i].sPairPoly == a || tasks[i].sPairPoly == b) {
+ tasks[i].sPairPoly = 0;
+ } else
+ MATHICGB_ASSERT(false);
+ }
+ }
+ }
+
+ void buildMatrixAndClear(std::vector<RowTask>& tasks, QuadMatrix& quadMatrix) {
+ MATHICGB_LOG_TIME(F4MatrixBuild2) <<
+ "\n***** Constructing matrix *****\n";
+
+ if (tasks.empty()) {
+ quadMatrix = QuadMatrix();
+ quadMatrix.ring = &ring();
+ return;
+ }
+
+ initializeRowsToReduce(tasks);
+
+ // The MonoRef's cannot be Mono's since enumerable_thread_specific
+ // apparently requires the stored data type to be copyable and
+ // Mono is not copyable.
+ struct ThreadData {
+ MonoRef tmp1;
+ MonoRef tmp2;
+ F4ProtoMatrix block;
+ };
+
+ mgb::mtbb::enumerable_thread_specific<ThreadData> threadData([&](){
+ // We need to grab a lock since monoid isn't internally synchronized.
+ mgb::mtbb::mutex::scoped_lock guard(mCreateColumnLock);
+ ThreadData data = {
+ *monoid().alloc().release(),
+ *monoid().alloc().release()
+ };
+ return data;
+ });
+
+ // Construct the matrix as pre-blocks
+ mgb::mtbb::parallel_do(tasks.begin(), tasks.end(),
+ [&](const RowTask& task, TaskFeeder& feeder)
+ {
+ auto& data = threadData.local();
+ const auto& poly = *task.poly;
+
+ // It is perfectly permissible for task.sPairPoly to be non-null. The
+ // assert is there because of an interaction between S-pair
+ // elimination/choice and the use of halves of S-pairs as reducers. The
+ // current effect of these is that *all* S-pairs have a component split
+ // off so that sPairPoly is always null (this is non-trivial to
+ // realize). So if this assert goes off, you've messed that interaction
+ // up somehow or you are using this class in some new way. So you can
+ // remove the assert if necessary.
+ MATHICGB_ASSERT(task.sPairPoly == 0);
+
+ if (task.sPairPoly != 0) {
+ monoid().colons(
+ poly.getLeadMonomial(),
+ task.sPairPoly->getLeadMonomial(),
+ data.tmp2,
+ data.tmp1
+ );
+ appendRowSPair
+ (poly, data.tmp1, *task.sPairPoly, data.tmp2, data.block, feeder);
+ return;
+ }
+ if (task.desiredLead == nullptr)
+ monoid().setIdentity(data.tmp1);
+ else
+ monoid().divide(poly.getLeadMonomial(), *task.desiredLead, data.tmp1);
+ appendRow(data.tmp1, *task.poly, data.block, feeder);
+ });
+ MATHICGB_ASSERT(!threadData.empty()); // as tasks empty causes early return
+
+ // Free the monomials from all the tasks
+ for (const auto& task : tasks)
+ if (task.desiredLead != nullptr)
+ monoid().freeRaw(*task.desiredLead.castAwayConst());
+ tasks.clear();
+
+ // Move the proto-matrices across all threads into the projection.
+ F4MatrixProjection projection
+ (ring(), static_cast<ColIndex>(mMap.entryCount()));
+ const auto end = threadData.end();
+ for (auto& data : threadData) {
+ monoid().freeRaw(data.tmp1);
+ monoid().freeRaw(data.tmp2);
+ projection.addProtoMatrix(std::move(data.block));
+ }
+
+ // Sort columns by monomial and tell the projection of the resulting order
+ MonomialMap<ColIndex>::Reader reader(mMap);
+ typedef std::pair<ColIndex, ConstMonoPtr> IndexMono;
+ auto toPtr = [](std::pair<ColIndex, ConstMonoRef> p) {
+ return std::make_pair(p.first, p.second.ptr());
+ };
+ std::vector<IndexMono> columns;
+ std::transform
+ (reader.begin(), reader.end(), std::back_inserter(columns), toPtr);
+
+ const auto cmp = [&](const IndexMono& a, const IndexMono b) {
+ return monoid().lessThan(*b.second, *a.second);
+ };
+ mgb::mtbb::parallel_sort(columns.begin(), columns.end(), cmp);
+
+ const auto colEnd = columns.end();
+ for (auto it = columns.begin(); it != colEnd; ++it) {
+ const auto p = *it;
+ projection.addColumn(p.first, Monoid::toOld(*p.second), mIsColumnToLeft[p.first]);
+ }
+
+ quadMatrix = projection.makeAndClear(mMemoryQuantum);
+
+ MATHICGB_LOG(F4MatrixSizes)
+ << "F4["
+ << mathic::ColumnPrinter::commafy(quadMatrix.rowCount())
+ << " by "
+ << mathic::ColumnPrinter::commafy(
+ quadMatrix.computeLeftColCount() + quadMatrix.computeRightColCount())
+ << "]" << std::endl;
+
+#ifdef MATHICGB_DEBUG
+ for (size_t side = 0; side < 2; ++side) {
+ auto& monos = side == 0 ?
+ quadMatrix.leftColumnMonomials : quadMatrix.rightColumnMonomials;
+ for (auto it = monos.begin(); it != monos.end(); ++it) {
+ MATHICGB_ASSERT(!it->isNull());
+ }
+ }
+ for (RowIndex row = 0; row < quadMatrix.topLeft.rowCount(); ++row) {
+ MATHICGB_ASSERT(quadMatrix.topLeft.entryCountInRow(row) > 0);
+ MATHICGB_ASSERT(quadMatrix.topLeft.leadCol(row) == row);
+ }
+ MATHICGB_ASSERT(quadMatrix.debugAssertValid());
+#endif
+ }
const PolyRing& ring() const {return mBasis.ring();}
- Builder(const PolyBasis& basis, const size_t memoryQuantum);
+ const Monoid& monoid() const {return mBasis.ring().monoid();}
+
+ Builder(const PolyBasis& basis, const size_t memoryQuantum):
+ mMemoryQuantum(memoryQuantum),
+ mTmp(basis.ring().monoid().alloc()),
+ mBasis(basis),
+ mMap(basis.ring())
+ {
+ // This assert has to be _NO_ASSUME since otherwise the compiler will
+ // assume that the error checking branch here cannot be taken and optimize
+ // it away.
+ const Scalar maxScalar = std::numeric_limits<Scalar>::max();
+ MATHICGB_ASSERT_NO_ASSUME(ring().charac() <= maxScalar);
+ if (ring().charac() > maxScalar)
+ mathic::reportInternalError("F4MatrixBuilder2: too large characteristic.");
+ }
typedef const Map::Reader ColReader;
typedef std::vector<monomial> Monomials;
@@ -42,51 +291,218 @@ public:
/// the column of interest is to grab a lock, and the lock being grabbed
/// is being grabbed inside createColumn.
MATHICGB_NO_INLINE
- std::pair<ColIndex, ConstMonomial> createColumn(
- const_monomial monoA,
- const_monomial monoB,
+ std::pair<ColIndex, ConstMonoRef> createColumn(
+ ConstMonoRef monoA,
+ ConstMonoRef monoB,
TaskFeeder& feeder
- );
+ ) {
+ mgb::mtbb::mutex::scoped_lock lock(mCreateColumnLock);
+ // see if the column exists now after we have synchronized
+ {
+ const auto found(ColReader(mMap).findProduct(monoA, monoB));
+ if (found.first != 0)
+ return std::make_pair(*found.first, *found.second);
+ }
+
+ // The column really does not exist, so we need to create it
+ monoid().multiply(monoA, monoB, mTmp);
+ if (!monoid().hasAmpleCapacity(mTmp))
+ mathic::reportError("Monomial exponent overflow in F4MatrixBuilder2.");
+
+ // look for a reducer of mTmp
+ const size_t reducerIndex = mBasis.classicReducer(Monoid::toOld(*mTmp.ptr()));
+ const bool insertLeft = (reducerIndex != static_cast<size_t>(-1));
+
+ // Create the new left or right column
+ if (mIsColumnToLeft.size() >= std::numeric_limits<ColIndex>::max())
+ throw std::overflow_error("Too many columns in QuadMatrix");
+ const auto newIndex = static_cast<ColIndex>(mIsColumnToLeft.size());
+ const auto inserted = mMap.insert(std::make_pair(mTmp.ptr(), newIndex));
+ mIsColumnToLeft.push_back(insertLeft);
+
+ // schedule new task if we found a reducer
+ if (insertLeft) {
+ RowTask task = {};
+ task.poly = &mBasis.poly(reducerIndex);
+ task.desiredLead = inserted.first.second;
+ feeder.add(task);
+ }
+
+ return std::make_pair(*inserted.first.first, *inserted.first.second);
+ }
+
/// Append multiple * poly to block, creating new columns as necessary.
void appendRow(
- const_monomial multiple,
+ ConstMonoRef multiple,
const Poly& poly,
F4ProtoMatrix& block,
TaskFeeder& feeder
- );
+ ) {
+ const auto begin = poly.begin();
+ const auto end = poly.end();
+ const auto count = static_cast<size_t>(std::distance(begin, end));
+ MATHICGB_ASSERT(count < std::numeric_limits<ColIndex>::max());
+ auto indices = block.makeRowWithTheseScalars(poly);
+
+ auto it = begin;
+ if ((count % 2) == 1) {
+ ColReader reader(mMap);
+ const auto col = findOrCreateColumn
+ (it.getMonomial(), multiple, reader, feeder);
+ MATHICGB_ASSERT(it.getCoefficient() < std::numeric_limits<Scalar>::max());
+ MATHICGB_ASSERT(it.getCoefficient());
+ *indices = col.first;
+ ++indices;
+ ++it;
+ }
+ updateReader:
+ ColReader colMap(mMap);
+ MATHICGB_ASSERT((std::distance(it, end) % 2) == 0);
+ while (it != end) {
+ MATHICGB_ASSERT(it.getCoefficient() < std::numeric_limits<Scalar>::max());
+ MATHICGB_ASSERT(it.getCoefficient() != 0);
+ const auto scalar1 = static_cast<Scalar>(it.getCoefficient());
+ const auto mono1 = it.mono();
+
+ auto it2 = it;
+ ++it2;
+ MATHICGB_ASSERT(it2.getCoefficient() < std::numeric_limits<Scalar>::max());
+ MATHICGB_ASSERT(it2.getCoefficient() != 0);
+ const auto scalar2 = static_cast<Scalar>(it2.getCoefficient());
+ const auto mono2 = it2.mono();
+
+ const auto colPair = colMap.findTwoProducts(Monoid::toOld(mono1), Monoid::toOld(mono2), Monoid::toOld(multiple));
+ if (colPair.first == 0 || colPair.second == 0) {
+ createColumn(mono1, multiple, feeder);
+ createColumn(mono2, multiple, feeder);
+ goto updateReader;
+ }
+
+ *indices = *colPair.first;
+ ++indices;
+ *indices = *colPair.second;
+ ++indices;
+
+ it = ++it2;
+ }
+ }
/// Append poly*multiply - sPairPoly*sPairMultiply to block, creating new
/// columns as necessary.
void appendRowSPair(
- const Poly* poly,
- monomial multiply,
- const Poly* sPairPoly,
- monomial sPairMultiply,
+ const Poly& poly,
+ ConstMonoRef multiply,
+ const Poly& sPairPoly,
+ ConstMonoRef sPairMultiply,
F4ProtoMatrix& block,
TaskFeeder& feeder
- );
+ ) {
+ MATHICGB_ASSERT(!poly.isZero());
+ auto itA = poly.begin();
+ const auto endA = poly.end();
+
+ MATHICGB_ASSERT(!sPairPoly.isZero());
+ auto itB = sPairPoly.begin();
+ const auto endB = sPairPoly.end();
+
+ // skip leading terms since they cancel
+ MATHICGB_ASSERT(itA.getCoefficient() == itB.getCoefficient());
+ ++itA;
+ ++itB;
+
+ // @todo: handle overflow of termCount addition here
+ MATHICGB_ASSERT(poly.termCount() + sPairPoly.termCount() - 2 <=
+ std::numeric_limits<ColIndex>::max());
+ const auto maxCols =
+ static_cast<ColIndex>(poly.termCount() + sPairPoly.termCount() - 2);
+ auto row = block.makeRow(maxCols);
+ const auto indicesBegin = row.first;
+
+ const ColReader colMap(mMap);
+
+ auto mulA = multiply;
+ auto mulB = sPairMultiply;
+ while (itB != endB && itA != endA) {
+ const auto colA = findOrCreateColumn
+ (itA.getMonomial(), mulA, colMap, feeder);
+ const auto colB = findOrCreateColumn
+ (itB.getMonomial(), mulB, colMap, feeder);
+ const auto cmp = monoid().compare(colA.second, colB.second);
+
+ coefficient coeff = 0;
+ ColIndex col;
+ if (cmp != LT) {
+ coeff = itA.getCoefficient();
+ col = colA.first;
+ ++itA;
+ }
+ if (cmp != GT) {
+ coeff = ring().coefficientSubtract(coeff, itB.getCoefficient());
+ col = colB.first;
+ ++itB;
+ }
+ MATHICGB_ASSERT(coeff < std::numeric_limits<Scalar>::max());
+ if (coeff != 0) {
+ *row.first++ = col;
+ *row.second++ = static_cast<Scalar>(coeff);
+ }
+ }
+
+ for (; itA != endA; ++itA) {
+ const auto colA = findOrCreateColumn
+ (itA.getMonomial(), mulA, colMap, feeder);
+ *row.first++ = colA.first;
+ *row.second++ = static_cast<Scalar>(itA.getCoefficient());
+ }
+
+ for (; itB != endB; ++itB) {
+ const auto colB = findOrCreateColumn
+ (itB.getMonomial(), mulB, colMap, feeder);
+ const auto negative = ring().coefficientNegate(itB.getCoefficient());
+ *row.first = colB.first;
+ ++row.first;
+ *row.second = static_cast<Scalar>(negative);
+ ++row.second;
+ }
+
+ const auto toRemove =
+ maxCols - static_cast<ColIndex>(row.first - indicesBegin);
+ block.removeLastEntries(block.rowCount() - 1, toRemove);
+ }
/// As createColumn, except with much better performance in the common
/// case that the column for monoA * monoB already exists. In particular,
/// no lock is grabbed in that case.
MATHICGB_NO_INLINE
- std::pair<ColIndex, ConstMonomial> findOrCreateColumn(
- const_monomial monoA,
- const_monomial monoB,
+ std::pair<ColIndex, ConstMonoRef> findOrCreateColumn(
+ ConstMonoRef monoA,
+ ConstMonoRef monoB,
TaskFeeder& feeder
- );
+ ) {
+ const auto col = ColReader(mMap).findProduct(monoA, monoB);
+ if (col.first != 0)
+ return std::make_pair(*col.first, *col.second);
+ return createColumn(monoA, monoB, feeder);
+ }
/// As the overload that does not take a ColReader parameter, except with
/// better performance in the common case that the column already exists
/// and colMap is up-to-date.
MATHICGB_INLINE
- std::pair<ColIndex, ConstMonomial> findOrCreateColumn(
- const_monomial monoA,
- const_monomial monoB,
+ std::pair<ColIndex, ConstMonoRef> findOrCreateColumn(
+ ConstMonoRef monoA,
+ ConstMonoRef monoB,
const ColReader& colMap,
TaskFeeder& feeder
- );
+ ) {
+ const auto col = colMap.findProduct(monoA, monoB);
+ if (col.first == 0) {
+ // The reader may be out of date, so try again with a fresh reader.
+ return findOrCreateColumn(monoA, monoB, feeder);
+ }
+ return std::make_pair(*col.first, *col.second);
+ }
/// The split into left and right columns is not done until the whole matrix
/// has been constructed. This vector keeps track of which side each column
@@ -103,7 +519,7 @@ public:
/// A monomial for temporary scratch calculations. Protected by
/// mCreateColumnLock.
- monomial mTmp;
+ Mono mTmp;
/// Mapping from monomials to column indices.
Map mMap;
@@ -112,39 +528,6 @@ public:
const PolyBasis& mBasis;
};
-MATHICGB_NO_INLINE
-std::pair<F4MatrixBuilder2::ColIndex, ConstMonomial>
-F4MatrixBuilder2::Builder::findOrCreateColumn(
- const const_monomial monoA,
- const const_monomial monoB,
- TaskFeeder& feeder
-) {
- MATHICGB_ASSERT(!monoA.isNull());
- MATHICGB_ASSERT(!monoB.isNull());
- const auto col = ColReader(mMap).findProduct(monoA, monoB);
- if (col.first != 0)
- return std::make_pair(*col.first, col.second);
- return createColumn(monoA, monoB, feeder);
-}
-
-MATHICGB_INLINE
-std::pair<F4MatrixBuilder2::ColIndex, ConstMonomial>
-F4MatrixBuilder2::Builder::findOrCreateColumn(
- const const_monomial monoA,
- const const_monomial monoB,
- const ColReader& colMap,
- TaskFeeder& feeder
-) {
- MATHICGB_ASSERT(!monoA.isNull());
- MATHICGB_ASSERT(!monoB.isNull());
- const auto col = colMap.findProduct(monoA, monoB);
- if (col.first == 0) {
- // The reader may be out of date, so try again with a fresh reader.
- return findOrCreateColumn(monoA, monoB, feeder);
- }
- return std::make_pair(*col.first, col.second);
-}
-
F4MatrixBuilder2::F4MatrixBuilder2(
const PolyBasis& basis,
const size_t memoryQuantum
@@ -153,23 +536,6 @@ F4MatrixBuilder2::F4MatrixBuilder2(
mMemoryQuantum(memoryQuantum)
{}
-F4MatrixBuilder2::Builder::Builder(
- const PolyBasis& basis,
- const size_t memoryQuantum
-):
- mMemoryQuantum(memoryQuantum),
- mTmp(basis.ring().allocMonomial()),
- mBasis(basis),
- mMap(basis.ring())
-{
- // This assert has to be _NO_ASSUME since otherwise the compiler will assume
- // that the error checking branch here cannot be taken and optimize it away.
- const Scalar maxScalar = std::numeric_limits<Scalar>::max();
- MATHICGB_ASSERT_NO_ASSUME(ring().charac() <= maxScalar);
- if (ring().charac() > maxScalar)
- mathic::reportInternalError("F4MatrixBuilder2: too large characteristic.");
-}
-
void F4MatrixBuilder2::addSPolynomialToMatrix(
const Poly& polyA,
const Poly& polyB
@@ -179,14 +545,9 @@ void F4MatrixBuilder2::addSPolynomialToMatrix(
MATHICGB_ASSERT(!polyB.isZero());
MATHICGB_ASSERT(polyB.isMonic());
- RowTask task;
- task.poly = &polyA;
- task.sPairPoly = &polyB;
- task.desiredLead = ring().allocMonomial();
- ring().monomialLeastCommonMultiple(
- task.poly->getLeadMonomial(),
- task.sPairPoly->getLeadMonomial(),
- task.desiredLead);
+ auto desiredLead = monoid().alloc();
+ monoid().lcm(polyA.getLeadMonomial(), polyB.getLeadMonomial(), *desiredLead);
+ RowTask task = {desiredLead.release(), &polyA, &polyB};
mTodo.push_back(task);
}
@@ -199,18 +560,16 @@ void F4MatrixBuilder2::addPolynomialToMatrix(const Poly& poly) {
mTodo.push_back(task);
}
-void F4MatrixBuilder2::addPolynomialToMatrix
-(const_monomial multiple, const Poly& poly) {
+void F4MatrixBuilder2::addPolynomialToMatrix(
+ ConstMonoRef multiple,
+ const Poly& poly
+) {
if (poly.isZero())
return;
- RowTask task = {};
- task.poly = &poly;
- task.desiredLead = ring().allocMonomial();
- ring().monomialMult(poly.getLeadMonomial(), multiple, task.desiredLead);
-
- MATHICGB_ASSERT(task.sPairPoly == 0);
-
+ auto desiredLead = monoid().alloc();
+ monoid().multiply(poly.getLeadMonomial(), multiple, desiredLead);
+ RowTask task = {desiredLead.release(), &poly, nullptr};
mTodo.push_back(task);
}
@@ -219,377 +578,4 @@ void F4MatrixBuilder2::buildMatrixAndClear(QuadMatrix& quadMatrix) {
builder.buildMatrixAndClear(mTodo, quadMatrix);
}
-void F4MatrixBuilder2::Builder::buildMatrixAndClear(
- std::vector<RowTask>& mTodo,
- QuadMatrix& quadMatrix
-) {
- MATHICGB_LOG_TIME(F4MatrixBuild2) <<
- "\n***** Constructing matrix *****\n";
-
- if (mTodo.empty()) {
- quadMatrix = QuadMatrix();
- quadMatrix.ring = &ring();
- return;
- }
-
- // Use halves of S-pairs as reducers to decrease the number of entries that
- // need to be looked up.
- mgb::mtbb::parallel_sort(mTodo.begin(), mTodo.end(),
- [&](const RowTask& a, const RowTask& b)
- {
- if (a.sPairPoly == 0)
- return b.sPairPoly != 0;
- else
- return ring().monomialLT(a.desiredLead, b.desiredLead);
- });
- const auto size = mTodo.size();
- for (size_t i = 0; i < size;) {
- if (mTodo[i].sPairPoly == 0) {
- ++i;
- continue;
- }
- MATHICGB_ASSERT(!mTodo[i].desiredLead.isNull());
- MATHICGB_ASSERT(ColReader(mMap).find(mTodo[i].desiredLead).first == 0);
-
- // Create column for the lead term that cancels in the S-pair
- if (mIsColumnToLeft.size() >= std::numeric_limits<ColIndex>::max())
- throw std::overflow_error("Too many columns in QuadMatrix");
- const auto newIndex = static_cast<ColIndex>(mIsColumnToLeft.size());
- const auto inserted =
- mMap.insert(std::make_pair(mTodo[i].desiredLead, newIndex));
- mIsColumnToLeft.push_back(true);
- const auto& mono = inserted.first.second;
-
- // Schedule the two parts of the S-pair as separate rows. This adds a row
- // while creating the column in the hash table without adding a reducer
- // removes a row, effectively making this operation equivalent to taking
- // half of the S-pair and using it as a reducer.
- RowTask newTask = {};
- newTask.sPairPoly = 0;
- newTask.poly = mTodo[i].sPairPoly;
- newTask.desiredLead = ring().allocMonomial();
- ring().monomialCopy(mono, newTask.desiredLead);
-
- // Now we can strip off any part of an S-pair with the same cancelling lead
- // term that equals a or b since those rows are in the matrix.
- const Poly* const a = mTodo[i].poly;
- const Poly* const b = mTodo[i].sPairPoly;
-
- mTodo[i].sPairPoly = 0;
-
- mTodo.push_back(newTask);
- for (++i; i < size; ++i) {
- if (mTodo[i].sPairPoly == 0)
- continue;
- MATHICGB_ASSERT(!mTodo[i].desiredLead.isNull());
- if (!ring().monomialEQ(mTodo[i].desiredLead, mono))
- break;
- if (mTodo[i].poly == a || mTodo[i].poly == b) {
- mTodo[i].poly = mTodo[i].sPairPoly;
- mTodo[i].sPairPoly = 0;
- } else if (mTodo[i].sPairPoly == a || mTodo[i].sPairPoly == b) {
- mTodo[i].sPairPoly = 0;
- } else
- MATHICGB_ASSERT(false);
- }
- }
-
- // Process pending rows until we are done. Note that the methods
- // we are calling here can add more pending items.
-
- struct ThreadData {
- F4ProtoMatrix block;
- monomial tmp1;
- monomial tmp2;
- };
-
- mgb::mtbb::enumerable_thread_specific<ThreadData> threadData([&](){
- ThreadData data;
- {
- mgb::mtbb::mutex::scoped_lock guard(mCreateColumnLock);
- data.tmp1 = ring().allocMonomial();
- data.tmp2 = ring().allocMonomial();
- }
- return std::move(data);
- });
-
- // Construct the matrix as pre-blocks
- mgb::mtbb::parallel_do(mTodo.begin(), mTodo.end(),
- [&](const RowTask& task, TaskFeeder& feeder)
- {
- auto& data = threadData.local();
- const auto& poly = *task.poly;
-
- // It is perfectly permissible for task.sPairPoly to be non-null. The
- // assert is there because of an interaction between S-pair
- // elimination/choice and the use of halves of S-pairs as reducers. The
- // current effect of these is that *all* S-pairs have a component split
- // off so that sPairPoly is always null (this is non-trivial to
- // realize). So if this assert goes off, you've messed that interaction
- // up somehow or you are using this class in some new way. So you can
- // remove the assert if necessary.
- MATHICGB_ASSERT(task.sPairPoly == 0);
-
- if (task.sPairPoly != 0) {
- ring().monomialColons(
- poly.getLeadMonomial(),
- task.sPairPoly->getLeadMonomial(),
- data.tmp2,
- data.tmp1
- );
- appendRowSPair
- (&poly, data.tmp1, task.sPairPoly, data.tmp2, data.block, feeder);
- return;
- }
- if (task.desiredLead.isNull())
- ring().monomialSetIdentity(data.tmp1);
- else
- ring().monomialDivide
- (task.desiredLead, poly.getLeadMonomial(), data.tmp1);
- appendRow(data.tmp1, *task.poly, data.block, feeder);
- });
- MATHICGB_ASSERT(!threadData.empty()); // as mTodo empty causes early return
-
- // Free the monomials from all the tasks
- const auto todoEnd = mTodo.end();
- for (auto it = mTodo.begin(); it != todoEnd; ++it)
- if (!it->desiredLead.isNull())
- ring().freeMonomial(it->desiredLead);
- mTodo.clear();
-
- // Move the proto-matrices across all threads into the projection.
- F4MatrixProjection projection
- (ring(), static_cast<ColIndex>(mMap.entryCount()));
- const auto end = threadData.end();
- for (auto it = threadData.begin(); it != end; ++it) {
- projection.addProtoMatrix(std::move(it->block));
- ring().freeMonomial(it->tmp1);
- ring().freeMonomial(it->tmp2);
- }
-
- // Sort columns by monomial and tell the projection of the resulting order
- MonomialMap<ColIndex>::Reader reader(mMap);
- typedef std::pair<ColIndex, ConstMonomial> IndexMono;
- std::vector<IndexMono> columns(reader.begin(), reader.end());
- const auto cmp = [&](const IndexMono& a, const IndexMono b) {
- return ring().monomialLT(b.second, a.second);
- };
- mgb::mtbb::parallel_sort(columns.begin(), columns.end(), cmp);
-
- const auto colEnd = columns.end();
- for (auto it = columns.begin(); it != colEnd; ++it) {
- const auto p = *it;
- projection.addColumn(p.first, p.second, mIsColumnToLeft[p.first]);
- }
-
- quadMatrix = projection.makeAndClear(mMemoryQuantum);
- threadData.clear();
-
- MATHICGB_LOG(F4MatrixSizes)
- << "F4["
- << mathic::ColumnPrinter::commafy(quadMatrix.rowCount())
- << " by "
- << mathic::ColumnPrinter::commafy(
- quadMatrix.computeLeftColCount() + quadMatrix.computeRightColCount())
- << "]" << std::endl;
-
-#ifdef MATHICGB_DEBUG
- for (size_t side = 0; side < 2; ++side) {
- auto& monos = side == 0 ?
- quadMatrix.leftColumnMonomials : quadMatrix.rightColumnMonomials;
- for (auto it = monos.begin(); it != monos.end(); ++it) {
- MATHICGB_ASSERT(!it->isNull());
- }
- }
- for (RowIndex row = 0; row < quadMatrix.topLeft.rowCount(); ++row) {
- MATHICGB_ASSERT(quadMatrix.topLeft.entryCountInRow(row) > 0);
- MATHICGB_ASSERT(quadMatrix.topLeft.leadCol(row) == row);
- }
- MATHICGB_ASSERT(quadMatrix.debugAssertValid());
-#endif
-}
-
-std::pair<F4MatrixBuilder2::ColIndex, ConstMonomial>
-F4MatrixBuilder2::Builder::createColumn(
- const const_monomial monoA,
- const const_monomial monoB,
- TaskFeeder& feeder
-) {
- MATHICGB_ASSERT(!monoA.isNull());
- MATHICGB_ASSERT(!monoB.isNull());
-
- mgb::mtbb::mutex::scoped_lock lock(mCreateColumnLock);
- // see if the column exists now after we have synchronized
- {
- const auto found(ColReader(mMap).findProduct(monoA, monoB));
- if (found.first != 0)
- return std::make_pair(*found.first, found.second);
- }
-
- // The column really does not exist, so we need to create it
- ring().monomialMult(monoA, monoB, mTmp);
- if (!ring().monomialHasAmpleCapacity(mTmp))
- mathic::reportError("Monomial exponent overflow in F4MatrixBuilder2.");
-
- // look for a reducer of mTmp
- const size_t reducerIndex = mBasis.classicReducer(mTmp);
- const bool insertLeft = (reducerIndex != static_cast<size_t>(-1));
-
- // Create the new left or right column
- if (mIsColumnToLeft.size() >= std::numeric_limits<ColIndex>::max())
- throw std::overflow_error("Too many columns in QuadMatrix");
- const auto newIndex = static_cast<ColIndex>(mIsColumnToLeft.size());
- const auto inserted = mMap.insert(std::make_pair(mTmp, newIndex));
- mIsColumnToLeft.push_back(insertLeft);
-
- // schedule new task if we found a reducer
- if (insertLeft) {
- RowTask task = {};
- task.poly = &mBasis.poly(reducerIndex);
- task.desiredLead = inserted.first.second.castAwayConst();
- feeder.add(task);
- }
-
- return std::make_pair(*inserted.first.first, inserted.first.second);
-}
-
-void F4MatrixBuilder2::Builder::appendRow(
- const const_monomial multiple,
- const Poly& poly,
- F4ProtoMatrix& block,
- TaskFeeder& feeder
-) {
- MATHICGB_ASSERT(!multiple.isNull());
-
- const auto begin = poly.begin();
- const auto end = poly.end();
- const auto count = static_cast<size_t>(std::distance(begin, end));
- MATHICGB_ASSERT(count < std::numeric_limits<ColIndex>::max());
- auto indices = block.makeRowWithTheseScalars(poly);
-
- auto it = begin;
- if ((count % 2) == 1) {
- ColReader reader(mMap);
- const auto col = findOrCreateColumn
- (it.getMonomial(), multiple, reader, feeder);
- MATHICGB_ASSERT(it.getCoefficient() < std::numeric_limits<Scalar>::max());
- MATHICGB_ASSERT(it.getCoefficient());
- *indices = col.first;
- ++indices;
- ++it;
- }
-updateReader:
- ColReader colMap(mMap);
- MATHICGB_ASSERT((std::distance(it, end) % 2) == 0);
- while (it != end) {
- MATHICGB_ASSERT(it.getCoefficient() < std::numeric_limits<Scalar>::max());
- MATHICGB_ASSERT(it.getCoefficient() != 0);
- const auto scalar1 = static_cast<Scalar>(it.getCoefficient());
- const const_monomial mono1 = it.getMonomial();
-
- auto it2 = it;
- ++it2;
- MATHICGB_ASSERT(it2.getCoefficient() < std::numeric_limits<Scalar>::max());
- MATHICGB_ASSERT(it2.getCoefficient() != 0);
- const auto scalar2 = static_cast<Scalar>(it2.getCoefficient());
- const const_monomial mono2 = it2.getMonomial();
-
- const auto colPair = colMap.findTwoProducts(mono1, mono2, multiple);
- if (colPair.first == 0 || colPair.second == 0) {
- createColumn(mono1, multiple, feeder);
- createColumn(mono2, multiple, feeder);
- goto updateReader;
- }
-
- *indices = *colPair.first;
- ++indices;
- *indices = *colPair.second;
- ++indices;
-
- it = ++it2;
- }
-}
-
-void F4MatrixBuilder2::Builder::appendRowSPair(
- const Poly* poly,
- monomial multiply,
- const Poly* sPairPoly,
- monomial sPairMultiply,
- F4ProtoMatrix& block,
- TaskFeeder& feeder
-) {
- MATHICGB_ASSERT(!poly->isZero());
- MATHICGB_ASSERT(!multiply.isNull());
- MATHICGB_ASSERT(sPairPoly != 0);
- Poly::const_iterator itA = poly->begin();
- const Poly::const_iterator endA = poly->end();
-
- MATHICGB_ASSERT(!sPairPoly->isZero());
- MATHICGB_ASSERT(!sPairMultiply.isNull());
- Poly::const_iterator itB = sPairPoly->begin();
- Poly::const_iterator endB = sPairPoly->end();
-
- // skip leading terms since they cancel
- MATHICGB_ASSERT(itA.getCoefficient() == itB.getCoefficient());
- ++itA;
- ++itB;
-
- // @todo: handle overflow of termCount addition here
- MATHICGB_ASSERT(poly->termCount() + sPairPoly->termCount() - 2 <=
- std::numeric_limits<ColIndex>::max());
- const auto maxCols =
- static_cast<ColIndex>(poly->termCount() + sPairPoly->termCount() - 2);
- auto row = block.makeRow(maxCols);
- const auto indicesBegin = row.first;
-
- const ColReader colMap(mMap);
-
- const const_monomial mulA = multiply;
- const const_monomial mulB = sPairMultiply;
- while (itB != endB && itA != endA) {
- const auto colA = findOrCreateColumn
- (itA.getMonomial(), mulA, colMap, feeder);
- const auto colB = findOrCreateColumn
- (itB.getMonomial(), mulB, colMap, feeder);
- const auto cmp = ring().monomialCompare(colA.second, colB.second);
-
- coefficient coeff = 0;
- ColIndex col;
- if (cmp != LT) {
- coeff = itA.getCoefficient();
- col = colA.first;
- ++itA;
- }
- if (cmp != GT) {
- coeff = ring().coefficientSubtract(coeff, itB.getCoefficient());
- col = colB.first;
- ++itB;
- }
- MATHICGB_ASSERT(coeff < std::numeric_limits<Scalar>::max());
- if (coeff != 0) {
- *row.first++ = col;
- *row.second++ = static_cast<Scalar>(coeff);
- }
- }
-
- for (; itA != endA; ++itA) {
- const auto colA = findOrCreateColumn
- (itA.getMonomial(), mulA, colMap, feeder);
- *row.first++ = colA.first;
- *row.second++ = static_cast<Scalar>(itA.getCoefficient());
- }
-
- for (; itB != endB; ++itB) {
- const auto colB = findOrCreateColumn
- (itB.getMonomial(), mulB, colMap, feeder);
- const auto negative = ring().coefficientNegate(itB.getCoefficient());
- *row.first++ = colB.first;
- *row.second++ = static_cast<Scalar>(negative);
- }
-
- const auto toRemove =
- maxCols - static_cast<ColIndex>(row.first - indicesBegin);
- block.removeLastEntries(block.rowCount() - 1, toRemove);
-}
-
MATHICGB_NAMESPACE_END
diff --git a/src/mathicgb/F4MatrixBuilder2.hpp b/src/mathicgb/F4MatrixBuilder2.hpp
index 3b21076..06ac20b 100755
--- a/src/mathicgb/F4MatrixBuilder2.hpp
+++ b/src/mathicgb/F4MatrixBuilder2.hpp
@@ -21,13 +21,14 @@ MATHICGB_NAMESPACE_BEGIN
/// very workable without an RAII monomial handle or a scope exit
/// functionality, so add one of those before fixing this.
class F4MatrixBuilder2 {
-private:
- typedef SparseMatrix::ColIndex ColIndex;
- typedef SparseMatrix::Scalar Scalar;
- typedef MonomialMap<ColIndex> Map;
- typedef SparseMatrix::RowIndex RowIndex;
-
public:
+ typedef PolyRing::Monoid Monoid;
+ typedef Monoid::Mono Mono;
+ typedef Monoid::MonoRef MonoRef;
+ typedef Monoid::ConstMonoRef ConstMonoRef;
+ typedef Monoid::MonoPtr MonoPtr;
+ typedef Monoid::ConstMonoPtr ConstMonoPtr;
+
/// memoryQuantum is how much to increase the memory size by each time the
/// current amount of memory is exhausted. A value of 0 indicates to start
/// small and double the quantum at each exhaustion.
@@ -46,7 +47,7 @@ public:
/// matrix. No ownership is taken, but poly must remain valid until
/// the matrix is constructed. multiple is copied, so it need not
/// remain valid.
- void addPolynomialToMatrix(const_monomial multiple, const Poly& poly);
+ void addPolynomialToMatrix(ConstMonoRef multiple, const Poly& poly);
/// As the overload with a multiple, where the multiple is 1.
void addPolynomialToMatrix(const Poly& poly);
@@ -71,6 +72,7 @@ public:
void buildMatrixAndClear(QuadMatrix& matrix);
const PolyRing& ring() const {return mBasis.ring();}
+ const Monoid& monoid() const {return mBasis.ring().monoid();}
private:
/// Represents the task of adding a row to the matrix. If sPairPoly is null
@@ -79,7 +81,7 @@ private:
/// where multiply and sPairMultiply are such that the leading terms become
/// desiredLead.
struct RowTask {
- monomial desiredLead; // multiply a monomial onto poly to get this lead
+ ConstMonoPtr desiredLead; // multiply a monomial onto poly to get this lead
const Poly* poly;
const Poly* sPairPoly;
};
diff --git a/src/mathicgb/FixedSizeMonomialMap.h b/src/mathicgb/FixedSizeMonomialMap.h
index 650ae10..5f484c9 100755
--- a/src/mathicgb/FixedSizeMonomialMap.h
+++ b/src/mathicgb/FixedSizeMonomialMap.h
@@ -26,8 +26,15 @@ MATHICGB_NAMESPACE_BEGIN
template<class T>
class FixedSizeMonomialMap {
public:
+ typedef PolyRing::Monoid Monoid;
+ typedef Monoid::Mono Mono;
+ typedef Monoid::MonoRef MonoRef;
+ typedef Monoid::ConstMonoRef ConstMonoRef;
+ typedef Monoid::MonoPtr MonoPtr;
+ typedef Monoid::ConstMonoPtr ConstMonoPtr;
+
typedef T mapped_type;
- typedef std::pair<const_monomial, mapped_type> value_type;
+ typedef std::pair<ConstMonoPtr, mapped_type> value_type;
/// Iterates through entries in the hash table.
class const_iterator;
@@ -44,7 +51,7 @@ public:
make_unique_array<Atomic<Node*>>(hashMaskToBucketCount(mHashToIndexMask))
),
mRing(ring),
- mNodeAlloc(sizeofNode(ring))
+ mNodeAlloc(Node::bytesPerNode(ring.monoid()))
{
// Calling new int[x] does not zero the array. std::atomic has a trivial
// constructor so the same thing is true of new atomic[x]. Calling
@@ -75,14 +82,16 @@ public:
mNodeAlloc(std::move(map.mNodeAlloc))
{
// We can store relaxed as the constructor does not run concurrently.
+ const auto relax = std::memory_order_relaxed;
+
setTableEntriesToNullRelaxed();
const auto tableEnd = map.mBuckets.get() + map.bucketCount();
for (auto tableIt = map.mBuckets.get(); tableIt != tableEnd; ++tableIt) {
for (Node* node = tableIt->load(); node != 0;) {
- const size_t index = hashToIndex(mRing.monomialHashValue(node->mono));
- Node* const next = node->next.load();
- node->next.store(mBuckets[index].load());
- mBuckets[index].store(node, std::memory_order_relaxed);
+ const size_t index = hashToIndex(monoid().hash(node->mono()));
+ const auto next = node->next(relax);
+ node->setNext(mBuckets[index].load(relax), relax);
+ mBuckets[index].store(node, relax);
node = next;
}
}
@@ -94,6 +103,7 @@ public:
}
const PolyRing& ring() const {return mRing;}
+ const Monoid& monoid() const {return mRing.monoid();}
/// The range [begin(), end()) contains all entries in the hash table.
/// Insertions invalidate all iterators. Beware that insertions can
@@ -111,57 +121,60 @@ public:
}
/// Returns the value associated to mono or null if there is no such value.
- /// Also returns an internal monomial that equals mono of such a monomial
+ /// Also returns an internal monomial that equals mono if such a monomial
/// exists.
- std::pair<const mapped_type*, ConstMonomial>
- find(const const_monomial mono) const {
- const HashValue monoHash = mRing.monomialHashValue(mono);
- const Node* node = bucketAtIndex(hashToIndex(monoHash));
- for (; node != 0; node = node->next.load(std::memory_order_consume)) {
+ std::pair<const mapped_type*, ConstMonoPtr>
+ find(ConstMonoRef mono) const {
+ const auto monoHash = monoid().hash(mono);
+ const auto* node = bucketAtIndex(hashToIndex(monoHash));
+ for (; node != 0; node = node->next(std::memory_order_consume)) {
// To my surprise, it seems to be faster to comment out this branch.
// I guess the hash table has too few collisions to make it worth it.
//if (monoHash != mRing.monomialHashValue(node->mono))
// continue;
- if (mRing.monomialEqualHintTrue(mono, node->mono))
- return std::make_pair(&node->value, node->mono);
+ if (monoid().equalHintTrue(mono, node->mono()))
+ return std::make_pair(&node->value, node->mono().ptr());
}
- return std::make_pair(static_cast<const mapped_type*>(0), ConstMonomial(0));
+ return std::make_pair(nullptr, ConstMonoPtr());
}
// As find on the product a*b but also returns the monomial that is the
// product.
MATHICGB_INLINE
- std::pair<const mapped_type*, ConstMonomial> findProduct(
- const const_monomial a,
- const const_monomial b
+ std::pair<const mapped_type*, ConstMonoPtr> findProduct(
+ ConstMonoRef a,
+ ConstMonoRef b
) const {
- const HashValue abHash = mRing.monomialHashOfProduct(a, b);
+ const HashValue abHash = monoid().hashOfProduct(a, b);
const Node* node = bucketAtIndex(hashToIndex(abHash));
- for (; node != 0; node = node->next.load(std::memory_order_consume)) {
+ for (; node != 0; node = node->next(std::memory_order_consume)) {
// To my surprise, it seems to be faster to comment out this branch.
// I guess the hash table has too few collisions to make it worth it.
//if (abHash != mRing.monomialHashValue(node->mono))
// continue;
- if (mRing.monomialIsProductOfHintTrue(a, b, node->mono))
- return std::make_pair(&node->value, node->mono);
+ if (monoid().isProductOfHintTrue(a, b, node->mono()))
+ return std::make_pair(&node->value, node->mono().ptr());
}
- return std::make_pair(static_cast<const mapped_type*>(0), ConstMonomial(0));
+ return std::make_pair(nullptr, nullptr);
}
/// As findProduct but looks for a1*b and a2*b at one time.
MATHICGB_INLINE
std::pair<const mapped_type*, const mapped_type*> findTwoProducts(
- const const_monomial a1,
- const const_monomial a2,
- const const_monomial b
+ ConstMonoRef a1,
+ ConstMonoRef a2,
+ ConstMonoRef b
) const {
- const HashValue a1bHash = mRing.monomialHashOfProduct(a1, b);
- const HashValue a2bHash = mRing.monomialHashOfProduct(a2, b);
+ const HashValue a1bHash = monoid().hashOfProduct(a1, b);
+ const HashValue a2bHash = monoid().hashOfProduct(a2, b);
const Node* const node1 = bucketAtIndex(hashToIndex(a1bHash));
const Node* const node2 = bucketAtIndex(hashToIndex(a2bHash));
- if (node1 != 0 && node2 != 0 && mRing.monomialIsTwoProductsOfHintTrue
- (a1, a2, b, node1->mono, node2->mono)
+ if (
+ node1 != 0 &&
+ node2 != 0 &&
+ monoid().isTwoProductsOfHintTrue
+ (a1, a2, b, node1->mono(), node2->mono())
)
return std::make_pair(&node1->value, &node2->value);
else
@@ -176,7 +189,7 @@ public:
/// inserted value equals the already present value.
///
/// p.first.second is a internal monomial that equals value.first.
- std::pair< std::pair<const mapped_type*, ConstMonomial>, bool>
+ std::pair< std::pair<const mapped_type*, ConstMonoPtr>, bool>
insert(const value_type& value) {
const mgb::mtbb::mutex::scoped_lock lockGuard(mInsertionMutex);
// find() loads buckets with memory_order_consume, so it may seem like
@@ -188,25 +201,27 @@ public:
// insertion mutex and by locking that mutex we have synchronized with
// all threads that previously did insertions.
{
- const auto found = find(value.first);
- if (found.first != 0)
- return std::make_pair(found, false); // key already present
+ const auto found = find(*value.first);
+ if (found.first != 0) {
+ auto p = std::make_pair(found.first, *found.second);
+ return std::make_pair(p, false); // key already present
+ }
}
const auto node = static_cast<Node*>(mNodeAlloc.alloc());
- const size_t index = hashToIndex(mRing.monomialHashValue(value.first));
+ const size_t index = hashToIndex(monoid().hash(*value.first));
// the constructor initializes the first field of node->mono, so
// it has to be called before copying the monomial.
new (node) Node(bucketAtIndex(index), value.second);
- {
- Monomial nodeTmp(node->mono);
- ring().monomialCopy(value.first, nodeTmp);
- }
+ monoid().copy(*value.first, node->mono());
+
// we cannot store with memory_order_relaxed here because unlocking the
// lock only synchronizes with threads who later grab the lock - it does
// not synchronize with reading threads since they do not grab the lock.
mBuckets[index].store(node, std::memory_order_release);
- return std::make_pair(std::make_pair(&node->value, node->mono), true); // successful insertion
+
+ auto p = std::make_pair(&node->value, node->constMono());
+ return std::make_pair(p, true); // successful insertion
}
/// This operation removes all entries from the table. This operation
@@ -238,12 +253,29 @@ private:
tableIt->store(0, std::memory_order_relaxed);
}
- struct Node {
- Node(Node* next, const mapped_type value): next(next), value(value) {}
+ class Node {
+ public:
+ Node(Node* next, const mapped_type value): mNext(next), value(value) {}
+
+ MonoRef mono() {return Monoid::toRef(mMono);}
+ ConstMonoRef mono() const {return Monoid::toRef(mMono);}
+ ConstMonoRef constMono() const {return Monoid::toRef(mMono);}
+
+ Node* next(std::memory_order order) {return mNext.load(order);}
+ const Node* next(std::memory_order order) const {return mNext.load(order);}
+ void setNext(Node* next, std::memory_order order) {
+ mNext.store(next, order);
+ }
- Atomic<Node*> next;
const mapped_type value;
- exponent mono[1];
+
+ static size_t bytesPerNode(const Monoid& monoid) {
+ return sizeof(Node) + sizeof(exponent) * (monoid.entryCount() - 1);
+ }
+
+ private:
+ Atomic<Node*> mNext;
+ exponent mMono[1];
};
static HashValue computeHashMask(const size_t requestedBucketCount) {
@@ -269,10 +301,6 @@ private:
return count;
}
- static size_t sizeofNode(const PolyRing& ring) {
- return sizeof(Node) - sizeof(exponent) + ring.maxMonomialByteSize();
- }
-
size_t hashToIndex(HashValue hash) const {
const auto index = hash & mHashToIndexMask;
MATHICGB_ASSERT(index == hash % bucketCount());
@@ -299,7 +327,7 @@ public:
class const_iterator {
public:
typedef std::forward_iterator_tag iterator_category;
- typedef std::pair<mapped_type, ConstMonomial> value_type;
+ typedef std::pair<mapped_type, ConstMonoRef> value_type;
typedef ptrdiff_t difference_type;
typedef size_t distance_type;
typedef value_type* pointer;
@@ -310,7 +338,7 @@ public:
const_iterator& operator++() {
MATHICGB_ASSERT(mNode != 0);
MATHICGB_ASSERT(mBucket != mBucketsEnd);
- const Node* const node = mNode->next.load(std::memory_order_consume);
+ const Node* const node = mNode->next(std::memory_order_consume);
if (node != 0)
mNode = node;
else
@@ -329,7 +357,7 @@ public:
const value_type operator*() const {
MATHICGB_ASSERT(mNode != 0);
- return std::make_pair(mNode->value, mNode->mono);
+ return std::make_pair(mNode->value, mNode->mono());
}
private:
diff --git a/src/mathicgb/MonoMonoid.hpp b/src/mathicgb/MonoMonoid.hpp
index def8cb3..ee9c08d 100755
--- a/src/mathicgb/MonoMonoid.hpp
+++ b/src/mathicgb/MonoMonoid.hpp
@@ -5,6 +5,7 @@
#include "MonoOrder.hpp"
#include "NonCopyable.hpp"
+#include <cstddef>
#include <vector>
#include <algorithm>
#include <memtailor.h>
@@ -372,6 +373,8 @@ public:
static ConstMonoRef toRef(const Exponent* e) {return ConstMonoRef(e);}
static Exponent* toOld(MonoRef e) {return rawPtr(e);}
static const Exponent* toOld(ConstMonoRef e) {return rawPtr(e);}
+ static Exponent* toOld(Mono& e) {return rawPtr(e);}
+ static const Exponent* toOld(const Mono& e) {return rawPtr(e);}
// *** Constructors and accessors
@@ -877,6 +880,9 @@ public:
/// Returns the number of gradings.
using Base::gradingCount;
+ /// Returns the number of entries per monomial. This includes entries
+ /// used fo r internal book-keeping, so it can be greater than varCount().
+ using Base::entryCount;
// *** Monomial mutating computations
@@ -969,7 +975,7 @@ public:
/// Sets mono to 1, which is the identity for multiplication.
void setIdentity(MonoRef mono) const {
- std::fill_n(rawPtr(mono), entryCount(), static_cast<Exponent>(0));
+ std::fill_n(rawPtr(mono), entryCount(), 0);
MATHICGB_ASSERT(debugValid(mono));
MATHICGB_ASSERT(isIdentity(mono));
@@ -1169,7 +1175,8 @@ public:
class ConstMonoPtr {
public:
- ConstMonoPtr(): mMono(0) {}
+ ConstMonoPtr(): mMono(nullptr) {}
+ ConstMonoPtr(std::nullptr_t): mMono(nullptr) {}
ConstMonoPtr(const ConstMonoPtr& mono): mMono(rawPtr(mono)) {}
ConstMonoPtr operator=(const ConstMonoPtr& mono) {
@@ -1183,8 +1190,15 @@ public:
/// to observe the ptr/ref distinction. Kill it with fire!
operator ConstMonoRef() const {return ConstMonoRef(*this);}
- bool isNull() const {return mMono == 0;}
- void toNull() {mMono = 0;}
+ bool isNull() const {return mMono == nullptr;}
+ void toNull() {mMono = nullptr;}
+
+ bool operator==(std::nullptr_t) const {return isNull();}
+ bool operator!=(std::nullptr_t) const {return !isNull();}
+
+ MonoPtr castAwayConst() const {
+ return MonoPtr(const_cast<Exponent*>(mMono));
+ }
private:
friend class MonoMonoid;
@@ -1197,7 +1211,8 @@ public:
class MonoPtr {
public:
- MonoPtr(): mMono(0) {}
+ MonoPtr(): mMono(nullptr) {}
+ MonoPtr(std::nullptr_t): mMono(nullptr) {}
MonoPtr(const MonoPtr& mono): mMono(mono.mMono) {}
MonoPtr operator=(const MonoPtr& mono) {
@@ -1211,13 +1226,17 @@ public:
/// to observe the ptr/ref distinction. Kill it with fire!
operator MonoRef() const {return MonoRef(*this);}
- bool isNull() const {return mMono == 0;}
- void toNull() {mMono = 0;}
+ bool isNull() const {return mMono == nullptr;}
+ void toNull() {mMono = nullptr;}
+
+ bool operator==(std::nullptr_t) const {return isNull();}
+ bool operator!=(std::nullptr_t) const {return !isNull();}
operator ConstMonoPtr() const {return ConstMonoPtr(mMono);}
private:
friend class MonoMonoid;
+ friend class ConstMonoPtr;
friend class PolyRing; // todo: remove
friend class Poly; // todo: remove
@@ -1229,7 +1248,8 @@ public:
class Mono : public NonCopyable<Mono> {
public:
- Mono(): mMono(), mPool(0) {}
+ Mono(): mMono(), mPool(nullptr) {}
+ Mono(std::nullptr_t): mMono(), mPool(nullptr) {}
/// Passes ownership of the resources of mono to this object. Mono must
/// have been allocated from pool and it must have no other owner.
@@ -1243,7 +1263,7 @@ public:
Mono(Mono&& mono): mMono(mono.mMono), mPool(mono.mPool) {
mono.mMono.toNull();
- mono.mPool = 0;
+ mono.mPool = nullptr;
}
~Mono() {toNull();}
@@ -1255,7 +1275,7 @@ public:
mono.mMono.toNull();
mPool = mono.mPool;
- mono.mPool = 0;
+ mono.mPool = nullptr;
}
/// Sets this object to null but does NOT free the resources previously
@@ -1264,20 +1284,31 @@ public:
/// already null then the returned MonoPtr is also null.
MonoPtr release() {
const auto oldPtr = ptr();
- mMono = 0;
- mPool = 0;
+ mMono.toNull();
+ mPool = nullptr;
return oldPtr;
}
bool isNull() const {return mMono.isNull();}
void toNull() {mPool->free(std::move(*this));}
- MonoPtr ptr() const {return mMono;}
+ MonoPtr ptr() {return mMono;}
+ ConstMonoPtr ptr() const {return mMono;}
+
+ MonoRef operator*() {
+ MATHICGB_ASSERT(!isNull());
+ return *mMono;
+ }
+
+ ConstMonoRef operator*() const {
+ MATHICGB_ASSERT(!isNull());
+ return *mMono;
+ }
- /// @todo: Get rid of this as soon as all code has been migrated
- /// to observe the ptr/ref distinction. Kill it with fire!
operator MonoPtr() const {return ptr();}
+ /// @todo: Get rid of this as soon as all code has been migrated
+ /// to observe the ptr/ref distinction. Kill it with fire!
operator MonoRef() const {
MATHICGB_ASSERT(!isNull());
return *mMono;
@@ -1307,6 +1338,7 @@ public:
private:
void operator=(const MonoRef&); // not available
friend class MonoMonoid;
+ friend class ConstMonoRef;
MonoRef(MonoPtr mono): mMono(mono) {}
Exponent* internalRawPtr() const {return rawPtr(mMono);}
@@ -1327,6 +1359,8 @@ public:
/// to observe the ptr/ref distinction. Kill it with fire!
operator ConstMonoPtr() const {return ptr();}
+ MonoRef castAwayConst() const {return MonoRef(mMono.castAwayConst());}
+
private:
void operator=(const MonoRef&); // not available
friend class MonoMonoid;
@@ -1363,8 +1397,8 @@ public:
if (mono.isNull())
return;
freeRaw(mono);
- mono.mMono = 0;
- mono.mPool = 0;
+ mono.mMono.toNull();
+ mono.mPool = nullptr;
}
void freeRaw(MonoRef mono) {mPool.free(rawPtr(mono));}
@@ -1812,12 +1846,6 @@ private:
// Layout in memory:
// [component] [exponents...] [order data...] [hash]
- /// Returns how many Exponents are necessary to store a
- /// monomial. This can include other data than the exponents, so
- /// this number can be larger than varCount().
- using Base::entryCount;
- //size_t entryCount() const {return mEntryCount;}
-
VarIndex componentIndex() const {
//static_assert(HasComponent, "");
return 0;
diff --git a/src/mathicgb/MonomialMap.hpp b/src/mathicgb/MonomialMap.hpp
index 5cd4af2..948c816 100755
--- a/src/mathicgb/MonomialMap.hpp
+++ b/src/mathicgb/MonomialMap.hpp
@@ -28,7 +28,7 @@ MATHICGB_NAMESPACE_BEGIN
///
/// 1) grab a reader X
/// 2) perform queries on X until done or there is a miss
-/// 3) replace X with a fresh reader Y
+/// 3) replace X with a fresh reader
/// 4) go to 2 if the miss is now a hit
/// 5) grab a lock shared with all writers
/// 6) perform the query again - a miss is now guaranteed to be accurate
@@ -41,6 +41,13 @@ MATHICGB_NAMESPACE_BEGIN
template<class T>
class MonomialMap {
public:
+ typedef PolyRing::Monoid Monoid;
+ typedef Monoid::Mono Mono;
+ typedef Monoid::MonoRef MonoRef;
+ typedef Monoid::ConstMonoRef ConstMonoRef;
+ typedef Monoid::MonoPtr MonoPtr;
+ typedef Monoid::ConstMonoPtr ConstMonoPtr;
+
typedef T mapped_type;
typedef FixedSizeMonomialMap<mapped_type> FixedSizeMap;
typedef typename FixedSizeMap::value_type value_type;
@@ -56,7 +63,9 @@ public:
}
~MonomialMap() {
- delete mMap.load();
+ // We can load with std::memory_order_relaxed because the destructor
+ // cannot run concurrently.
+ delete mMap.load(std::memory_order_relaxed);
}
const PolyRing& ring() const {return mRing;}
@@ -92,16 +101,16 @@ public:
/// inserted. Also returns the internal monomial that matches mono if such
/// a monomial exists. Misses can be spurious! Read the comments on the parent
/// class.
- std::pair<const mapped_type*, ConstMonomial>
- find(const_monomial mono) const {
+ std::pair<const mapped_type*, ConstMonoPtr>
+ find(ConstMonoRef mono) const {
return mMap.find(mono);
}
// As find but looks for the product of a and b and also returns the
// monomal that is the product.
- std::pair<const mapped_type*, ConstMonomial> findProduct(
- const const_monomial a,
- const const_monomial b
+ std::pair<const mapped_type*, ConstMonoPtr> findProduct(
+ ConstMonoRef a,
+ ConstMonoRef b
) const {
return mMap.findProduct(a, b);
}
@@ -132,7 +141,11 @@ public:
/// Removes all entries from the hash table. This requires mutual exclusion
/// from and synchronization with all readers and writers.
- void clearNonConcurrent() {mMap.load()->clearNonConcurrent();}
+ void clearNonConcurrent() {
+ // We can load with std::memory_order_relaxed because this method
+ // requires external synchronization.
+ mMap.load(std::memory_order_relaxed)->clearNonConcurrent();
+ }
/// Makes value.first map to value.second unless value.first is already
/// present in the map - in that case nothing is done. If p is the returned
@@ -141,7 +154,7 @@ public:
/// equal value.second if an insertion was not performed - unless the
/// inserted value equals the already present value. p.first.second is an
/// internal monomial that equals value.first.
- std::pair< std::pair<const mapped_type*, ConstMonomial>, bool>
+ std::pair<std::pair<const mapped_type*, ConstMonoPtr>, bool>
insert(const value_type& value) {
const mgb::mtbb::mutex::scoped_lock lockGuard(mInsertionMutex);
@@ -180,8 +193,8 @@ public:
return p;
}
- /// Return the number of entries. This method requires locking so do not
- /// call too much. The count may have
+ /// Return the number of entries. This method uses internal synchronization
+ /// so do not call too much or you'll get degraded performance.
size_t entryCount() const {
const mgb::mtbb::mutex::scoped_lock lockGuard(mInsertionMutex);
// We can load with std::memory_order_relaxed because we are holding the
@@ -211,7 +224,7 @@ private:
/// keep these around as we have no way to determine if there are still
/// readers looking at them. This could be changed at the cost of
/// more overhead in the Reader constructor and destructor.
- std::vector< std::unique_ptr<FixedSizeMap>> mOldMaps;
+ std::vector<std::unique_ptr<FixedSizeMap>> mOldMaps;
};
MATHICGB_NAMESPACE_END
diff --git a/src/mathicgb/Poly.hpp b/src/mathicgb/Poly.hpp
index 01c22c0..5041499 100755
--- a/src/mathicgb/Poly.hpp
+++ b/src/mathicgb/Poly.hpp
@@ -14,6 +14,10 @@ MATHICGB_NAMESPACE_BEGIN
class Poly {
public:
+ typedef PolyRing::Monoid Monoid;
+ typedef Monoid::ConstMonoRef ConstMonoRef;
+ typedef Monoid::MonoRef MonoRef;
+
Poly(const PolyRing& ring) : R(&ring) {MATHICGB_ASSERT(R != 0);}
void parse(std::istream &i); // reads into this, sorts terms
@@ -42,6 +46,7 @@ public:
iterator operator++() { ++ic; im += monsize; return *this; }
coefficient &getCoefficient() const { return *ic; }
monomial getMonomial() const { return &*im; }
+ MonoRef mono() const {return getMonomial();}
size_t operator-(const iterator &b) const { return ic - b.ic; }
friend bool operator==(const iterator &a, const iterator &b);
friend bool operator!=(const iterator &a, const iterator &b);
@@ -76,6 +81,7 @@ public:
const_iterator operator++() { ++ic; im += monsize; return *this; }
coefficient getCoefficient() const { return *ic; }
const_monomial getMonomial() const { return &*im; }
+ ConstMonoRef mono() const {return getMonomial();}
size_t operator-(const const_iterator &b) const { return ic - b.ic; }
friend bool operator==(const const_iterator &a, const const_iterator &b);
friend bool operator!=(const const_iterator &a, const const_iterator &b);
diff --git a/src/mathicgb/QuadMatrixBuilder.cpp b/src/mathicgb/QuadMatrixBuilder.cpp
index 490ef5f..225b9f6 100755
--- a/src/mathicgb/QuadMatrixBuilder.cpp
+++ b/src/mathicgb/QuadMatrixBuilder.cpp
@@ -72,9 +72,11 @@ namespace {
auto p(inserted.first);
toMonomial.back() = copied;
- MATHICGB_ASSERT(ring.monomialEqualHintTrue(copied, p.second));
+ MATHICGB_ASSERT(ring.monoid().equalHintTrue(copied, *p.second));
MATHICGB_ASSERT(*p.first == QuadMatrixBuilder::LeftRightColIndex(colCount, left));
- return std::make_pair(*p.first, p.second);
+
+ auto ptr = const_cast<exponent*>(Monoid::toOld(*p.second));
+ return std::make_pair(*p.first, ptr);
} catch (...) {
toMonomial.pop_back();
ring.freeMonomial(copied);
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/mathicgb.git
More information about the debian-science-commits
mailing list