[mathicgb] 204/393: Added MonoMonoid along with tests. The purpose of MonoMonoid is to take over the monomial functionality from PolyRing.
Doug Torrance
dtorrance-guest at moszumanska.debian.org
Fri Apr 3 15:59:00 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 7849dfdc987adf97bb61772a2722c350fb4967e6
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date: Wed Mar 27 22:03:31 2013 +0100
Added MonoMonoid along with tests. The purpose of MonoMonoid is to take over the monomial functionality from PolyRing.
---
Makefile.am | 4 +-
src/mathicgb/MonoMonoid.hpp | 452 ++++++++++++++++++++++++++++++++++++++++++++
src/mathicgb/PrimeField.hpp | 290 ++++++++++++++--------------
src/test/MonoMonoid.cpp | 233 +++++++++++++++++++++++
4 files changed, 832 insertions(+), 147 deletions(-)
diff --git a/Makefile.am b/Makefile.am
index 79337a5..6bfe5e1 100755
--- a/Makefile.am
+++ b/Makefile.am
@@ -68,7 +68,7 @@ libmathicgb_la_SOURCES = src/mathicgb/BjarkeGeobucket2.cpp \
src/mathicgb/F4ProtoMatrix.cpp src/mathicgb/F4MatrixProject.hpp \
src/mathicgb/F4MatrixProjection.cpp src/mathicgb/ScopeExit.hpp \
src/mathicgb.cpp src/mathicgb.h src/mathicgb/mtbb.hpp \
- src/mathicgb/PrimeField.hpp
+ src/mathicgb/PrimeField.hpp src/mathicgb/MonoMonoid.hpp
# The headers that libmathicgb installs.
@@ -115,7 +115,7 @@ unittest_SOURCES=src/test/FreeModuleOrderTest.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/mathicgb.cpp src/test/PrimeField.cpp
+ src/test/mathicgb.cpp src/test/PrimeField.cpp src/test/MonoMonoid.cpp
else
diff --git a/src/mathicgb/MonoMonoid.hpp b/src/mathicgb/MonoMonoid.hpp
new file mode 100755
index 0000000..f284ca7
--- /dev/null
+++ b/src/mathicgb/MonoMonoid.hpp
@@ -0,0 +1,452 @@
+#ifndef MATHICGB_MONO_MONOID_GUARD
+#define MATHICGB_MONO_MONOID_GUARD
+
+#include <vector>
+#include <algorithm>
+#include <memtailor.h>
+
+/// Implements the monoid of (monic) monomials with integer
+/// non-negative exponents. T must be an unsigned integer type that is
+/// used to store each exponent of a monomial.
+///
+/// TODO: support grading and comparison.
+template<class E>
+class MonoMonoid {
+public:
+ // *** Types
+
+ // Integer index representing a variable. Indices start at 0 and go
+ // up to varCount() - 1 where varCount() is the number of variables.
+ typedef size_t VarIndex;
+
+ /// The type of each exponent of a monomial.
+ typedef E Exponent;
+
+ /// Iterator for the exponents in a monomial.
+ typedef const Exponent* const_iterator;
+
+ /// Represents a monomial and manages the memory underlying it. To
+ /// refer to a non-owned monomial or to refer to a Mono, use MonoRef
+ /// or ConstMonoRef. Do not use Mono& or Mono* if you do not have
+ /// to, since that implies a double indirection when accessing the
+ /// monomial.
+ class Mono;
+
+ /// A reference to a non-const monomial. Cannot be null, cannot be
+ /// reassigned to refer to a different monomial and does not connote
+ /// ownership - the same semantics as C++ references.
+ class MonoRef;
+
+ /// A reference to a monomial. As MonoRef, but you cannot change the
+ /// monomial through this reference. Prefer this class over the
+ /// other reference/pointer classes unless there is a reason not to.
+ class ConstMonoRef;
+
+ /// A pointer to a non-const monomial. Can be null and can be
+ /// reassigned to refer to a different monomial - the same semantics
+ /// as C++ pointers. Does not connote ownership.
+ class MonoPtr;
+
+ /// A pointer to a monomial. As MonoPtr, but you cannot change the
+ /// monomial through this pointer.
+ class ConstMonoPtr;
+
+ /// A pool of memory for monomials.
+ ///
+ /// @todo: This approach is a poor fit for variable-sized
+ /// monomials. So prefer other solutions where reasonable.
+ class MonoPool;
+
+ /// A vector of monomials. The interface is a subset of
+ /// std::vector. Monomials can be appended (push_back). Only the
+ /// last monomial can be mutated and monomials cannot be reordered
+ /// or removed. These restrictions should make it easier to support
+ /// variable-sized monomials in future. Change it if you need to
+ /// break these restrictions, but first try to find an alternative.
+ class MonoVector;
+
+
+ // *** Constructors and accessors
+
+ MonoMonoid(const VarIndex varCount): mVarCount(varCount) {}
+
+ bool operator==(const MonoMonoid& monoid) const {
+ return this == &monoid;
+ }
+
+ bool operator!=(const MonoMonoid& monoid) const {
+ return !(*this == monoid);
+ }
+
+ /// Returns the number of variables. This is also the number of
+ /// exponents in the exponent vector of a monomial.
+ VarIndex varCount() const {return mVarCount;}
+
+
+ // *** Monomial accessors and queries
+
+ /// Returns iterator to the first exponent.
+ const_iterator begin(ConstMonoRef mono) const {return mono.rawPtr();}
+
+ /// Returns iterator to one-past-the-end of the range of exponents.
+ const_iterator end(ConstMonoRef mono) const {
+ return mono.rawPtr() + entriesPerMono();
+ }
+
+ /// Returns the exponent of var in mono.
+ Exponent exponent(ConstMonoRef mono, const VarIndex var) const {
+ MATHICGB_ASSERT(var < varCount());
+ return mono.rawPtr()[var];
+ }
+
+ bool equal(ConstMonoRef a, ConstMonoRef b) const {
+ return std::equal(begin(a), end(a), begin(b));
+ }
+
+ /// Returns true if all the exponents of mono are zero. In other
+ /// words, returns true if mono is the identity for multiplication
+ /// of monomials.
+ bool isIdentity(ConstMonoRef mono) const {
+ return std::all_of(begin(mono), end(mono), [](Exponent e) {return e == 0;});
+ }
+
+
+ // *** Monomial mutating computations
+
+ void copy(ConstMonoRef from, MonoRef to) const {
+ std::copy_n(from.rawPtr(), entriesPerMono(), to.rawPtr());
+ }
+
+ void setExponent(const VarIndex var, const Exponent e, MonoRef mono) const {
+ MATHICGB_ASSERT(var < varCount());
+ mono.rawPtr()[var] = e;
+ }
+
+ void setIdentity(MonoRef mono) const {
+ std::fill_n(mono.rawPtr(), entriesPerMono(), static_cast<Exponent>(0));
+ MATHICGB_ASSERT(isIdentity(mono));
+ }
+
+
+ // *** Classes for holding and referring to monomials
+
+ class ConstMonoPtr {
+ public:
+ ConstMonoPtr(): mMono(0) {}
+ ConstMonoPtr(const ConstMonoPtr& mono): mMono(mono.rawPtr()) {}
+
+ ConstMonoPtr operator=(const ConstMonoPtr& mono) {
+ mMono = mono.mMono;
+ return *this;
+ }
+
+ ConstMonoRef operator*() const {return *this;}
+
+ bool isNull() const {return mMono == 0;}
+ void toNull() {mMono = 0;}
+
+ private:
+ friend class MonoVector;
+ friend class MonoPtr;
+ friend class ConstMonoRef;
+
+ const Exponent* rawPtr() const {return mMono;}
+ ConstMonoPtr(const Exponent* mono): mMono(mono) {}
+
+ const Exponent* mMono;
+ };
+
+ class MonoPtr {
+ public:
+ MonoPtr(): mMono(0) {}
+ MonoPtr(const MonoPtr& mono): mMono(mono.rawPtr()) {}
+
+ MonoPtr operator=(const MonoPtr& mono) {
+ mMono = mono.mMono;
+ return *this;
+ }
+
+ MonoRef operator*() const {return *this;}
+
+ bool isNull() const {return mMono == 0;}
+ void toNull() {mMono = 0;}
+
+ operator ConstMonoPtr() const {return ConstMonoPtr(mMono);}
+
+ private:
+ friend class Mono;
+ friend class MonoRef;
+ friend class MonoVector;
+ friend class MonoPool;
+
+ Exponent* rawPtr() const {return mMono;}
+ MonoPtr(Exponent* mono): mMono(mono) {}
+
+ Exponent* mMono;
+ };
+
+ class Mono {
+ public:
+ Mono(): mMono(), mPool(0) {}
+
+ Mono(Mono&& mono): mMono(mono.mMono), mPool(mono.mPool) {
+ mono.mMono.toNull();
+ mono.mPool = 0;
+ }
+
+ ~Mono() {toNull();}
+
+ void operator=(Mono&& mono) {
+ toNull();
+
+ mMono = mono.mMono;
+ mono.mMono.toNull();
+
+ mPool = mono.mPool;
+ mono.mPool = 0;
+ }
+
+ bool isNull() const {return mMono.isNull();}
+ void toNull() {mPool->free(*this);}
+
+ MonoPtr ptr() const {return mMono;}
+
+ operator MonoRef() const {
+ MATHICGB_ASSERT(!isNull());
+ return *mMono;
+ }
+
+ private:
+ Mono(const Mono&); // not available
+ void operator=(const Mono&); // not available
+ friend class MonoPool;
+
+ Mono(const MonoPtr mono, MonoPool& pool):
+ mMono(mono), mPool(&pool) {}
+
+ Exponent* rawPtr() const {return mMono.rawPtr();}
+
+ MonoPtr mMono;
+ MonoPool* mPool;
+ };
+
+ class MonoRef {
+ public:
+ MonoPtr ptr() const {return mMono;}
+
+ operator ConstMonoRef() const {return *static_cast<ConstMonoPtr>(mMono);}
+
+ private:
+ void operator=(const MonoRef&); // not available
+ friend class MonoMonoid;
+ friend class MonoPtr;
+
+ MonoRef(MonoPtr mono): mMono(mono) {}
+ Exponent* rawPtr() const {return mMono.rawPtr();}
+
+ const MonoPtr mMono;
+ };
+
+ class ConstMonoRef {
+ public:
+ ConstMonoRef(const Mono& mono): mMono(mono.ptr()) {
+ MATHICGB_ASSERT(!mono.isNull());
+ }
+
+ ConstMonoPtr ptr() const {return mMono;}
+
+ private:
+ void operator=(const MonoRef&); // not available
+ friend class MonoMonoid;
+ friend class ConstMonoPtr;
+
+ ConstMonoRef(ConstMonoPtr mono): mMono(mono) {}
+ const Exponent* rawPtr() const {return mMono.rawPtr();}
+
+ const ConstMonoPtr mMono;
+ };
+
+
+ // *** Classes that provide memory resources for monomials
+
+ class MonoPool {
+ public:
+ MonoPool(const MonoMonoid& monoid):
+ mMonoid(monoid),
+ mPool(sizeof(Exponent) * mMonoid.entriesPerMono()) {}
+
+ Mono alloc() {
+ const auto ptr = static_cast<Exponent*>(mPool.alloc());
+ Mono mono(ptr, *this);
+ monoid().setIdentity(mono);
+ return mono;
+ }
+
+ void free(Mono& mono) {free(std::move(mono));}
+ void free(Mono&& mono) {
+ if (mono.isNull())
+ return;
+ mPool.free(mono.rawPtr());
+ mono.mMono = 0;
+ mono.mPool = 0;
+ }
+
+ const MonoMonoid& monoid() const {return mMonoid;}
+
+ private:
+ const MonoMonoid& mMonoid;
+ memt::BufferPool mPool;
+ };
+
+ class MonoVector {
+ private:
+ typedef std::vector<Exponent> RawVector;
+
+ public:
+ /// Class for iterating through the monomials in a MonoVector.
+ ///
+ /// There is no operator->() since MonoRef does not have any
+ /// relevant methods to call. Implement it if you need it.
+ ///
+ /// There are no postfix increment operator as prefix is
+ /// better. Add it if you y need it (you probably do not).
+ ///
+ /// We could make this a random access iterator, but that would
+ /// make it tricky to support variable-sized exponent vectors
+ /// (e.g. sparse) in future and so far we have not needed random
+ /// access.
+ class const_iterator {
+ public:
+ typedef std::forward_iterator_tag iterator_category;
+ typedef ConstMonoPtr value_type;
+
+ const_iterator(): mIt(), mEntriesPerMono(0) {}
+ const_iterator(const const_iterator& it):
+ mIt(it.mIt), mEntriesPerMono(it.mEntriesPerMono) {}
+
+ bool operator==(const const_iterator& it) const {return mIt == it.mIt;}
+ bool operator!=(const const_iterator& it) const {return mIt != it.mIt;}
+
+ ConstMonoRef operator*() {
+ MATHICGB_ASSERT(debugValid());
+ return *ConstMonoPtr(&*mIt);
+ }
+
+ const_iterator operator++() {
+ MATHICGB_ASSERT(debugValid());
+ mIt += mEntriesPerMono;
+ return *this;
+ }
+
+ private:
+ friend class MonoVector;
+ bool debugValid() {return mEntriesPerMono > 0;}
+
+ const_iterator(
+ typename RawVector::const_iterator it,
+ size_t entriesPerMono
+ ): mIt(it), mEntriesPerMono(entriesPerMono) {}
+
+ typename RawVector::const_iterator mIt;
+ size_t mEntriesPerMono;
+ };
+
+ // ** Constructors and assignment
+ MonoVector(const MonoMonoid& monoid): mMonoid(monoid) {}
+ MonoVector(const MonoVector& v): mMonos(v.mMonos), mMonoid(v.monoid()) {}
+ MonoVector(MonoVector&& v):
+ mMonos(std::move(v.mMonos)), mMonoid(v.monoid()) {}
+
+ MonoVector& operator=(const MonoVector& v) {
+ MATHICGB_ASSERT(monoid() == v.monoid());
+ mMonos = v.mMonos;
+ return *this;
+ }
+
+ MonoVector& operator=(MonoVector&& v) {
+ MATHICGB_ASSERT(monoid() == v.monoid());
+ mMonos = std::move(v.mMonos);
+ return *this;
+ }
+
+ // ** Iterators
+ const_iterator begin() const {
+ return const_iterator(mMonos.begin(), mMonoid.entriesPerMono());
+ }
+
+ const_iterator end() const {
+ return const_iterator(mMonos.end(), mMonoid.entriesPerMono());
+ }
+
+ const_iterator cbegin() const {return begin();}
+ const_iterator cend() const {return end();}
+
+ // ** Capacity
+ size_t size() const {return mMonos.size() / monoid().entriesPerMono();}
+ bool empty() const {return mMonos.empty();}
+
+ // ** Element access
+ ConstMonoRef front() const {
+ MATHICGB_ASSERT(!empty());
+ return *begin();
+ }
+
+ MonoRef back() {
+ MATHICGB_ASSERT(!empty());
+ const auto offset = mMonos.size() - monoid().entriesPerMono();
+ return *MonoPtr(mMonos.data() + offset);
+ }
+
+ ConstMonoRef back() const {
+ MATHICGB_ASSERT(!empty());
+ const auto offset = mMonos.size() - monoid().entriesPerMono();
+ return *ConstMonoPtr(mMonos.data() + offset);
+ }
+
+ // ** Modifiers
+ void push_back(ConstMonoRef mono) {
+ const auto offset = mMonos.size();
+ mMonos.resize(offset + monoid().entriesPerMono());
+ monoid().copy(mono, *MonoPtr(mMonos.data() + offset));
+ }
+
+ /// Appends the identity.
+ void push_back() {
+ const auto offset = mMonos.size();
+ mMonos.resize(offset + monoid().entriesPerMono());
+ }
+
+ void swap(MonoVector& v) {
+ MATHICGB_ASSERT(&monoid() == &v.monoid());
+ mMonos.swap(v.mMonos);
+ }
+
+ void clear() {mMonos.clear();}
+
+ // ** Relational operators
+ bool operator==(const MonoVector& v) const {
+ MATHICGB_ASSERT(monoid() == v.monoid());
+ return mMonos == v.mMonos;
+ }
+ bool operator!=(const MonoVector& v) const {return !(*this == v);}
+
+ // ** Other
+ size_t memoryBytesUsed() const {
+ return mMonos.capacity() * sizeof(mMonos[0]);
+ }
+
+ const MonoMonoid& monoid() const {return mMonoid;}
+
+ private:
+ RawVector mMonos;
+ const MonoMonoid& mMonoid;
+ };
+
+
+private:
+ size_t entriesPerMono() const {return varCount();}
+
+ const VarIndex mVarCount;
+};
+
+#endif
diff --git a/src/mathicgb/PrimeField.hpp b/src/mathicgb/PrimeField.hpp
index 7a90c54..b00c5e4 100644
--- a/src/mathicgb/PrimeField.hpp
+++ b/src/mathicgb/PrimeField.hpp
@@ -1,146 +1,146 @@
-#ifndef MATHICGB_PRIME_FIELD_GUARD
-#define MATHICGB_PRIME_FIELD_GUARD
-
-#include <limits>
-#include <type_traits>
-#include <ostream>
-
-/// Implements arithmetic in a prime field. T must be an unsigned integer type
-/// that is used to store the elements of the field. The characteristic of the
-/// field must be a prime not exceeding std::numeric_limits<T>::max().
-template<class T>
-class PrimeField {
-public:
- class Element {
- public:
- static_assert(!std::numeric_limits<T>::is_signed, "");
- static_assert(std::numeric_limits<T>::is_integer, "");
-
- Element(const Element& e): mValue(e.value()) {}
-
- Element& operator=(const Element& e) {
- mValue = e.value();
- return *this;
- }
-
- bool operator==(const Element e) const {return value() == e.value();}
- bool operator!=(const Element e) const {return !(*this == e);}
- T value() const {return mValue;}
-
- private:
- friend class PrimeField;
- Element(const T value): mValue(value) {}
-
- friend class PrimeFile;
- T mValue;
- };
-
- PrimeField(const T primeCharacteristic): mCharac(primeCharacteristic) {}
-
- Element zero() const {return Element(0);}
- Element one() const {return Element(1);}
-
- bool isZero(const Element a) const {return a == zero();}
- bool isOne(const Element a) const {return a == one();}
-
- T charac() const {return mCharac;}
-
- template<class Integer>
- Element toElement(Integer&& i) const {
- typedef typename std::remove_reference<Integer>::type NoRefInteger;
- static_assert(std::numeric_limits<NoRefInteger>::is_integer, "");
-
- // We need to take the modulus of i to put it into the range [0;charac()).
- // That is more tricky to get right than it might seem.
- //
- // The sign of a % b is implementation defined in C++ if either of a or b
- // are negative. We need the positive remainder so the operands have to be
- // positive. Another reason for this is that we could not allow b to be
- // converted to a signed integer since it might not be representable that
- // way.
- //
- // If Integer is signed and i is std::numeric_limits<Integer>::min() then
- // it is undefined behavior to evaluate the expression -i since -i is not
- // representable, leading to a signed overflow. So we have to cast to
- // unsigned before doing the minus.
- typedef typename std::make_unsigned<NoRefInteger>::type Unsigned;
- if (i < 0) {
- // Negate i to get a positive number, then negate again to cancel out
- // the first negation. The first cast to unsigned is to avoid
- // undefined behavior from -i. The second is there because apparently,
- // at least on GCC, the unary - re-introduces signed and we need to
- // be unsigned to avoid sign extension. We need zero extension.
- const auto unsignedNegative =
- static_cast<Unsigned>(-static_cast<Unsigned>(i));
- return negative(Element(unsignedNegative % charac()));
- } else
- return Element(static_cast<Unsigned>(i) % charac());
- }
-
- T toValue(Element e) const {return e.value;}
-
- Element sum(const Element a, const Element b) const {
- const auto s = a.value() + b.value();
- // The sum overflowed if and only if a.value() > s. In that case
- // subtraction of charac() will overflow again in the other direction,
- // leaving us with the correct result.
- if (a.value() > s || s >= charac()) {
- MATHICGB_ASSERT(s - charac() < charac());
- return Element(s - charac());
- } else
- return Element(s);
- }
-
- Element difference(const Element a, const Element b) const {
- if (a.value() < b.value()) {
- MATHICGB_ASSERT(a.value() - b.value() + charac() < charac());
- return Element(a.value() - b.value() + charac());
- } else
- return Element(a.value() - b.value());
- }
-
- Element negative(const Element a) const {
- if (a.value() == 0)
- return a;
- else
- return negativeNonZero(a);
- }
-
- Element negativeNonZero(const Element a) const {
- MATHICGB_ASSERT(!isZero(a));
- return Element(charac() - a.value());
- }
-
- Element product(const Element a, const Element b) const;
-
+#ifndef MATHICGB_PRIME_FIELD_GUARD
+#define MATHICGB_PRIME_FIELD_GUARD
+
+#include <limits>
+#include <type_traits>
+#include <ostream>
+
+/// Implements arithmetic in a prime field. T must be an unsigned integer type
+/// that is used to store the elements of the field. The characteristic of the
+/// field must be a prime not exceeding std::numeric_limits<T>::max().
+template<class T>
+class PrimeField {
+public:
+ class Element {
+ public:
+ static_assert(!std::numeric_limits<T>::is_signed, "");
+ static_assert(std::numeric_limits<T>::is_integer, "");
+
+ Element(const Element& e): mValue(e.value()) {}
+
+ Element& operator=(const Element& e) {
+ mValue = e.value();
+ return *this;
+ }
+
+ bool operator==(const Element e) const {return value() == e.value();}
+ bool operator!=(const Element e) const {return !(*this == e);}
+ T value() const {return mValue;}
+
+ private:
+ friend class PrimeField;
+ Element(const T value): mValue(value) {}
+
+ friend class PrimeFile;
+ T mValue;
+ };
+
+ PrimeField(const T primeCharacteristic): mCharac(primeCharacteristic) {}
+
+ Element zero() const {return Element(0);}
+ Element one() const {return Element(1);}
+
+ bool isZero(const Element a) const {return a == zero();}
+ bool isOne(const Element a) const {return a == one();}
+
+ T charac() const {return mCharac;}
+
+ template<class Integer>
+ Element toElement(Integer&& i) const {
+ typedef typename std::remove_reference<Integer>::type NoRefInteger;
+ static_assert(std::numeric_limits<NoRefInteger>::is_integer, "");
+
+ // We need to take the modulus of i to put it into the range [0;charac()).
+ // That is more tricky to get right than it might seem.
+ //
+ // The sign of a % b is implementation defined in C++ if either of a or b
+ // are negative. We need the positive remainder so the operands have to be
+ // positive. Another reason for this is that we could not allow b to be
+ // converted to a signed integer since it might not be representable that
+ // way.
+ //
+ // If Integer is signed and i is std::numeric_limits<Integer>::min() then
+ // it is undefined behavior to evaluate the expression -i since -i is not
+ // representable, leading to a signed overflow. So we have to cast to
+ // unsigned before doing the minus.
+ typedef typename std::make_unsigned<NoRefInteger>::type Unsigned;
+ if (i < 0) {
+ // Negate i to get a positive number, then negate again to cancel out
+ // the first negation. The first cast to unsigned is to avoid
+ // undefined behavior from -i. The second is there because apparently,
+ // at least on GCC, the unary - re-introduces signed and we need to
+ // be unsigned to avoid sign extension. We need zero extension.
+ const auto unsignedNegative =
+ static_cast<Unsigned>(-static_cast<Unsigned>(i));
+ return negative(Element(unsignedNegative % charac()));
+ } else
+ return Element(static_cast<Unsigned>(i) % charac());
+ }
+
+ T toValue(Element e) const {return e.value;}
+
+ Element sum(const Element a, const Element b) const {
+ const auto s = a.value() + b.value();
+ // The sum overflowed if and only if a.value() > s. In that case
+ // subtraction of charac() will overflow again in the other direction,
+ // leaving us with the correct result.
+ if (a.value() > s || s >= charac()) {
+ MATHICGB_ASSERT(s - charac() < charac());
+ return Element(s - charac());
+ } else
+ return Element(s);
+ }
+
+ Element difference(const Element a, const Element b) const {
+ if (a.value() < b.value()) {
+ MATHICGB_ASSERT(a.value() - b.value() + charac() < charac());
+ return Element(a.value() - b.value() + charac());
+ } else
+ return Element(a.value() - b.value());
+ }
+
+ Element negative(const Element a) const {
+ if (a.value() == 0)
+ return a;
+ else
+ return negativeNonZero(a);
+ }
+
+ Element negativeNonZero(const Element a) const {
+ MATHICGB_ASSERT(!isZero(a));
+ return Element(charac() - a.value());
+ }
+
+ Element product(const Element a, const Element b) const;
+
/// Returns the multiplicative inverse a^-1 mod charac(). a must not be zero.
- Element inverse(const Element a) const;
-
-private:
- const T mCharac;
-};
-
-namespace PrimeFieldInternal {
+ Element inverse(const Element a) const;
+
+private:
+ const T mCharac;
+};
+
+namespace PrimeFieldInternal {
template<class T>
struct ModularProdType {};
template<> struct ModularProdType<uint8> {typedef uint16 type;};
template<> struct ModularProdType<uint16> {typedef uint32 type;};
template<> struct ModularProdType<uint32> {typedef uint64 type;};
-}
-
-template<class T>
-auto PrimeField<T>::product(
- const Element a,
- const Element b
-) const -> Element {
+}
+
+template<class T>
+auto PrimeField<T>::product(
+ const Element a,
+ const Element b
+) const -> Element {
typedef typename PrimeFieldInternal::ModularProdType<T>::type BigT;
BigT bigProd = static_cast<BigT>(a.value()) * b.value();
MATHICGB_ASSERT(a.value() == 0 || bigProd / a.value() == b.value());
return Element(static_cast<T>(bigProd % charac()));
-}
-
-template<class T>
-auto PrimeField<T>::inverse(const Element elementA) const -> Element {
+}
+
+template<class T>
+auto PrimeField<T>::inverse(const Element elementA) const -> Element {
// We do two turns of the extended Euclidian algorithm per
// loop. Usually the sign of x changes each time through the loop,
// but we avoid that by representing every other x as its negative,
@@ -185,15 +185,15 @@ auto PrimeField<T>::inverse(const Element elementA) const -> Element {
MATHICGB_ASSERT(isOne(product(elementA, inverseElement)));
return inverseElement;
}
-
-
-template<class T>
-std::ostream& operator<<(
- std::ostream& out,
- typename PrimeField<T>::Element e
-) {
- out << e.value();
- return out;
-}
-
-#endif
+
+
+template<class T>
+std::ostream& operator<<(
+ std::ostream& out,
+ typename PrimeField<T>::Element e
+) {
+ out << e.value();
+ return out;
+}
+
+#endif
diff --git a/src/test/MonoMonoid.cpp b/src/test/MonoMonoid.cpp
new file mode 100755
index 0000000..0ae7e77
--- /dev/null
+++ b/src/test/MonoMonoid.cpp
@@ -0,0 +1,233 @@
+#include "mathicgb/stdinc.h"
+
+#include "mathicgb/MonoMonoid.hpp"
+#include <gtest/gtest.h>
+#include <sstream>
+
+// expect(i,j) encodes a matrix with interesting bit patterns that
+// are supposed to be likely to surface errors in how monomials are
+// stored inside a vector.
+uint32 expect(size_t mono, size_t var, size_t varCount) {
+ const auto unique = var + varCount * mono + 1;
+
+ while (true) {
+ // 000
+ if (mono == 0)
+ return 0;
+ --mono;
+
+ // 100
+ // 010
+ // 001
+ if (mono < varCount)
+ return var == mono ? unique : 0;
+ mono -= varCount;
+
+ // 000
+ // 100
+ // 110
+ // 111
+ if (mono < varCount + 1)
+ return var < mono ? unique : 0;
+ mono -= varCount + 1;
+
+ // 111
+ // 011
+ // 001
+ // 000
+ if (mono < varCount + 1)
+ return var >= mono ? unique : 0;
+ mono -= varCount + 1;
+
+ // 101
+ // 010
+ if (mono < 4)
+ return (var % 2) == (mono % 2) ? unique : 0;
+ mono -= 4;
+
+ // 100100
+ // 010010
+ // 001001
+ if (mono < 6)
+ return (var % 3) == (mono % 3) ? unique : 0;
+ mono -= 6;
+
+ // mix the patterns
+ mono += var % 17;
+ }
+};
+
+TEST(MonoMonoid, VarCount) {
+ ASSERT_EQ(0, MonoMonoid<uint8>(0).varCount());
+ ASSERT_EQ(1000 * 1000, MonoMonoid<uint8>(1000 * 1000).varCount());
+ ASSERT_EQ(1, MonoMonoid<uint16>(1).varCount());
+ ASSERT_EQ(2, MonoMonoid<uint32>(2).varCount());
+ ASSERT_EQ(12, MonoMonoid<uint64>(12).varCount());
+}
+
+TEST(MonoMonoid, MonoVector) {
+ typedef MonoMonoid<uint32> Monoid;
+ typedef Monoid::VarIndex VarIndex;
+ typedef Monoid::MonoVector MonoVector;
+
+ Monoid monoid(13);
+ MonoVector v(monoid);
+ MonoVector v2(monoid);
+ ASSERT_EQ(v2.monoid(), monoid);
+ const auto varCount = monoid.varCount();
+
+ ASSERT_TRUE(v.empty());
+ size_t count = 1000;
+
+
+ // Not a correctness error, but empty vectors should preferably not
+ // use any memory.
+ ASSERT_EQ(0, v.memoryBytesUsed());
+
+ for (size_t i = 0; i < count; ++i) {
+ ASSERT_EQ(i, v.size());
+ v.push_back(); // push_back, no param
+ ASSERT_GT(v.memoryBytesUsed(), 0);
+ ASSERT_FALSE(v.empty()); // empty
+ ASSERT_EQ(i + 1, v.size()); // size
+
+ ASSERT_TRUE(monoid.isIdentity(v.back())); // isIdentity true, back non-const
+ bool allZero = true;
+ for (VarIndex var = 0; var < varCount; ++var) {
+ const auto exponent = expect(i, var, varCount);
+ if (exponent != 0) {
+ allZero = false;
+ monoid.setExponent(var, exponent, v.back());
+ }
+ }
+ ASSERT_EQ(allZero, monoid.isIdentity(v.back())); // isIdentity false
+ v2.push_back(v.back()); // push_back with param
+ ASSERT_TRUE(monoid.equal(v.back(), v2.back()));
+ }
+ auto it = v.begin();
+ ASSERT_EQ(it, v.cbegin());
+ for (size_t i = 0; i < count; ++i, ++it) {
+ MonoVector::const_iterator tmp;
+ ASSERT_TRUE(it != tmp);
+ tmp = it;
+ ASSERT_EQ(tmp, it);
+ ASSERT_TRUE(v.end() != it);
+
+ for (VarIndex var = 0; var < monoid.varCount(); ++var) {
+ ASSERT_EQ(expect(i, var, varCount), monoid.exponent(*it, var));
+ }
+ }
+ ASSERT_EQ(v.end(), it);
+ ASSERT_EQ(v.cend(), it);
+
+ ASSERT_EQ(v, v2); // operator== true
+ monoid.setExponent(0, 1 + monoid.exponent(v2.back(), 0), v2.back());
+ ASSERT_TRUE(v != v2); // operator!=, true, same length
+
+ auto& vc = const_cast<const MonoVector&>(v);
+ ASSERT_TRUE(monoid.equal(v.front(), *v2.begin())); // front, non-const
+ ASSERT_TRUE(monoid.equal(vc.front(), *v2.begin())); // front, const
+ ASSERT_TRUE(monoid.equal(vc.back(), v.back())); // back, non-const
+
+ auto v3(v2); // copy constructor
+ ASSERT_EQ(v3.monoid(), monoid);
+ ASSERT_TRUE(v != v3 && v2 == v3);
+ v2.swap(v); // member swap
+ ASSERT_TRUE(v == v3 && v2 != v3);
+ std::swap(v, v2); // std::swap
+ ASSERT_TRUE(v != v3 && v2 == v3);
+ using std::swap;
+ swap(v, v2); // let compiler decide which swap to use
+ ASSERT_TRUE(v == v3 && v2 != v3);
+ swap(v, v2); // get back to original state
+ ASSERT_TRUE(v != v3 && v2 == v3);
+
+ ASSERT_FALSE(v3 != v2); // operator!=, false, same length
+ v3.push_back();
+ ASSERT_TRUE(v3 != v2); // operator!=, true, different length
+
+
+ ASSERT_FALSE(v3 == v);
+ v3 = v; // copy assignment
+ ASSERT_EQ(v3.monoid(), monoid);
+ ASSERT_EQ(v3, v);
+
+ ASSERT_FALSE(v3.empty());
+ v2 = std::move(v3); // move assignment
+ ASSERT_EQ(v2.monoid(), monoid);
+ ASSERT_EQ(v2, v);
+ ASSERT_TRUE(v3.empty());
+
+ ASSERT_FALSE(v2.empty());
+ auto v4(std::move(v2)); // move constructor
+ ASSERT_EQ(v4.monoid(), monoid);
+ ASSERT_TRUE(v2.empty());
+ ASSERT_EQ(v4, v);
+
+ ASSERT_FALSE(v.empty());
+ v.clear();
+ ASSERT_TRUE(v.empty());
+}
+
+TEST(MonoMonoid, MonoPool) {
+ typedef MonoMonoid<uint32> Monoid;
+ typedef Monoid::VarIndex VarIndex;
+ typedef Monoid::Mono Mono;
+
+ for (int q = 0; q < 2; ++q) {
+ Monoid monoid(13);
+ Monoid::MonoPool pool(monoid);
+ const auto varCount = monoid.varCount();
+
+ const auto count = 1000;
+ std::vector<Mono> monos;
+ for (int i = 0; i < count; ++i) {
+ pool.alloc();
+ pool.free(pool.alloc());
+ auto m1 = pool.alloc();
+ ASSERT_TRUE(monoid.isIdentity(m1));
+ auto m2 = pool.alloc();
+ ASSERT_TRUE(monoid.isIdentity(m2));
+ for (VarIndex var = 0; var < varCount; ++var) {
+ monoid.setExponent(var, 1, m1);
+ monoid.setExponent(var, 1, m2);
+ }
+ if (i > 10) {
+ using std::swap;
+ swap(m2, monos[i - 10]);
+ }
+ monos.push_back(std::move(m1));
+ }
+
+ // This ensures that we get to each entry in monos exactly once.
+ MATHICGB_ASSERT((count % 17) != 0);
+ int i = 0;
+ do {
+ MATHICGB_ASSERT(!monos[i].isNull());
+ ASSERT_FALSE(monoid.isIdentity(monos[i]));
+ pool.free(monos[i]);
+ ASSERT_TRUE(monos[i].isNull());
+ pool.free(monos[i]);
+ ASSERT_TRUE(monos[i].isNull());
+ i = (i + 17) % count;
+ } while (i != 0);
+
+ // If the ordering of monomials inside the pool has anything to do with
+ // allocation and deallocation order, then the monomials inside the
+ // pool are at this point all jumbled around. All the entries were also
+ // non-zero before, so we test that new allocations are the identity.
+
+ for (int i = 0; i < count; ++i) {
+ monos[i] = pool.alloc();
+ ASSERT_TRUE(monoid.isIdentity(monos[i]));
+ for (VarIndex var = 0; var < varCount; ++var)
+ monoid.setExponent(var, expect(i, var, varCount), monos[i]);
+ }
+ for (int i = 0; i < count; ++i) {
+ for (VarIndex var = 0; var < varCount; ++var) {
+ ASSERT_EQ(expect(i, var, varCount), monoid.exponent(monos[i], var));
+ }
+ }
+ // everything should be free'd now. Let's do all that again.
+ }
+}
--
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