[mathicgb] 188/393: Added a library interface to MathicGB along with tests.

Doug Torrance dtorrance-guest at moszumanska.debian.org
Fri Apr 3 15:58:58 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 b7027fcdccb97371745d4585c2e52b76f069999e
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date:   Fri Mar 15 18:59:47 2013 +0100

    Added a library interface to MathicGB along with tests.
---
 Makefile.am                                        |   6 +-
 build/vs12/mathicgb-test/mathicgb-test.vcxproj     |   1 +
 .../mathicgb-test/mathicgb-test.vcxproj.filters    |   3 +
 src/mathicgb.cpp                                   | 560 +++++++++++++++++++++
 src/mathicgb.h                                     | 523 +++++++++++++++++++
 src/mathicgb/PolyBasis.cpp                         |  18 +
 src/mathicgb/PolyBasis.hpp                         |  16 +-
 src/mathicgb/PolyRing.cpp                          |   6 +-
 src/mathicgb/PolyRing.hpp                          |   8 +-
 src/test/mathicgb.cpp                              | 249 +++++++++
 10 files changed, 1373 insertions(+), 17 deletions(-)

diff --git a/Makefile.am b/Makefile.am
index 566d38c..8baf096 100755
--- a/Makefile.am
+++ b/Makefile.am
@@ -68,7 +68,8 @@ libmathicgb_la_SOURCES = src/mathicgb/BjarkeGeobucket2.cpp		\
   src/mathicgb/F4MatrixBuilder2.hpp src/mathicgb/F4MatrixBuilder2.cpp \
   src/mathicgb/LogDomainSet.cpp src/mathicgb/F4ProtoMatrix.hpp \
   src/mathicgb/F4ProtoMatrix.cpp src/mathicgb/F4MatrixProject.hpp \
-  src/mathicgb/F4MatrixProjection.cpp src/mathicgb/ScopeExit.hpp
+  src/mathicgb/F4MatrixProjection.cpp src/mathicgb/ScopeExit.hpp \
+  src/mathicgb.cpp src/mathicgb.h
 
 # When making a distribution file, Automake knows to include all files
 # that are necessary to build the project. EXTRA_DIST specifies files
@@ -104,4 +105,5 @@ unittest_SOURCES=src/test/FreeModuleOrderTest.cpp						\
   src/test/gtestInclude.cpp src/test/testMain.cpp src/test/gb-test.cpp	\
   src/test/ideals.cpp src/test/poly-test.cpp src/test/ideals.hpp		\
   src/test/SparseMatrix.cpp src/test/QuadMatrixBuilder.cpp				\
-  src/test/F4MatrixBuilder.cpp src/test/F4MatrixReducer.cpp
+  src/test/F4MatrixBuilder.cpp src/test/F4MatrixReducer.cpp             \
+  src/test/mathicgb.cpp
diff --git a/build/vs12/mathicgb-test/mathicgb-test.vcxproj b/build/vs12/mathicgb-test/mathicgb-test.vcxproj
index 3a6e5e5..6c53744 100755
--- a/build/vs12/mathicgb-test/mathicgb-test.vcxproj
+++ b/build/vs12/mathicgb-test/mathicgb-test.vcxproj
@@ -402,6 +402,7 @@
     <ClCompile Include="..\..\..\src\test\gb-test.cpp" />
     <ClCompile Include="..\..\..\src\test\gtestInclude.cpp" />
     <ClCompile Include="..\..\..\src\test\ideals.cpp" />
+    <ClCompile Include="..\..\..\src\test\mathicgb.cpp" />
     <ClCompile Include="..\..\..\src\test\poly-test.cpp" />
     <ClCompile Include="..\..\..\src\test\QuadMatrixBuilder.cpp" />
     <ClCompile Include="..\..\..\src\test\SparseMatrix.cpp" />
diff --git a/build/vs12/mathicgb-test/mathicgb-test.vcxproj.filters b/build/vs12/mathicgb-test/mathicgb-test.vcxproj.filters
index 1a7a20f..eb741c1 100755
--- a/build/vs12/mathicgb-test/mathicgb-test.vcxproj.filters
+++ b/build/vs12/mathicgb-test/mathicgb-test.vcxproj.filters
@@ -45,6 +45,9 @@
     <ClCompile Include="..\..\..\src\test\testMain.cpp">
       <Filter>Source Files</Filter>
     </ClCompile>
