[mathicgb] 62/393: Matrix construction is now a bit faster due to hash table allowing look-up of a product without storing the product anywhere.
Doug Torrance
dtorrance-guest at moszumanska.debian.org
Fri Apr 3 15:58:33 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 2261565ed91772fddc80acd33f52fa15e3bfb39f
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date: Mon Oct 15 21:30:24 2012 +0200
Matrix construction is now a bit faster due to hash table allowing look-up of a product without storing the product anywhere.
---
build/vs12/mathicgb.sln | 3 +++
src/mathicgb/F4MatrixBuilder.cpp | 51 +++++++++++++++-----------------------
src/mathicgb/F4MatrixBuilder.hpp | 9 +++----
src/mathicgb/MonomialMap.hpp | 47 +++++++++++++++++++++++++++++++++--
src/mathicgb/PairTriangle.cpp | 9 ++++---
src/mathicgb/PolyRing.hpp | 32 +++++++++++++++++++++---
src/mathicgb/QuadMatrixBuilder.hpp | 16 +++++++++---
7 files changed, 118 insertions(+), 49 deletions(-)
diff --git a/build/vs12/mathicgb.sln b/build/vs12/mathicgb.sln
index 8a6a242..ad9225a 100755
--- a/build/vs12/mathicgb.sln
+++ b/build/vs12/mathicgb.sln
@@ -139,4 +139,7 @@ Global
GlobalSection(SolutionProperties) = preSolution
HideSolutionNode = FALSE
EndGlobalSection
+ GlobalSection(Performance) = preSolution
+ HasPerformanceSessions = true
+ EndGlobalSection
EndGlobal
diff --git a/src/mathicgb/F4MatrixBuilder.cpp b/src/mathicgb/F4MatrixBuilder.cpp
index 1c4a629..c7041ea 100755
--- a/src/mathicgb/F4MatrixBuilder.cpp
+++ b/src/mathicgb/F4MatrixBuilder.cpp
@@ -2,7 +2,7 @@
#include "F4MatrixBuilder.hpp"
F4MatrixBuilder::F4MatrixBuilder(const PolyBasis& basis):
- mBasis(basis), mBuilder(basis.ring()) {}
+ mBasis(basis), mBuilder(basis.ring()), tmp(basis.ring().allocMonomial()) {}
void F4MatrixBuilder::addSPolynomialToMatrix
(const Poly& polyA, const Poly& polyB) {
@@ -139,7 +139,11 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
}
if (c != 0) {
MATHICGB_ASSERT(c < std::numeric_limits<QuadMatrixBuilder::Scalar>::max());
- mBuilder.appendEntryBottom(createOrFindColumnOf(mono), static_cast<QuadMatrixBuilder::Scalar>(c));
+ auto col = mBuilder.findColumn(mono);
+ if (!col.valid())
+ col = this->createColumn(mono);
+ mBuilder.appendEntryBottom(col,
+ static_cast<QuadMatrixBuilder::Scalar>(c));
}
}
ring().freeMonomial(monoA);
@@ -153,10 +157,7 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
RowTask task = mTodo.back();
MATHICGB_ASSERT(ring().hashValid(task.multiple));
mTodo.pop_back();
- if (task.useAsReducer)
- appendRowTop(task.multiple, *task.poly);
- else
- appendRowBottom(task.multiple, *task.poly);
+ appendRowTop(task.multiple, *task.poly);
}
mBuilder.sortColumnsLeft(mBasis.order());
@@ -165,13 +166,11 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
}
F4MatrixBuilder::LeftRightColIndex
-F4MatrixBuilder::createOrFindColumnOf(const_monomial mono) {
+F4MatrixBuilder::createColumn(const_monomial mono) {
MATHICGB_ASSERT(ring().hashValid(mono));
- LeftRightColIndex colIndex = mBuilder.findColumn(mono);
- if (colIndex.valid())
- return colIndex;
+ MATHICGB_ASSERT(!mBuilder.findColumn(mono).valid());
- // mono did not already have a column so look for a reducer
+ // look for a reducer of mono
size_t reducerIndex = mBasis.divisor(mono);
if (reducerIndex == static_cast<size_t>(-1))
return mBuilder.createColumnRight(mono);
@@ -179,7 +178,6 @@ F4MatrixBuilder::createOrFindColumnOf(const_monomial mono) {
// schedule the reducer to be added as a row
RowTask task;
task.poly = &mBasis.poly(reducerIndex);
- task.useAsReducer = true;
task.multiple = ring().allocMonomial();
ring().monomialDivideToNegative
(mono, task.poly->getLeadMonomial(), task.multiple);
@@ -190,27 +188,18 @@ F4MatrixBuilder::createOrFindColumnOf(const_monomial mono) {
}
void F4MatrixBuilder::appendRowTop(const_monomial multiple, const Poly& poly) {
- monomial mono = ring().allocMonomial();
Poly::const_iterator end = poly.end();
for (Poly::const_iterator it = poly.begin(); it != end; ++it) {
- ring().monomialMult(it.getMonomial(), multiple, mono);
- MATHICGB_ASSERT(it.getCoefficient() < std::numeric_limits<QuadMatrixBuilder::Scalar>::max());
- mBuilder.appendEntryTop(createOrFindColumnOf(mono), static_cast<QuadMatrixBuilder::Scalar>(it.getCoefficient()));
+ MATHICGB_ASSERT(it.getCoefficient() <
+ std::numeric_limits<QuadMatrixBuilder::Scalar>::max());
+ auto col = mBuilder.findColumnProduct(it.getMonomial(), multiple);
+ if (!col.valid()) {
+ ring().monomialMult(it.getMonomial(), multiple, tmp);
+ col = createColumn(tmp);
+ }
+ const auto scalar =
+ static_cast<QuadMatrixBuilder::Scalar>(it.getCoefficient());
+ mBuilder.appendEntryTop(col, scalar);
}
mBuilder.rowDoneTopLeftAndRight();
- ring().freeMonomial(mono);
-}
-
-void F4MatrixBuilder::appendRowBottom
-(const_monomial multiple, const Poly& poly) {
- monomial mono = ring().allocMonomial();
- Poly::const_iterator end = poly.end();
- for (Poly::const_iterator it = poly.begin(); it != end; ++it) {
- ring().monomialMult(it.getMonomial(), multiple, mono);
- MATHICGB_ASSERT(it.getCoefficient() < std::numeric_limits<QuadMatrixBuilder::Scalar>::max());
- mBuilder.appendEntryBottom
- (createOrFindColumnOf(mono), static_cast<QuadMatrixBuilder::Scalar>(it.getCoefficient()));
- }
- mBuilder.rowDoneBottomLeftAndRight();
- ring().freeMonomial(mono);
}
diff --git a/src/mathicgb/F4MatrixBuilder.hpp b/src/mathicgb/F4MatrixBuilder.hpp
index 8855b63..ca60dae 100755
--- a/src/mathicgb/F4MatrixBuilder.hpp
+++ b/src/mathicgb/F4MatrixBuilder.hpp
@@ -61,12 +61,11 @@ public:
const PolyRing& ring() const {return mBuilder.ring();}
private:
- /** Returns a left or right column that represents mono. Creates a
- new column and schedules a new row to reduce that column if necessary. */
- LeftRightColIndex createOrFindColumnOf(const_monomial mono);
+ /** Creates a column with monomial label mono and schedules a new row to
+ reduce that column if possible. */
+ LeftRightColIndex createColumn(const_monomial mono);
void appendRowTop(const_monomial multiple, const Poly& poly);
- void appendRowBottom(const_monomial multiple, const Poly& poly);
/// 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
@@ -80,11 +79,11 @@ private:
/// Represents the task of adding a row representing poly*multiple.
struct RowTask {
- bool useAsReducer; // if true: put in top part of matrix
const Poly* poly;
monomial multiple;
};
+ monomial tmp;
std::vector<SPairTask> mSPairTodo;
std::vector<RowTask> mTodo;
const PolyBasis& mBasis;
diff --git a/src/mathicgb/MonomialMap.hpp b/src/mathicgb/MonomialMap.hpp
index 22f9678..36b9922 100755
--- a/src/mathicgb/MonomialMap.hpp
+++ b/src/mathicgb/MonomialMap.hpp
@@ -57,7 +57,6 @@ namespace MonomialMapInternal {
public:
Hash(const PolyRing& ring): mRing(ring) {}
size_t operator()(const_monomial m) const {
- MATHICGB_ASSERT(mRing.hashValid(m));
return mRing.monomialHashValue(m);
}
const PolyRing& ring() const {return mRing;}
@@ -124,21 +123,59 @@ namespace MonomialMapInternal {
HashAllocator;
typedef std::unordered_map<const_monomial, MTT, Hash, Equal, HashAllocator>
Map;
+ typedef typename Map::iterator iterator;
+ typedef typename Map::const_iterator const_iterator;
MapClass(const PolyRing& ring):
mArena(),
- mMap(100, Hash(ring), Equal(ring), HashAllocator(mArena))
+ mMap(100, Hash(ring), Equal(ring), HashAllocator(mArena)),
+ mTmp(ring.allocMonomial())
{
mMap.max_load_factor(0.3f);
}
+ ~MapClass() {
+ ring().freeMonomial(mTmp);
+ }
+
Map& map() {return mMap;}
const Map& map() const {return mMap;}
const PolyRing& ring() const {return mMap.key_eq().ring();}
+ iterator findProduct(const_monomial a, const_monomial b) {
+ ring().setGivenHash(mTmp, ring().monomialHashOfProduct(a, b));
+ size_t bucket = mMap.bucket(mTmp);
+ auto stop = mMap.end(bucket);
+ for (auto it = mMap.begin(bucket); it != stop; ++it)
+ if (ring().monomialIsProductOf(a, b, it->first))
+ return it;
+ return mMap.end();
+ }
+
+ const_iterator findProduct(const_monomial a, const_monomial b) const {
+ iterator it = const_cast<MapClass<MTT>*>(this)->findProduct(a, b);
+ return const_iterator(it);
+ }
+
+/*
+ size_t bucket = mMonomialToCol.
+ iterator lower_bound(const key_type& _Keyval)
+ { // find leftmost not less than _Keyval in mutable hash table
+ size_type _Bucket = _Hashval(_Keyval);
+ for (_Unchecked_iterator _Where = _Begin(_Bucket);
+ _Where != _End(_Bucket); ++_Where)
+ if (!((_Traits&)*this)(this->_Kfn(*_Where), _Keyval))
+ return (((_Traits&)*this)(_Keyval,
+ this->_Kfn(*_Where)) ? end() : _Make_iter(_Where));
+ return (end());
+ }
+
+ */
+
private:
memt::Arena mArena;
Map mMap;
+ monomial mTmp;
};
#endif
}
@@ -161,6 +198,12 @@ public:
const_iterator end() const {return map().end();}
iterator find(const_monomial m) {return map().find(m);}
const_iterator find(const_monomial m) const {return map().find(m);}
+ iterator findProduct(const_monomial a, const_monomial b) {
+ return mMap.findProduct(a, b);
+ }
+ const_iterator findProduct(const_monomial a, const_monomial b) const {
+ return mMap.findProduct(a, b);
+ }
size_t size() const {return map().size();}
void clear() {map().clear();}
diff --git a/src/mathicgb/PairTriangle.cpp b/src/mathicgb/PairTriangle.cpp
index eb2f17a..39326d5 100755
--- a/src/mathicgb/PairTriangle.cpp
+++ b/src/mathicgb/PairTriangle.cpp
@@ -36,18 +36,19 @@ void PairTriangle::addPair(size_t index, monomial orderBy) {
}
namespace {
+ // Iterator that accesses the field i based on a passed-in iterator.
template<class PairIterator>
class IndexIterator {
public:
typedef typename PairIterator::iterator_category iterator_category;
- typedef typename PairIterator::value_type value_type;
+ typedef decltype(reinterpret_cast<typename PairIterator::value_type*>(0)->i) value_type;
typedef typename PairIterator::difference_type difference_type;
- typedef typename PairIterator::pointer pointer;
- typedef typename PairIterator::reference reference;
+ typedef value_type* pointer;
+ typedef value_type& reference;
IndexIterator(PairIterator pairIterator): mIterator(pairIterator) {}
IndexIterator& operator++() {++mIterator; return *this;}
- size_t operator*() const {return mIterator->i;}
+ const value_type operator*() const {return mIterator->i;}
difference_type operator-(const IndexIterator<PairIterator>& it) const {
return mIterator - it.mIterator;
}
diff --git a/src/mathicgb/PolyRing.hpp b/src/mathicgb/PolyRing.hpp
index 023a429..3eb6641 100755
--- a/src/mathicgb/PolyRing.hpp
+++ b/src/mathicgb/PolyRing.hpp
@@ -321,7 +321,7 @@ public:
// Monomial Routines //////////////////////
///////////////////////////////////////////
- size_t monomialHashValue(ConstMonomial m) const { return static_cast<size_t>(m[mHashIndex]); }
+ exponent monomialHashValue(ConstMonomial m) const {return m[mHashIndex];}
exponent monomialExponent(ConstMonomial m, size_t var) const {
return m[var+1];
@@ -345,6 +345,8 @@ public:
inline void setHashOnly(Monomial& a) const;
+ void setGivenHash(Monomial& a, exponent hash) const {a[mHashIndex] = hash;}
+
bool hashValid(const_monomial m) const;
bool weightsCorrect(ConstMonomial a) const;
@@ -389,15 +391,23 @@ public:
void monomialMultTo(Monomial &a, ConstMonomial b) const; // a *= b
+ /// returns true if b divides a, in this case, result is set to b//a.
bool monomialDivide(ConstMonomial a, ConstMonomial b, Monomial &result) const;
- // returns true if b divides a, in this case, result is set to b//a.
+ /// sets result to a/b, even if that produces negative exponents.
void monomialDivideToNegative(ConstMonomial a, ConstMonomial b, Monomial &result) const;
- // sets result to a/b, even if that produces negative exponents.
+ /// returns true if b divides a. Components are ignored.
bool monomialIsDivisibleBy(ConstMonomial a, ConstMonomial b) const;
- // returns true if b divides a. Components are ignored.
+ /// Returns true if ab is the product of a and b.
+ bool monomialIsProductOf
+ (ConstMonomial a, ConstMonomial b, ConstMonomial ab) const;
+
+ /// Returns the hash of the product of a and b.
+ exponent monomialHashOfProduct(ConstMonomial a, ConstMonomial b) const {
+ return a[mHashIndex] + b[mHashIndex];
+ }
void monomialCopy(ConstMonomial a, Monomial &result) const;
@@ -587,6 +597,20 @@ inline bool PolyRing::monomialEQ(ConstMonomial a, ConstMonomial b) const
return true;
}
+inline bool PolyRing::monomialIsProductOf(
+ ConstMonomial a,
+ ConstMonomial b,
+ ConstMonomial ab
+) const {
+ MATHICGB_ASSERT(hashValid(a));
+ MATHICGB_ASSERT(hashValid(b));
+ MATHICGB_ASSERT(hashValid(ab));
+ for (size_t i = mHashIndex; i != static_cast<size_t>(-1); --i)
+ if (ab[i] != a[i] + b[i])
+ return false;
+ return true;
+}
+
inline void PolyRing::monomialMult(ConstMonomial a,
ConstMonomial b,
Monomial &result) const
diff --git a/src/mathicgb/QuadMatrixBuilder.hpp b/src/mathicgb/QuadMatrixBuilder.hpp
index b36aa8e..7a425c9 100755
--- a/src/mathicgb/QuadMatrixBuilder.hpp
+++ b/src/mathicgb/QuadMatrixBuilder.hpp
@@ -164,9 +164,9 @@ class QuadMatrixBuilder {
// *** Querying columns
- // Returns a column for the findThis monomial. Searches on both the
- // left and right side. Returns an invalid index if no such column
- // exists.
+ /** Returns a column for the findThis monomial. Searches on both the
+ left and right side. Returns an invalid index if no such column
+ exists. */
LeftRightColIndex findColumn(const_monomial findThis) const {
auto it = mMonomialToCol.find(findThis);
if (it != mMonomialToCol.end())
@@ -175,6 +175,16 @@ class QuadMatrixBuilder {
return LeftRightColIndex();
}
+ /// As findColumn, but looks for the monomial that is the product of a and b.
+ LeftRightColIndex findColumnProduct(const_monomial a, const_monomial b) const
+ {
+ auto it = mMonomialToCol.findProduct(a, b);
+ if (it != mMonomialToCol.end())
+ return it->second;
+ else
+ return LeftRightColIndex();
+ }
+
const_monomial monomialOfLeftCol(ColIndex col) const {
MATHICGB_ASSERT(col < mMonomialsLeft.size());
return mMonomialsLeft[col];
--
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