[mathicgb] 45/393: Fixed a bug: we were inappropriately assuming that all input polynomials have their terms already sorted. Also added asserts in places checking that a polynomial is sorted.
Doug Torrance
dtorrance-guest at moszumanska.debian.org
Fri Apr 3 15:58:30 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 bfc0b259da0b245b992e2aa89e38f26d0f9c6354
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date: Thu Oct 4 20:57:14 2012 +0200
Fixed a bug: we were inappropriately assuming that all input polynomials have their terms already sorted. Also added asserts in places checking that a polynomial is sorted.
---
src/mathicgb/BjarkeGeobucket2.cpp | 9 +--
src/mathicgb/F4MatrixReducer.cpp | 15 ++++
src/mathicgb/F4Reducer.cpp | 21 +++++
src/mathicgb/Ideal.cpp | 23 +++---
src/mathicgb/Poly.cpp | 157 +++++++++++++++++++++++++++-----------
src/mathicgb/Poly.hpp | 13 +++-
src/mathicgb/PolyRing.cpp | 3 +-
src/test/F4MatrixReducer.cpp | 8 +-
src/test/QuadMatrixBuilder.cpp | 10 +--
src/test/poly-test.cpp | 6 +-
10 files changed, 192 insertions(+), 73 deletions(-)
diff --git a/src/mathicgb/BjarkeGeobucket2.cpp b/src/mathicgb/BjarkeGeobucket2.cpp
index ef384a7..36ecf3a 100755
--- a/src/mathicgb/BjarkeGeobucket2.cpp
+++ b/src/mathicgb/BjarkeGeobucket2.cpp
@@ -59,14 +59,13 @@ void BjarkeGeobucket2::insertTail(const_term multiplier, const Poly *g1)
void BjarkeGeobucket2::insert(monomial multiplier, const Poly *g1)
{
- HashPoly M;
+ MATHICGB_ASSERT(g1 != 0);
+ MATHICGB_ASSERT(g1->termsAreInDescendingOrder());
+ HashPoly M;
mHashTableOLD.insert(multiplier, g1->begin(), g1->end(), M);
-
if (!M.empty())
- {
- mHeap.push(M.begin(),M.end());
- }
+ mHeap.push(M.begin(),M.end());
stats_n_inserts++;
stats_n_compares += mHeap.getConfiguration().getComparisons();
diff --git a/src/mathicgb/F4MatrixReducer.cpp b/src/mathicgb/F4MatrixReducer.cpp
index aea7a76..1bdfe05 100755
--- a/src/mathicgb/F4MatrixReducer.cpp
+++ b/src/mathicgb/F4MatrixReducer.cpp
@@ -10,11 +10,14 @@
#include <map>
#include <string>
#include <cstdio>
+#include <iostream>
#ifdef _OPENMP
#include <omp.h>
#endif
+extern int tracingLevel;
+
template<class T>
class DenseRow {
public:
@@ -762,6 +765,18 @@ void concatenateMatricesHorizontal
void F4MatrixReducer::reduce
(const PolyRing& ring, QuadMatrix& matrix, SparseMatrix& newPivots) {
+ if (tracingLevel >= 3)
+ std::cerr << "Row reducing (" << matrix.topLeft.rowCount()
+ << " + " << matrix.bottomLeft.rowCount()
+ << ") x (" << matrix.topLeft.colCount()
+ << " + " << matrix.topRight.colCount()
+ << ") matrix.\n#nz entry: top-left "
+ << matrix.topLeft.entryCount()
+ << ", top right " << matrix.topRight.entryCount()
+ << ", bot left " << matrix.bottomLeft.entryCount()
+ << ", bot right " << matrix.bottomRight.entryCount()
+ << '\n';
+
if (ring.charac() > std::numeric_limits<SparseMatrix::Scalar>::max())
throw std::overflow_error("Too large modulus in F4 matrix computation.");
diff --git a/src/mathicgb/F4Reducer.cpp b/src/mathicgb/F4Reducer.cpp
index 0b4b41a..fe5b7cc 100755
--- a/src/mathicgb/F4Reducer.cpp
+++ b/src/mathicgb/F4Reducer.cpp
@@ -3,6 +3,9 @@
#include "F4MatrixBuilder.hpp"
#include "F4MatrixReducer.hpp"
+#include <iostream>
+
+extern int tracingLevel;
F4Reducer::F4Reducer(
const PolyRing& ring,
@@ -14,6 +17,10 @@ F4Reducer::F4Reducer(
std::unique_ptr<Poly> F4Reducer::classicReduce
(const Poly& poly, const PolyBasis& basis) {
+ if (tracingLevel >= 2)
+ std::cerr <<
+ "F4Reducer: Using fall-back reducer for single classic reduction\n";
+
std::unique_ptr<Poly> p;
p = mFallback->classicReduce(poly, basis);
mSigStats = mFallback->sigStats();
@@ -24,6 +31,10 @@ std::unique_ptr<Poly> F4Reducer::classicReduce
std::unique_ptr<Poly> F4Reducer::classicTailReduce
(const Poly& poly, const PolyBasis& basis) {
std::unique_ptr<Poly> p;
+ if (tracingLevel >= 2)
+ std::cerr <<
+ "F4Reducer: Using fall-back reducer for single classic tail reduction\n";
+
p = mFallback->classicTailReduce(poly, basis);
mSigStats = mFallback->sigStats();
mClassicStats = mFallback->classicStats();
@@ -32,6 +43,9 @@ std::unique_ptr<Poly> F4Reducer::classicTailReduce
std::unique_ptr<Poly> F4Reducer::classicReduceSPoly
(const Poly& a, const Poly& b, const PolyBasis& basis) {
+ if (tracingLevel >= 2)
+ std::cerr << "F4Reducer: Reducing single S-pair.\n";
+
QuadMatrix qm;
{
F4MatrixBuilder builder(basis);
@@ -65,6 +79,9 @@ void F4Reducer::classicReduceSPolyGroup
if (spairs.empty())
return;
+ if (tracingLevel >= 2)
+ std::cerr << "F4Reducer: Reducing " << spairs.size() << " S-pair.\n";
+
SparseMatrix reduced;
std::vector<monomial> monomials;
{
@@ -84,6 +101,10 @@ void F4Reducer::classicReduceSPolyGroup
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);
diff --git a/src/mathicgb/Ideal.cpp b/src/mathicgb/Ideal.cpp
old mode 100644
new mode 100755
index 25b4ed7..1f388e4
--- a/src/mathicgb/Ideal.cpp
+++ b/src/mathicgb/Ideal.cpp
@@ -17,6 +17,8 @@ Ideal::~Ideal()
}
void Ideal::insert(std::unique_ptr<Poly> p) {
+ MATHICGB_ASSERT(p.get() != 0);
+ MATHICGB_ASSERT(p->termsAreInDescendingOrder());
mGenerators.reserve(mGenerators.size() + 1);
mGenerators.push_back(p.release());
}
@@ -40,21 +42,22 @@ void Ideal::sort(FreeModuleOrder& order) {
std::sort(mGenerators.begin(), mGenerators.end(), cmp);
}
-Ideal *Ideal::parse(std::istream &i)
+Ideal *Ideal::parse(std::istream& in)
{
- PolyRing *R = PolyRing::read(i);
+ PolyRing *R = PolyRing::read(in);
size_t npolys;
- i >> npolys;
+ in >> npolys;
Ideal *result = new Ideal(*R);
- for (size_t j= 0; j<npolys; j++)
- {
- Poly *g = new Poly(R);
- while (std::isspace(i.peek())) i.get();
- g->parse(i);
- result->mGenerators.push_back(g);
- }
+ for (size_t j= 0; j < npolys; j++) {
+ auto g = make_unique<Poly>(R);
+ while (std::isspace(in.peek()))
+ in.get();
+ g->parse(in);
+ result->insert(std::move(g));
+ }
return result;
}
+
void Ideal::display(std::ostream &o, bool print_comp) const
{
mRing.write(o);
diff --git a/src/mathicgb/Poly.cpp b/src/mathicgb/Poly.cpp
index 18243e1..c3cb019 100755
--- a/src/mathicgb/Poly.cpp
+++ b/src/mathicgb/Poly.cpp
@@ -3,8 +3,8 @@
#include "stdinc.h"
#include "Poly.hpp"
#include <ostream>
-#include <istream>
#include <iostream>
+#include <algorithm>
// Format for input/output:
// #terms term1 term2 ...
@@ -23,6 +23,56 @@ void Poly::copy(Poly &result) const
std::copy(monoms.begin(), monoms.end(), result.monoms.begin());
}
+void Poly::sortTermsDescending() {
+ struct Cmp {
+ public:
+ Cmp(const Poly& poly): mPoly(poly) {}
+
+ bool operator()(size_t a, size_t b) {
+ MATHICGB_ASSERT(a < mPoly.nTerms());
+ MATHICGB_ASSERT(b < mPoly.nTerms());
+ return mPoly.R->monomialLT(mPoly.monomialAt(b), mPoly.monomialAt(a));
+ }
+
+ private:
+ const Poly& mPoly;
+ };
+
+ const size_t count = nTerms();
+ std::vector<size_t> ordered(count);
+ for (size_t i = 0; i < count; ++i)
+ ordered[i] = i;
+
+ std::sort(ordered.begin(), ordered.end(), Cmp(*this));
+
+ Poly poly(R);
+ for (size_t i = 0; i < count; ++i)
+ poly.appendTerm(coefficientAt(ordered[i]), monomialAt(ordered[i]));
+ *this = std::move(poly);
+
+ MATHICGB_ASSERT(termsAreInDescendingOrder());
+}
+
+monomial Poly::monomialAt(size_t index) {
+ MATHICGB_ASSERT(index < nTerms());
+ return &monoms[index * R->maxMonomialSize()];
+}
+
+const_monomial Poly::monomialAt(size_t index) const {
+ MATHICGB_ASSERT(index < nTerms());
+ return &monoms[index * R->maxMonomialSize()];
+}
+
+coefficient& Poly::coefficientAt(size_t index) {
+ MATHICGB_ASSERT(index < nTerms());
+ return coeffs[index];
+}
+
+const coefficient Poly::coefficientAt(size_t index) const {
+ MATHICGB_ASSERT(index < nTerms());
+ return coeffs[index];
+}
+
void Poly::appendTerm(coefficient a, const_monomial m)
{
// the monomial will be copied on.
@@ -192,54 +242,55 @@ void Poly::dump() const
std::cout << std::endl;
}
-void Poly::parse(std::istream &i)
+void Poly::parseDoNotOrder(std::istream& i)
{
char next = i.peek();
- if (next == '0')
- {
+ if (next == '0') {
+ i.get();
+ return;
+ }
+
+ for (;;) {
+ bool is_neg = false;
+ next = i.peek();
+ if (next == '+') {
i.get();
- return;
- }
- for (;;)
- {
- bool is_neg = false;
next = i.peek();
- if (next == '+')
- {
- i.get();
- next = i.peek();
- }
- else if (next == '-')
- {
- is_neg = true;
+ } else if (next == '-') {
+ is_neg = true;
+ i.get();
+ next = i.peek();
+ } if (isdigit(next) || isalpha(next) || next == '<') {
+ int a = 1;
+ coefficient b;
+ if (isdigit(next))
+ {
+ i >> a;
+ next = i.peek();
+ }
+ if (is_neg) a = -a;
+ size_t firstloc = monoms.size();
+ monoms.resize(firstloc + R->maxMonomialSize());
+ monomial m = &monoms[firstloc];
+ R->coefficientFromInt(b,a);
+ coeffs.push_back(b);
+ if (isalpha(next) || next == '<')
+ R->monomialParse(i, m);
+ else
+ R->monomialSetIdentity(m); // have to do this to set hash value
+ MATHICGB_ASSERT(ring().hashValid(m));
+ next = i.peek();
+ if (next == '>')
i.get();
- next = i.peek();
- }
- if (isdigit(next) || isalpha(next) || next == '<')
- {
- int a = 1;
- coefficient b;
- if (isdigit(next))
- {
- i >> a;
- next = i.peek();
- }
- if (is_neg) a = -a;
- size_t firstloc = monoms.size();
- monoms.resize(firstloc + R->maxMonomialSize());
- monomial m = &monoms[firstloc];
- R->coefficientFromInt(b,a);
- coeffs.push_back(b);
- if (isalpha(next) || next == '<')
- R->monomialParse(i, m);
- else
- R->monomialSetIdentity(m); // have to do this to set hash value
- MATHICGB_ASSERT(ring().hashValid(m));
- next = i.peek();
- if (next == '>') i.get();
- }
- else return;
- }
+ }
+ else
+ break;
+ }
+}
+
+void Poly::parse(std::istream& in) {
+ parseDoNotOrder(in);
+ sortTermsDescending();
}
void Poly::display(std::ostream &o, bool print_comp) const
@@ -292,6 +343,24 @@ std::ostream& operator<<(std::ostream& out, const Poly& p) {
p.see(false);
}
+bool Poly::termsAreInDescendingOrder() const {
+ if (isZero())
+ return true;
+
+ auto stop = end();
+ auto it = begin();
+ auto previous = it;
+ ++it;
+ while (it != stop) {
+ if (R->monomialCompare(previous.getMonomial(), it.getMonomial()) == LT)
+ return false;
+ previous = it;
+ ++it;
+ }
+ return true;
+}
+
+
// Local Variables:
// compile-command: "make -C .. "
// indent-tabs-mode: nil
diff --git a/src/mathicgb/Poly.hpp b/src/mathicgb/Poly.hpp
index 674f13d..b64cb9e 100755
--- a/src/mathicgb/Poly.hpp
+++ b/src/mathicgb/Poly.hpp
@@ -12,7 +12,8 @@ class Poly {
public:
Poly(const PolyRing *R0) : R(R0) {MATHICGB_ASSERT(R != 0);}
- void parse(std::istream &i); // reads into this.
+ void parse(std::istream &i); // reads into this, sorts terms
+ void parseDoNotOrder(std::istream &i); // reads into this, does not sort terms
void display(std::ostream &o, bool print_comp=true) const;
void see(bool print_comp) const;
@@ -61,12 +62,20 @@ public:
}
};
+ void sortTermsDescending();
+
void append(iterator &first, iterator &last);
// Insert [first, last) onto the end of this.
// This invalidates iterators on this.
Poly *copy() const;
+ // Cannot call this monomial() since that is already a type :-(
+ monomial monomialAt(size_t index);
+ const_monomial monomialAt(size_t index) const;
+ coefficient& coefficientAt(size_t index);
+ const coefficient coefficientAt(size_t index) const;
+
void appendTerm(coefficient a, const_monomial m);
// all iterators are invalid after this
@@ -112,6 +121,8 @@ public:
void dump() const; // used for debugging
+ bool termsAreInDescendingOrder() const;
+
private:
const PolyRing *R;
std::vector<coefficient> coeffs;
diff --git a/src/mathicgb/PolyRing.cpp b/src/mathicgb/PolyRing.cpp
index 0f29bbd..594399c 100755
--- a/src/mathicgb/PolyRing.cpp
+++ b/src/mathicgb/PolyRing.cpp
@@ -96,7 +96,8 @@ void PolyRing::monomialSetIdentity(Monomial& result) const
void PolyRing::monomialEi(size_t i, Monomial &result) const
{
- for (size_t j=mTopIndex; j != static_cast<size_t>(-1); --j) result[j] = 0;
+ for (size_t j=mTopIndex; j != static_cast<size_t>(-1); --j)
+ result[j] = 0;
*result = static_cast<int>(i); // todo: handle overflow or change representation
result[mHashIndex] = static_cast<int>(i); // todo: handle overflow or change representation
}
diff --git a/src/test/F4MatrixReducer.cpp b/src/test/F4MatrixReducer.cpp
index 781df0a..d86e36d 100755
--- a/src/test/F4MatrixReducer.cpp
+++ b/src/test/F4MatrixReducer.cpp
@@ -9,13 +9,13 @@
#include <sstream>
TEST(F4MatrixReducer, Reduce) {
- auto ring = ringFromString("101 6 1\n1 1 1 1 1 1");
+ auto ring = ringFromString("101 6 1\n10 1 1 1 1 1");
QuadMatrix m;
m.ring = ring.get();
Poly p(ring.get());
- std::istringstream in("a1+a2+a3+a4+b1+b2+b3+b4+b5");
+ std::istringstream in("a4+a3+a2+a1+b5+b4+b3+b2+b1");
p.parse(in);
size_t count = 0;
for (Poly::iterator it = p.begin(); it != p.end(); ++it) {
@@ -100,8 +100,8 @@ TEST(F4MatrixReducer, Reduce) {
MATHICGB_ASSERT(m.debugAssertValid());
const char* origStr =
- "Left columns: a a2 a3 a4\n"
- "Right columns: b b2 b3 b4 b5\n"
+ "Left columns: a4 a3 a2 a\n"
+ "Right columns: b5 b4 b3 b2 b\n"
"0: 0#1 1#2 3#3 | 0: 2#8 \n"
"1: 1#1 2#3 | 1: 3#9 \n"
"2: 2#1 3#7 | 2: 4#10 \n"
diff --git a/src/test/QuadMatrixBuilder.cpp b/src/test/QuadMatrixBuilder.cpp
index 1ad2536..4c0543a 100755
--- a/src/test/QuadMatrixBuilder.cpp
+++ b/src/test/QuadMatrixBuilder.cpp
@@ -22,7 +22,7 @@ namespace {
{
Poly p(&b.ring());
std::istringstream in(left);
- p.parse(in);
+ p.parseDoNotOrder(in);
size_t colCount = 0;
for (Poly::iterator it = p.begin(); it != p.end(); ++it) {
QuadMatrixBuilder::ColIndex col = b.createColumnLeft(it.getMonomial());
@@ -39,7 +39,7 @@ namespace {
{
Poly p(&b.ring());
std::istringstream in(right);
- p.parse(in);
+ p.parseDoNotOrder(in);
size_t colCount = 0;
for (Poly::iterator it = p.begin(); it != p.end(); ++it) {
QuadMatrixBuilder::ColIndex col = b.createColumnRight(it.getMonomial());
@@ -71,7 +71,7 @@ TEST(QuadMatrixBuilder, Empty) {
TEST(QuadMatrixBuilder, Construction) {
std::unique_ptr<PolyRing> ring(ringFromString("32003 6 1\n1 1 1 1 1 1"));
QuadMatrixBuilder b(*ring);
- createColumns("a<1>+<0>", "b<0>+c<0>+bc<0>", b);
+ createColumns("a<1>+<0>", "bc<0>+b<0>+c<0>", b);
// top row: nothing, nothing
b.rowDoneTopLeftAndRight();
@@ -95,7 +95,7 @@ TEST(QuadMatrixBuilder, Construction) {
const char* matrixStr =
"Left columns: a 1\n"
- "Right columns: b c bc\n"
+ "Right columns: bc b c\n"
"0: | 0: \n"
"1: 0#1 1#2 | 1: 2#3\n"
" | \n"
@@ -114,7 +114,7 @@ TEST(QuadMatrixBuilder, ColumnQuery) {
// coefficient 1X=left, 2X=right, 30=not there, % 10 = column index
std::istringstream in
("10a<1>+11<0>+20b<0>+21c<0>+22bc<0>+30ab<0>+30e<0>+10a<1>");
- p.parse(in);
+ p.parseDoNotOrder(in);
for (Poly::iterator it = p.begin(); it != p.end(); ++it) {
QuadMatrixBuilder::LeftRightColIndex col = b.findColumn(it.getMonomial());
if (it.getCoefficient() / 10 == 3)
diff --git a/src/test/poly-test.cpp b/src/test/poly-test.cpp
index 2edc069..953fa86 100755
--- a/src/test/poly-test.cpp
+++ b/src/test/poly-test.cpp
@@ -57,7 +57,7 @@ TEST(Poly,readwrite) {
std::unique_ptr<PolyRing> R(ringFromString("32003 6 1\n1 1 1 1 1 1"));
Poly f(R.get());
std::stringstream ifil(f1);
- f.parse(ifil);
+ f.parseDoNotOrder(ifil);
std::ostringstream o;
f.display(o,true);
EXPECT_EQ(o.str(), f1);
@@ -68,7 +68,7 @@ bool testPolyParse(PolyRing* R, std::string s)
// parse poly, then see if it matches the orig string
Poly f(R);
std::istringstream i(s);
- f.parse(i);
+ f.parseDoNotOrder(i);
std::ostringstream o;
f.display(o);
// std::cout << "orig = " << s << std::endl;
@@ -80,7 +80,7 @@ bool testPolyParse2(PolyRing* R, std::string s, std::string answer)
// parse poly, then see if it matches the orig string
Poly f(R);
std::istringstream i(s);
- f.parse(i);
+ f.parseDoNotOrder(i);
std::ostringstream o;
f.display(o);
// std::cout << "orig = " << s << std::endl;
--
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