+    <ClCompile Include="..\..\..\src\test\mathicgb.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="..\..\..\src\test\ideals.hpp">
diff --git a/src/mathicgb.cpp b/src/mathicgb.cpp
new file mode 100644
index 0000000..f5e03ad
--- /dev/null
+++ b/src/mathicgb.cpp
@@ -0,0 +1,560 @@
+#include "mathicgb/stdinc.h"
+#include "mathicgb.h"
+
+#include "mathicgb/Ideal.hpp"
+#include "mathicgb/PolyRing.hpp"
+#include "mathicgb/Poly.hpp"
+#include "mathicgb/Reducer.hpp"
+#include "mathicgb/BuchbergerAlg.hpp"
+#include <mathic.h>
+
+namespace {
+  bool isPrime(unsigned int n) {
+    if (n == 0 || n == 1)
+      return false;
+    if (n == 2 || n == 3)
+      return true;
+    return true; // todo: make better test
+  }
+};
+
+#ifndef MATHICGB_ASSERT
+#ifdef MATHICGB_DEBUG
+#include <cassert>
+#define MATHICGB_ASSERT(X) assert(X)
+#else
+#define MATHICGB_ASSERT(X)
+#endif
+#endif
+
+#define MATHICGB_STREAM_CHECK(X, MSG) \
+  do { \
+    const bool value = (X); \
+    if (!value) { \
+      const bool ignoreMe = false; \
+      MATHICGB_ASSERT(( \
+        "MathicGB stream protocol error: "#MSG \
+        "\nAssert expression: "#X"\n", \
+        false \
+      )); \
+      throw std::invalid_argument( \
+        "MathicGB stream protocol error: "#MSG \
+        "\nAssert expression: "#X"\n" \
+      ); \
+    } \
+  } while (false)
+
+namespace mgbi {
+  struct StreamStateChecker::Pimpl {
+    Pimpl(Coefficient modulus, VarIndex varCount):
+      modulus(modulus),
+      varCount(varCount),
+      state(Initial),
+      claimedPolyCount(0),
+      seenPolyCount(0),
+      claimedTermCount(0),
+      seenTermCount(0)
+    {}
+
+    bool debugAssertValid() const;
+
+    enum State {
+      Initial,
+      MakingIdeal,
+      MakingPoly,
+      MakingTerm,
+      HasIdeal
+    };
+
+    const Coefficient modulus;
+    const VarIndex varCount;
+
+    State state;
+    size_t claimedPolyCount;
+    size_t seenPolyCount;
+    size_t claimedTermCount;
+    size_t seenTermCount;
+    VarIndex lastVar;
+  };
+
+  bool StreamStateChecker::Pimpl::debugAssertValid() const {
+#ifdef MATHICGB_DEBUG
+    MATHICGB_ASSERT(this != 0);
+    switch (state) {
+    case Initial:
+    case MakingIdeal:
+    case MakingPoly:
+    case MakingTerm:
+    case HasIdeal:
+      break;
+
+    default:
+      MATHICGB_ASSERT(false);
+      return false;
+    }
+#endif
+    return true;
+  };
+
+  StreamStateChecker::StreamStateChecker(const Coefficient modulus, const VarIndex varCount):
+    mPimpl(new Pimpl(modulus, varCount))
+  {
+    try {
+      MATHICGB_STREAM_CHECK(isPrime(modulus), "The modulus must be prime");
+      MATHICGB_ASSERT(mPimpl->debugAssertValid());
+    } catch (...) {
+      delete mPimpl;
+    }
+  }
+
+  StreamStateChecker::~StreamStateChecker() {
+    MATHICGB_ASSERT(mPimpl->debugAssertValid());
+    delete mPimpl;
+  }
+
+  void StreamStateChecker::idealBegin(size_t polyCount) {
+    MATHICGB_ASSERT(mPimpl->debugAssertValid());
+
+    MATHICGB_STREAM_CHECK(
+      mPimpl->state == Pimpl::Initial || mPimpl->state == Pimpl::HasIdeal,
+      "idealBegin() must not be called twice "
+      "without an intervening call to idealDone()."
+    );
+    mPimpl->state = Pimpl::MakingIdeal;
+    mPimpl->claimedPolyCount = polyCount;
+    mPimpl->seenPolyCount = 0;
+
+    MATHICGB_ASSERT(mPimpl->debugAssertValid());
+  }
+
+  void StreamStateChecker::appendPolynomialBegin(size_t termCount) {
+    MATHICGB_ASSERT(mPimpl->debugAssertValid());
+
+    MATHICGB_STREAM_CHECK(
+      mPimpl->state != Pimpl::Initial && mPimpl->state != Pimpl::HasIdeal,
+      "appendPolynomialBegin() must only be called after idealBegin() and "
+      "before idealEnd()."
+    );
+    MATHICGB_STREAM_CHECK(
+      mPimpl->state == Pimpl::MakingIdeal,
+     "appendPolynomialBegin() must not be called twice without "
+     "an intervening call to appendPolynomialDone()."
+    );
+    MATHICGB_STREAM_CHECK(
+      mPimpl->seenPolyCount < mPimpl->claimedPolyCount,
+      "The number of polynomials in an ideal must not exceed the amount "
+      "passed to idealBegin()."
+    );
+    mPimpl->state = Pimpl::MakingPoly;
+    mPimpl->seenPolyCount += 1;
+    mPimpl->claimedTermCount = termCount;
+    mPimpl->seenTermCount = 0;
+
+    MATHICGB_ASSERT(mPimpl->debugAssertValid());
+  }
+
+  void StreamStateChecker::appendTermBegin() {
+    MATHICGB_ASSERT(mPimpl->debugAssertValid());
+
+    MATHICGB_STREAM_CHECK(
+      mPimpl->state != Pimpl::Initial &&
+        mPimpl->state != Pimpl::HasIdeal &&
+        mPimpl->state != Pimpl::MakingIdeal,
+      "appendTermBegin() must only be called after appendPolynomialBegin() "
+      "and before appendPolynomialDone()."
+    );
+    MATHICGB_STREAM_CHECK(
+      mPimpl->state == Pimpl::MakingPoly,
+      "appendTermBegin() must not be called twice without an intervening "
+      "call to appendTermDone()."
+    );
+    MATHICGB_STREAM_CHECK(
+      mPimpl->seenTermCount < mPimpl->claimedTermCount,
+      "The number of terms in a polynomial must not exceed the amount "
+      "passed to appendPolynomialBegin()."
+    );
+     
+    mPimpl->state = Pimpl::MakingTerm;
+    mPimpl->seenTermCount += 1;
+    mPimpl->lastVar = std::numeric_limits<decltype(mPimpl->lastVar)>::max();
+
+    MATHICGB_ASSERT(mPimpl->debugAssertValid());
+  }
+
+  void StreamStateChecker::appendExponent(VarIndex index, Exponent exponent) {
+    MATHICGB_ASSERT(mPimpl->debugAssertValid());
+
+    MATHICGB_STREAM_CHECK(
+      mPimpl->state == Pimpl::MakingTerm,
+      "appendExponent must only be called after appendTermBegin() and before "
+      "appendTermDone()."
+    );
+    MATHICGB_STREAM_CHECK(
+      index < mPimpl->varCount,
+      "The index passed to appendExponent must be strictly less than "
+      "the number of variables."
+    );
+    MATHICGB_STREAM_CHECK(
+      mPimpl->lastVar ==
+        std::numeric_limits<decltype(mPimpl->lastVar)>::max() ||
+      mPimpl->lastVar < index,
+      "The variable indices passed to appendExponent must be in strictly "
+      "increasing order within each monomial."
+    );
+
+    mPimpl->lastVar = index;
+    
+    MATHICGB_ASSERT(mPimpl->debugAssertValid());
+  }
+
+  void StreamStateChecker::appendTermDone(Coefficient coefficient) {
+    MATHICGB_ASSERT(mPimpl->debugAssertValid());
+
+    MATHICGB_STREAM_CHECK(
+      coefficient > 0,
+      "The coefficient passed to appendTermDone() must be strictly positive."
+    );
+    MATHICGB_STREAM_CHECK(
+      coefficient < mPimpl->modulus,
+      "The coefficient passed to appendTermDone() must be strictly less "
+      "then the modulus."
+    );
+    MATHICGB_STREAM_CHECK(
+      mPimpl->state == Pimpl::MakingTerm,
+      "appendTermDone() must only be called after appendTermBegin()."
+    );
+    mPimpl->state = Pimpl::MakingPoly;
+    
+    MATHICGB_ASSERT(mPimpl->debugAssertValid());
+  }
+
+  void StreamStateChecker::appendPolynomialDone() {
+    MATHICGB_ASSERT(mPimpl->debugAssertValid());
+
+    MATHICGB_STREAM_CHECK(
+      mPimpl->state == Pimpl::MakingPoly,
+      "appendPolynomialDone() must only be called after appendPolynomialBegin()."
+    );
+    MATHICGB_STREAM_CHECK(
+      mPimpl->seenTermCount == mPimpl->claimedTermCount,
+      "The number of terms in a polynomial must match the amount "
+      "passed to appendPolynomialBegin()."
+    );
+    mPimpl->state = Pimpl::MakingIdeal;
+
+    MATHICGB_ASSERT(mPimpl->debugAssertValid());
+  }
+
+  void StreamStateChecker::idealDone() {
+    MATHICGB_ASSERT(mPimpl->debugAssertValid());
+
+    MATHICGB_STREAM_CHECK(
+      mPimpl->state == Pimpl::MakingIdeal,
+      "idealDone() must only be called after idealBegin()."
+    );
+    MATHICGB_STREAM_CHECK(
+      mPimpl->seenPolyCount == mPimpl->claimedPolyCount,
+      "The number of polynomials in an ideal must match the amount "
+      "passed to idealBegin()."
+    );
+
+    mPimpl->state = Pimpl::HasIdeal;
+
+    MATHICGB_ASSERT(mPimpl->debugAssertValid());
+  }
+
+  bool StreamStateChecker::hasIdeal() const {
+    MATHICGB_ASSERT(mPimpl->debugAssertValid());
+    return mPimpl->state == Pimpl::HasIdeal;
+  }
+}
+
+namespace mgb {
+  // ** Implementation of the class GroebnerConfiguration
+
+  struct GroebnerConfiguration::Pimpl {
+    Pimpl(Coefficient modulus, VarIndex varCount):
+#ifdef MATHICGB_DEBUG
+      hasBeenDestroyed(false),
+#endif
+      modulus(modulus),
+      varCount(varCount)
+    {
+    }
+
+    MATHICGB_IF_DEBUG(bool hasBeenDestroyed);
+    const Coefficient modulus;
+    const VarIndex varCount;
+  };
+
+  GroebnerConfiguration::GroebnerConfiguration(
+    Coefficient modulus,
+    VarIndex varCount
+  ):
+    mPimpl(new Pimpl(modulus, varCount))
+  {
+    if (modulus > std::numeric_limits<unsigned short>::max()) {
+      MATHICGB_ASSERT_NO_ASSUME(false);
+      std::ostringstream str;
+      str << "Modulus " << modulus
+        << " is too large. MathicGB only supports 16 bit moduli.";
+      mathic::reportError(str.str());
+    }
+    if (!isPrime(modulus)) {
+      MATHICGB_ASSERT_NO_ASSUME(false);
+      std::ostringstream str;
+      str << "Modulus " << modulus
+        << " is not prime. MathicGB only supports prime fields.";
+      mathic::reportError(str.str());
+    }
+  }
+
+  GroebnerConfiguration::GroebnerConfiguration(
+    const GroebnerConfiguration& conf
+  ):
+    mPimpl(new Pimpl(*conf.mPimpl))
+  {
+    MATHICGB_ASSERT(conf.mPimpl->modulus != 0);
+    MATHICGB_ASSERT_NO_ASSUME(!conf.mPimpl->hasBeenDestroyed);
+  }
+
+  GroebnerConfiguration::~GroebnerConfiguration() {
+    MATHICGB_ASSERT(mPimpl != 0);
+    MATHICGB_ASSERT_NO_ASSUME(!mPimpl->hasBeenDestroyed);
+    MATHICGB_IF_DEBUG(mPimpl->hasBeenDestroyed = true);
+
+    delete mPimpl;
+  }
+
+  auto GroebnerConfiguration::modulus() const -> Coefficient {
+    return mPimpl->modulus;
+  }
+ 
+  auto GroebnerConfiguration::varCount() const -> VarIndex {
+    return mPimpl->varCount;
+  }
+}
+
+// ** Implementation of class GroebnerInputIdealStream
+namespace mgb {
+  struct GroebnerInputIdealStream::Pimpl {
+    Pimpl(const GroebnerConfiguration& conf):
+#ifdef MATHICGB_DEBUG
+      hasBeenDestroyed(false),
+      checker(conf.modulus(), conf.varCount()),
+#endif MATHICGB_DEBUG
+      // @todo: varCount should not be int. Fix PolyRing constructor,
+      // then remove this static_cast.
+      ring(conf.modulus(), static_cast<int>(conf.varCount()), 1),
+      ideal(ring),
+      poly(ring),
+      monomial(ring.allocMonomial()),
+      conf(conf)
+    {}
+
+    ~Pimpl() {
+      ring.freeMonomial(monomial);
+    }
+
+    MATHICGB_IF_DEBUG(bool hasBeenDestroyed);
+    MATHICGB_IF_DEBUG(::mgbi::StreamStateChecker checker);
+    const PolyRing ring;
+    Ideal ideal;
+    Poly poly;
+    Monomial monomial;
+    const GroebnerConfiguration conf;
+  };
+
+  GroebnerInputIdealStream::GroebnerInputIdealStream(
+    const GroebnerConfiguration& conf
+  ):
+    mExponents(new Exponent[conf.varCount()]),
+    mPimpl(new Pimpl(conf))
+  {
+    MATHICGB_ASSERT(debugAssertValid());
+  }
+
+  GroebnerInputIdealStream::~GroebnerInputIdealStream() {
+    MATHICGB_ASSERT(debugAssertValid());
+    MATHICGB_ASSERT(mExponents != 0);
+    MATHICGB_ASSERT(mPimpl != 0);
+    MATHICGB_ASSERT(!mPimpl->hasBeenDestroyed);
+    MATHICGB_IF_DEBUG(mPimpl->hasBeenDestroyed = true);
+    delete mPimpl;
+    delete[] mExponents;
+  }
+
+  const GroebnerConfiguration& GroebnerInputIdealStream::configuration() {
+    return mPimpl->conf;
+  }
+
+  auto GroebnerInputIdealStream::modulus() const -> Coefficient {
+    return mPimpl->conf.modulus();
+  }
+
+  auto GroebnerInputIdealStream::varCount() const -> VarIndex {
+    return mPimpl->conf.varCount();
+  }
+
+  void GroebnerInputIdealStream::idealBegin(size_t polyCount) {
+    MATHICGB_ASSERT(debugAssertValid());
+    MATHICGB_IF_DEBUG(mPimpl->checker.idealBegin(polyCount));
+    MATHICGB_ASSERT(mPimpl->poly.isZero());
+    MATHICGB_ASSERT(mPimpl->ideal.empty());
+
+    mPimpl->ideal.reserve(polyCount);
+
+    MATHICGB_ASSERT(debugAssertValid());
+  }
+
+  void GroebnerInputIdealStream::appendPolynomialBegin(size_t termCount) {
+    MATHICGB_ASSERT(debugAssertValid());
+    MATHICGB_IF_DEBUG(mPimpl->checker.appendPolynomialBegin(termCount));
+    MATHICGB_ASSERT(mPimpl->poly.isZero());
+
+    mPimpl->poly.reserve(termCount);
+
+    MATHICGB_ASSERT(debugAssertValid());
+  }
+
+  void GroebnerInputIdealStream::appendTermBegin() {
+    MATHICGB_ASSERT(debugAssertValid());
+    MATHICGB_IF_DEBUG(mPimpl->checker.appendTermBegin());
+
+    std::fill_n(mExponents, varCount(), 0);
+
+    MATHICGB_ASSERT(debugAssertValid());
+  }
+
+  void GroebnerInputIdealStream::appendTermDone(Coefficient coefficient) {
+    MATHICGB_ASSERT(debugAssertValid());
+    MATHICGB_IF_DEBUG(mPimpl->checker.appendTermDone(coefficient));
+
+    // @todo: do this directly into the polynomial instead of copying a second
+    // time.
+    mPimpl->ring.monomialSetExponents(mPimpl->monomial, mExponents);
+    mPimpl->poly.appendTerm(coefficient, mPimpl->monomial);
+
+    MATHICGB_ASSERT(debugAssertValid());
+  }
+
+  void GroebnerInputIdealStream::appendPolynomialDone() {
+    MATHICGB_ASSERT(debugAssertValid());
+    MATHICGB_IF_DEBUG(mPimpl->checker.appendPolynomialDone());
+
+    // todo: avoid copy here by the following changes
+    // todo: give Ideal a Poly&& insert
+    // todo: give Poly a Poly&& constructor
+    auto poly = make_unique<Poly>(std::move(mPimpl->poly));
+    if (!poly->termsAreInDescendingOrder())
+      poly->sortTermsDescending();
+    mPimpl->ideal.insert(std::move(poly));
+    mPimpl->poly.setToZero();
+
+    MATHICGB_ASSERT(debugAssertValid());
+  }
+
+  void GroebnerInputIdealStream::idealDone() {
+    MATHICGB_ASSERT(debugAssertValid());
+    MATHICGB_IF_DEBUG(mPimpl->checker.idealDone());
+  }
+
+  bool GroebnerInputIdealStream::debugAssertValid() const {
+    MATHICGB_ASSERT(this != 0);
+    MATHICGB_ASSERT(mExponents != 0);
+    MATHICGB_ASSERT(mPimpl != 0);
+    MATHICGB_ASSERT(!mPimpl->hasBeenDestroyed);
+    MATHICGB_ASSERT(!mPimpl->monomial.isNull());
+    MATHICGB_ASSERT(&mPimpl->ideal.ring() == &mPimpl->ring);
+    MATHICGB_ASSERT(&mPimpl->poly.ring() == &mPimpl->ring);
+    MATHICGB_ASSERT(mPimpl->ring.getNumVars() == mPimpl->conf.varCount());
+    MATHICGB_ASSERT(mPimpl->ring.charac() == mPimpl->conf.modulus());
+    return true;
+  }
+}
+
+// ** Implementation of class mgbi::PimplOf
+namespace mgbi {
+  class PimplOf {
+  public:
+    template<class T>
+    typename T::Pimpl& operator()(T& t) {
+      MATHICGB_ASSERT(t.mPimpl != 0);
+      return *t.mPimpl;
+    }
+  };
+}
+
+// ** Implementation of mgbi::IdealAdapter
+namespace mgbi {
+  struct IdealAdapter::Pimpl {
+    std::unique_ptr<Ideal> ideal;
+  };
+
+  IdealAdapter::IdealAdapter():
+    mPimpl(new Pimpl()) {
+  }
+
+  IdealAdapter::~IdealAdapter() {
+    MATHICGB_ASSERT(mPimpl != 0);
+    delete mPimpl;
+  }
+
+  auto IdealAdapter::varCount() const -> VarIndex {
+    MATHICGB_ASSERT(mPimpl->ideal.get() != 0);
+    return mPimpl->ideal->ring().getNumVars();
+  }
+
+  size_t IdealAdapter::polyCount() const {
+    MATHICGB_ASSERT(mPimpl->ideal.get() != 0);
+    return mPimpl->ideal->size();
+  }
+
+  size_t IdealAdapter::termCount(PolyIndex poly) const {
+    MATHICGB_ASSERT(mPimpl->ideal.get() != 0);
+    MATHICGB_ASSERT(poly < mPimpl->ideal->size());
+    return mPimpl->ideal->getPoly(poly)->nTerms();
+  }
+
+  auto IdealAdapter::term(
+    PolyIndex poly,
+    TermIndex term
+  ) const -> std::pair<Coefficient, const Exponent*> {
+    MATHICGB_ASSERT(mPimpl->ideal.get() != 0);
+    MATHICGB_ASSERT(poly < mPimpl->ideal->size());
+    const auto& p = *mPimpl->ideal->getPoly(poly);
+
+    MATHICGB_ASSERT(term < p.nTerms());
+    return std::make_pair(
+      p.coefficientAt(term),
+      p.monomialAt(term).unsafeGetRepresentation() + 1
+    );
+  }
+}
+
+// ** Implementation of function mgbi::internalComputeGreobnerBasis
+namespace mgbi {
+  void internalComputeGroebnerBasis(
+    GroebnerInputIdealStream& inputWhichWillBeCleared,
+    IdealAdapter& output
+  ) {
+    /// @todo: make a scheme where the output Groebner basis is freed
+    /// polynomial-by-polynomial as data is transferred to out. Also
+    /// make it so that ideal is not copied.
+
+    auto& ideal = PimplOf()(inputWhichWillBeCleared).ideal;
+    const auto& ring = ideal.ring();
+    const auto varCount = ring.getNumVars();
+
+    // Compute Groebner basis
+    const auto reducer = Reducer::makeReducer(Reducer::Reducer_F4_New, ring);
+    BuchbergerAlg alg(ideal, 4, *reducer, 2, true, 0);
+    alg.setSPairGroupSize(10000);
+    alg.setReducerMemoryQuantum(100 * 1024);
+    alg.setUseAutoTopReduction(true);
+    alg.setUseAutoTailReduction(false);
+    alg.computeGrobnerBasis();
+
+    PimplOf()(output).ideal = alg.basis().toIdealAndRetireAll();
+  }
+}
diff --git a/src/mathicgb.h b/src/mathicgb.h
new file mode 100644
index 0000000..e12f9ec
--- /dev/null
+++ b/src/mathicgb.h
@@ -0,0 +1,523 @@
+#ifndef MATHICGB_MATHICGB_GUARD
+#define MATHICGB_MATHICGB_GUARD
+
+#include <ostream>
+
+// The main function in this file is computeGroebnerBasis. See the comment
+// preceding that function for an example of how to use this library
+// interface.
+
+/// Code in this namespace is not part of the public interface of MathicGB.
+/// If you take a dependency on code in this namespace, except indirectly
+/// through using code in the mgb namespace, then you are not using the
+/// library interface as intended.
+namespace mgbi { // Not part of the public interface of MathicGB
+  class PimplOf;
+}
+
+/// The classes and functions in this namespace make up the public interface
+/// of MathicGB. You should not have to update your code when upgrading to a
+/// newer minor revision of MathicGB if you only use the public interface.
+namespace mgb { // Part of the public interface of MathicGB
+  /// Use this class to describe a configuration of a Groebner basis algorithm
+  /// that you want to run.
+  ///
+  /// @todo: expose more of the available functionality.
+  class GroebnerConfiguration {
+  public:
+    typedef unsigned int Coefficient;
+    typedef size_t VarIndex;
+    typedef int Exponent;
+
+    GroebnerConfiguration(Coefficient modulus, VarIndex varCount);
+    GroebnerConfiguration(const GroebnerConfiguration& conf);
+    ~GroebnerConfiguration();
+
+    Coefficient modulus() const;
+    VarIndex varCount() const;
+
+  private:
+    struct Pimpl;
+    Pimpl* mPimpl;
+  };
+
+  /// After making a configuration, use this class to communicate the input
+  /// basis you want want to run a Groebner basis algorithm on.
+  class GroebnerInputIdealStream {
+  public:
+    GroebnerInputIdealStream(const GroebnerConfiguration& conf);
+    ~GroebnerInputIdealStream();
+
+    typedef GroebnerConfiguration::Coefficient Coefficient;
+    typedef GroebnerConfiguration::VarIndex VarIndex;
+    typedef GroebnerConfiguration::Exponent Exponent;
+
+    const GroebnerConfiguration& configuration();
+    Coefficient modulus() const;
+    VarIndex varCount() const;
+
+    void idealBegin(size_t polyCount);
+    void appendPolynomialBegin(size_t termCount);
+    void appendTermBegin();
+
+    /// The sequence of indices appended to a term must be in strictly
+    /// ascending order.
+    void appendExponent(VarIndex index, Exponent exponent);
+    void appendTermDone(Coefficient coefficient);
+    void appendPolynomialDone();
+    void idealDone();
+
+  private:
+    bool debugAssertValid() const;
+    Exponent* const mExponents;
+
+    struct Pimpl;
+    friend class mgbi::PimplOf;
+    Pimpl* const mPimpl;
+  };
+
+  /// After making a configuration and an ideal, use this function to compute
+  /// a Groebner basis. The output basis is constructed on output, which must
+  /// resemble GroebnerInputIdealStream by having the following functions.
+  ///
+  ///   - modulus() const;
+  ///   - varCount() const;
+  ///   - idealBegin(size_t polyCount);
+  ///   - void appendPolynomialBegin(size_t termCount);
+  ///   - void appendTermBegin();
+  ///   - void appendExponent(VarIndex index, Exponent exponent);
+  ///   - void appendTermDone(Coefficient coefficient);
+  ///   - void appendPolynomialDone();
+  ///   - void idealDone();
+  ///
+  /// ** Example
+  ///
+  /// This example uses a default configuration, constructs an ideal and then
+  /// outputs the Groebner basis to a NullIdealStream which does not do
+  /// anything with the output. However, we wrap the NullIdealStream in a
+  /// IdealStreamLog, which prints out the method calls done on the stream
+  /// to std::cerr. We also wrap the GroebnerInputIdealStream in a
+  /// IdealStreamChecker which checks that we are correctly following the
+  /// protocol of GroebnerInputIdealStream - that's only recommended when
+  /// debugging as it is slow.
+  ///
+  /// GroebnerConfiguration configuration(101, 4); // mod 101, 4 variables
+  /// GroebnerInputIdealStream input(configuration);
+  /// IdealStreamChecker<GroebnerInputIdealStream> checked(input);
+  /// checked.idealBegin(2); // describe ideal with 2 basis elements
+  ///   checked.appendPolynomial(2); // describe generator with 2 terms
+  ///     checked.appendTermBegin();
+  ///       checked.appendExponent(0, 40); // x0^40
+  ///     checked.appendTermDone(3); // 3 * x0^40
+  ///     checked.appendTermBegin();
+  ///       checked.appendExponent(1, 5); // x1^5
+  ///       checked.appendExponent(2, 7); // x2^7
+  ///     checked.appendTermDone(11); // 11 * x1^5 * x2^7
+  ///   checked.appendPolynomialDone(); // 3 * x0^40 + 11 * x1^5 * x2^7
+  ///   checked.appendPolynomialBegin(1);
+  ///     checked.appendTermBegin();
+  ///     checked.appendTermDone(13); // 13
+  ///   checked.appendPolynomialDone(); // 13
+  /// checked.idealDone(); // the generators are 3*x0^40 + 11*x1^5*x2^7 and 13
+  /// NullIdealStream nullStream;
+  /// IdealStreamLog<NullIdealStream> logStream(nullStream, std::cerr);
+  /// computeGroebnerBasis(input, logStream);
+  ///
+  /// The ideal constructed on the passed-in GroebnerInputIdealStream will
+  /// be cleared. If you need it again for something else, you will have
+  /// to re-construct it.
+  template<class OutputStream>
+  void computeGroebnerBasis(
+    GroebnerInputIdealStream& inputWhichWillBeCleared,
+    OutputStream& output
+  );
+
+  /// Passes on all method calls to an inner ideal stream while printing out
+  /// what methods get called to an std::ostream.
+  template<class Stream = NullIdealStream>
+  class IdealStreamLog {
+  public:
+    typedef GroebnerConfiguration::Coefficient Coefficient;
+    typedef GroebnerConfiguration::VarIndex VarIndex;
+    typedef GroebnerConfiguration::Exponent Exponent;
+
+    /// All calls are written to log and then passed on to stream.
+    IdealStreamLog(std::ostream& log, Stream& stream);
+
+    /// All calls are written to log.
+    IdealStreamLog(std::ostream& log, Coefficient modulus, VarIndex varCount);
+
+    ~IdealStreamLog();
+
+    Coefficient modulus() const;
+    VarIndex varCount() const;
+
+    void idealBegin(size_t polyCount);
+    void appendPolynomialBegin(size_t termCount);
+    void appendTermBegin();
+    void appendExponent(VarIndex index, Exponent exponent);
+    void appendTermDone(Coefficient coefficient);
+    void appendPolynomialDone();
+    void idealDone();
+
+  private:
+    const Coefficient mModulus;
+    const VarIndex mVarCount;
+    Stream* const mStream;
+    std::ostream& mLog;
+  };
+
+  /// An ideal stream that simply ignores all the method calls it receives.
+  /// Can be handy in combination with IdealStreamLog or as a temporary
+  /// stand-in if you have not yet written your own ideal stream.
+  class NullIdealStream {
+  public:
+    typedef GroebnerConfiguration::Coefficient Coefficient;
+    typedef GroebnerConfiguration::VarIndex VarIndex;
+    typedef GroebnerConfiguration::Exponent Exponent;
+
+    NullIdealStream(Coefficient modulus, VarIndex varCount);
+
+    Coefficient modulus() const {return mModulus;}
+    VarIndex varCount() const {return mVarCount;}
+
+    void idealBegin(size_t polyCount) {}
+    void appendPolynomialBegin(size_t termCount) {}
+    void appendTermBegin() {}
+    void appendExponent(VarIndex index, Exponent exponent) {}
+    void appendTermDone(Coefficient coefficient) {}
+    void appendPolynomialDone() {}
+    void idealDone() {}
+
+  private:
+    const Coefficient mModulus;
+    const VarIndex mVarCount;
+  };
+
+  template<class OutputStream>
+  void streamSimpleIdeal(OutputStream& output);
+}
+
+namespace mgbi { // Not part of the public interface of MathicGB
+  using namespace mgb;
+
+  /// Class that checks to see if the protocol for an ideal stream is followed.
+  /// A method will return false if that is not the case. This class is not
+  /// itself an ideal stream - it is intended to be used with ideal streams for
+  /// debugging.
+  class StreamStateChecker {
+  public:
+    typedef GroebnerConfiguration::Coefficient Coefficient;
+    typedef GroebnerConfiguration::VarIndex VarIndex;
+    typedef GroebnerConfiguration::Exponent Exponent;
+
+    StreamStateChecker(const Coefficient modulus, const VarIndex varCount);
+    ~StreamStateChecker();
+
+    void idealBegin(size_t polyCount);
+    void appendPolynomialBegin(size_t termCount);
+    void appendTermBegin();
+    void appendExponent(VarIndex index, Exponent exponent);
+    void appendTermDone(Coefficient coefficient);
+    void appendPolynomialDone();
+    void idealDone();
+
+    bool hasIdeal() const;
+
+  private:
+    struct Pimpl;
+    Pimpl* const mPimpl;
+  };
+}
+
+namespace mgb { // Part of the public interface of MathicGB
+
+  /// Use this class to check that you are following the correct protocol
+  /// for calling methods on an ideal stream. This has significant overhead
+  /// so it is not recommended for production use. If you have built the
+  /// MathicGB library in debug mode then this is already automatically
+  /// used for GroebnerInputIdealStream.
+  template<class Stream>
+  class IdealStreamChecker {
+  public:
+    typedef GroebnerConfiguration::Coefficient Coefficient;
+    typedef GroebnerConfiguration::VarIndex VarIndex;
+    typedef GroebnerConfiguration::Exponent Exponent;
+
+    IdealStreamChecker(Stream& stream);
+
+    Coefficient modulus() const;
+    VarIndex varCount() const;
+
+    void idealBegin(size_t polyCount);
+    void appendPolynomialBegin(size_t termCount);
+    void appendTermBegin();
+    void appendExponent(VarIndex index, Exponent exponent);
+    void appendTermDone(Coefficient coefficient);
+    void appendPolynomialDone();
+    void idealDone();
+
+  private:
+    Stream& mStream;
+    ::mgbi::StreamStateChecker mChecker;
+  };
+}
+
+// ********************************************************
+// Nothing below this line is part of the public interface of MathicGB.
+
+#ifdef MATHICGB_DEBUG
+#include <cassert>
+#endif
+
+namespace mgb {
+  // ** Functions
+
+  // This method is made inline to avoid the overhead from calling a function
+  // for every exponent. This is also why mExponents is not inside the pimpl -
+  // otherwise we couldn't access it fro here. That then explains why mExponents
+  // is a raw pointer instead of a std::vector - the compiler for the caller
+  // and the library must agree on the memory layout of the object and that is
+  // less likely to introduce problems for a raw pointer than for a
+  // std::vector. In particular, doing it this way allows the library and
+  // the caller to use different implementations of the STL.
+  inline void GroebnerInputIdealStream::appendExponent(
+    const VarIndex index,
+    Exponent exponent
+  ) {
+#ifdef MATHICGB_DEBUG
+    assert(index < varCount());
+#endif
+    mExponents[index] = exponent;
+  }
+
+
+  // ** Implementation of the class IdealStreamLog
+  // This class has to be inline as it is a template.
+
+  template<class Stream> 
+  IdealStreamLog<Stream>::IdealStreamLog(std::ostream& log, Stream& stream):
+    mModulus(stream.modulus()),
+    mVarCount(stream.varCount()),
+    mStream(&stream),
+    mLog(log)
+  {
+    mLog << "IdealStreamLog s(stream, log); // modulus=" << mModulus
+      << ", varCount=" << mVarCount << '\n';
+  }
+
+  template<class Stream> 
+  IdealStreamLog<Stream>::IdealStreamLog(
+    std::ostream& log,
+    Coefficient modulus,
+    VarIndex varCount
+  ):
+    mModulus(modulus),
+    mVarCount(varCount),
+    mStream(0),
+    mLog(log)
+  {
+    mLog << "IdealStreamLog s(stream, " << mModulus << ", " << mVarCount << ");\n";
+  }
+
+  template<class Stream> 
+  IdealStreamLog<Stream>::~IdealStreamLog() {
+    mLog << "// s.~IdealStreamLog();\n";
+  }
+
+  template<class Stream> 
+  typename IdealStreamLog<Stream>::Coefficient
+  IdealStreamLog<Stream>::modulus() const {
+    return mModulus;
+  }
+
+  template<class Stream> 
+  typename IdealStreamLog<Stream>::VarIndex
+  IdealStreamLog<Stream>::varCount() const {
+    return mVarCount;
+  }
+
+  template<class Stream> 
+  void IdealStreamLog<Stream>::idealBegin(size_t polyCount) {
+    mLog << "s.idealBegin(" << polyCount << "); // polyCount\n";
+    if (mStream != 0)
+      mStream->idealBegin(polyCount);
+  }
+
+  template<class Stream> 
+  void IdealStreamLog<Stream>::appendPolynomialBegin(size_t termCount) {
+    mLog << "s.appendPolynomialBegin(" << termCount << ");\n";
+    if (mStream != 0)
+      mStream->appendPolynomialBegin(termCount);
+  }
+
+  template<class Stream> 
+  void IdealStreamLog<Stream>::appendTermBegin() {
+    mLog << "s.appendTermBegin();\n";
+    if (mStream != 0)
+      mStream->appendTermBegin();
+  }
+
+  template<class Stream> 
+  void IdealStreamLog<Stream>::appendExponent(VarIndex index, Exponent exponent) {
+    mLog << "s.appendExponent(" << index << ", " << exponent <<
+      "); // index, exponent\n";
+    if (mStream != 0)
+      mStream->appendExponent(index, exponent);
+  }
+
+  template<class Stream> 
+  void IdealStreamLog<Stream>::appendTermDone(Coefficient coefficient) {
+    mLog << "s.appendTermDone(" << coefficient << "); // coefficient\n";
+    if (mStream != 0)
+      mStream->appendTermDone(coefficient);
+  }
+
+  template<class Stream> 
+  void IdealStreamLog<Stream>::appendPolynomialDone() {
+    mLog << "s.appendPolynomialDone();\n";
+    if (mStream != 0)
+      mStream->appendPolynomialDone();
+  }
+
+  template<class Stream> 
+  void IdealStreamLog<Stream>::idealDone() {
+    mLog << "s.idealDone();\n";
+    if (mStream != 0)
+      mStream->idealDone();
+  }
+
+
+  // ** Implementation of the class IdealStreamChecker
+  // This class has to be inline as it is a template.
+
+  template<class Stream>
+  IdealStreamChecker<Stream>::IdealStreamChecker(Stream& stream):
+    mStream(stream),
+    mChecker(stream.modulus(), stream.varCount())
+  {}
+
+  template<class Stream>
+  typename IdealStreamChecker<Stream>::Coefficient
+  IdealStreamChecker<Stream>::modulus() const {
+    return mStream.modulus();
+  }
+
+  template<class Stream>
+  typename IdealStreamChecker<Stream>::VarIndex
+  IdealStreamChecker<Stream>::varCount() const {
+    return mStream.varCount();
+  }
+
+  template<class Stream>
+  void IdealStreamChecker<Stream>::idealBegin(size_t polyCount) {
+    mChecker.idealBegin(polyCount);
+    mStream.idealBegin(polyCount);
+  }
+
+  template<class Stream>
+  void IdealStreamChecker<Stream>::appendPolynomialBegin(size_t termCount) {
+    mChecker.appendPolynomialBegin(termCount);
+    mStream.appendPolynomialBegin(termCount);
+  }
+
+  template<class Stream>
+  void IdealStreamChecker<Stream>::appendTermBegin() {
+    mChecker.appendTermBegin();
+    mStream.appendTermBegin();
+  }
+
+  template<class Stream>
+  void IdealStreamChecker<Stream>::appendExponent(VarIndex index, Exponent exponent) {
+    mChecker.appendExponent(index, exponent);
+    mStream.appendExponent(index, exponent);
+  }
+
+  template<class Stream>
+  void IdealStreamChecker<Stream>::appendTermDone(Coefficient coefficient) {
+    mChecker.appendTermDone(coefficient);
+    mStream.appendTermDone(coefficient);
+  }
+
+  template<class Stream>
+  void IdealStreamChecker<Stream>::appendPolynomialDone() {
+    mChecker.appendPolynomialDone();
+    mStream.appendPolynomialDone();
+  }
+
+  template<class Stream>
+  void IdealStreamChecker<Stream>::idealDone() {
+    mChecker.idealDone();
+    mStream.idealDone();
+  }
+
+  // ** Implementation of the class NullIdealStream
+  // This class has to be inline as it is a template.
+  inline NullIdealStream::NullIdealStream(
+    Coefficient modulus,
+    VarIndex varCount
+  ):
+    mModulus(modulus), mVarCount(varCount) {}
+}
+
+namespace mgbi {
+  /// Used to read an internal MathicGB ideal without exposing the type of
+  /// the ideal.
+  class IdealAdapter {
+  public:
+    typedef GroebnerConfiguration::Coefficient Coefficient;
+    typedef GroebnerConfiguration::VarIndex VarIndex;
+    typedef GroebnerConfiguration::Exponent Exponent;
+    typedef size_t PolyIndex;
+    typedef size_t TermIndex;
+
+    IdealAdapter();
+    ~IdealAdapter();
+
+    VarIndex varCount() const;
+    size_t polyCount() const;
+    size_t termCount(PolyIndex poly) const;
+    std::pair<Coefficient, const Exponent*> term
+      (PolyIndex poly, TermIndex term) const;
+
+  private:
+    friend class ::mgbi::PimplOf;
+    struct Pimpl;
+    Pimpl* mPimpl;
+  };
+
+  void internalComputeGroebnerBasis(
+    GroebnerInputIdealStream& inputWhichWillBeCleared,
+    IdealAdapter& output
+  );
+}
+
+namespace mgb {
+  template<class OutputStream>
+  void computeGroebnerBasis(
+    GroebnerInputIdealStream& inputWhichWillBeCleared,
+    OutputStream& output
+  ) {
+    mgbi::IdealAdapter ideal;
+    mgbi::internalComputeGroebnerBasis(inputWhichWillBeCleared, ideal);
+
+    const auto varCount = ideal.varCount();
+    const auto polyCount = ideal.polyCount();
+    output.idealBegin(polyCount);
+    for (size_t polyIndex = 0; polyIndex < polyCount; ++polyIndex) {
+      const auto termCount = ideal.termCount(polyIndex);
+      output.appendPolynomialBegin(termCount);
+      for (size_t termIndex = 0; termIndex < termCount; ++termIndex) {
+        output.appendTermBegin();
+        const auto term = ideal.term(polyIndex, termIndex);
+        for (size_t var = 0; var < varCount; ++var)
+          output.appendExponent(var, term.second[var]);
+        output.appendTermDone(term.first);
+      }
+      output.appendPolynomialDone();
+    }
+    output.idealDone();
+  }
+}
+
+#endif
diff --git a/src/mathicgb/PolyBasis.cpp b/src/mathicgb/PolyBasis.cpp
index 761cf6c..73cd31b 100755
--- a/src/mathicgb/PolyBasis.cpp
+++ b/src/mathicgb/PolyBasis.cpp
@@ -87,6 +87,24 @@ void PolyBasis::insert(std::unique_ptr<Poly> poly) {
   MATHICGB_ASSERT(mEntries.back().poly != 0);
 }
 
