[mathicgb] 83/393: Parallelized also the construction of the bottom matrix of S-pairs and removed some memory leaks of monomials.
Doug Torrance
dtorrance-guest at moszumanska.debian.org
Fri Apr 3 15:58:36 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 d3ad560b3707c50af2f2569fbbba64b41b078811
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date: Fri Oct 26 16:20:34 2012 +0200
Parallelized also the construction of the bottom matrix of S-pairs and removed some memory leaks of monomials.
---
src/mathicgb/F4MatrixBuilder.cpp | 342 ++++++++++++++++++++++-----------------
src/mathicgb/F4MatrixBuilder.hpp | 51 +++---
src/mathicgb/F4Reducer.cpp | 17 ++
src/mathicgb/PolyRing.hpp | 65 +++++---
src/mathicgb/QuadMatrix.cpp | 66 ++++++--
src/test/F4MatrixBuilder.cpp | 3 +-
6 files changed, 337 insertions(+), 207 deletions(-)
diff --git a/src/mathicgb/F4MatrixBuilder.cpp b/src/mathicgb/F4MatrixBuilder.cpp
index 116a9ee..4ae2668 100755
--- a/src/mathicgb/F4MatrixBuilder.cpp
+++ b/src/mathicgb/F4MatrixBuilder.cpp
@@ -5,6 +5,21 @@
#include <omp.h>
#endif
+MATHIC_INLINE QuadMatrixBuilder::LeftRightColIndex
+ F4MatrixBuilder::findOrCreateColumn
+(
+ const const_monomial monoA,
+ const const_monomial monoB,
+ QuadMatrixBuilder& builder
+) {
+ MATHICGB_ASSERT(!monoA.isNull());
+ MATHICGB_ASSERT(!monoB.isNull());
+ const auto col = builder.findColumnProduct(monoA, monoB);
+ if (col.valid())
+ return col;
+ return createColumn(builder, monoA, monoB);
+}
+
F4MatrixBuilder::F4MatrixBuilder(const PolyBasis& basis, const int threadCount):
mThreadCount(threadCount),
mBasis(basis),
@@ -12,6 +27,13 @@ F4MatrixBuilder::F4MatrixBuilder(const PolyBasis& basis, const int threadCount):
mTmp(basis.ring().allocMonomial())
{
MATHICGB_ASSERT(threadCount >= 1);
+
+ // This assert 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("F4MatrixBuilder: too large characteristic.");
}
void F4MatrixBuilder::addSPolynomialToMatrix
@@ -25,17 +47,18 @@ void F4MatrixBuilder::addSPolynomialToMatrix
ring().monomialLeastCommonMultiple
(polyA.getLeadMonomial(), polyB.getLeadMonomial(), lcm);
- SPairTask task = {};
+ RowTask task;
+ task.addToTop = false;
- task.polyA = &polyA;
- task.multiplyA = ring().allocMonomial();
- ring().monomialDivide(lcm, polyA.getLeadMonomial(), task.multiplyA);
+ task.poly = &polyA;
+ task.multiply = ring().allocMonomial();
+ ring().monomialDivide(lcm, polyA.getLeadMonomial(), task.multiply);
- task.polyB = &polyB;
- task.multiplyB = ring().allocMonomial();
- ring().monomialDivide(lcm, polyB.getLeadMonomial(), task.multiplyB);
+ task.sPairPoly = &polyB;
+ task.sPairMultiply = ring().allocMonomial();
+ ring().monomialDivide(lcm, polyB.getLeadMonomial(), task.sPairMultiply);
- mSPairTodo.push_back(task);
+ mTodo.push_back(task);
ring().freeMonomial(lcm);
}
@@ -43,15 +66,16 @@ void F4MatrixBuilder::addPolynomialToMatrix(const Poly& poly) {
if (poly.isZero())
return;
- SPairTask task = {};
+ RowTask task = {};
+ task.addToTop = false;
- task.polyA = &poly;
- task.multiplyA = ring().allocMonomial();
- ring().monomialSetIdentity(task.multiplyA);
+ task.poly = &poly;
+ task.multiply = ring().allocMonomial();
+ ring().monomialSetIdentity(task.multiply);
- MATHICGB_ASSERT(task.polyB == 0);
- MATHICGB_ASSERT(task.multiplyB.unsafeGetRepresentation() == 0);
- mSPairTodo.push_back(task);
+ MATHICGB_ASSERT(task.sPairPoly == 0);
+ MATHICGB_ASSERT(task.sPairMultiply.isNull());
+ mTodo.push_back(task);
}
void F4MatrixBuilder::addPolynomialToMatrix
@@ -60,106 +84,20 @@ void F4MatrixBuilder::addPolynomialToMatrix
if (poly.isZero())
return;
- SPairTask task = {};
+ RowTask task = {};
+ task.addToTop = false;
- task.polyA = &poly;
- task.multiplyA = ring().allocMonomial();
- ring().monomialCopy(multiple, task.multiplyA);
+ task.poly = &poly;
+ task.multiply = ring().allocMonomial();
+ ring().monomialCopy(multiple, task.multiply);
- MATHICGB_ASSERT(task.polyB == 0);
- MATHICGB_ASSERT(task.multiplyB.unsafeGetRepresentation() == 0);
- mSPairTodo.push_back(task);
+ MATHICGB_ASSERT(task.sPairPoly == 0);
+ MATHICGB_ASSERT(task.sPairMultiply.isNull());
+ mTodo.push_back(task);
}
void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
- // todo: detect and remove duplicate input rows.
// todo: prefer sparse/old reducers among the inputs.
- monomial mono = ring().allocMonomial();
-
- for (auto it = mSPairTodo.begin(); it != mSPairTodo.end(); ++it) {
- Poly::const_iterator itA = it->polyA->begin();
- Poly::const_iterator endA = it->polyA->end();
- MATHICGB_ASSERT(itA != endA);
- MATHICGB_ASSERT(it->multiplyA.unsafeGetRepresentation() != 0);
- MATHICGB_ASSERT(ring().hashValid(it->multiplyA));
-
- // make itB==endB if there is no polyB
- Poly::const_iterator itB;
- Poly::const_iterator endB;
- if (it->polyB != 0) {
- MATHICGB_ASSERT(it->multiplyB.unsafeGetRepresentation() != 0);
- MATHICGB_ASSERT(ring().hashValid(it->multiplyB));
- itB = it->polyB->begin();
- endB = it->polyB->end();
- MATHICGB_ASSERT(itB != endB);
-
- // skip leading terms since they cancel
- ++itA;
- ++itB;
- } else {
- // set polyB to an empty range if null
- itB = endA;
- endB = endA;
- MATHICGB_ASSERT(itB == endB);
- }
-
- monomial monoA = ring().allocMonomial();
- monomial monoB = ring().allocMonomial();
- monomial mono = ring().allocMonomial();
-
- if (itA != endA)
- ring().monomialMult(itA.getMonomial(), it->multiplyA, monoA);
- if (itB != endB)
- ring().monomialMult(itB.getMonomial(), it->multiplyB, monoB);
-
- while (true) {
- bool popA = false;
- bool popB = false;
- if (itA == endA) {
- if (itB == endB)
- break;
- popB = true;
- } else if (itB == endB)
- popA = true;
- else {
- int cmp = ring().monomialCompare(monoA, monoB);
- if (cmp != GT)
- popB = true;
- if (cmp != LT)
- popA = true;
- }
-
- MATHICGB_ASSERT(popA || popB);
- coefficient c = 0;
- if (popA) {
- std::swap(mono, monoA);
- ring().coefficientAddTo(c, itA.getCoefficient());
- ++itA;
- if (itA != endA)
- ring().monomialMult(itA.getMonomial(), it->multiplyA, monoA);
- }
- if (popB) {
- std::swap(mono, monoB);
- coefficient cB = itB.getCoefficient();
- ring().coefficientNegateTo(cB);
- ring().coefficientAddTo(c, cB);
- ++itB;
- if (itB != endB)
- ring().monomialMult(itB.getMonomial(), it->multiplyB, monoB);
- }
- if (c != 0) {
- MATHICGB_ASSERT(c < std::numeric_limits<QuadMatrixBuilder::Scalar>::max());
- auto col = mBuilder.findColumn(mono);
- if (!col.valid())
- col = this->createColumn(mono, mBuilder);
- mBuilder.appendEntryBottom(col,
- static_cast<QuadMatrixBuilder::Scalar>(c));
- }
- }
- ring().freeMonomial(monoA);
- ring().freeMonomial(monoB);
- mBuilder.rowDoneBottomLeftAndRight();
- }
// Process pending rows until we are done. Note that the methods
// we are calling here can add more items to mTodo.
@@ -184,6 +122,13 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
decltype(mTodo) currentTasks;
while (!mTodo.empty()) {
+ for (auto it = currentTasks.begin(); it != currentTasks.end(); ++it) {
+ MATHICGB_ASSERT(!it->multiply.isNull());
+ ring().freeMonomial(it->multiply);
+ MATHICGB_ASSERT(it->sPairMultiply.isNull() == (it->sPairPoly == 0));
+ if (!it->sPairMultiply.isNull())
+ ring().freeMonomial(it->sPairMultiply);
+ }
currentTasks.clear();
mTodo.swap(currentTasks);
const auto taskCountOMP = static_cast<OMPIndex>(currentTasks.size());
@@ -196,15 +141,18 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
MATHICGB_ASSERT(omp_get_thread_num() < threadData.size());
const ThreadData& td = threadData[omp_get_thread_num()];
QuadMatrixBuilder& builder = *td.builder;
- const monomial tmp = td.tmp;
#else
- const monomial tmp = mTmp;
QuadMatrixBuilder& builder = mBuilder;
#endif
+
const RowTask task = currentTasks[taskIndex];
- MATHICGB_ASSERT(ring().hashValid(task.multiple));
+ MATHICGB_ASSERT(ring().hashValid(task.multiply));
- appendRowTop(task.multiple, *task.poly, builder, tmp);
+ if (task.addToTop) {
+ MATHICGB_ASSERT(task.sPairPoly == 0);
+ appendRowTop(task.multiply, *task.poly, builder);
+ } else
+ appendRowBottom(task, builder);
}
}
@@ -217,7 +165,6 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
QuadMatrix qm(it->builder->buildMatrixAndClear());
mBuilder.takeRowsFrom(std::move(qm));
delete it->builder;
- delete[] it->tmp.unsafeGetRepresentation();
}
}
threadData.clear();
@@ -229,53 +176,61 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
F4MatrixBuilder::LeftRightColIndex
F4MatrixBuilder::createColumn(
- const const_monomial mono,
- QuadMatrixBuilder& builder
+ QuadMatrixBuilder& builder,
+ const const_monomial monoA,
+ const const_monomial monoB
) {
- QuadMatrixBuilder::LeftRightColIndex newCol;
+ MATHICGB_ASSERT(!monoA.isNull());
+ MATHICGB_ASSERT(!monoB.isNull());
+ // Must be declared up here since we cannot return out of an omp critical
+ // section.
+ QuadMatrixBuilder::LeftRightColIndex newCol;
#pragma omp critical
{
- MATHICGB_ASSERT(ring().hashValid(mono));
- MATHICGB_ASSERT(!builder.findColumn(mono).valid());
+ ring().monomialMult(monoA, monoB, mTmp);
+
+ MATHICGB_ASSERT(ring().hashValid(mTmp));
+ MATHICGB_ASSERT(!builder.findColumn(mTmp).valid());
#ifndef _OPENMP
MATHICGB_ASSERT(&builder == &mBuilder);
#endif
#ifdef _OPENMP
- // builder does not have a column for mono but mBuilder might
+ // builder does not have a column for mTmp but mBuilder might
if (&mBuilder != &builder)
- newCol = mBuilder.findColumn(mono);
+ newCol = mBuilder.findColumn(mTmp);
if (!newCol.valid())
#endif
{
// we need to add a new column to mBuilder
- // look for a reducer of mono
- size_t reducerIndex = mBasis.divisor(mono);
+ // look for a reducer of mTmp
+ size_t reducerIndex = mBasis.divisor(mTmp);
if (reducerIndex == static_cast<size_t>(-1)) {
- newCol = mBuilder.createColumnRight(mono);
+ newCol = mBuilder.createColumnRight(mTmp);
} else {
- newCol = mBuilder.createColumnLeft(mono);
+ newCol = mBuilder.createColumnLeft(mTmp);
// schedule the reducer to be added as a row
- RowTask task;
+ RowTask task = {};
+ task.addToTop = true;
task.poly = &mBasis.poly(reducerIndex);
- task.multiple = ring().allocMonomial();
+ task.multiply = ring().allocMonomial();
ring().monomialDivideToNegative
- (mono, task.poly->getLeadMonomial(), task.multiple);
- MATHICGB_ASSERT(ring().hashValid(task.multiple));
+ (mTmp, task.poly->getLeadMonomial(), task.multiply);
+ MATHICGB_ASSERT(ring().hashValid(task.multiply));
mTodo.push_back(task);
}
}
MATHICGB_ASSERT(newCol.valid());
- MATHICGB_ASSERT(ring().monomialEQ(mBuilder.monomialOfCol(newCol), mono));
- MATHICGB_ASSERT(newCol == mBuilder.findColumn(mono));
+ MATHICGB_ASSERT(ring().monomialEQ(mBuilder.monomialOfCol(newCol), mTmp));
+ MATHICGB_ASSERT(newCol == mBuilder.findColumn(mTmp));
#ifdef _OPENMP
// If we are running in parallel we need to catch builder up to date with
// with respect to the columns in mBuilder. We have just ensured that
- // mBuilder has a column for mono so now builder will also get that column.
+ // mBuilder has a column for mTmp so now builder will also get that column.
if (&builder != &mBuilder) {
// add missing left columns
@@ -324,31 +279,120 @@ F4MatrixBuilder::createColumn(
}
#endif
- MATHICGB_ASSERT(ring().monomialEQ(mBuilder.monomialOfCol(mBuilder.findColumn(mono)), mono));
- MATHICGB_ASSERT(newCol == mBuilder.findColumn(mono));
+ MATHICGB_ASSERT(ring().monomialEQ(mBuilder.monomialOfCol(mBuilder.findColumn(mTmp)), mTmp));
+ MATHICGB_ASSERT(newCol == mBuilder.findColumn(mTmp));
- MATHICGB_ASSERT(builder.findColumn(mono).valid());
- MATHICGB_ASSERT(ring().monomialEQ(builder.monomialOfCol(builder.findColumn(mono)), mono));
- MATHICGB_ASSERT(newCol == builder.findColumn(mono));
+ MATHICGB_ASSERT(builder.findColumn(mTmp).valid());
+ MATHICGB_ASSERT(ring().monomialEQ(builder.monomialOfCol(builder.findColumn(mTmp)), mTmp));
+ MATHICGB_ASSERT(newCol == builder.findColumn(mTmp));
}
#endif
}
return newCol;
}
-void F4MatrixBuilder::appendRowTop(const const_monomial multiple, const Poly& poly, QuadMatrixBuilder& builder, monomial tmp) {
- Poly::const_iterator end = poly.end();
+void F4MatrixBuilder::appendRowBottom(
+ const_monomial multiple,
+ const bool negate,
+ const Poly::const_iterator begin,
+ const Poly::const_iterator end,
+ QuadMatrixBuilder& builder
+) {
+ // todo: eliminate the code-duplication between here and appendRowTop.
+ MATHICGB_ASSERT(!multiple.isNull());
+ MATHICGB_ASSERT(&builder != 0);
+
+ for (auto it = begin; it != end; ++it) {
+ const auto col = findOrCreateColumn(it.getMonomial(), multiple, builder);
+ const auto origScalar = it.getCoefficient();
+ MATHICGB_ASSERT(origScalar != 0);
+ const auto possiblyNegated =
+ negate ? ring().coefficientNegateNonZero(origScalar) : origScalar;
+ MATHICGB_ASSERT(possiblyNegated < std::numeric_limits<Scalar>::max());
+ builder.appendEntryBottom(col, static_cast<Scalar>(possiblyNegated));
+ }
+ builder.rowDoneBottomLeftAndRight();
+}
+
+void F4MatrixBuilder::appendRowTop(
+ const const_monomial multiple,
+ const Poly& poly,
+ QuadMatrixBuilder& builder
+) {
+ MATHICGB_ASSERT(!multiple.isNull());
+ MATHICGB_ASSERT(&poly != 0);
+ MATHICGB_ASSERT(&builder != 0);
+
+ const Poly::const_iterator end = poly.end();
for (Poly::const_iterator it = poly.begin(); it != end; ++it) {
- MATHICGB_ASSERT(it.getCoefficient() <
- std::numeric_limits<QuadMatrixBuilder::Scalar>::max());
- auto col = builder.findColumnProduct(it.getMonomial(), multiple);
- if (!col.valid()) {
- ring().monomialMult(it.getMonomial(), multiple, tmp);
- col = createColumn(tmp, builder);
- }
- const auto scalar =
- static_cast<QuadMatrixBuilder::Scalar>(it.getCoefficient());
- builder.appendEntryTop(col, scalar);
+ const auto col = findOrCreateColumn(it.getMonomial(), multiple, builder);
+ MATHICGB_ASSERT(it.getCoefficient() < std::numeric_limits<Scalar>::max());
+ MATHICGB_ASSERT(it.getCoefficient());
+ builder.appendEntryTop(col, static_cast<Scalar>(it.getCoefficient()));
}
builder.rowDoneTopLeftAndRight();
}
+
+void F4MatrixBuilder::appendRowBottom(
+ const RowTask& task,
+ QuadMatrixBuilder& builder
+) {
+ MATHICGB_ASSERT(!task.addToTop);
+ MATHICGB_ASSERT(!task.poly->isZero());
+ MATHICGB_ASSERT(!task.multiply.isNull());
+ MATHICGB_ASSERT(ring().hashValid(task.multiply));
+ Poly::const_iterator itA = task.poly->begin();
+ const Poly::const_iterator endA = task.poly->end();
+
+ if (task.sPairPoly == 0) {
+ appendRowBottom(task.multiply, false, itA, endA, builder);
+ return;
+ }
+
+ MATHICGB_ASSERT(!task.sPairPoly->isZero());
+ MATHICGB_ASSERT(!task.sPairMultiply.isNull());
+ MATHICGB_ASSERT(ring().hashValid(task.sPairMultiply));
+ Poly::const_iterator itB = task.sPairPoly->begin();
+ Poly::const_iterator endB = task.sPairPoly->end();
+
+ // skip leading terms since they cancel
+ MATHICGB_ASSERT(itA.getCoefficient() == itB.getCoefficient());
+ ++itA;
+ ++itB;
+
+ const const_monomial mulA = task.multiply;
+ const const_monomial mulB = task.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
+ // (itA == itA && itB == endB) then you need to do that yourself.
+ if (itB == endB) {
+ appendRowBottom(mulA, false, itA, endA, builder);
+ break;
+ }
+ if (itA == endA) {
+ appendRowBottom(mulB, true, itB, endB, builder);
+ break;
+ }
+
+ coefficient coeff = 0;
+ LeftRightColIndex col;
+ const auto colA = findOrCreateColumn(itA.getMonomial(), mulA, builder);
+ const auto colB = findOrCreateColumn(itB.getMonomial(), mulB, builder);
+ const auto cmp = ring().monomialCompare
+ (builder.monomialOfCol(colA), builder.monomialOfCol(colB));
+ if (cmp != LT) {
+ coeff = itA.getCoefficient();
+ col = colA;
+ ++itA;
+ }
+ if (cmp != GT) {
+ coeff = ring().coefficientSubtract(coeff, itB.getCoefficient());
+ col = colB;
+ ++itB;
+ }
+ MATHICGB_ASSERT(coeff < std::numeric_limits<Scalar>::max());
+ if (coeff != 0)
+ builder.appendEntryBottom(col, static_cast<Scalar>(coeff));
+ }
+}
diff --git a/src/mathicgb/F4MatrixBuilder.hpp b/src/mathicgb/F4MatrixBuilder.hpp
index 42f6757..9732092 100755
--- a/src/mathicgb/F4MatrixBuilder.hpp
+++ b/src/mathicgb/F4MatrixBuilder.hpp
@@ -19,6 +19,7 @@ class F4MatrixBuilder {
private:
typedef QuadMatrixBuilder::ColIndex ColIndex;
typedef QuadMatrixBuilder::LeftRightColIndex LeftRightColIndex;
+ typedef QuadMatrixBuilder::Scalar Scalar;
public:
F4MatrixBuilder(const PolyBasis& basis, int threadCount);
@@ -61,31 +62,43 @@ public:
const PolyRing& ring() const {return mBuilder.ring();}
private:
- /** Creates a column with monomial label mono and schedules a new row to
- reduce that column if possible. */
- LeftRightColIndex createColumn(const_monomial mono, QuadMatrixBuilder& builder);
-
- void appendRowTop(const_monomial multiple, const Poly& poly, QuadMatrixBuilder& builder, monomial tmp);
-
- /// Represents an S-pair that was added to the matrix for reduction
- /// or, if polyB is null, a polynomial that was added to the matrix
- /// for reduction.
- struct SPairTask {
- const Poly* polyA;
- monomial multiplyA;
- const Poly* polyB;
- monomial multiplyB;
- };
-
- /// Represents the task of adding a row representing poly*multiple.
+ /** 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. */
+ MATHIC_NO_INLINE LeftRightColIndex createColumn
+ (QuadMatrixBuilder& builder, const_monomial monoA, const_monomial monoB);
+
+ /// Represents the task of adding a row to the matrix. If sPairPoly is null
+ /// then the row to add is multiply * poly. Otherwise, the row to add is
+ /// multiply * poly - sPairMultiply * sPairPoly
+ /// where sPairMultiply makes the lead terms cancel.
struct RowTask {
+ bool addToTop; // add the row to the bottom if false
const Poly* poly;
- monomial multiple;
+ monomial multiply;
+ const Poly* sPairPoly;
+ monomial sPairMultiply;
};
+ void appendRowTop(
+ const_monomial multiple,
+ const Poly& poly,
+ QuadMatrixBuilder& builder);
+ void appendRowBottom(const RowTask& task, QuadMatrixBuilder& builder);
+ void appendRowBottom(
+ const_monomial multiple,
+ bool negate,
+ Poly::const_iterator begin,
+ Poly::const_iterator end,
+ QuadMatrixBuilder& builder
+ );
+
+ MATHIC_INLINE LeftRightColIndex findOrCreateColumn
+ (const_monomial monoA, const_monomial monoB, QuadMatrixBuilder& builder);
+
+
const int mThreadCount;
monomial mTmp;
- std::vector<SPairTask> mSPairTodo;
std::vector<RowTask> mTodo;
const PolyBasis& mBasis;
QuadMatrixBuilder mBuilder;
diff --git a/src/mathicgb/F4Reducer.cpp b/src/mathicgb/F4Reducer.cpp
index 1d29ae9..dcd9cea 100755
--- a/src/mathicgb/F4Reducer.cpp
+++ b/src/mathicgb/F4Reducer.cpp
@@ -67,6 +67,13 @@ std::unique_ptr<Poly> F4Reducer::classicReduceSPoly
MATHICGB_ASSERT(reduced.rowCount() == 1);
reduced.rowToPolynomial(0, qm.rightColumnMonomials, *p);
}
+
+ for (auto it = qm.leftColumnMonomials.begin();
+ it != qm.leftColumnMonomials.end(); ++it)
+ mRing.freeMonomial(*it);
+ for (auto it = qm.rightColumnMonomials.begin();
+ it != qm.rightColumnMonomials.end(); ++it)
+ mRing.freeMonomial(*it);
return p;
}
@@ -107,6 +114,9 @@ void F4Reducer::classicReduceSPolySet
}
F4MatrixReducer(mThreadCount).reduce(basis.ring(), qm, reduced);
monomials = std::move(qm.rightColumnMonomials);
+ for (auto it = qm.leftColumnMonomials.begin();
+ it != qm.leftColumnMonomials.end(); ++it)
+ mRing.freeMonomial(*it);
}
if (tracingLevel >= 2)
@@ -118,6 +128,8 @@ void F4Reducer::classicReduceSPolySet
reduced.rowToPolynomial(row, monomials, *p);
reducedOut.push_back(std::move(p));
}
+ for (auto it = monomials.begin(); it != monomials.end(); ++it)
+ mRing.freeMonomial(*it);
}
void F4Reducer::classicReducePolySet
@@ -155,6 +167,9 @@ void F4Reducer::classicReducePolySet
}
F4MatrixReducer(mThreadCount).reduce(basis.ring(), qm, reduced);
monomials = std::move(qm.rightColumnMonomials);
+ for (auto it = qm.leftColumnMonomials.begin();
+ it != qm.leftColumnMonomials.end(); ++it)
+ mRing.freeMonomial(*it);
}
if (tracingLevel >= 2)
@@ -166,6 +181,8 @@ void F4Reducer::classicReducePolySet
reduced.rowToPolynomial(row, monomials, *p);
reducedOut.push_back(std::move(p));
}
+ for (auto it = monomials.begin(); it != monomials.end(); ++it)
+ mRing.freeMonomial(*it);
}
Poly* F4Reducer::regularReduce
diff --git a/src/mathicgb/PolyRing.hpp b/src/mathicgb/PolyRing.hpp
index b81be5c..9756f39 100755
--- a/src/mathicgb/PolyRing.hpp
+++ b/src/mathicgb/PolyRing.hpp
@@ -291,8 +291,10 @@ public:
- coefficient toCoefficient(int64 a) const;
- coefficient coefficientNegate(coefficient result) const;
+ coefficient toCoefficient(int64 value) const;
+ coefficient coefficientNegate(coefficient coeff) const;
+ coefficient coefficientNegateNonZero(coefficient coeff) const;
+ coefficient coefficientSubtract(coefficient a, coefficient b) const;
void coefficientFromInt(coefficient &result, int a) const;
void coefficientSetOne(coefficient &result) const { result = 1; }
@@ -984,16 +986,6 @@ inline void PolyRing::coefficientDivide(coefficient a, coefficient b, coefficien
MATHICGB_ASSERT(result < mCharac);
}
-inline coefficient PolyRing::toCoefficient(int64 a) const {
- auto modLong = a % mCharac;
- if (modLong < 0)
- modLong += mCharac;
- const coefficient mod = static_cast<coefficient>(modLong);
- MATHICGB_ASSERT(0 <= mod);
- MATHICGB_ASSERT(mod < mCharac);
- return mod;
-}
-
inline void PolyRing::coefficientFromInt(coefficient &result, int a) const
{
result = toCoefficient(a);
@@ -1006,20 +998,47 @@ inline void PolyRing::coefficientAddOneTo(coefficient &result) const
result = 0;
}
-inline void PolyRing::coefficientNegateTo(coefficient &result) const
- // result = -result
-{
+inline void PolyRing::coefficientNegateTo(coefficient& result) const {
+ MATHICGB_ASSERT(result < mCharac);
if (result != 0)
- result = mCharac - result;
+ result = coefficientNegateNonZero(result);
}
-inline coefficient PolyRing::coefficientNegate(coefficient coef) const
- // result = -result
-{
- if (coef == 0)
- return 0;
- else
- return mCharac - coef;
+inline coefficient PolyRing::toCoefficient(const int64 value) const {
+ auto modLong = value % mCharac;
+ if (modLong < 0)
+ modLong += mCharac;
+ MATHICGB_ASSERT(0 <= modLong);
+ MATHICGB_ASSERT(modLong < mCharac);
+ const auto mod = static_cast<coefficient>(modLong);
+ MATHICGB_ASSERT(0 <= mod);
+ MATHICGB_ASSERT(mod < mCharac);
+ return mod;
+}
+
+inline coefficient PolyRing::coefficientNegate(const coefficient coeff) const {
+ MATHICGB_ASSERT(coeff < mCharac);
+ return coeff == 0 ? 0 : coefficientNegateNonZero(coeff);
+}
+
+inline coefficient PolyRing::coefficientNegateNonZero(
+ const coefficient coeff
+) const {
+ MATHICGB_ASSERT(coeff != 0);
+ MATHICGB_ASSERT(coeff < mCharac);
+ return mCharac - coeff;
+}
+
+inline coefficient PolyRing::coefficientSubtract(
+ const coefficient a,
+ const coefficient b
+) const {
+ MATHICGB_ASSERT(a < mCharac);
+ MATHICGB_ASSERT(b < mCharac);
+ const auto diff = a < b ? a + (mCharac - b) : a - b;
+ MATHICGB_ASSERT(diff < mCharac);
+ MATHICGB_ASSERT((diff + b) % mCharac == a);
+ return diff;
}
inline void PolyRing::coefficientAddTo(coefficient &result, coefficient a, coefficient b) const
diff --git a/src/mathicgb/QuadMatrix.cpp b/src/mathicgb/QuadMatrix.cpp
index a820c35..7413c74 100755
--- a/src/mathicgb/QuadMatrix.cpp
+++ b/src/mathicgb/QuadMatrix.cpp
@@ -65,9 +65,6 @@ std::string QuadMatrix::toString() const {
}
QuadMatrix QuadMatrix::toCanonical() const {
- std::vector<size_t> rows;
- for (size_t row = 0; row < topLeft.rowCount(); ++row)
- rows.push_back(row);
class RowComparer {
public:
RowComparer(const SparseMatrix& matrix): mMatrix(matrix) {}
@@ -76,27 +73,66 @@ QuadMatrix QuadMatrix::toCanonical() const {
// then update this code.
MATHICGB_ASSERT(!mMatrix.emptyRow(a));
MATHICGB_ASSERT(!mMatrix.emptyRow(b));
- return mMatrix.leadCol(a) > mMatrix.leadCol(b);
+ auto itA = mMatrix.rowBegin(a);
+ const auto endA = mMatrix.rowEnd(a);
+ auto itB = mMatrix.rowBegin(b);
+ const auto endB = mMatrix.rowEnd(b);
+ for (; itA != endA; ++itA, ++itB) {
+ if (itB == endB)
+ return true;
+
+ if (itA.index() > itB.index())
+ return true;
+ if (itA.index() < itB.index())
+ return false;
+
+ if (itA.scalar() > itB.scalar())
+ return false;
+ if (itA.scalar() < itB.scalar())
+ return true;
+ }
+ return false;
}
private:
const SparseMatrix& mMatrix;
};
- {
- RowComparer comparer(topLeft);
- std::sort(rows.begin(), rows.end(), comparer);
- }
+ // todo: eliminate left/right code duplication here
QuadMatrix matrix;
- matrix.topLeft.clear(topLeft.colCount());
- matrix.topRight.clear(topRight.colCount());
- for (size_t i = 0; i < rows.size(); ++i) {
- matrix.topLeft.appendRow(topLeft, rows[i]);
- matrix.topRight.appendRow(topRight, rows[i]);
+ { // left side
+ std::vector<size_t> rows;
+ for (size_t row = 0; row < topLeft.rowCount(); ++row)
+ rows.push_back(row);
+ {
+ RowComparer comparer(topLeft);
+ std::sort(rows.begin(), rows.end(), comparer);
+ }
+
+ matrix.topLeft.clear(topLeft.colCount());
+ matrix.topRight.clear(topRight.colCount());
+ for (size_t i = 0; i < rows.size(); ++i) {
+ matrix.topLeft.appendRow(topLeft, rows[i]);
+ matrix.topRight.appendRow(topRight, rows[i]);
+ }
+ }
+ { // right side
+ std::vector<size_t> rows;
+ for (size_t row = 0; row < bottomLeft.rowCount(); ++row)
+ rows.push_back(row);
+ {
+ RowComparer comparer(bottomLeft);
+ std::sort(rows.begin(), rows.end(), comparer);
+ }
+
+ matrix.bottomLeft.clear(bottomLeft.colCount());
+ matrix.bottomRight.clear(bottomRight.colCount());
+ for (size_t i = 0; i < rows.size(); ++i) {
+ matrix.bottomLeft.appendRow(bottomLeft, rows[i]);
+ matrix.bottomRight.appendRow(bottomRight, rows[i]);
+ }
}
- matrix.bottomLeft = bottomLeft;
- matrix.bottomRight = bottomRight;
matrix.leftColumnMonomials = leftColumnMonomials;
matrix.rightColumnMonomials = rightColumnMonomials;
matrix.ring = ring;
diff --git a/src/test/F4MatrixBuilder.cpp b/src/test/F4MatrixBuilder.cpp
index 5661a18..c5d8191 100755
--- a/src/test/F4MatrixBuilder.cpp
+++ b/src/test/F4MatrixBuilder.cpp
@@ -142,7 +142,8 @@ TEST(F4MatrixBuilder, DirectReducers) {
" | \n"
"0: 0#1 1#1 2#1 | 0: 0#1\n"
"1: 0#1 1#2 2#3 | 1: 0#4\n";
- ASSERT_EQ(str, qm.toCanonical().toString()) << "** qm:\n" << qm;
+ qm = qm.toCanonical();
+ ASSERT_EQ(str, qm.toString()) << "** qm:\n" << qm;
}
}
--
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