[mathicgb] 49/393: F4Reducer now uses F4 reduction to reduce sets of polynomials (as opposed to just sets of S-polynomials). Also did a few microreductions that increased speed on hcyclic8 by about 10%.
Doug Torrance
dtorrance-guest at moszumanska.debian.org
Fri Apr 3 15:58:31 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 b6a273c51964989c17fa210ea6003517b9d08861
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date: Fri Oct 5 18:47:53 2012 +0200
F4Reducer now uses F4 reduction to reduce sets of polynomials (as opposed to just sets of S-polynomials). Also did a few microreductions that increased speed on hcyclic8 by about 10%.
---
src/mathicgb/BuchbergerAlg.cpp | 4 +-
src/mathicgb/F4MatrixBuilder.cpp | 25 +++++++++++--
src/mathicgb/F4MatrixBuilder.hpp | 4 ++
src/mathicgb/F4Reducer.cpp | 55 +++++++++++++++++++++++++---
src/mathicgb/Poly.cpp | 16 +++-----
src/mathicgb/Poly.hpp | 10 +++++
src/mathicgb/PolyRing.cpp | 72 +-----------------------------------
src/mathicgb/PolyRing.hpp | 75 ++++++++++++++++++++++++++++++++++++++
src/mathicgb/QuadMatrixBuilder.cpp | 5 ++-
9 files changed, 172 insertions(+), 94 deletions(-)
diff --git a/src/mathicgb/BuchbergerAlg.cpp b/src/mathicgb/BuchbergerAlg.cpp
index 98a7329..d51d429 100755
--- a/src/mathicgb/BuchbergerAlg.cpp
+++ b/src/mathicgb/BuchbergerAlg.cpp
@@ -32,8 +32,10 @@ BuchbergerAlg::BuchbergerAlg(
{
// Reduce and insert the generators of the ideal into the starting basis
size_t const idealSize = ideal.size();
+ std::vector<std::unique_ptr<Poly> > polys;
for (size_t gen = 0; gen != idealSize; ++gen)
- insertReducedPoly(mReducer->classicReduce(*ideal.getPoly(gen), mBasis));
+ polys.push_back(make_unique<Poly>(*ideal.getPoly(gen)));
+ insertPolys(polys);
}
void BuchbergerAlg::insertPolys
diff --git a/src/mathicgb/F4MatrixBuilder.cpp b/src/mathicgb/F4MatrixBuilder.cpp
index a680c11..98fda78 100755
--- a/src/mathicgb/F4MatrixBuilder.cpp
+++ b/src/mathicgb/F4MatrixBuilder.cpp
@@ -7,7 +7,9 @@ F4MatrixBuilder::F4MatrixBuilder(const PolyBasis& basis):
void F4MatrixBuilder::addSPolynomialToMatrix
(const Poly& polyA, const Poly& polyB) {
MATHICGB_ASSERT(!polyA.isZero());
+ MATHICGB_ASSERT(polyA.isMonic());
MATHICGB_ASSERT(!polyB.isZero());
+ MATHICGB_ASSERT(polyB.isMonic());
monomial lcm = ring().allocMonomial();
ring().monomialLeastCommonMultiple
@@ -27,6 +29,21 @@ void F4MatrixBuilder::addSPolynomialToMatrix
ring().freeMonomial(lcm);
}
+void F4MatrixBuilder::addPolynomialToMatrix(const Poly& poly) {
+ if (poly.isZero())
+ return;
+
+ SPairTask task = {};
+
+ task.polyA = &poly;
+ task.multiplyA = ring().allocMonomial();
+ ring().monomialSetIdentity(task.multiplyA);
+
+ MATHICGB_ASSERT(task.polyB == 0);
+ MATHICGB_ASSERT(task.multiplyB.unsafeGetRepresentation() == 0);
+ mSPairTodo.push_back(task);
+}
+
void F4MatrixBuilder::addPolynomialToMatrix
(const_monomial multiple, const Poly& poly) {
MATHICGB_ASSERT(ring().hashValid(multiple));
@@ -65,6 +82,10 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
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;
@@ -72,10 +93,6 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
MATHICGB_ASSERT(itB == endB);
}
- // skip leading terms since they cancel
- //++itA;
- //++itB;
-
monomial monoA = ring().allocMonomial();
monomial monoB = ring().allocMonomial();
monomial mono = ring().allocMonomial();
diff --git a/src/mathicgb/F4MatrixBuilder.hpp b/src/mathicgb/F4MatrixBuilder.hpp
index c2dc425..8855b63 100755
--- a/src/mathicgb/F4MatrixBuilder.hpp
+++ b/src/mathicgb/F4MatrixBuilder.hpp
@@ -38,6 +38,10 @@ public:
remain valid. */
void addPolynomialToMatrix(const_monomial multiple, const Poly& poly);
+ /// as the overload with a multiple, just letting multiple be the
+ /// identity.
+ void addPolynomialToMatrix(const Poly& poly);
+
/** Builds an F4 matrix to the specifications given. Also clears the
information in this object.
diff --git a/src/mathicgb/F4Reducer.cpp b/src/mathicgb/F4Reducer.cpp
index eca815a..02b8032 100755
--- a/src/mathicgb/F4Reducer.cpp
+++ b/src/mathicgb/F4Reducer.cpp
@@ -75,12 +75,20 @@ void F4Reducer::classicReduceSPolySet
const PolyBasis& basis,
std::vector<std::unique_ptr<Poly> >& reducedOut)
{
+ if (spairs.size() <= 1) {
+ if (tracingLevel >= 2)
+ std::cerr << "F4Reducer: Using fall-back reducer for "
+ << spairs.size() << " S-pairs.\n";
+ mFallback->classicReduceSPolySet(spairs, basis, reducedOut);
+ return;
+ }
+
reducedOut.clear();
if (spairs.empty())
return;
if (tracingLevel >= 2)
- std::cerr << "F4Reducer: Reducing " << spairs.size() << " S-pair.\n";
+ std::cerr << "F4Reducer: Reducing " << spairs.size() << " S-polynomials.\n";
SparseMatrix reduced;
std::vector<monomial> monomials;
@@ -108,7 +116,6 @@ void F4Reducer::classicReduceSPolySet
for (SparseMatrix::RowIndex row = 0; row < reduced.rowCount(); ++row) {
auto p = make_unique<Poly>(&basis.ring());
reduced.rowToPolynomial(row, monomials, *p);
-
reducedOut.push_back(std::move(p));
}
}
@@ -118,10 +125,46 @@ void F4Reducer::classicReducePolySet
const PolyBasis& basis,
std::vector<std::unique_ptr<Poly> >& reducedOut)
{
- for (auto it = polys.begin(); it != polys.end(); ++it) {
- auto reducedPoly = classicReduce(**it, basis);
- if (!reducedPoly->isZero())
- reducedOut.push_back(std::move(reducedPoly));
+ if (polys.size() <= 1) {
+ if (tracingLevel >= 2)
+ std::cerr << "F4Reducer: Using fall-back reducer for "
+ << polys.size() << " polynomials.\n";
+ mFallback->classicReducePolySet(polys, basis, reducedOut);
+ return;
+ }
+
+ reducedOut.clear();
+ if (polys.empty())
+ return;
+
+ if (tracingLevel >= 2)
+ std::cerr << "F4Reducer: Reducing " << polys.size() << " polynomials.\n";
+
+ SparseMatrix reduced;
+ std::vector<monomial> monomials;
+ {
+ QuadMatrix qm;
+ {
+ F4MatrixBuilder builder(basis);
+ for (auto it = polys.begin(); it != polys.end(); ++it)
+ builder.addPolynomialToMatrix(**it);
+ builder.buildMatrixAndClear(qm);
+
+ // there has to be something to reduce
+ MATHICGB_ASSERT(qm.bottomLeft.rowCount() > 0);
+ }
+ F4MatrixReducer(mThreadCount).reduce(basis.ring(), qm, reduced);
+ monomials = std::move(qm.rightColumnMonomials);
+ }
+
+ if (tracingLevel >= 2)
+ std::cerr << "F4Reducer: Extracted " << reduced.rowCount()
+ << " non-zero rows\n";
+
+ for (SparseMatrix::RowIndex row = 0; row < reduced.rowCount(); ++row) {
+ auto p = make_unique<Poly>(&basis.ring());
+ reduced.rowToPolynomial(row, monomials, *p);
+ reducedOut.push_back(std::move(p));
}
}
diff --git a/src/mathicgb/Poly.cpp b/src/mathicgb/Poly.cpp
index c187cd6..708cec6 100755
--- a/src/mathicgb/Poly.cpp
+++ b/src/mathicgb/Poly.cpp
@@ -73,15 +73,6 @@ const coefficient Poly::coefficientAt(size_t index) const {
return coeffs[index];
}
-void Poly::appendTerm(coefficient a, const_monomial m)
-{
- // the monomial will be copied on.
- coeffs.push_back(a);
- size_t len = R->monomialSize(m);
- exponent const * e = m.unsafeGetRepresentation();
- monoms.insert(monoms.end(), e, e + len);
-}
-
void Poly::append(iterator &first, iterator &last)
{
for ( ; first != last; ++first)
@@ -104,8 +95,11 @@ void Poly::multByCoefficient(coefficient c)
R->coefficientMultTo(*i, c);
}
-void Poly::makeMonic()
-{
+bool Poly::isMonic() const {
+ return !isZero() && R->coefficientIsOne(getLeadCoefficient());
+}
+
+void Poly::makeMonic() {
if (isZero())
return;
coefficient c = getLeadCoefficient();
diff --git a/src/mathicgb/Poly.hpp b/src/mathicgb/Poly.hpp
index b64cb9e..3fbaa77 100755
--- a/src/mathicgb/Poly.hpp
+++ b/src/mathicgb/Poly.hpp
@@ -99,6 +99,7 @@ public:
void multByCoefficient(coefficient a);
void makeMonic();
+ bool isMonic() const;
const_monomial getLeadMonomial() const { return &(monoms[0]); }
const_coefficient getLeadCoefficient() const { return coeffs[0]; }
@@ -149,6 +150,15 @@ inline bool operator!=(const Poly::const_iterator &a, const Poly::const_iterator
return a.ic != b.ic;
}
+inline void Poly::appendTerm(coefficient a, const_monomial m)
+{
+ // the monomial will be copied on.
+ coeffs.push_back(a);
+ size_t len = R->monomialSize(m);
+ exponent const * e = m.unsafeGetRepresentation();
+ monoms.insert(monoms.end(), e, e + len);
+}
+
#endif
// Local Variables:
diff --git a/src/mathicgb/PolyRing.cpp b/src/mathicgb/PolyRing.cpp
index 594399c..d13bd02 100755
--- a/src/mathicgb/PolyRing.cpp
+++ b/src/mathicgb/PolyRing.cpp
@@ -60,7 +60,7 @@ void PolyRing::setWeightsAndHash(Monomial& a1) const
setHashOnly(a1);
}
-inline bool PolyRing::weightsCorrect(ConstMonomial a1) const
+bool PolyRing::weightsCorrect(ConstMonomial a1) const
{
exponent const *a = a1.mValue;
++a;
@@ -184,14 +184,6 @@ void PolyRing::monomialGreatestCommonDivisor(ConstMonomial a,
setWeightsOnly(g);
}
-bool PolyRing::monomialIsLeastCommonMultiple(
- ConstMonomial a,
- ConstMonomial b,
- ConstMonomial l) const
-{
- return monomialIsLeastCommonMultipleNoWeights(a, b, l) && weightsCorrect(l);
-}
-
bool PolyRing::monomialIsLeastCommonMultipleNoWeights(
ConstMonomial a,
ConstMonomial b,
@@ -336,7 +328,7 @@ bool PolyRing::monomialIsLeastCommonMultiple(
return monomialIsLeastCommonMultipleNoWeights(a, b, l) && weightsCorrect(l);
}
-inline bool PolyRing::weightsCorrect(const_monomial a) const
+bool PolyRing::weightsCorrect(const_monomial a) const
{
++a;
const int *wts = &mWeights[0];
@@ -665,49 +657,6 @@ void PolyRing::write(std::ostream &o) const
}
}
-void PolyRing::coefficientFromInt(coefficient &result, int a) const
-{
- result = a % mCharac;
- if (result < 0) result += mCharac;
-}
-
-void PolyRing::coefficientAddOneTo(coefficient &result) const
-{
- result++;
- if (result == mCharac) result = 0;
-}
-void PolyRing::coefficientNegateTo(coefficient &result) const
- // result = -result
-{
- result = mCharac - result;
-}
-void PolyRing::coefficientAddTo(coefficient &result, coefficient a, coefficient b) const
-// result += a*b
-{
- mStats.n_addmult++;
- long c = a * b + result;
- result = c % mCharac;
-}
-void PolyRing::coefficientAddTo(coefficient &result, coefficient a) const
- // result += a
-{
- mStats.n_add++;
- result += a;
- if (result >= mCharac) result -= mCharac;
-}
-void PolyRing::coefficientMultTo(coefficient &result, coefficient a) const
- // result *= a
-{
- mStats.n_mult++;
- long b = result * a;
- result = b % mCharac;
-}
-void PolyRing::coefficientMult(coefficient a, coefficient b, coefficient &result) const
-{
- mStats.n_mult++;
- long c = b * a;
- result = c % mCharac;
-}
void gcd_extended(long a, long b, long &u, long &v, long &g)
{
@@ -728,23 +677,6 @@ void gcd_extended(long a, long b, long &u, long &v, long &g)
}
}
-void PolyRing::coefficientReciprocalTo(coefficient& result) const
-{
- MATHICGB_ASSERT(result != 0);
- mStats.n_recip++;
- result = modularInverse(result, mCharac);
-}
-
-void PolyRing::coefficientDivide(coefficient a, coefficient b, coefficient &result) const
- // result = a/b
-{
- mStats.n_divide++;
- result = (a * modularInverse(b, mCharac)) % mCharac;
- MATHICGB_ASSERT((result * b) % mCharac == a);
- MATHICGB_ASSERT(result >= 0);
- MATHICGB_ASSERT(result < mCharac);
-}
-
// Local Variables:
// compile-command: "make -C .. "
// indent-tabs-mode: nil
diff --git a/src/mathicgb/PolyRing.hpp b/src/mathicgb/PolyRing.hpp
index 89dfbc9..a1bb68b 100755
--- a/src/mathicgb/PolyRing.hpp
+++ b/src/mathicgb/PolyRing.hpp
@@ -853,8 +853,83 @@ inline bool PolyRing::monomialRelativelyPrime(const_monomial a, const_monomial b
return false;
return true;
}
+
#endif
+inline bool PolyRing::monomialIsLeastCommonMultiple(
+ ConstMonomial a,
+ ConstMonomial b,
+ ConstMonomial l) const
+{
+ return monomialIsLeastCommonMultipleNoWeights(a, b, l) && weightsCorrect(l);
+}
+
+inline void PolyRing::coefficientReciprocalTo(coefficient& result) const
+{
+ MATHICGB_ASSERT(result != 0);
+ mStats.n_recip++;
+ result = modularInverse(result, mCharac);
+}
+
+inline void PolyRing::coefficientDivide(coefficient a, coefficient b, coefficient &result) const
+ // result = a/b
+{
+ mStats.n_divide++;
+ result = (a * modularInverse(b, mCharac)) % mCharac;
+ MATHICGB_ASSERT((result * b) % mCharac == a);
+ MATHICGB_ASSERT(result >= 0);
+ MATHICGB_ASSERT(result < mCharac);
+}
+
+inline void PolyRing::coefficientFromInt(coefficient &result, int a) const
+{
+ result = a % mCharac;
+ if (result < 0) result += mCharac;
+}
+
+inline void PolyRing::coefficientAddOneTo(coefficient &result) const
+{
+ result++;
+ if (result == mCharac) result = 0;
+}
+
+inline void PolyRing::coefficientNegateTo(coefficient &result) const
+ // result = -result
+{
+ result = mCharac - result;
+}
+
+inline void PolyRing::coefficientAddTo(coefficient &result, coefficient a, coefficient b) const
+// result += a*b
+{
+ mStats.n_addmult++;
+ long c = a * b + result;
+ result = c % mCharac;
+}
+
+inline void PolyRing::coefficientAddTo(coefficient &result, coefficient a) const
+ // result += a
+{
+ mStats.n_add++;
+ result += a;
+ if (result >= mCharac) result -= mCharac;
+}
+
+inline void PolyRing::coefficientMultTo(coefficient &result, coefficient a) const
+ // result *= a
+{
+ mStats.n_mult++;
+ long b = result * a;
+ result = b % mCharac;
+}
+
+inline void PolyRing::coefficientMult(coefficient a, coefficient b, coefficient &result) const
+{
+ mStats.n_mult++;
+ long c = b * a;
+ result = c % mCharac;
+}
+
#endif
// Local Variables:
diff --git a/src/mathicgb/QuadMatrixBuilder.cpp b/src/mathicgb/QuadMatrixBuilder.cpp
index 24974e2..209cea7 100755
--- a/src/mathicgb/QuadMatrixBuilder.cpp
+++ b/src/mathicgb/QuadMatrixBuilder.cpp
@@ -30,7 +30,7 @@ namespace {
if (colCount == std::numeric_limits<QuadMatrixBuilder::ColIndex>::max())
throw std::overflow_error("Too many columns in QuadMatrixBuilder");
- toMonomial.reserve(toMonomial.size() + 1);
+ toMonomial.push_back(0); // allocate memory ahead of time to avoid bad_alloc
monomial copied = ring.allocMonomial();
ring.monomialCopy(mono, copied);
try {
@@ -38,10 +38,11 @@ namespace {
(copied,
QuadMatrixBuilder::LeftRightColIndex(colCount, left)));
} catch (...) {
+ toMonomial.pop_back();
ring.freeMonomial(copied);
throw;
}
- toMonomial.push_back(copied); // no throw due to reserve
+ toMonomial.back() = copied;
top.ensureAtLeastThisManyColumns(colCount + 1);
bottom.ensureAtLeastThisManyColumns(colCount + 1);
--
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