[mathicgb] 48/393: Now only add 1 bottom row to F4 matrices per S-pair.

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 27f3d7169988a51b1628f4d97fe7c6ce2f68c609
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date:   Fri Oct 5 16:53:43 2012 +0200

    Now only add 1 bottom row to F4 matrices per S-pair.
---
 src/mathicgb/F4MatrixBuilder.cpp | 129 ++++++++++++++++++++++++++++++---------
 src/mathicgb/F4MatrixBuilder.hpp |  50 +++++++++------
 src/mathicgb/F4Reducer.cpp       |   4 +-
 src/test/F4MatrixBuilder.cpp     |  62 ++++++++++++-------
 4 files changed, 174 insertions(+), 71 deletions(-)

diff --git a/src/mathicgb/F4MatrixBuilder.cpp b/src/mathicgb/F4MatrixBuilder.cpp
index ad33cf8..a680c11 100755
--- a/src/mathicgb/F4MatrixBuilder.cpp
+++ b/src/mathicgb/F4MatrixBuilder.cpp
@@ -4,7 +4,8 @@
 F4MatrixBuilder::F4MatrixBuilder(const PolyBasis& basis):
   mBasis(basis), mBuilder(basis.ring()) {}
 
-void F4MatrixBuilder::addTwoRowsForSPairToMatrix(const Poly& polyA, const Poly& polyB) {
+void F4MatrixBuilder::addSPolynomialToMatrix
+(const Poly& polyA, const Poly& polyB) {
   MATHICGB_ASSERT(!polyA.isZero());
   MATHICGB_ASSERT(!polyB.isZero());
 
@@ -12,28 +13,35 @@ void F4MatrixBuilder::addTwoRowsForSPairToMatrix(const Poly& polyA, const Poly&
   ring().monomialLeastCommonMultiple
     (polyA.getLeadMonomial(), polyB.getLeadMonomial(), lcm);
 
-  monomial multiple = ring().allocMonomial();
+  SPairTask task = {};
 
-  ring().monomialDivide(lcm, polyA.getLeadMonomial(), multiple);
-  addRowToMatrix(multiple, polyA);
+  task.polyA = &polyA;
+  task.multiplyA = ring().allocMonomial();
+  ring().monomialDivide(lcm, polyA.getLeadMonomial(), task.multiplyA);
 
-  ring().monomialDivide(lcm, polyB.getLeadMonomial(), multiple);
-  addRowToMatrix(multiple, polyB);
+  task.polyB = &polyB;
+  task.multiplyB = ring().allocMonomial();
+  ring().monomialDivide(lcm, polyB.getLeadMonomial(), task.multiplyB);
 
+  mSPairTodo.push_back(task);
   ring().freeMonomial(lcm);
-  ring().freeMonomial(multiple);
 }
 
-void F4MatrixBuilder::addRowToMatrix
+void F4MatrixBuilder::addPolynomialToMatrix
 (const_monomial multiple, const Poly& poly) {
-  RowTask task;
-  task.useAsReducer = false; // to be updated later
-  task.poly = &poly;
+  MATHICGB_ASSERT(ring().hashValid(multiple));
+  if (poly.isZero())
+    return;
 
-  task.multiple = ring().allocMonomial();
-  ring().monomialCopy(multiple, task.multiple);
-  MATHICGB_ASSERT(ring().hashValid(task.multiple));
-  mTodo.push_back(task);
+  SPairTask task = {};
+
+  task.polyA = &poly;
+  task.multiplyA = ring().allocMonomial();
+  ring().monomialCopy(multiple, task.multiplyA);
+
+  MATHICGB_ASSERT(task.polyB == 0);
+  MATHICGB_ASSERT(task.multiplyB.unsafeGetRepresentation() == 0);
+  mSPairTodo.push_back(task);
 }
 
 void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
@@ -41,21 +49,83 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
   // todo: prefer sparse/old reducers among the inputs.
   monomial mono = ring().allocMonomial();
 
-  // Decide which input rows are going to be used as reducers and
-  // create pivot columns ahead of time for those reducers so that
-  // we do not add a reducer for those columns later on.
-  typedef std::vector<RowTask>::iterator TaskIter;
-  TaskIter end = mTodo.end();
-  for (TaskIter it = mTodo.begin(); it != end; ++it) {
-    MATHICGB_ASSERT(ring().hashValid(it->multiple));
-    ring().monomialMult(it->multiple, it->poly->getLeadMonomial(), mono);
-    LeftRightColIndex leadCol = mBuilder.findColumn(mono);
-    it->useAsReducer = !leadCol.valid();
-    if (it->useAsReducer) {
-      // create column so we know later on that we already have a
-      // reducer for this column.
-      mBuilder.createColumnLeft(mono);
+  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);
+    } else {
+      // set polyB to an empty range if null
+      itB = endA;
+      endB = endA;
+      MATHICGB_ASSERT(itB == endB);
+    }
+
+    // skip leading terms since they cancel
+    //++itA;
+    //++itB;
+
+    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)
+        mBuilder.appendEntryBottom(createOrFindColumnOf(mono), c);
     }