+std::unique_ptr<Poly> PolyBasis::retire(size_t index) {
+  MATHICGB_ASSERT(index < size());
+  MATHICGB_ASSERT(!retired(index));
+  mDivisorLookup->remove(leadMonomial(index));
+  std::unique_ptr<Poly> poly(mEntries[index].poly);
+  mEntries[index].poly = 0;
+  mEntries[index].retired = true;
+  return poly;
+}
+
+std::unique_ptr<Ideal> PolyBasis::toIdealAndRetireAll() {
+  auto ideal = make_unique<Ideal>(ring());
+  for (size_t i = 0; i < size(); ++i)
+    if (!retired(i))
+      ideal->insert(retire(i));
+  return ideal;
+}
+
 size_t PolyBasis::divisor(const_monomial mon) const {
   size_t index = divisorLookup().divisor(mon);
   MATHICGB_ASSERT((index == static_cast<size_t>(-1)) ==
diff --git a/src/mathicgb/PolyBasis.hpp b/src/mathicgb/PolyBasis.hpp
index 4b38816..4cdebe8 100755
--- a/src/mathicgb/PolyBasis.hpp
+++ b/src/mathicgb/PolyBasis.hpp
@@ -114,16 +114,12 @@ public:
 
   // Retires the basis element at index, which frees the memory associated
   // to it, including the basis element polynomial, and marks it as retired. 
-  // todo: implement
-  std::unique_ptr<Poly> retire(size_t index) {
-    MATHICGB_ASSERT(index < size());
-    MATHICGB_ASSERT(!retired(index));
-    mDivisorLookup->remove(leadMonomial(index));
-    std::unique_ptr<Poly> poly(mEntries[index].poly);
-    mEntries[index].poly = 0;
-    mEntries[index].retired = true;
-    return poly;
-  }
+  std::unique_ptr<Poly> retire(size_t index);
+
+  /// Returns an ideal containing all non-retired basis elements and
+  /// retires all those basis elements. The point of the simultaneous
+  /// retirement is that this way no polynomials need be copied.
+  std::unique_ptr<Ideal> toIdealAndRetireAll();
 
   // Returns true of the basis element at index has been retired.
   bool retired(size_t index) const {
diff --git a/src/mathicgb/PolyRing.cpp b/src/mathicgb/PolyRing.cpp
index ae5c45e..8a38518 100755
--- a/src/mathicgb/PolyRing.cpp
+++ b/src/mathicgb/PolyRing.cpp
@@ -26,11 +26,11 @@ PolyRing::PolyRing(coefficient p0,
     mMaxMonomialSize(nvars + nweights + 2),
     mMaxMonomialByteSize(mMaxMonomialSize * sizeof(exponent)),
     mMonomialPool(mMaxMonomialByteSize),
-    mTotalDegreeGradedOnly(false)
+    mTotalDegreeGradedOnly(nweights == 1)
 {
-  // set weights to the default value 1
+  // set weights to the default value -1
   mWeights.resize(mNumVars * mNumWeights);
-  std::fill(mWeights.begin(), mWeights.end(), 1);
+  std::fill(mWeights.begin(), mWeights.end(), static_cast<exponent>(-1));
 
   resetCoefficientStats();
   srand(0);
diff --git a/src/mathicgb/PolyRing.hpp b/src/mathicgb/PolyRing.hpp
index c3b9638..6d65d0f 100755
--- a/src/mathicgb/PolyRing.hpp
+++ b/src/mathicgb/PolyRing.hpp
@@ -348,6 +348,7 @@ public:
   }
 
   void monomialSetExponents(Monomial m, exponent* exponents) const {
+    *m = 0;
     std::memcpy(
       m.unsafeGetRepresentation() + 1,
       exponents,
@@ -701,6 +702,9 @@ inline bool PolyRing::monomialIsProductOfHintTrue(
 
   uint64 orOfXor = 0;
   for (size_t i = mNumVars / 2; i != static_cast<size_t>(-1); --i) {
+    MATHICGB_ASSERT(a[i*2] >= 0);
+    MATHICGB_ASSERT(i == mNumVars / 2 || a[i*2+1] >= 0);
+
     uint64 A, B, AB;
     // We have to use std::memcpy here because just casting to a int64 breaks
     // the strict aliasing rule which implies undefined behavior. Both MSVC and
@@ -781,7 +785,7 @@ inline void PolyRing::setWeightsOnly(Monomial& a1) const
 {
   exponent *a = a1.unsafeGetRepresentation();
   a++;
-  auto wts = &mWeights[0];
+  auto wts = mWeights.data();
   for (size_t i = 0; i < mNumWeights; ++i)
     {
       exponent result = 0;
@@ -1034,7 +1038,7 @@ inline bool PolyRing::monomialHasStrictlyLargerExponent(
 inline void PolyRing::setWeightsOnly(monomial a) const
 {
   a++;
-  auto wts = &mWeights[0];
+  auto wts = mWeights.data();
   for (size i = 0; i < mNumWeights; ++i) {
     exponent result = 0;
     for (size_t j=0; j<mNumVars; ++j)
diff --git a/src/test/mathicgb.cpp b/src/test/mathicgb.cpp
new file mode 100644
index 0000000..5178eca
--- /dev/null
+++ b/src/test/mathicgb.cpp
@@ -0,0 +1,249 @@
+#include "mathicgb/stdinc.h"
+
+#include <gtest/gtest.h>
+#include "mathicgb.h"
+
+namespace {
+  template<class Stream>
+  void makeBasis(Stream& s) {
+    s.idealBegin(2);
+      s.appendPolynomialBegin(2); // x^2 - y
+        s.appendTermBegin();
+          s.appendExponent(0,2);
+        s.appendTermDone(1);
+        s.appendTermBegin();
+          s.appendExponent(1,1);
+        s.appendTermDone(s.modulus() - 1);
+      s.appendPolynomialDone();
+      s.appendPolynomialBegin(2); // x^3-z
+        s.appendTermBegin();
+          s.appendExponent(0,3);
+        s.appendTermDone(1);
+        s.appendTermBegin();
+          s.appendExponent(2,1);
+        s.appendTermDone(s.modulus() - 1);
+      s.appendPolynomialDone();
+    s.idealDone();
+  }
+
+  template<class Stream>
+  void makeGroebnerBasis(Stream& s) {
+    s.idealBegin(3);
+      s.appendPolynomialBegin(2); // x^2 - y
+        s.appendTermBegin();
+          s.appendExponent(0, 2);
+          s.appendExponent(1, 0);
+          s.appendExponent(2, 0);
+        s.appendTermDone(1);
+        s.appendTermBegin();
+          s.appendExponent(0, 0);
+          s.appendExponent(1, 1);
+          s.appendExponent(2, 0);
+        s.appendTermDone(s.modulus() - 1);
+      s.appendPolynomialDone();
+      s.appendPolynomialBegin(2); // xy - z
+        s.appendTermBegin();
+          s.appendExponent(0, 1);
+          s.appendExponent(1, 1);
+          s.appendExponent(2, 0);
+        s.appendTermDone(1);
+        s.appendTermBegin();
+          s.appendExponent(0, 0);
+          s.appendExponent(1, 0);
+          s.appendExponent(2, 1);
+        s.appendTermDone(s.modulus() - 1);
+      s.appendPolynomialDone();
+      s.appendPolynomialBegin(2); // y^2 - xy
+        s.appendTermBegin();
+          s.appendExponent(0, 0);
+          s.appendExponent(1, 2);
+          s.appendExponent(2, 0);
+        s.appendTermDone(1);
+        s.appendTermBegin();
+          s.appendExponent(0, 1);
+          s.appendExponent(1, 0);
+          s.appendExponent(2, 1);
+        s.appendTermDone(s.modulus() - 1);
+      s.appendPolynomialDone();
+    s.idealDone();
+  }
+}
+
+TEST(MathicGBLib, NullIdealStream) {
+  {
+    mgb::NullIdealStream stream(2, 3);
+    ASSERT_EQ(2, stream.modulus());
+    ASSERT_EQ(3, stream.varCount());
+    makeBasis(stream);
+  }
+
+  {
+    mgb::NullIdealStream stream(101, 0);
+    ASSERT_EQ(101, stream.modulus());
+    ASSERT_EQ(0, stream.varCount());
+  }
+}
+
+TEST(MathicGBLib, IdealStreamLog) {
+  {
+    const char* const idealStr = 
+      "s.idealBegin(2); // polyCount\n"
+      "s.appendPolynomialBegin(2);\n"
+      "s.appendTermBegin();\n"
+      "s.appendExponent(0, 2); // index, exponent\n"
+      "s.appendTermDone(1); // coefficient\n"
+      "s.appendTermBegin();\n"
+      "s.appendExponent(1, 1); // index, exponent\n"
+      "s.appendTermDone(6); // coefficient\n"
+      "s.appendPolynomialDone();\n"
+      "s.appendPolynomialBegin(2);\n"
+      "s.appendTermBegin();\n"
+      "s.appendExponent(0, 3); // index, exponent\n"
+      "s.appendTermDone(1); // coefficient\n"
+      "s.appendTermBegin();\n"
+      "s.appendExponent(2, 1); // index, exponent\n"
+      "s.appendTermDone(6); // coefficient\n"
+      "s.appendPolynomialDone();\n"
+      "s.idealDone();\n";
+
+    std::ostringstream out1;
+    mgb::IdealStreamLog<> stream1(out1, 7, 3);
+
+    mgb::IdealStreamChecker<decltype(stream1)> checker(stream1);
+
+    std::ostringstream out2;
+    mgb::IdealStreamLog<decltype(checker)> stream2(out2, checker);
+
+    std::ostringstream out3;
+    mgb::IdealStreamLog<decltype(stream2)> stream3(out3, stream2);
+
+    ASSERT_EQ(7, stream1.modulus());
+    ASSERT_EQ(3, stream1.varCount());
+    ASSERT_EQ(7, checker.modulus());
+    ASSERT_EQ(3, checker.varCount());
+    ASSERT_EQ(7, stream2.modulus());
+    ASSERT_EQ(3, stream2.varCount());
+    ASSERT_EQ(7, stream3.modulus());
+    ASSERT_EQ(3, stream3.varCount());
+
+    makeBasis(stream3);
+    const auto str1 = std::string(
+      "IdealStreamLog s(stream, 7, 3);\n"
+    ) + idealStr;
+    ASSERT_EQ(str1, out1.str()) << "Displayed expected:\n" << out1.str();
+
+    const auto str2 = std::string(
+      "IdealStreamLog s(stream, log); // modulus=7, varCount=3\n"
+    ) + idealStr;
+    ASSERT_EQ(str2, out2.str()) << "Displayed expected:\n" << out2.str();
+    ASSERT_EQ(str2, out3.str()) << "Displayed expected:\n" << out3.str();
+  }
+
+  // The ideal <> in no variables
+  {
+    std::ostringstream out;
+    mgb::IdealStreamLog<> stream(out, 101, 0);
+    ASSERT_EQ(101, stream.modulus());
+    ASSERT_EQ(0, stream.varCount());
+    stream.idealBegin(0);
+    stream.idealDone();
+  }
+
+  // The ideal <1, 0> in no variables
+  {
+    std::ostringstream out;
+    mgb::IdealStreamLog<> stream(out, 101, 0);
+    ASSERT_EQ(101, stream.modulus());
+    ASSERT_EQ(0, stream.varCount());
+    stream.idealBegin(2);
+      stream.appendPolynomialBegin(0); // 1
+        stream.appendTermBegin();
+        stream.appendTermDone(1);
+      stream.appendPolynomialDone();
+      stream.appendPolynomialBegin(0); // 0
+      stream.appendPolynomialDone();
+    stream.idealDone();
+  }
+}
+
+TEST(MathicGBLib, ZeroIdealGB) {
+  mgb::GroebnerConfiguration configuration(2, 0);
+  mgb::GroebnerInputIdealStream input(configuration);
+  std::ostringstream out;
+  mgb::IdealStreamLog<> logStream(out, 2, 0);
+
+  input.idealBegin(0);
+  input.idealDone();
+  mgb::computeGroebnerBasis(input, logStream);
+
+  const auto msg =
+    "IdealStreamLog s(stream, 2, 0);\n"
+    "s.idealBegin(0); // polyCount\n"
+    "s.idealDone();\n";
+  EXPECT_EQ(msg, out.str());
+}
+
+TEST(MathicGBLib, OneIdealGB) {
+  mgb::GroebnerConfiguration configuration(2, 0);
+  mgb::GroebnerInputIdealStream input(configuration);
+  std::ostringstream out;
+  mgb::IdealStreamLog<> logStream(out, 2, 0);
+
+  input.idealBegin(1);
+  input.appendPolynomialBegin(1);
+  input.appendTermBegin();
+  input.appendTermDone(1);
+  input.appendPolynomialDone();
+  input.idealDone();
+  mgb::computeGroebnerBasis(input, logStream);
+
+  const auto msg =
+     "IdealStreamLog s(stream, 2, 0);\n"
+     "s.idealBegin(1); // polyCount\n"
+     "s.appendPolynomialBegin(1);\n"
+     "s.appendTermBegin();\n"
+     "s.appendTermDone(1); // coefficient\n"
+     "s.appendPolynomialDone();\n"
+     "s.idealDone();\n";
+  EXPECT_EQ(msg, out.str());
+}
+
+TEST(MathicGBLib, EasyGB) {
+  mgb::GroebnerConfiguration configuration(101, 3);
+  mgb::GroebnerInputIdealStream input(configuration);
+  std::ostringstream computedStr;
+  mgb::IdealStreamLog<> computed(computedStr, 101, 3);
+  mgb::IdealStreamChecker<decltype(computed)> checked(computed);
+
+  makeBasis(input);
+  mgb::computeGroebnerBasis(input, checked);
+
+  std::ostringstream correctStr;
+  mgb::IdealStreamLog<> correct(correctStr, 101, 3);
+  mgb::IdealStreamChecker<decltype(correct)> correctChecked(correct);
+  makeGroebnerBasis(correctChecked);
+
+  EXPECT_EQ(correctStr.str(), computedStr.str())
+    << "\nDisplayed expected:\n" << correctStr.str()
+    << "\nDisplayed computed:\n" << computedStr.str();
+}
+
+TEST(MathicGBLib, EasyReGB) {
+  mgb::GroebnerConfiguration configuration(101, 3);
+  mgb::GroebnerInputIdealStream input(configuration);
+  std::ostringstream computedStr;
+  mgb::IdealStreamLog<> computed(computedStr, 101, 3);
+  mgb::IdealStreamChecker<decltype(computed)> checked(computed);
+
+  makeGroebnerBasis(input);
+  mgb::computeGroebnerBasis(input, checked);
+
+  std::ostringstream correctStr;
+  mgb::IdealStreamLog<> correct(correctStr, 101, 3);
+  mgb::IdealStreamChecker<decltype(correct)> correctChecked(correct);
+  makeGroebnerBasis(correctChecked);
+
+  EXPECT_EQ(correctStr.str(), computedStr.str())
+    << "\nDisplayed expected:\n" << correctStr.str()
+    << "\nDisplayed computed:\n" << computedStr.str();
+}

-- 
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