[mathicgb] 366/393: Migrated F4MatrixBuilder, QuadMatrix and a few related classes to use MonoMonoid.
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 bc394f44722fe8a3c4ff97cfab366b9366183d05
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date: Tue Sep 17 15:28:15 2013 +0200
Migrated F4MatrixBuilder, QuadMatrix and a few related classes to use MonoMonoid.
---
src/mathicgb/F4MatrixBuilder.cpp | 149 +++++++++++++++++-------------------
src/mathicgb/F4MatrixBuilder.hpp | 47 ++++++------
src/mathicgb/F4MatrixBuilder2.cpp | 8 +-
src/mathicgb/F4MatrixProjection.cpp | 6 +-
src/mathicgb/F4MatrixProjection.hpp | 11 ++-
src/mathicgb/F4MatrixReducer.cpp | 45 +----------
src/mathicgb/F4MatrixReducer.hpp | 2 -
src/mathicgb/F4Reducer.cpp | 75 +++++++++---------
src/mathicgb/MonoMonoid.hpp | 8 ++
src/mathicgb/MonomialMap.hpp | 6 +-
src/mathicgb/PolyRing.hpp | 2 +-
src/mathicgb/QuadMatrix.cpp | 47 ++++++------
src/mathicgb/QuadMatrix.hpp | 32 ++++++--
src/mathicgb/QuadMatrixBuilder.cpp | 12 +--
src/mathicgb/QuadMatrixBuilder.hpp | 17 ++--
src/mathicgb/SparseMatrix.cpp | 4 +-
src/mathicgb/SparseMatrix.hpp | 2 +-
src/test/F4MatrixBuilder.cpp | 10 +--
src/test/F4MatrixReducer.cpp | 4 +-
src/test/QuadMatrixBuilder.cpp | 41 ++++++----
src/test/SparseMatrix.cpp | 4 +-
21 files changed, 263 insertions(+), 269 deletions(-)
diff --git a/src/mathicgb/F4MatrixBuilder.cpp b/src/mathicgb/F4MatrixBuilder.cpp
index 9052631..fd0796d 100755
--- a/src/mathicgb/F4MatrixBuilder.cpp
+++ b/src/mathicgb/F4MatrixBuilder.cpp
@@ -14,12 +14,10 @@ MATHICGB_NAMESPACE_BEGIN
MATHICGB_NO_INLINE
auto F4MatrixBuilder::findOrCreateColumn(
- const const_monomial monoA,
- const const_monomial monoB,
+ ConstMonoRef monoA,
+ ConstMonoRef 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);
@@ -28,13 +26,11 @@ auto F4MatrixBuilder::findOrCreateColumn(
MATHICGB_INLINE
auto F4MatrixBuilder::findOrCreateColumn(
- const const_monomial monoA,
- const const_monomial monoB,
+ ConstMonoRef monoA,
+ ConstMonoRef 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);
@@ -43,9 +39,9 @@ auto F4MatrixBuilder::findOrCreateColumn(
MATHICGB_NO_INLINE
void F4MatrixBuilder::createTwoColumns(
- const const_monomial monoA1,
- const const_monomial monoA2,
- const const_monomial monoB,
+ ConstMonoRef monoA1,
+ ConstMonoRef monoA2,
+ ConstMonoRef monoB,
TaskFeeder& feeder
) {
createColumn(monoA1, monoB, feeder);
@@ -56,7 +52,7 @@ F4MatrixBuilder::F4MatrixBuilder(
const PolyBasis& basis,
const size_t memoryQuantum
):
- mTmp(basis.ring().allocMonomial()),
+ mTmp(basis.ring().monoid().alloc()),
mBasis(basis),
mMap(basis.ring()),
mMonomialsLeft(),
@@ -100,27 +96,29 @@ void F4MatrixBuilder::addPolynomialToMatrix(const Poly& poly) {
}
void F4MatrixBuilder::addPolynomialToMatrix
-(const_monomial multiple, const Poly& poly) {
+ (ConstMonoRef multiple, const Poly& poly)
+{
if (poly.isZero())
return;
+ auto desiredLead = monoid().alloc();
+ monoid().multiply(poly.getLeadMonomial(), multiple, desiredLead);
RowTask task = {};
task.addToTop = false;
task.poly = &poly;
- task.desiredLead = ring().allocMonomial();
- ring().monomialMult(poly.getLeadMonomial(), multiple, task.desiredLead);
+ task.desiredLead = desiredLead.release();
MATHICGB_ASSERT(task.sPairPoly == 0);
mTodo.push_back(task);
}
void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
+ MATHICGB_ASSERT(&matrix.ring() == &ring());
MATHICGB_LOG_TIME(F4MatrixBuild) <<
"\n***** Constructing matrix *****\n";
if (mTodo.empty()) {
- matrix = QuadMatrix();
- matrix.ring = &ring();
+ matrix.clear();
return;
}
@@ -131,19 +129,23 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
struct ThreadData {
QuadMatrixBuilder builder;
- monomial tmp1;
- monomial tmp2;
+ MonoRef tmp1;
+ MonoRef tmp2;
};
mgb::mtbb::enumerable_thread_specific<ThreadData> threadData([&](){
- ThreadData data = {QuadMatrixBuilder(
- ring(), mMap, mMonomialsLeft, mMonomialsRight, mBuilder.memoryQuantum()
- )};
- {
- mgb::mtbb::mutex::scoped_lock guard(mCreateColumnLock);
- data.tmp1 = ring().allocMonomial();
- data.tmp2 = ring().allocMonomial();
- }
+ mgb::mtbb::mutex::scoped_lock guard(mCreateColumnLock);
+ ThreadData data = {
+ QuadMatrixBuilder(
+ ring(),
+ mMap,
+ mMonomialsLeft,
+ mMonomialsRight,
+ mBuilder.memoryQuantum()
+ ),
+ *monoid().alloc().release(),
+ *monoid().alloc().release()
+ };
return std::move(data);
});
@@ -156,21 +158,20 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
if (task.sPairPoly != 0) {
MATHICGB_ASSERT(!task.addToTop);
- ring().monomialColons(
+ monoid().colons(
poly.getLeadMonomial(),
task.sPairPoly->getLeadMonomial(),
data.tmp2,
data.tmp1
);
appendRowBottom
- (&poly, data.tmp1, task.sPairPoly, data.tmp2, data.builder, feeder);
+ (poly, data.tmp1, *task.sPairPoly, data.tmp2, data.builder, feeder);
return;
}
- if (task.desiredLead.isNull())
- ring().monomialSetIdentity(data.tmp1);
+ if (task.desiredLead == nullptr)
+ monoid().setIdentity(data.tmp1);
else
- ring().monomialDivide
- (task.desiredLead, poly.getLeadMonomial(), data.tmp1);
+ monoid().divide(poly.getLeadMonomial(), *task.desiredLead, data.tmp1);
if (task.addToTop)
appendRowTop(data.tmp1, *task.poly, builder, feeder);
else
@@ -180,18 +181,17 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
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);
+ for (auto& mono : mTodo)
+ if (!mono.desiredLead.isNull())
+ monoid().freeRaw(*mono.desiredLead.castAwayConst());
mTodo.clear();
auto& builder = threadData.begin()->builder;
- const auto end = threadData.end();
- for (auto it = threadData.begin(); it != end; ++it) {
- if (&it->builder != &builder)
- builder.takeRowsFrom(it->builder.buildMatrixAndClear());
- ring().freeMonomial(it->tmp1);
- ring().freeMonomial(it->tmp2);
+ for (auto& data : threadData) {
+ if (&data.builder != &builder)
+ builder.takeRowsFrom(data.builder.buildMatrixAndClear());
+ monoid().freeRaw(data.tmp1);
+ monoid().freeRaw(data.tmp2);
}
matrix = builder.buildMatrixAndClear();
threadData.clear();
@@ -203,15 +203,15 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
const auto end = reader.end();
for (auto it = reader.begin(); it != end; ++it) {
const auto p = *it;
- monomial copy = ring().allocMonomial();
- ring().monoid().copy(p.second, copy);
+ auto copy = monoid().alloc();
+ monoid().copy(p.second, copy);
auto& monos = p.first.left() ?
matrix.leftColumnMonomials : matrix.rightColumnMonomials;
const auto index = p.first.index();
if (monos.size() <= index)
monos.resize(index + 1);
MATHICGB_ASSERT(monos[index].isNull());
- monos[index] = copy;
+ monos[index] = copy.release();
}
}
#ifdef MATHICGB_DEBUG
@@ -228,13 +228,10 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
}
auto F4MatrixBuilder::createColumn(
- const const_monomial monoA,
- const const_monomial monoB,
+ ConstMonoRef monoA,
+ ConstMonoRef monoB,
TaskFeeder& feeder
) -> std::pair<F4MatrixBuilder::LeftRightColIndex, ConstMonoRef> {
- 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
{
@@ -244,12 +241,12 @@ auto F4MatrixBuilder::createColumn(
}
// The column really does not exist, so we need to create it
- ring().monomialMult(monoA, monoB, mTmp);
- if (!ring().monomialHasAmpleCapacity(mTmp))
+ monoid().multiply(monoA, monoB, mTmp);
+ if (!monoid().hasAmpleCapacity(mTmp))
mathic::reportError("Monomial exponent overflow in F4MatrixBuilder.");
// look for a reducer of mTmp
- const size_t reducerIndex = mBasis.classicReducer(mTmp);
+ const size_t reducerIndex = mBasis.classicReducer(Monoid::toOld(*mTmp));
const bool insertLeft = (reducerIndex != static_cast<size_t>(-1));
// Create the new left or right column
@@ -257,26 +254,25 @@ auto F4MatrixBuilder::createColumn(
if (colCount == std::numeric_limits<ColIndex>::max())
throw std::overflow_error("Too many columns in QuadMatrix");
const auto inserted = mMap.insert
- (std::make_pair(mTmp, LeftRightColIndex(colCount, insertLeft)));
+ (std::make_pair(mTmp.ptr(), LeftRightColIndex(colCount, insertLeft)));
++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 = ptr;
+ task.desiredLead = inserted.first.second;
feeder.add(task);
}
- return std::make_pair(*inserted.first.first, Monoid::toRef(ptr));
+ return std::make_pair(*inserted.first.first, *inserted.first.second);
}
void F4MatrixBuilder::appendRowBottom(
- const_monomial multiple,
+ ConstMonoRef multiple,
const bool negate,
const Poly::const_iterator begin,
const Poly::const_iterator end,
@@ -284,7 +280,6 @@ void F4MatrixBuilder::appendRowBottom(
TaskFeeder& feeder
) {
// todo: eliminate the code-duplication between here and appendRowTop.
- MATHICGB_ASSERT(!multiple.isNull());
MATHICGB_ASSERT(&builder != 0);
auto it = begin;
@@ -311,12 +306,11 @@ updateReader:
}
void F4MatrixBuilder::appendRowTop(
- const const_monomial multiple,
+ ConstMonoRef multiple,
const Poly& poly,
QuadMatrixBuilder& builder,
TaskFeeder& feeder
) {
- MATHICGB_ASSERT(!multiple.isNull());
MATHICGB_ASSERT(&poly != 0);
MATHICGB_ASSERT(&builder != 0);
@@ -339,14 +333,14 @@ updateReader:
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();
+ 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 const_monomial mono2 = it2.getMonomial();
+ const auto mono2 = it2.mono();
const auto colPair = colMap.findTwoProducts(mono1, mono2, multiple);
if (colPair.first == 0 || colPair.second == 0) {
@@ -362,23 +356,20 @@ updateReader:
}
void F4MatrixBuilder::appendRowBottom(
- const Poly* poly,
- monomial multiply,
- const Poly* sPairPoly,
- monomial sPairMultiply,
+ const Poly& poly,
+ ConstMonoRef multiply,
+ const Poly& sPairPoly,
+ ConstMonoRef sPairMultiply,
QuadMatrixBuilder& builder,
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(!poly.isZero());
+ 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();
+ MATHICGB_ASSERT(!sPairPoly.isZero());
+ Poly::const_iterator itB = sPairPoly.begin();
+ Poly::const_iterator endB = sPairPoly.end();
// skip leading terms since they cancel
MATHICGB_ASSERT(itA.getCoefficient() == itB.getCoefficient());
@@ -387,8 +378,8 @@ void F4MatrixBuilder::appendRowBottom(
const ColReader colMap(mMap);
- const const_monomial mulA = multiply;
- const const_monomial mulB = sPairMultiply;
+ const auto mulA = multiply;
+ const auto mulB = sPairMultiply;
while (true) {
// Watch out: we are depending on appendRowBottom to finish the row, so
// if you decide not to call that function in case
@@ -408,7 +399,7 @@ void F4MatrixBuilder::appendRowBottom(
(itA.getMonomial(), mulA, colMap, feeder);
const auto colB = findOrCreateColumn
(itB.getMonomial(), mulB, colMap, feeder);
- const auto cmp = ring().monoid().compare(colA.second, colB.second);
+ const auto cmp = 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 6e6ac90..a1ab4d7 100755
--- a/src/mathicgb/F4MatrixBuilder.hpp
+++ b/src/mathicgb/F4MatrixBuilder.hpp
@@ -26,7 +26,7 @@ private:
typedef QuadMatrixBuilder::LeftRightColIndex LeftRightColIndex;
typedef QuadMatrixBuilder::Scalar Scalar;
typedef QuadMatrixBuilder::Map Map;
- typedef QuadMatrixBuilder::MonomialsType Monomials;
+ typedef QuadMatrixBuilder::Monomials Monomials;
public:
typedef PolyRing::Monoid Monoid;
@@ -55,7 +55,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, just letting multiple be the
/// identity.
@@ -78,6 +78,7 @@ public:
void buildMatrixAndClear(QuadMatrix& matrix);
const PolyRing& ring() const {return mBuilder.ring();}
+ const Monoid& monoid() const {return mBuilder.ring().monoid();}
private:
typedef const MonomialMap<LeftRightColIndex>::Reader ColReader;
@@ -88,40 +89,38 @@ private:
/// where sPairMultiply makes the lead terms cancel.
struct RowTask {
bool addToTop; // add the row to the bottom if false
- monomial desiredLead; // multiply monomial onto poly to get this lead
+ ConstMonoPtr desiredLead; // multiply monomial onto poly to get this lead
const Poly* poly;
const Poly* sPairPoly;
- monomial sPairMultiply;
+ ConstMonoPtr sPairMultiply;
};
- typedef ::mgb::mtbb::parallel_do_feeder<RowTask> TaskFeeder;
+ typedef mtbb::parallel_do_feeder<RowTask> TaskFeeder;
/// Creates a column with monomial label x and schedules a new row to
/// 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, ConstMonoRef>
- createColumn(
- const_monomial monoA,
- const_monomial monoB,
- TaskFeeder& feeder
- );
+ createColumn(ConstMonoRef monoA, ConstMonoRef monoB, TaskFeeder& feeder);
void appendRowTop(
- const_monomial multiple,
+ ConstMonoRef multiple,
const Poly& poly,
QuadMatrixBuilder& builder,
TaskFeeder& feeder
);
+
void appendRowBottom(
- const Poly* poly,
- monomial multiply,
- const Poly* sPairPoly,
- monomial sPairMultiply,
+ const Poly& poly,
+ ConstMonoRef multiply,
+ const Poly& sPairPoly,
+ ConstMonoRef sPairMultiply,
QuadMatrixBuilder& builder,
TaskFeeder& feeder
);
+
void appendRowBottom(
- const_monomial multiple,
+ ConstMonoRef multiple,
bool negate,
Poly::const_iterator begin,
Poly::const_iterator end,
@@ -132,32 +131,32 @@ private:
MATHICGB_NO_INLINE
std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonoRef>
findOrCreateColumn(
- const_monomial monoA,
- const_monomial monoB,
+ ConstMonoRef monoA,
+ ConstMonoRef monoB,
TaskFeeder& feeder
);
MATHICGB_INLINE
std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonoRef>
findOrCreateColumn(
- const_monomial monoA,
- const_monomial monoB,
+ ConstMonoRef monoA,
+ ConstMonoRef monoB,
const ColReader& colMap,
TaskFeeder& feeder
);
MATHICGB_NO_INLINE
void createTwoColumns(
- const const_monomial monoA1,
- const const_monomial monoA2,
- const const_monomial monoB,
+ ConstMonoRef monoA1,
+ ConstMonoRef monoA2,
+ ConstMonoRef monoB,
TaskFeeder& feeder
);
mgb::mtbb::mutex mCreateColumnLock;
ColIndex mLeftColCount;
ColIndex mRightColCount;
- monomial mTmp;
+ Mono mTmp;
const PolyBasis& mBasis;
Monomials mMonomialsLeft;
Monomials mMonomialsRight;
diff --git a/src/mathicgb/F4MatrixBuilder2.cpp b/src/mathicgb/F4MatrixBuilder2.cpp
index 55e3a41..1588abf 100755
--- a/src/mathicgb/F4MatrixBuilder2.cpp
+++ b/src/mathicgb/F4MatrixBuilder2.cpp
@@ -127,12 +127,12 @@ public:
}
void buildMatrixAndClear(std::vector<RowTask>& tasks, QuadMatrix& quadMatrix) {
+ MATHICGB_ASSERT(&quadMatrix.ring() == &ring());
MATHICGB_LOG_TIME(F4MatrixBuild2) <<
"\n***** Constructing matrix *****\n";
if (tasks.empty()) {
- quadMatrix = QuadMatrix();
- quadMatrix.ring = &ring();
+ quadMatrix.clear();
return;
}
@@ -259,7 +259,7 @@ public:
const PolyRing& ring() const {return mBasis.ring();}
const Monoid& monoid() const {return mBasis.ring().monoid();}
- Builder(const PolyBasis& basis, const size_t memoryQuantum):
+ Builder(const PolyBasis& basis, const size_t memoryQuantum):
mMemoryQuantum(memoryQuantum),
mTmp(basis.ring().monoid().alloc()),
mBasis(basis),
@@ -372,7 +372,7 @@ public:
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));
+ const auto colPair = colMap.findTwoProducts(mono1, mono2, multiple);
if (colPair.first == 0 || colPair.second == 0) {
createColumn(mono1, multiple, feeder);
createColumn(mono2, multiple, feeder);
diff --git a/src/mathicgb/F4MatrixProjection.cpp b/src/mathicgb/F4MatrixProjection.cpp
index a59321d..5e186a5 100755
--- a/src/mathicgb/F4MatrixProjection.cpp
+++ b/src/mathicgb/F4MatrixProjection.cpp
@@ -269,8 +269,7 @@ done:;
bottom.appendRowsPermuted(tb.moveBottom());
// Move the data into place
- QuadMatrix qm;
- qm.ring = &ring();
+ QuadMatrix qm(ring());
qm.leftColumnMonomials = std::move(mLeftMonomials);
qm.rightColumnMonomials = std::move(mRightMonomials);
@@ -346,11 +345,10 @@ QuadMatrix F4MatrixProjection::makeAndClearTwoStep(const size_t quantum) {
}
MATHICGB_ASSERT(tb.debugAssertValid());
- QuadMatrix qm;
+ QuadMatrix qm(ring());
auto left = projectRows(tb, quantum, lr.moveLeft());
auto right = projectRows(tb, quantum, lr.moveRight());
- qm.ring = &ring();
qm.topLeft = std::move(left.first);
qm.bottomLeft = std::move(left.second);
qm.topRight = std::move(right.first);
diff --git a/src/mathicgb/F4MatrixProjection.hpp b/src/mathicgb/F4MatrixProjection.hpp
index 16eebed..cddcd2d 100755
--- a/src/mathicgb/F4MatrixProjection.hpp
+++ b/src/mathicgb/F4MatrixProjection.hpp
@@ -14,6 +14,13 @@ MATHICGB_NAMESPACE_BEGIN
class F4MatrixProjection {
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 SparseMatrix::RowIndex RowIndex;
typedef SparseMatrix::ColIndex ColIndex;
typedef SparseMatrix::Scalar Scalar;
@@ -48,8 +55,8 @@ private:
std::vector<ColProjectTo> mColProjectTo;
std::vector<F4ProtoMatrix*> mMatrices;
- std::vector<monomial> mLeftMonomials;
- std::vector<monomial> mRightMonomials;
+ std::vector<ConstMonoPtr> mLeftMonomials;
+ std::vector<ConstMonoPtr> mRightMonomials;
const PolyRing& mRing;
};
diff --git a/src/mathicgb/F4MatrixReducer.cpp b/src/mathicgb/F4MatrixReducer.cpp
index 8a8cc27..977bd48 100755
--- a/src/mathicgb/F4MatrixReducer.cpp
+++ b/src/mathicgb/F4MatrixReducer.cpp
@@ -201,10 +201,8 @@ namespace {
const SparseMatrix& reduceByRight = qm.topRight;
const auto leftColCount = qm.computeLeftColCount();
-// static_cast<SparseMatrix::ColIndex>(qm.leftColumnMonomials.size());
const auto rightColCount =
static_cast<SparseMatrix::ColIndex>(qm.computeRightColCount());
-// static_cast<SparseMatrix::ColIndex>(qm.rightColumnMonomials.size());
MATHICGB_ASSERT(leftColCount == reduceByLeft.rowCount());
const auto pivotCount = leftColCount;
const auto rowCount = toReduceLeft.rowCount();
@@ -488,7 +486,6 @@ namespace {
denseRow.makeUnitary(modulus, col);
}
}});
- //std::cout << "done reducing that batch" << std::endl;
reduced.clear();
std::sort(nextReducers.begin(), nextReducers.end());
@@ -553,36 +550,6 @@ namespace {
}
}
-// todo: use auto instead of these typedefs where possible/reasonable
-// todo: do SparseMatrix::Scalar instead of Scalar and remove this typedef :: DONE
-//typedef SparseMatrix::Scalar Scalar; :: DONE
-//typedef SparseMatrix::RowIndex RowIndex; // todo: same :: DONE
-//typedef SparseMatrix::ColIndex ColIndex; // todo: same :: DONE
-
-//Scalar modPrime = 11; // todo: remove this variable :: DONE
-/*
-class SharwanMatrix {
-public:
-// typedefs for scalar, row index and col index
-// typedef Row to be your representation of a row
-
- const Scalar& operator[](RowIndex row, ColIndex col) const {return mMatrix[row][col];}
- Scalar& operator[](RowIndex row, ColIndex col) {return mMatrix[row][col];}
-
- Row& operator[](RowIndex) {}
- const Row& operator[](RowIndex) const {}
-
- // example of setter. Do not make a setter for modulus, row index or col index. No setters, except for entries of the matrix.
- void setX(int value) {mX = value;}
-
-
-// store matrix, modulus, rowCount and colCount
-// accessor for getting modulus: modulus()
-private:
- int mX; // todo: remove, just example
- // all member variables go here. member x is written mX.
-};
-*/
void addRowMultipleInplace(
std::vector< std::vector<SparseMatrix::Scalar> >& matrix,
const SparseMatrix::RowIndex addRow,
@@ -621,11 +588,8 @@ void makeRowUnitary(
auto multiply = modularInverse(leadingScalar, modulus);
for(SparseMatrix::ColIndex col = leadingCol; col < colCount; ++col)
matrix[row][col] = modularProduct(matrix[row][col], multiply, modulus);
-
- // todo: use modularProduct on above line ::DONE
}
-// todo: make this take a parameter startAtCol ::DONE
SparseMatrix::ColIndex leadingColumn(
const std::vector< std::vector<SparseMatrix::Scalar>>& matrix,
const SparseMatrix::RowIndex row,
@@ -707,9 +671,7 @@ SparseMatrix reduceToEchelonFormShrawan(
}
}
- // todo: make modPrime a parameter and rename it to modulus. :: DONE
- // modPrime = modulus; :: DONE
- rowReducedEchelonMatrix(matrix, colCount, modulus);
+ rowReducedEchelonMatrix(matrix, colCount, modulus);
// convert reduced matrix to SparseMatrix.
SparseMatrix reduced;
@@ -731,7 +693,6 @@ SparseMatrix reduceToEchelonFormShrawanDelayedModulus(
const SparseMatrix& toReduce,
SparseMatrix::Scalar modulus
) {
- // todo: implement delayed modulus
const SparseMatrix::RowIndex rowCount = toReduce.rowCount();
const auto colCount = toReduce.computeColCount();
@@ -793,8 +754,8 @@ SparseMatrix F4MatrixReducer::reducedRowEchelonForm(
else
return reduceToEchelonFormShrawan(matrix, mModulus);
} else {
- // todo: actually do some work to determine a good way to determine
- // when to use the sparse method, or alternatively make some some
+ // todo: actually do some work to find a good way to determine
+ // when to use the sparse method, or alternatively make some
// sort of hybrid.
if (matrix.computeDensity() < 0.02)
return reduceToEchelonFormSparse(matrix, mModulus);
diff --git a/src/mathicgb/F4MatrixReducer.hpp b/src/mathicgb/F4MatrixReducer.hpp
index 29e2d42..904ea81 100755
--- a/src/mathicgb/F4MatrixReducer.hpp
+++ b/src/mathicgb/F4MatrixReducer.hpp
@@ -20,8 +20,6 @@ class PolyRing;
class F4MatrixReducer {
public:
/// The ring used is Z/pZ where modulus is the prime p.
- ///
- ///
F4MatrixReducer(coefficient modulus);
/// Reduces the bottom rows by the top rows and returns the bottom right
diff --git a/src/mathicgb/F4Reducer.cpp b/src/mathicgb/F4Reducer.cpp
index 491c662..93e26b3 100755
--- a/src/mathicgb/F4Reducer.cpp
+++ b/src/mathicgb/F4Reducer.cpp
@@ -8,6 +8,7 @@
#include "F4MatrixReducer.hpp"
#include "QuadMatrix.hpp"
#include "LogDomain.hpp"
+#include "CFile.hpp"
#include <iostream>
#include <limits>
@@ -101,6 +102,9 @@ public:
virtual std::string description() const;
virtual size_t getMemoryUse() const;
+ const PolyRing& ring() const {return mRing;}
+ const Monoid& monoid() const {return mRing.monoid();}
+
private:
void saveMatrix(const QuadMatrix& matrix);
@@ -140,8 +144,7 @@ std::unique_ptr<Poly> F4Reducer::classicReduce
std::cerr <<
"F4Reducer: Using fall-back reducer for single classic reduction\n";
- auto p = mFallback->classicReduce(poly, basis);
- return std::move(p);
+ return mFallback->classicReduce(poly, basis);
}
std::unique_ptr<Poly> F4Reducer::classicTailReduce
@@ -151,8 +154,7 @@ std::unique_ptr<Poly> F4Reducer::classicTailReduce
std::cerr <<
"F4Reducer: Using fall-back reducer for single classic tail reduction\n";
- auto p = mFallback->classicTailReduce(poly, basis);
- return std::move(p);
+ return mFallback->classicTailReduce(poly, basis);
}
std::unique_ptr<Poly> F4Reducer::classicReduceSPoly(
@@ -163,8 +165,7 @@ std::unique_ptr<Poly> F4Reducer::classicReduceSPoly(
if (tracingLevel >= 2)
std::cerr << "F4Reducer: "
"Using fall-back reducer for single classic S-pair reduction\n";
- auto p = mFallback->classicReduceSPoly(a, b, basis);
- return std::move(p);
+ return mFallback->classicReduceSPoly(a, b, basis);
}
void F4Reducer::classicReduceSPolySet(
@@ -172,7 +173,7 @@ void F4Reducer::classicReduceSPolySet(
const PolyBasis& basis,
std::vector<std::unique_ptr<Poly>>& reducedOut
) {
- if (spairs.size() <= 1 && false) {
+ if (spairs.size() <= 1) {
if (tracingLevel >= 2)
std::cerr << "F4Reducer: Using fall-back reducer for "
<< spairs.size() << " S-pairs.\n";
@@ -182,27 +183,25 @@ void F4Reducer::classicReduceSPolySet(
reducedOut.clear();
MATHICGB_ASSERT(!spairs.empty());
- if (tracingLevel >= 2 && false)
+ if (tracingLevel >= 2)
std::cerr << "F4Reducer: Reducing " << spairs.size() << " S-polynomials.\n";
SparseMatrix reduced;
- std::vector<monomial> monomials;
+ QuadMatrix::Monomials monomials;
{
- QuadMatrix qm;
+ QuadMatrix qm(basis.ring());
{
if (mType == OldType) {
F4MatrixBuilder builder(basis, mMemoryQuantum);
- for (auto it = spairs.begin(); it != spairs.end(); ++it) {
+ for (const auto& spair : spairs)
builder.addSPolynomialToMatrix
- (basis.poly(it->first), basis.poly(it->second));
- }
+ (basis.poly(spair.first), basis.poly(spair.second));
builder.buildMatrixAndClear(qm);
} else {
F4MatrixBuilder2 builder(basis, mMemoryQuantum);
- for (auto it = spairs.begin(); it != spairs.end(); ++it) {
+ for (const auto& spair : spairs)
builder.addSPolynomialToMatrix
- (basis.poly(it->first), basis.poly(it->second));
- }
+ (basis.poly(spair.first), basis.poly(spair.second));
builder.buildMatrixAndClear(qm);
}
}
@@ -214,9 +213,8 @@ void F4Reducer::classicReduceSPolySet(
reduced = F4MatrixReducer(basis.ring().charac()).
reducedRowEchelonFormBottomRight(qm);
monomials = std::move(qm.rightColumnMonomials);
- const auto end = qm.leftColumnMonomials.end();
- for (auto it = qm.leftColumnMonomials.begin(); it != end; ++it)
- mRing.freeMonomial(*it);
+ for (auto& mono : qm.leftColumnMonomials)
+ monoid().freeRaw(mono.castAwayConst());
}
for (SparseMatrix::RowIndex row = 0; row < reduced.rowCount(); ++row) {
@@ -224,9 +222,8 @@ void F4Reducer::classicReduceSPolySet(
reduced.rowToPolynomial(row, monomials, *p);
reducedOut.push_back(std::move(p));
}
- const auto end = monomials.end();
- for (auto it = monomials.begin(); it != end; ++it)
- mRing.freeMonomial(*it);
+ for (auto& mono : monomials)
+ monoid().freeRaw(mono.castAwayConst());
}
void F4Reducer::classicReducePolySet(
@@ -234,7 +231,7 @@ void F4Reducer::classicReducePolySet(
const PolyBasis& basis,
std::vector< std::unique_ptr<Poly> >& reducedOut
) {
- if (polys.size() <= 1 && false) {
+ if (polys.size() <= 1) {
if (tracingLevel >= 2)
std::cerr << "F4Reducer: Using fall-back reducer for "
<< polys.size() << " polynomials.\n";
@@ -245,23 +242,23 @@ void F4Reducer::classicReducePolySet(
reducedOut.clear();
MATHICGB_ASSERT(!polys.empty());
- if (tracingLevel >= 2 && false)
+ if (tracingLevel >= 2)
std::cerr << "F4Reducer: Reducing " << polys.size() << " polynomials.\n";
SparseMatrix reduced;
- std::vector<monomial> monomials;
+ QuadMatrix::Monomials monomials;
{
- QuadMatrix qm;
+ QuadMatrix qm(ring());
{
if (mType == OldType) {
F4MatrixBuilder builder(basis, mMemoryQuantum);
- for (auto it = polys.begin(); it != polys.end(); ++it)
- builder.addPolynomialToMatrix(**it);
+ for (const auto& poly : polys)
+ builder.addPolynomialToMatrix(*poly);
builder.buildMatrixAndClear(qm);
} else {
F4MatrixBuilder2 builder(basis, mMemoryQuantum);
- for (auto it = polys.begin(); it != polys.end(); ++it)
- builder.addPolynomialToMatrix(**it);
+ for (const auto& poly : polys)
+ builder.addPolynomialToMatrix(*poly);
builder.buildMatrixAndClear(qm);
}
}
@@ -273,9 +270,8 @@ void F4Reducer::classicReducePolySet(
reduced = F4MatrixReducer(basis.ring().charac()).
reducedRowEchelonFormBottomRight(qm);
monomials = std::move(qm.rightColumnMonomials);
- for (auto it = qm.leftColumnMonomials.begin();
- it != qm.leftColumnMonomials.end(); ++it)
- mRing.freeMonomial(*it);
+ for (auto& mono : qm.leftColumnMonomials)
+ monoid().freeRaw(mono.castAwayConst());
}
if (tracingLevel >= 2 && false)
@@ -287,9 +283,8 @@ void F4Reducer::classicReducePolySet(
reduced.rowToPolynomial(row, monomials, *p);
reducedOut.push_back(std::move(p));
}
- const auto end = monomials.end();
- for (auto it = monomials.begin(); it != end; ++it)
- mRing.freeMonomial(*it);
+ for (auto& mono : monomials)
+ monoid().freeRaw(mono.castAwayConst());
}
std::unique_ptr<Poly> F4Reducer::regularReduce(
@@ -328,10 +323,10 @@ void F4Reducer::saveMatrix(const QuadMatrix& matrix) {
fileName << mStoreToFile << '-' << mMatrixSaveCount << ".qmat";
if (tracingLevel > 2)
std::cerr << "F4Reducer: Saving matrix to " << fileName.str() << '\n';
- FILE* file = fopen(fileName.str().c_str(), "wb");
- // @todo: fix leak of file on exception
- matrix.write(static_cast<SparseMatrix::Scalar>(mRing.charac()), file);
- fclose(file);
+
+ CFile file(fileName.str(), "wb");
+ matrix.write
+ (static_cast<SparseMatrix::Scalar>(mRing.charac()), file.handle());
}
std::unique_ptr<Reducer> makeF4Reducer(
diff --git a/src/mathicgb/MonoMonoid.hpp b/src/mathicgb/MonoMonoid.hpp
index ee9c08d..a1b65d0 100755
--- a/src/mathicgb/MonoMonoid.hpp
+++ b/src/mathicgb/MonoMonoid.hpp
@@ -1167,8 +1167,16 @@ public:
}
Mono alloc() const {return mPool.alloc();}
+
void free(Mono&& mono) const {mPool.free(std::move(mono));}
+
void freeRaw(MonoRef mono) const {mPool.freeRaw(mono);}
+
+ void freeRaw(MonoPtr mono) const {
+ if (mono != nullptr)
+ freeRaw(*mono);
+ }
+
bool fromPool(ConstMonoRef mono) const {mPool.fromPool(mono);}
// *** Classes for holding and referring to monomials
diff --git a/src/mathicgb/MonomialMap.hpp b/src/mathicgb/MonomialMap.hpp
index 948c816..d50e689 100755
--- a/src/mathicgb/MonomialMap.hpp
+++ b/src/mathicgb/MonomialMap.hpp
@@ -120,9 +120,9 @@ public:
/// loop.
MATHICGB_INLINE
std::pair<const mapped_type*, const mapped_type*> findTwoProducts(
- const const_monomial a1,
- const const_monomial a2,
- const const_monomial b
+ const ConstMonoRef a1,
+ const ConstMonoRef a2,
+ const ConstMonoRef b
) const {
return mMap.findTwoProducts(a1, a2, b);
}
diff --git a/src/mathicgb/PolyRing.hpp b/src/mathicgb/PolyRing.hpp
index c3781f3..5fe70de 100755
--- a/src/mathicgb/PolyRing.hpp
+++ b/src/mathicgb/PolyRing.hpp
@@ -299,7 +299,7 @@ public:
// Only call this method for monomials returned by allocMonomial().
void freeMonomial(Monomial m) const {
- monoid().freeRaw(m);
+ monoid().freeRaw(Monoid::MonoPtr(m));
}
// Free monomials allocated here by calling freeMonomial().
diff --git a/src/mathicgb/QuadMatrix.cpp b/src/mathicgb/QuadMatrix.cpp
index ddc8ff2..9c45a74 100755
--- a/src/mathicgb/QuadMatrix.cpp
+++ b/src/mathicgb/QuadMatrix.cpp
@@ -3,6 +3,7 @@
#include "stdinc.h"
#include "QuadMatrix.hpp"
+#include "MathicIO.hpp"
#include "mtbb.hpp"
#include <mathic.h>
#include <ostream>
@@ -43,17 +44,15 @@ void QuadMatrix::print(std::ostream& out) const {
// column monomials
out << "Left columns:";
- const auto leftColCount = leftColumnMonomials.size();
- for (ColIndex leftCol = 0; leftCol < leftColCount; ++leftCol) {
+ for (const auto& mono : leftColumnMonomials) {
out << ' ';
- ring->monomialDisplay(out, leftColumnMonomials[leftCol], false, true);
+ MathicIO<>().writeMonomial(monoid(), false, *mono, out);
}
out << "\nRight columns:";
- const auto rightColCount = rightColumnMonomials.size();
- for (ColIndex rightCol = 0; rightCol < rightColCount; ++rightCol) {
+ for (const auto& mono : rightColumnMonomials) {
out << ' ';
- ring->monomialDisplay(out, rightColumnMonomials[rightCol], false, true);
+ MathicIO<>().writeMonomial(monoid(), false, *mono, out);
}
out << '\n';
@@ -211,7 +210,7 @@ QuadMatrix QuadMatrix::toCanonical() const {
const auto rightColCount = rightColumnMonomials.size();
// todo: eliminate left/right code duplication here
- QuadMatrix matrix;
+ QuadMatrix matrix(ring());
{ // left side
std::vector<SparseMatrix::RowIndex> rows;
for (SparseMatrix::RowIndex row = 0; row < topLeft.rowCount(); ++row)
@@ -247,7 +246,6 @@ QuadMatrix QuadMatrix::toCanonical() const {
matrix.leftColumnMonomials = leftColumnMonomials;
matrix.rightColumnMonomials = rightColumnMonomials;
- matrix.ring = ring;
return std::move(matrix);
}
@@ -258,34 +256,38 @@ std::ostream& operator<<(std::ostream& out, const QuadMatrix& qm) {
}
namespace {
+ template<class Monoid>
class ColumnComparer {
public:
- ColumnComparer(const PolyRing& ring): mRing(ring) {}
+ ColumnComparer(const Monoid& monoid): mMonoid(monoid) {}
- typedef SparseMatrix::ColIndex ColIndex;
- typedef std::pair<monomial, ColIndex> Pair;
+ template<class Pair>
bool operator()(const Pair& a, const Pair b) const {
- return mRing.monomialLT(b.first, a.first);
+ return mMonoid.lessThan(*b.first, *a.first);
}
private:
- const PolyRing& mRing;
+ const Monoid& mMonoid;
};
std::vector<SparseMatrix::ColIndex> sortColumnMonomialsAndMakePermutation(
- std::vector<monomial>& monomials,
- const PolyRing& ring
+ QuadMatrix::Monomials& monomials,
+ const QuadMatrix::Monoid& monoid
) {
typedef SparseMatrix::ColIndex ColIndex;
MATHICGB_ASSERT(monomials.size() <= std::numeric_limits<ColIndex>::max());
const ColIndex colCount = static_cast<ColIndex>(monomials.size());
// Monomial needs to be non-const as we are going to put these
// monomials back into the vector of monomials which is not const.
- std::vector< std::pair<monomial, ColIndex>> columns;
+ std::vector<std::pair<QuadMatrix::ConstMonoPtr, ColIndex>> columns;
columns.reserve(colCount);
for (ColIndex col = 0; col < colCount; ++col)
columns.push_back(std::make_pair(monomials[col], col));
- std::sort(columns.begin(), columns.end(), ColumnComparer(ring));
+ std::sort(
+ columns.begin(),
+ columns.end(),
+ ColumnComparer<QuadMatrix::Monoid>(monoid)
+ );
// Apply sorting permutation to monomials. This is why it is necessary to
// copy the values in monomial out of there: in-place application of a
@@ -293,8 +295,10 @@ namespace {
MATHICGB_ASSERT(columns.size() == colCount);
MATHICGB_ASSERT(monomials.size() == colCount);
for (size_t col = 0; col < colCount; ++col) {
- MATHICGB_ASSERT(col == 0 ||
- ring.monomialLT(columns[col].first, columns[col - 1].first));
+ MATHICGB_ASSERT(
+ col == 0 ||
+ monoid.lessThan(*columns[col].first, *columns[col - 1].first)
+ );
monomials[col] = columns[col].first;
}
@@ -318,10 +322,10 @@ void QuadMatrix::sortColumnsLeftRightParallel() {
mgb::mtbb::parallel_for(0, 2, 1, [&](int i) {
if (i == 0)
leftPermutation =
- sortColumnMonomialsAndMakePermutation(leftColumnMonomials, *ring);
+ sortColumnMonomialsAndMakePermutation(leftColumnMonomials, monoid());
else
rightPermutation =
- sortColumnMonomialsAndMakePermutation(rightColumnMonomials, *ring);
+ sortColumnMonomialsAndMakePermutation(rightColumnMonomials, monoid());
});
mgb::mtbb::parallel_for(0, 4, 1, [&](int i) {
@@ -354,7 +358,6 @@ SparseMatrix::Scalar QuadMatrix::read(FILE* file) {
leftColumnMonomials.clear();
rightColumnMonomials.clear();
- ring = 0;
const auto topLeftModulus = topLeft.read(file);
const auto topRightModulus = topRight.read(file);
diff --git a/src/mathicgb/QuadMatrix.hpp b/src/mathicgb/QuadMatrix.hpp
index 91485e2..d784064 100755
--- a/src/mathicgb/QuadMatrix.hpp
+++ b/src/mathicgb/QuadMatrix.hpp
@@ -15,12 +15,17 @@ MATHICGB_NAMESPACE_BEGIN
/// into one matrix divided into top left, top right, bottom left and
/// bottom right. This is a convenient representation of the matrices
/// encountered in the F4 polynomial reduction algorithm.
-///
-/// So far this is just a data class used as the output of a
-/// QuadMatrixBuilder or F4MatrixBuilder.
class QuadMatrix {
public:
- QuadMatrix() {}
+ typedef PolyRing::Monoid Monoid;
+ typedef Monoid::Mono Mono;
+ typedef Monoid::MonoRef MonoRef;
+ typedef Monoid::ConstMonoRef ConstMonoRef;
+ typedef Monoid::MonoPtr MonoPtr;
+ typedef Monoid::ConstMonoPtr ConstMonoPtr;
+
+ QuadMatrix(): mRing(nullptr) {}
+ QuadMatrix(const PolyRing& ring): mRing(&ring) {}
QuadMatrix(QuadMatrix&& matrix):
topLeft(std::move(matrix.topLeft)),
@@ -29,22 +34,28 @@ public:
bottomRight(std::move(matrix.bottomRight)),
leftColumnMonomials(std::move(matrix.leftColumnMonomials)),
rightColumnMonomials(std::move(matrix.rightColumnMonomials)),
- ring(matrix.ring)
+ mRing(&matrix.ring())
{}
QuadMatrix& operator=(QuadMatrix&& matrix) {
+ MATHICGB_ASSERT(mRing == matrix.mRing);
this->~QuadMatrix();
new (this) QuadMatrix(std::move(matrix));
return *this;
}
+ void clear() {
+ *this = QuadMatrix(ring());
+ }
+
+ typedef std::vector<ConstMonoPtr> Monomials;
+
SparseMatrix topLeft;
SparseMatrix topRight;
SparseMatrix bottomLeft;
SparseMatrix bottomRight;
- std::vector<monomial> leftColumnMonomials;
- std::vector<monomial> rightColumnMonomials;
- const PolyRing* ring;
+ Monomials leftColumnMonomials;
+ Monomials rightColumnMonomials;
/// Prints whole matrix to out in human-readable format. Useful for
/// debugging.
@@ -87,9 +98,14 @@ public:
/// Asserts internal invariants if asserts are turned on.
bool debugAssertValid() const;
+ const PolyRing& ring() const {return *mRing;}
+ const Monoid& monoid() const {return ring().monoid();}
+
private:
QuadMatrix(const QuadMatrix&); // not available
void operator=(const QuadMatrix&); // not available
+
+ const PolyRing* const mRing;
};
std::ostream& operator<<(std::ostream& out, const QuadMatrix& qm);
diff --git a/src/mathicgb/QuadMatrixBuilder.cpp b/src/mathicgb/QuadMatrixBuilder.cpp
index 225b9f6..805e868 100755
--- a/src/mathicgb/QuadMatrixBuilder.cpp
+++ b/src/mathicgb/QuadMatrixBuilder.cpp
@@ -12,8 +12,8 @@ MATHICGB_NAMESPACE_BEGIN
QuadMatrixBuilder::QuadMatrixBuilder(
const PolyRing& ring,
Map& map,
- MonomialsType& monomialsLeft,
- MonomialsType& monomialsRight,
+ Monomials& monomialsLeft,
+ Monomials& monomialsRight,
const size_t memoryQuantum
):
mMonomialToCol(map),
@@ -26,7 +26,7 @@ QuadMatrixBuilder::QuadMatrixBuilder(
{}
void QuadMatrixBuilder::takeRowsFrom(QuadMatrix&& matrix) {
- MATHICGB_ASSERT(&ring() == matrix.ring);
+ MATHICGB_ASSERT(&ring() == &matrix.ring());
MATHICGB_ASSERT(matrix.debugAssertValid());
mTopLeft.takeRowsFrom(std::move(matrix.topLeft));
@@ -59,7 +59,8 @@ namespace {
if (colCount == std::numeric_limits<QuadMatrixBuilder::ColIndex>::max())
throw std::overflow_error("Too many columns in QuadMatrixBuilder");
- toMonomial.push_back(0); // allocate memory now to avoid bad_alloc later
+ // allocate memory in toMonomial now to avoid bad_alloc later
+ toMonomial.emplace_back(nullptr);
monomial copied = ring.allocMonomial();
ring.monomialCopy(mono, copied);
std::pair<QuadMatrixBuilder::LeftRightColIndex, ConstMonomial> p;
@@ -114,13 +115,12 @@ QuadMatrixBuilder::createColumnRight(
}
QuadMatrix QuadMatrixBuilder::buildMatrixAndClear() {
- QuadMatrix out;
+ QuadMatrix out(ring());
mTopLeft.swap(out.topLeft);
mTopRight.swap(out.topRight);
mBottomLeft.swap(out.bottomLeft);
mBottomRight.swap(out.bottomRight);
- out.ring = &ring();
mTopLeft.clear();
mTopRight.clear();
diff --git a/src/mathicgb/QuadMatrixBuilder.hpp b/src/mathicgb/QuadMatrixBuilder.hpp
index 2d44d9d..4d3c634 100755
--- a/src/mathicgb/QuadMatrixBuilder.hpp
+++ b/src/mathicgb/QuadMatrixBuilder.hpp
@@ -23,6 +23,13 @@ class QuadMatrix;
/// class that allows step-wise construction of a final product.
class QuadMatrixBuilder {
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 SparseMatrix::RowIndex RowIndex;
typedef SparseMatrix::ColIndex ColIndex;
typedef SparseMatrix::Scalar Scalar;
@@ -84,13 +91,13 @@ class QuadMatrixBuilder {
};
typedef MonomialMap<LeftRightColIndex> Map;
- typedef std::vector<monomial> MonomialsType;
+ typedef std::vector<ConstMonoPtr> Monomials;
QuadMatrixBuilder(
const PolyRing& ring,
Map& map,
- MonomialsType& monomialsLeft,
- MonomialsType& monomialsRight,
+ Monomials& monomialsLeft,
+ Monomials& monomialsRight,
size_t memoryQuantum = 0
);
@@ -188,8 +195,8 @@ class QuadMatrixBuilder {
QuadMatrix buildMatrixAndClear();
private:
- MonomialsType& mMonomialsLeft; /// stores one monomial per left column
- MonomialsType& mMonomialsRight; /// stores one monomial per right column
+ Monomials& mMonomialsLeft; /// stores one monomial per left column
+ Monomials& mMonomialsRight; /// stores one monomial per right column
/// Used for fast determination of which column has a given monomial.
Map& mMonomialToCol;
diff --git a/src/mathicgb/SparseMatrix.cpp b/src/mathicgb/SparseMatrix.cpp
index fb7824f..7801475 100755
--- a/src/mathicgb/SparseMatrix.cpp
+++ b/src/mathicgb/SparseMatrix.cpp
@@ -33,7 +33,7 @@ void SparseMatrix::takeRowsFrom(SparseMatrix&& matrix) {
void SparseMatrix::rowToPolynomial(
const RowIndex row,
- const std::vector<monomial>& colMonomials,
+ const std::vector<PolyRing::Monoid::ConstMonoPtr>& colMonomials,
Poly& poly
) {
poly.setToZero();
@@ -42,7 +42,7 @@ void SparseMatrix::rowToPolynomial(
for (auto it = rowBegin(row); it != end; ++it) {
MATHICGB_ASSERT(it.index() < colMonomials.size());
if (it.scalar() != 0)
- poly.appendTerm(it.scalar(), colMonomials[it.index()]);
+ poly.appendTerm(it.scalar(), Monoid::toOld(*colMonomials[it.index()]));
}
MATHICGB_ASSERT(poly.termsAreInDescendingOrder());
}
diff --git a/src/mathicgb/SparseMatrix.hpp b/src/mathicgb/SparseMatrix.hpp
index beeea73..c38d4a3 100755
--- a/src/mathicgb/SparseMatrix.hpp
+++ b/src/mathicgb/SparseMatrix.hpp
@@ -241,7 +241,7 @@ public:
/// Let poly be the dot product of colMonomials and the given row.
void rowToPolynomial(
RowIndex row,
- const std::vector<monomial>& colMonomials,
+ const std::vector<PolyRing::Monoid::ConstMonoPtr>& colMonomials,
Poly& poly);
/// Reorders the rows so that the index of the leading column in
diff --git a/src/test/F4MatrixBuilder.cpp b/src/test/F4MatrixBuilder.cpp
index 5e292af..9b24100 100755
--- a/src/test/F4MatrixBuilder.cpp
+++ b/src/test/F4MatrixBuilder.cpp
@@ -61,7 +61,7 @@ TEST(F4MatrixBuilder, Empty) {
BuilderMaker maker;
F4MatrixBuilder& builder = maker.create();
- QuadMatrix matrix;
+ QuadMatrix matrix(maker.ring());
builder.buildMatrixAndClear(matrix);
ASSERT_EQ(0, matrix.topLeft.rowCount());
ASSERT_EQ(0, matrix.bottomLeft.rowCount());
@@ -82,7 +82,7 @@ TEST(F4MatrixBuilder, SPair) {
const Poly& p3 = maker.addBasisElement("c2d+3");
F4MatrixBuilder& builder = maker.create();
builder.addSPolynomialToMatrix(p1, p2);
- QuadMatrix qm;
+ QuadMatrix qm(builder.ring());
builder.buildMatrixAndClear(qm);
const char* const str1 =
"Left columns: c2d\n"
@@ -109,7 +109,7 @@ TEST(F4MatrixBuilder, OneByOne) {
const Poly& p = maker.addBasisElement("a");
F4MatrixBuilder& builder = maker.create();
builder.addPolynomialToMatrix(p.getLeadMonomial(), p);
- QuadMatrix qm;
+ QuadMatrix qm(builder.ring());
builder.buildMatrixAndClear(qm);
const char* str =
"Left columns: a2\n"
@@ -144,7 +144,7 @@ TEST(F4MatrixBuilder, DirectReducers) {
builder.addPolynomialToMatrix(p2.getLeadMonomial(), p2);
}
- QuadMatrix qm;
+ QuadMatrix qm(builder.ring());
builder.buildMatrixAndClear(qm);
const char* str =
@@ -168,7 +168,7 @@ TEST(F4MatrixBuilder, IteratedReducer) {
const Poly& p2 = maker.addBasisElement("a-1");
F4MatrixBuilder& builder = maker.create();
builder.addPolynomialToMatrix(p1.getLeadMonomial(), p2);
- QuadMatrix qm;
+ QuadMatrix qm(builder.ring());
builder.buildMatrixAndClear(qm);
const char* str =
"Left columns: a5 a4 a3 a2 a\n"
diff --git a/src/test/F4MatrixReducer.cpp b/src/test/F4MatrixReducer.cpp
index a5bd60c..00b18b7 100755
--- a/src/test/F4MatrixReducer.cpp
+++ b/src/test/F4MatrixReducer.cpp
@@ -14,9 +14,7 @@ using namespace mgb;
TEST(F4MatrixReducer, Reduce) {
auto ring = ringFromString("101 6 1\n10 1 1 1 1 1");
-
- QuadMatrix m;
- m.ring = ring.get();
+ QuadMatrix m(*ring);
Poly p(*ring);
std::istringstream in("a4+a3+a2+a1+b5+b4+b3+b2+b1");
diff --git a/src/test/QuadMatrixBuilder.cpp b/src/test/QuadMatrixBuilder.cpp
index fdc2aac..5c7f146 100755
--- a/src/test/QuadMatrixBuilder.cpp
+++ b/src/test/QuadMatrixBuilder.cpp
@@ -8,6 +8,7 @@
#include "mathicgb/io-util.hpp"
#include "mathicgb/Basis.hpp"
#include "mathicgb/QuadMatrix.hpp"
+#include "mathicgb/MathicIO.hpp"
#include <gtest/gtest.h>
using namespace mgb;
@@ -19,6 +20,16 @@ namespace {
return out.str();
}
+ template<class Monoid>
+ std::string monoToStr(
+ const Monoid& monoid,
+ typename Monoid::ConstMonoRef a
+ ) {
+ std::ostringstream out;
+ MathicIO<Monoid>().writeMonomial(monoid, false, a, out);
+ return out.str();
+ }
+
void createColumns(const char* left, const char* right, QuadMatrixBuilder& b)
{
const PolyRing& ring = b.ring();
@@ -61,8 +72,8 @@ TEST(QuadMatrixBuilder, Empty) {
// test a builder with no rows and no columns
PolyRing ring(2, PolyRing::Monoid(0));
QuadMatrixBuilder::Map map(ring);
- QuadMatrixBuilder::MonomialsType monoLeft;
- QuadMatrixBuilder::MonomialsType monoRight;
+ QuadMatrixBuilder::Monomials monoLeft;
+ QuadMatrixBuilder::Monomials monoRight;
QuadMatrixBuilder b(ring, map, monoLeft, monoRight);
const char* matrixStr =
"Left columns:\n"
@@ -79,8 +90,8 @@ TEST(QuadMatrixBuilder, Empty) {
TEST(QuadMatrixBuilder, Construction) {
std::unique_ptr<PolyRing> ring(ringFromString("32003 6 1\n1 1 1 1 1 1"));
QuadMatrixBuilder::Map map(*ring);
- QuadMatrixBuilder::MonomialsType monoLeft;
- QuadMatrixBuilder::MonomialsType monoRight;
+ QuadMatrixBuilder::Monomials monoLeft;
+ QuadMatrixBuilder::Monomials monoRight;
QuadMatrixBuilder b(*ring, map, monoLeft, monoRight);
createColumns("a<1>+<0>", "bc<0>+b<0>+c<0>", b);
@@ -122,8 +133,8 @@ TEST(QuadMatrixBuilder, Construction) {
TEST(QuadMatrixBuilder, ColumnQuery) {
std::unique_ptr<PolyRing> ring(ringFromString("32003 6 1\n1 1 1 1 1 1"));
QuadMatrixBuilder::Map map(*ring);
- QuadMatrixBuilder::MonomialsType monoLeft;
- QuadMatrixBuilder::MonomialsType monoRight;
+ QuadMatrixBuilder::Monomials monoLeft;
+ QuadMatrixBuilder::Monomials monoRight;
QuadMatrixBuilder b(*ring, map, monoLeft, monoRight);
createColumns("a<1>+<0>", "b<0>+c<0>+bc<0>", b);
@@ -156,8 +167,8 @@ TEST(QuadMatrixBuilder, SortColumns) {
// one row top, no rows bottom, no columns
{
QuadMatrixBuilder::Map map(*ring);
- QuadMatrixBuilder::MonomialsType monoLeft;
- QuadMatrixBuilder::MonomialsType monoRight;
+ QuadMatrixBuilder::Monomials monoLeft;
+ QuadMatrixBuilder::Monomials monoRight;
QuadMatrixBuilder b(*ring, map, monoLeft, monoRight);
b.rowDoneTopLeftAndRight();
auto matrix = b.buildMatrixAndClear();
@@ -173,8 +184,8 @@ TEST(QuadMatrixBuilder, SortColumns) {
{
QuadMatrixBuilder::Map map(*ring);
- QuadMatrixBuilder::MonomialsType monoLeft;
- QuadMatrixBuilder::MonomialsType monoRight;
+ QuadMatrixBuilder::Monomials monoLeft;
+ QuadMatrixBuilder::Monomials monoRight;
QuadMatrixBuilder b(*ring, map, monoLeft, monoRight);
createColumns("<0>+a<0>", "b<0>+bcd<0>+bc<0>", b);
b.appendEntryTopLeft(0,1);
@@ -219,8 +230,8 @@ TEST(QuadMatrixBuilder, SortColumns) {
TEST(QuadMatrixBuilder, BuildAndClear) {
std::unique_ptr<PolyRing> ring(ringFromString("32003 6 1\n1 1 1 1 1 1"));
QuadMatrixBuilder::Map map(*ring);
- QuadMatrixBuilder::MonomialsType monoLeft;
- QuadMatrixBuilder::MonomialsType monoRight;
+ QuadMatrixBuilder::Monomials monoLeft;
+ QuadMatrixBuilder::Monomials monoRight;
QuadMatrixBuilder b(*ring, map, monoLeft, monoRight);
createColumns("a<1>+<0>", "b<0>+c<0>+bc<0>", b);
@@ -243,6 +254,8 @@ TEST(QuadMatrixBuilder, BuildAndClear) {
ASSERT_EQ("0: 2#4\n", qm.bottomRight.toString());
ASSERT_EQ(2, qm.leftColumnMonomials.size());
ASSERT_EQ(3, qm.rightColumnMonomials.size());
- ASSERT_EQ("a", monToStr(*ring, qm.leftColumnMonomials[0]));
- ASSERT_EQ("b", monToStr(*ring, qm.rightColumnMonomials[0]));
+
+ const auto& monoid = ring->monoid();
+ ASSERT_EQ("a", monoToStr(monoid, *qm.leftColumnMonomials[0]));
+ ASSERT_EQ("b", monoToStr(monoid, *qm.rightColumnMonomials[0]));
}
diff --git a/src/test/SparseMatrix.cpp b/src/test/SparseMatrix.cpp
index 23e3ef1..b0b961f 100755
--- a/src/test/SparseMatrix.cpp
+++ b/src/test/SparseMatrix.cpp
@@ -66,9 +66,9 @@ TEST(SparseMatrix, Simple) {
TEST(SparseMatrix, toRow) {
auto ring = ringFromString("32003 6 1\n1 1 1 1 1 1");
auto polyForMonomials = parsePoly(*ring, "a5+a4+a3+a2+a1+a0");
- std::vector<monomial> monomials;
+ std::vector<PolyRing::Monoid::ConstMonoPtr> monomials;
for (auto it = polyForMonomials->begin(); it != polyForMonomials->end(); ++it)
- monomials.push_back(it.getMonomial());
+ monomials.push_back(it.mono().ptr());
SparseMatrix mat(5);
mat.clear();
--
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