+    ring().freeMonomial(monoA);
+    ring().freeMonomial(monoB);
+    mBuilder.rowDoneBottomLeftAndRight();
   }
 
   // Process pending rows until we are done. Note that the methods
@@ -77,6 +147,7 @@ void F4MatrixBuilder::buildMatrixAndClear(QuadMatrix& matrix) {
 
 F4MatrixBuilder::LeftRightColIndex
 F4MatrixBuilder::createOrFindColumnOf(const_monomial mono) {
+  MATHICGB_ASSERT(ring().hashValid(mono));
   LeftRightColIndex colIndex = mBuilder.findColumn(mono);
   if (colIndex.valid())
     return colIndex;
diff --git a/src/mathicgb/F4MatrixBuilder.hpp b/src/mathicgb/F4MatrixBuilder.hpp
index 86632f3..c2dc425 100755
--- a/src/mathicgb/F4MatrixBuilder.hpp
+++ b/src/mathicgb/F4MatrixBuilder.hpp
@@ -23,35 +23,35 @@ private:
 public:
   F4MatrixBuilder(const PolyBasis& basis);
 
-  /** Schedules two rows to be added to the matrix whose linear span
-    includes the S-polynomial between polyA and polyB. More precisely,
-    the two rows represent (B:A)*polyA and (A:B)*polyB where
-    A=lead(polyA) and B=lead(polyB). */
-  void addTwoRowsForSPairToMatrix(const Poly& polyA, const Poly& polyB);
+  /** Schedules a row representing the S-polynomial between polyA and
+    polyB to be added to the matrix. No ownership is taken, but polyA
+    and polyB must remain valid until the matrix is constructed.
+
+    Currently, the two monomials must be monic, though this is just
+    because they always happen to be monic so there was no reason to
+    support the non-monic case. */
+  void addSPolynomialToMatrix(const Poly& polyA, const Poly& polyB);
 
   /** Schedules a row representing multiple*poly to be added to the
-      matrix. No ownership is taken, but poly must remain valid until
-      the matrix is constructed. multiple is copied so there is no
-      requirement there. */
-  void addRowToMatrix(const_monomial multiple, const Poly& poly);
+    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);
 
   /** Builds an F4 matrix to the specifications given. Also clears the
     information in this object.
 
     The right columns are ordered by decreasing monomial of that
     column according to the order from the basis. The left columns are
-    order in some way so that the first entry in each row (the pivot)
-    has a lower index than any other entries in that row.
+    ordered in some way so that the first entry in each top row (the
+    pivot) has a lower index than any other entries in that row.
 
     The matrix contains a reducer/pivot for every monomial that can be
-    reduced by the basis and that is present in the matrix. Note that
-    the client-added rows also count as reducers so their lead terms
-    will not get another reducer added automatically -- specifically,
-    adding an S-polynomial will not do what you want because its lead
-    term may have no reducer in the matrix other than itself. Instead,
-    add the two polynomials that you would have subtracted from each
-    other to form the S-polynomial - or even better call the method
-    that adds an S-pair for you. */
+    reduced by the basis and that is present in the matrix. There is
+    no guarantee that the bottom part of the matrix contains rows that
+    exactly correspond to the polynomials that have been scheduled to
+    be added to the matrix. It is only guaranteed that the matrix has
+    the same row-space as though that had been the case. */
   void buildMatrixAndClear(QuadMatrix& matrix);
 
   const PolyRing& ring() const {return mBuilder.ring();}
@@ -64,12 +64,24 @@ private:
   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
+  /// for reduction.
+  struct SPairTask {
+    const Poly* polyA;
+    monomial multiplyA;
+    const Poly* polyB;
+    monomial multiplyB;
+  };
+
   /// 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;
   };
+
+  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 e03c21c..eca815a 100755
--- a/src/mathicgb/F4Reducer.cpp
+++ b/src/mathicgb/F4Reducer.cpp
@@ -49,7 +49,7 @@ std::unique_ptr<Poly> F4Reducer::classicReduceSPoly
   QuadMatrix qm;
   {
     F4MatrixBuilder builder(basis);
-    builder.addTwoRowsForSPairToMatrix(a, b);
+    builder.addSPolynomialToMatrix(a, b);
     builder.buildMatrixAndClear(qm);
 
     // there has to be something to reduce
@@ -89,7 +89,7 @@ void F4Reducer::classicReduceSPolySet
     {
       F4MatrixBuilder builder(basis);
       for (auto it = spairs.begin(); it != spairs.end(); ++it) {
-        builder.addTwoRowsForSPairToMatrix
+        builder.addSPolynomialToMatrix
           (basis.poly(it->first), basis.poly(it->second));
       }
       builder.buildMatrixAndClear(qm);
diff --git a/src/test/F4MatrixBuilder.cpp b/src/test/F4MatrixBuilder.cpp
index d475592..86a9eb2 100755
--- a/src/test/F4MatrixBuilder.cpp
+++ b/src/test/F4MatrixBuilder.cpp
@@ -64,26 +64,45 @@ TEST(F4MatrixBuilder, Empty) {
   ASSERT_EQ(0, matrix.rightColumnMonomials.size());
 }
 
+TEST(F4MatrixBuilder, SPair) {
+  BuilderMaker maker;
+  const Poly& p1 = maker.addBasisElement("a4c2-d");
+  const Poly& p2 = maker.addBasisElement("a4b+d");
+  // S-pair of p1 and p2 is -c2d-bd
+  const Poly& p3 = maker.addBasisElement("c2d+3");
+  F4MatrixBuilder& builder = maker.create();
+  builder.addSPolynomialToMatrix(p1, p2);
+  QuadMatrix qm;
+  builder.buildMatrixAndClear(qm);
+  const char* str = 
+    "Left columns: c2d\n"
+    "Right columns: bd 1\n"
+    "0: 0#1   | 0: 1#3  \n"
+    "         |         \n"
+    "0: 0#100 | 0: 0#100\n";
+  ASSERT_EQ(str, qm.toString());
+}
+
 TEST(F4MatrixBuilder, OneByOne) {
   BuilderMaker maker;
   const Poly& p = maker.addBasisElement("a");
   F4MatrixBuilder& builder = maker.create();
-  builder.addRowToMatrix(p.getLeadMonomial(), p);
+  builder.addPolynomialToMatrix(p.getLeadMonomial(), p);
   QuadMatrix qm;
   builder.buildMatrixAndClear(qm);
   const char* str = 
     "Left columns: a2\n"
     "Right columns:\n"
-    "0: 0#1              | 0:                 \n"
-    "                    |                    \n"
-    "matrix with no rows | matrix with no rows\n";
-  ASSERT_EQ(str, qm.toString());
+    "0: 0#1 | 0:\n"
+    "       |   \n"
+    "0: 0#1 | 0:\n";
+  ASSERT_EQ(str, qm.toString()) << "** qm:\n" << qm;
 }
 
 TEST(F4MatrixBuilder, DirectReducers) {
   BuilderMaker maker;
-  maker.addBasisElement("a6<0>"); // reducer ==, but won't be used as not added
-  maker.addBasisElement("a3b2<0>+a3c"); // reducer ==
+  maker.addBasisElement("a6<0>"); // reducer == to lead term
+  maker.addBasisElement("a3b2<0>+a3c"); // reducer == to lower order term
   maker.addBasisElement("c<0>"); // reducer divides
   maker.addBasisElement("d2<0>"); // does not divide
   F4MatrixBuilder& builder = maker.create();
@@ -92,14 +111,14 @@ TEST(F4MatrixBuilder, DirectReducers) {
   { 
     std::istringstream in("a3<0>+b2+c+d");
     p1.parse(in);
-    builder.addRowToMatrix(p1.getLeadMonomial(), p1);
+    builder.addPolynomialToMatrix(p1.getLeadMonomial(), p1);
   }
 
   Poly p2(&builder.ring());
   {
     std::istringstream in("a3<0>+2b2+3c+4d");
     p2.parse(in);
-    builder.addRowToMatrix(p2.getLeadMonomial(), p2);
+    builder.addPolynomialToMatrix(p2.getLeadMonomial(), p2);
   }
 
   QuadMatrix qm;
@@ -110,12 +129,13 @@ TEST(F4MatrixBuilder, DirectReducers) {
     "Right columns: a3d\n"
     "0: 2#1         | 0:    \n"
     "1: 1#1 2#1     | 1:    \n"
-    "2: 0#1 1#1 2#1 | 2: 0#1\n"
+    "2: 0#1         | 2:    \n"
     "               |       \n"
-    "0: 0#1 1#2 2#3 | 0: 0#4\n";
+    "0: 0#1 1#1 2#1 | 0: 0#1\n"
+    "1: 0#1 1#2 2#3 | 1: 0#4\n";
   // This quest is currently fragile because of the possibility of
   // reordering of the rows.
-  ASSERT_EQ(str, qm.toString());
+  ASSERT_EQ(str, qm.toString()) << "** qm:\n" << qm;
 }
 
 TEST(F4MatrixBuilder, IteratedReducer) {
@@ -123,18 +143,18 @@ TEST(F4MatrixBuilder, IteratedReducer) {
   const Poly& p1 = maker.addBasisElement("a4-a3");
   const Poly& p2 = maker.addBasisElement("a-1");
   F4MatrixBuilder& builder = maker.create();
-  builder.addRowToMatrix(p1.getLeadMonomial(), p2);
+  builder.addPolynomialToMatrix(p1.getLeadMonomial(), p2);
   QuadMatrix qm;
   builder.buildMatrixAndClear(qm);
   const char* str = 
     "Left columns: a5 a4 a3 a2 a\n"
     "Right columns: 1\n"
-    "0: 0#1 1#100        | 0:                 \n"
-    "1: 1#1 2#100        | 1:                 \n"
-    "2: 2#1 3#100        | 2:                 \n"
-    "3: 3#1 4#100        | 3:                 \n"
-    "4: 4#1              | 4: 0#100           \n"
-    "                    |                    \n"
-    "matrix with no rows | matrix with no rows\n";
-  ASSERT_EQ(str, qm.toString());
+    "0: 1#1 2#100 | 0:      \n"
+    "1: 2#1 3#100 | 1:      \n"
+    "2: 3#1 4#100 | 2:      \n"
+    "3: 4#1       | 3: 0#100\n"
+    "4: 0#1 1#100 | 4:      \n"
+    "             |         \n"
+    "0: 0#1 1#100 | 0:      \n";
+  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