[mathicgb] 225/393: Moving SPairs towards proper use of MonoMonoid.
Doug Torrance
dtorrance-guest at moszumanska.debian.org
Fri Apr 3 15:59:07 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 129b86c5be5c17387343464904c0b285c65522a3
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date: Fri Apr 5 17:18:58 2013 +0200
Moving SPairs towards proper use of MonoMonoid.
---
src/mathicgb/MonoMonoid.hpp | 243 ++++++++++++++++++++++++++++++++++++++++----
src/mathicgb/PolyBasis.cpp | 182 ---------------------------------
src/mathicgb/PolyBasis.hpp | 24 -----
src/mathicgb/PolyRing.cpp | 14 ---
src/mathicgb/PolyRing.hpp | 27 ++---
src/mathicgb/SPairs.cpp | 221 ++++++++++++++++++++++------------------
src/mathicgb/SPairs.hpp | 98 +++++++++++++-----
src/test/MonoMonoid.cpp | 130 ++++++++++++++++++++++--
8 files changed, 553 insertions(+), 386 deletions(-)
diff --git a/src/mathicgb/MonoMonoid.hpp b/src/mathicgb/MonoMonoid.hpp
index ec54cf9..8f79365 100755
--- a/src/mathicgb/MonoMonoid.hpp
+++ b/src/mathicgb/MonoMonoid.hpp
@@ -14,8 +14,6 @@
/// Implements the monoid of (monic) monomials with integer
/// non-negative exponents. Exponent must be an unsigned integer type that is
/// used to store each exponent of a monomial.
-///
-/// TODO: support grading and comparison.
template<
class Exponent,
bool HasComponent = true,
@@ -115,6 +113,8 @@ public:
friend class PolyRing;
static MonoRef toRef(Exponent* e) {return MonoRef(e);}
static ConstMonoRef toRef(const Exponent* e) {return ConstMonoRef(e);}
+ static Exponent* toOld(MonoRef e) {return rawPtr(e);}
+ static const Exponent* toOld(ConstMonoRef e) {return rawPtr(e);}
@@ -134,7 +134,8 @@ public:
return true;
}()
),
- mGrading(mGradingIsTotalDegree ? std::vector<Exponent>() : grading)
+ mGrading(mGradingIsTotalDegree ? std::vector<Exponent>() : grading),
+ mPool(*this)
{
if (!mGradingIsTotalDegree) {
// Take negative values since reverse lex makes bigger values
@@ -158,13 +159,41 @@ public:
mOrderIndexEnd(HasComponent + mVarCount + StoreOrder),
mHashCoefficients(mVarCount),
mGradingIsTotalDegree(true),
- mGrading()
+ mGrading(),
+ mPool(*this)
{
std::srand(0); // To use the same hash coefficients every time.
for (VarIndex var = 0; var < varCount; ++var)
mHashCoefficients[var] = static_cast<HashValue>(std::rand());
}
+ /// Copies the ordering and varCount from monoid.
+ template<class E2, bool HC2, bool SH2, bool SO2>
+ MonoMonoid(const MonoMonoid<E2, HC2, SH2, SO2>& monoid):
+ mVarCount(monoid.varCount()),
+ mOrderEntryCount(StoreOrder),
+ mOrderIndexBegin(HasComponent + mVarCount),
+ mOrderIndexEnd(HasComponent + mVarCount + StoreOrder),
+ mHashCoefficients(mVarCount),
+ mGradingIsTotalDegree(monoid.mGradingIsTotalDegree),
+ mGrading(monoid.mGrading),
+ mPool(*this)
+ {}
+ // the compiler-generated copy constructor is a better match for overload
+ // resolution than the template constructor above, so we need to provide
+ // our own copy constructor to prevent the compiler-generated one from
+ // being used.
+ MonoMonoid(const MonoMonoid& monoid):
+ mVarCount(monoid.varCount()),
+ mOrderEntryCount(StoreOrder),
+ mOrderIndexBegin(HasComponent + mVarCount),
+ mOrderIndexEnd(HasComponent + mVarCount + StoreOrder),
+ mHashCoefficients(mVarCount),
+ mGradingIsTotalDegree(monoid.mGradingIsTotalDegree),
+ mGrading(monoid.mGrading),
+ mPool(*this)
+ {}
+
bool operator==(const MonoMonoid& monoid) const {
return this == &monoid;
}
@@ -223,6 +252,19 @@ public:
return true;
}
+ template<class MonoidA>
+ bool equal(
+ const MonoidA& monoidA,
+ typename MonoidA::ConstMonoRef a,
+ ConstMonoRef b
+ ) const {
+ // todo: assert compatible
+ for (VarIndex var = 0; var < varCount(); ++var)
+ if (monoidA.exponent(a, var) != exponent(b, var))
+ return false;
+ return true;
+ }
+
/// As equal(), but optimized for the case where true is returned.
bool equalHintTrue(ConstMonoRef a, ConstMonoRef b) const {
// if a[i] != b[i] then a[i] ^ b[i] != 0, so the or of all xors is zero
@@ -332,16 +374,46 @@ public:
}
/// Returns true if a divides b. Equal monomials divide each other.
- bool divides(ConstMonoRef a, ConstMonoRef b) const {
- for (auto i = entriesIndexBegin(); i < exponentsIndexEnd(); ++i)
- if (access(a, i) > access(b, i))
+ bool divides(ConstMonoRef div, ConstMonoRef into) const {
+ // todo: enable this when the code works with it
+ //if (HasComponent && component(div) != component(into))
+ // return false;
+ for (auto i = exponentsIndexBegin(); i < exponentsIndexEnd(); ++i)
+ if (access(div, i) > access(into, i))
+ return false;
+ return true;
+ }
+
+ template<class MonoidA>
+ bool divides(
+ const MonoidA& monoidA,
+ typename MonoidA::ConstMonoRef a,
+ ConstMonoRef b
+ ) const {
+ // todo: fix other divisibility functions to work properly for component too.
+ MATHICGB_ASSERT(monoidA.varCount() == varCount());
+ MATHICGB_ASSERT(!MonoidA::HasComponent || HasComponent);
+ MATHICGB_ASSERT(monoidA.debugValid(a));
+ MATHICGB_ASSERT(debugValid(b));
+ // todo: enable this when the code works with it
+ //if (HasComponent && component(div) != component(into))
+ // return false;
+ //if (
+ // MonoidA::HasComponent &&
+ // HasComponent &&
+ // monoidA.component(a) != component(b)
+ //)
+ // return false;
+ for (VarIndex var = 0; var < varCount(); ++var)
+ if (monoidA.exponent(a, var) > exponent(b, var))
return false;
return true;
}
/// Returns true if div divides lcm(a, b).
bool dividesLcm(ConstMonoRef div, ConstMonoRef a, ConstMonoRef b) const {
- MATHICGB_ASSERT(!HasComponent || component(a) == component(b));
+ MATHICGB_ASSERT(debugLcmCheck(*this, a, *this, b));
+ MATHICGB_ASSERT(debugValid(div));
for (auto i = exponentsIndexBegin(); i != exponentsIndexEnd(); ++i) {
const auto dive = access(div, i);
@@ -351,9 +423,29 @@ public:
return true;
}
+ template<class MonoidDiv, class MonoidA>
+ bool dividesLcm(
+ const MonoidDiv& monoidDiv,
+ typename MonoidDiv::ConstMonoRef div,
+ const MonoidA& monoidA,
+ typename MonoidA::ConstMonoRef a,
+ ConstMonoRef b
+ ) const {
+ MATHICGB_ASSERT(monoidDiv.debugLcmCheck(monoidA, a, *this, b));
+ MATHICGB_ASSERT(monoidDiv.debugValid(div));
+
+ for (VarIndex var = 0; var < varCount(); ++var) {
+ const auto e = monoidDiv.exponent(div, var);
+ if (e > monoidA.exponent(a, var) && e > exponent(b, var))
+ return false;
+ }
+ return true;
+ }
+
/// Returns true if lcm(a,b) == lcmAB.
bool isLcm(ConstMonoRef a, ConstMonoRef b, ConstMonoRef lcmAB) const {
- MATHICGB_ASSERT(!HasComponent || component(a) == component(b));
+ MATHICGB_ASSERT(debugLcmCheck(*this, a, *this, b));
+ MATHICGB_ASSERT(debugValid(lcmAB));
for (auto i = exponentsIndexBegin(); i != exponentsIndexEnd(); ++i)
if (access(lcmAB, i) != std::max(access(a, i), access(b, i)))
@@ -361,6 +453,38 @@ public:
return true;
}
+ template<class MonoidA, class MonoidB>
+ bool isLcm(
+ const MonoidA& monoidA,
+ typename MonoidA::ConstMonoRef a,
+ const MonoidB& monoidB,
+ typename MonoidB::ConstMonoRef b,
+ ConstMonoRef lcmAB
+ ) const {
+ MATHICGB_ASSERT(debugLcmCheck(monoidA, a, monoidB, b));
+ MATHICGB_ASSERT(debugValid(lcmAB));
+
+ if (HasComponent) {
+ if (MonoidA::HasComponent) {
+ if (monoidA.component(a) != component(lcmAB))
+ return false;
+ } else {
+ MATHICGB_ASSERT(MonoidB::HasComponent);
+ if (monoidB.component(b) != component(lcmAB))
+ return false;
+ }
+ }
+
+ for (VarIndex var = 0; var < varCount(); ++var) {
+ if (
+ ptr(lcmAB, exponentsIndexBegin())[var] !=
+ std::max(monoidA.exponent(a, var), monoidB.exponent(b, var))
+ )
+ return false;
+ }
+ return true;
+ }
+
// Graded reverse lexicographic order. The grading is total degree.
CompareResult compare(ConstMonoRef a, ConstMonoRef b) const {
MATHICGB_ASSERT(debugOrderValid(a));
@@ -430,6 +554,32 @@ public:
MATHICGB_ASSERT(debugValid(to));
}
+ template<class MonoidFrom>
+ void copy(
+ const MonoidFrom& monoidFrom,
+ typename MonoidFrom::ConstMonoRef from,
+ MonoRef to
+ ) const {
+ // todo: extract this in checker method
+ MATHICGB_ASSERT(HasComponent == MonoidFrom::HasComponent);
+ MATHICGB_ASSERT(monoidFrom.debugValid(from));
+ MATHICGB_ASSERT(monoidFrom.varCount() == varCount());
+ MATHICGB_ASSERT
+ ((std::is_same<Exponent, typename MonoidFrom::Exponent>::value));
+
+ if (HasComponent)
+ access(to, componentIndex()) = monoidFrom.component(from);
+ for (VarIndex var = 0; var < varCount(); ++var)
+ ptr(to, exponentsIndexBegin())[var] = monoidFrom.exponent(from, var);
+ if (StoreOrder)
+ access(to, orderIndexBegin()) = monoidFrom.degree(from);
+ if (StoreHash)
+ access(to, hashIndex()) = monoidFrom.hash(from);
+
+ MATHICGB_ASSERT(debugValid(to));
+ // todo: check equal
+ }
+
/// Set the exponent of var to newExponent in mono.
void setExponent(
const VarIndex var,
@@ -577,6 +727,7 @@ public:
MATHICGB_ASSERT(debugValid(bColonA));
}
+ /// Sets lcmAB to the lcm of a and b.
void lcm(ConstMonoRef a, ConstMonoRef b, MonoRef lcmAB) const {
if (HasComponent) {
MATHICGB_ASSERT(component(a) == component(b));
@@ -588,22 +739,39 @@ public:
setHash(lcmAB);
MATHICGB_ASSERT(debugValid(lcmAB));
+ MATHICGB_ASSERT(isLcm(a, b, lcmAB));
}
- /// todo: get rid of this
- void lcmRaw(ConstMonoRef a, ConstMonoRef b, MonoRef lcmAB) const {
- MATHICGB_ASSERT(component(a) == component(b));
+ template<class MonoidA, class MonoidB>
+ void lcm(
+ const MonoidA& monoidA,
+ typename MonoidA::ConstMonoRef a,
+ const MonoidB& monoidB,
+ typename MonoidB::ConstMonoRef b,
+ MonoRef lcmAB
+ ) const {
+ MATHICGB_ASSERT(debugLcmCheck(monoidA, a, monoidB, b));
- // Loop also sets component.
- access(lcmAB, componentIndex()) = 0;
- for (auto i = exponentsIndexBegin(); i != exponentsIndexEnd(); ++i)
- access(lcmAB, i) = std::max(access(a, i), access(b, i));
- //setOrderData(lcmAB);
- //setHash(lcmAB);
+ if (HasComponent) {
+ access(lcmAB, componentIndex()) =
+ MonoidA::HasComponent ? monoidA.component(a) : monoidB.component(b);
+ }
- //MATHICGB_ASSERT(debugValid(lcmAB));
+ for (VarIndex var = 0; var < varCount(); ++var) {
+ ptr(lcmAB, exponentsIndexBegin())[var] =
+ std::max(monoidA.exponent(a, var), monoidB.exponent(b, var));
+ }
+
+ setOrderData(lcmAB);
+ setHash(lcmAB);
+
+ MATHICGB_ASSERT(debugValid(lcmAB));
+ MATHICGB_ASSERT(isLcm(monoidA, a, monoidB, b, lcmAB));
}
+ Mono alloc() const {return mPool.alloc();}
+ void free(Mono&& mono) const {mPool.free(std::move(mono));}
+
/// Parses a monomial out of a string. Valid examples: 1 abc a2bc
/// aA. Variable names are case sensitive. Whitespace terminates the
/// parse as does any other character that is not a letter or a
@@ -777,6 +945,9 @@ public:
const MonoMonoid& monoid() const {return mMonoid;}
private:
+ MonoPool(const MonoPool&); // not available
+ void operator=(const MonoPool&); // not available
+
const MonoMonoid& mMonoid;
memt::BufferPool mPool;
};
@@ -953,6 +1124,11 @@ public:
private:
+ // Grants access to other template instantiations.
+ template<class E2, bool HC2, bool SH2, bool SO2>
+ friend class MonoMonoid;
+
+ // The main point here is to grant access to rawPtr().
friend class Mono;
friend class MonoRef;
friend class ConstMonoRef;
@@ -967,6 +1143,33 @@ private:
return true;
}
+ template<class MonoidA, class MonoidB>
+ bool debugLcmCheck(
+ const MonoidA& monoidA,
+ typename MonoidA::ConstMonoRef a,
+ const MonoidB& monoidB,
+ typename MonoidB::ConstMonoRef b
+ ) const {
+ MATHICGB_ASSERT(monoidA.varCount() == varCount());
+ MATHICGB_ASSERT(monoidB.varCount() == varCount());
+ MATHICGB_ASSERT
+ ((std::is_same<Exponent, typename MonoidA::Exponent>::value));
+ MATHICGB_ASSERT
+ ((std::is_same<Exponent, typename MonoidB::Exponent>::value));
+ MATHICGB_ASSERT
+ (HasComponent == (MonoidA::HasComponent || MonoidB::HasComponent));
+ MATHICGB_ASSERT(monoidA.debugValid(a));
+ MATHICGB_ASSERT(monoidB.debugValid(b));
+ MATHICGB_ASSERT(
+ !HasComponent ||
+ !MonoidA::HasComponent ||
+ !MonoidB::HasComponent ||
+ monoidA.component(a) == monoidB.component(b)
+ );
+
+ return true;
+ }
+
// *** Accessing fields of a monomial
template<class M>
static auto rawPtr(M&& m) -> decltype(m.internalRawPtr()) {
@@ -1178,6 +1381,8 @@ private:
/// vector with this vector. If mGradingIsTotalDegree is true then
/// mGrading is empty but implicitly it is a vector of ones.
std::vector<Exponent> mGrading;
+
+ mutable MonoPool mPool;
};
namespace MonoMonoidHelper {
diff --git a/src/mathicgb/PolyBasis.cpp b/src/mathicgb/PolyBasis.cpp
index 73cd31b..3a5190f 100755
--- a/src/mathicgb/PolyBasis.cpp
+++ b/src/mathicgb/PolyBasis.cpp
@@ -185,188 +185,6 @@ size_t PolyBasis::getMemoryUse() const {
return sum;
}
-bool PolyBasis::buchbergerLcmCriterion(size_t a, size_t b) const {
- MATHICGB_ASSERT(a != b);
- MATHICGB_ASSERT(!retired(a));
- MATHICGB_ASSERT(!retired(b));
- // We don't need to set the weights on the lcm since we will only use it
- // for testing divisibility, not for determining order.
- monomial lcmAB = mRing.allocMonomial();
- mRing.monomialLeastCommonMultipleNoWeights
- (leadMonomial(a), leadMonomial(b), lcmAB);
- bool value = buchbergerLcmCriterion(a, b, lcmAB);
- mRing.freeMonomial(lcmAB);
- return value;
-}
-
-bool PolyBasis::buchbergerLcmCriterion
- (size_t a, size_t b, const_monomial lcmAB) const
-{
- MATHICGB_ASSERT(a != b);
- MATHICGB_ASSERT(!retired(a));
- MATHICGB_ASSERT(!retired(b));
- MATHICGB_ASSERT(mRing.monomialIsLeastCommonMultipleNoWeights
- (leadMonomial(a), leadMonomial(b), lcmAB));
-
- class Criterion : public DivisorLookup::EntryOutput {
- public:
- Criterion(size_t a, size_t b, const_monomial lcmAB, const PolyBasis& basis):
- mA(a), mB(b),
- mLcmAB(lcmAB),
- mRing(basis.ring()),
- mBasis(basis),
- mHit(static_cast<size_t>(-1)),
- mAlmostApplies(false) {}
-
- virtual bool proceed(size_t index) {
- MATHICGB_ASSERT(index < mBasis.size());
- MATHICGB_ASSERT(!applies()); // should have stopped search in this case
- MATHICGB_ASSERT(
- mRing.monomialIsDivisibleBy(mLcmAB, mBasis.leadMonomial(index)));
- if (!mBasis.leadMinimal(index))
- return true;
- if (index == mA || index == mB)
- return true;
- mAlmostApplies = true;
-
- if (index < mA || index < mB) {
- // check lcm(a,index) != lcm(a,b) <=>
- // exists i such that max(a[i], c[i]) != max(a[i],b[i]) <=>
- // exists i such that b[i] > a[i] && b[i] > c[i]
- const_monomial leadA = mBasis.leadMonomial(mA);
- const_monomial leadB = mBasis.leadMonomial(mB);
- const_monomial leadC = mBasis.leadMonomial(index);
- if (//index > mA &&
- !mRing.monomialHasStrictlyLargerExponent(leadB, leadC, leadA))
- return true;
-
- // check lcm(b,index) != lcm(a,b)
- if (//index > mB &&
- !mRing.monomialHasStrictlyLargerExponent(leadA, leadC, leadB))
- return true;
- }
-
- mHit = index;
- return false; // stop search
- }
-
- const_monomial lcmAB() const {return mLcmAB;}
- bool almostApplies() const {return mAlmostApplies;}
- bool applies() const {return mHit != static_cast<size_t>(-1);}
- size_t hit() const {return mHit;}
-
- private:
- size_t mA;
- size_t mB;
- const_monomial mLcmAB;
- const PolyRing& mRing;
- const PolyBasis& mBasis;
- size_t mHit; // the divisor that made the criterion apply
- bool mAlmostApplies; // applies ignoring lcm(a,b)=lcm(a,c) complication
- };
-
- ++mStats.buchbergerLcmQueries;
- bool applies = false;
- bool almostApplies = false;
- {
- Criterion criterion(a, b, lcmAB, *this);
- if (mUseBuchbergerLcmHitCache) {
- // Check cacheB first since when I tried this there was a higher hit rate
- // for cacheB than cacheA. Might not be a persistent phenomenon, but
- // there's no downside to trying out cacheB first so I'm going for that.
- //
- // I update the cache if the second check is a hit but not if the first
- // check is a hit. In the one test I did, the worst hit rate was from
- // updating the cache every time, the second best hit rate was from
- // not updating the cache (from cache hits) and the best hit rate was
- // from doing this.
- //
- // The idea is that when the first cache check is a hit,
- // the second cache member might have been a hit too, and updating it
- // might replace a high hit rate element with a low hit rate element,
- // which would be bad. When the second cache check is a hit, we know
- // that the first one wasn't (or we would have taken an early exit),
- // so we have reason to suspect that the first cache element is not
- // a high hit rate element. So it should be better to replace it.
- // That idea seems to be right since it worked better in the one
- // test I did.
- size_t cacheB = mBuchbergerLcmHitCache[b];
- if (!applies && !retired(cacheB) &&
- ring().monomialIsDivisibleBy
- (criterion.lcmAB(), leadMonomial(cacheB)))
- applies = !criterion.Criterion::proceed(cacheB);
-
- size_t cacheA = mBuchbergerLcmHitCache[a];
- if (!applies && !retired(cacheA) &&
- ring().monomialIsDivisibleBy
- (criterion.lcmAB(), leadMonomial(cacheA))) {
- applies = !criterion.Criterion::proceed(cacheA);
- if (applies)
- mBuchbergerLcmHitCache[b] = cacheA;
- }
- }
- if (applies)
- ++mStats.buchbergerLcmCacheHits;
- else {
- MATHICGB_ASSERT(!criterion.applies());
- mDivisorLookup->divisors(criterion.lcmAB(), criterion);
- applies = criterion.applies();
-
- if (mUseBuchbergerLcmHitCache && applies) {
- MATHICGB_ASSERT(criterion.hit() < size());
- mBuchbergerLcmHitCache[a] = criterion.hit();
- mBuchbergerLcmHitCache[b] = criterion.hit();
- }
- }
- if (!applies)
- almostApplies = criterion.almostApplies();
- }
-
- if (applies)
- ++mStats.buchbergerLcmHits;
- else if (almostApplies)
- ++mStats.buchbergerLcmNearHits;
-
- MATHICGB_ASSERT(applies == buchbergerLcmCriterionSlow(a, b));
- return applies;
-}
-
-bool PolyBasis::buchbergerLcmCriterionSlow(size_t a, size_t b) const {
- MATHICGB_ASSERT(a != b);
- MATHICGB_ASSERT(!retired(a));
- MATHICGB_ASSERT(!retired(b));
-
- monomial lcmAB = ring().allocMonomial();
- monomial lcm = ring().allocMonomial();
- ring().monomialLeastCommonMultiple
- (leadMonomial(a), leadMonomial(b), lcmAB);
- size_t stop = size();
- size_t i = 0;
- for (; i < stop; ++i) {
- if (retired(i) || !leadMinimal(i))
- continue;
- if (!ring().monomialIsDivisibleBy(lcmAB, leadMonomial(i)))
- continue;
- if (i == a || i == b)
- continue;
- if (i < a || i < b) {
- ring().monomialLeastCommonMultiple
- (leadMonomial(a), leadMonomial(i), lcm);
- if (ring().monomialEQ(lcmAB, lcm))
- continue;
-
- ring().monomialLeastCommonMultiple
- (leadMonomial(b), leadMonomial(i), lcm);
- if (ring().monomialEQ(lcmAB, lcm))
- continue;
- }
- break;
- }
- ring().freeMonomial(lcmAB);
- ring().freeMonomial(lcm);
- return i != stop;
-}
-
PolyBasis::Entry::Entry():
poly(0),
leadMinimal(0),
diff --git a/src/mathicgb/PolyBasis.hpp b/src/mathicgb/PolyBasis.hpp
index 4cdebe8..ddf71e9 100755
--- a/src/mathicgb/PolyBasis.hpp
+++ b/src/mathicgb/PolyBasis.hpp
@@ -72,30 +72,6 @@ public:
};
Stats stats() const {return mStats;}
- // Returns true if Buchberger's second criterion for eliminating useless
- // S-pairs applies to the pair (a,b). Let
- // l(a,b) = lcm(lead(a), lead(b)).
- // The criterion says that if there is some other basis element c such that
- // lead(c)|l(a,b)
- // and
- // l(a,c) reduces to zero, and
- // l(b,c) reduces to zero
- // then (a,b) will reduce to zero (using classic non-signature reduction).
- //
- // This criterion is less straight forward to apply in case for example
- // l(a,b) = l(a,c) = l(b,c)
- // since then there is the potential to erroneously eliminate all the three
- // pairs among a,b,c on the assumption that the other two pairs will reduce
- // to zero. In such cases, we eliminate the pair with the lowest indexes.
- // This allows removing generators that get non-minimal lead term without
- // problems.
- bool buchbergerLcmCriterion(size_t a, size_t b) const;
-
- // As the overload with an lcmAb parameter, except lcmAB must be the lcm of
- // the lead monomials of a and b. Then this quantity does not have to be
- // computed for the criterion.
- bool buchbergerLcmCriterion(size_t a, size_t b, const_monomial lcmAB) const;
-
// As the non-slow version, but uses simpler and slower code.
bool buchbergerLcmCriterionSlow(size_t a, size_t b) const;
diff --git a/src/mathicgb/PolyRing.cpp b/src/mathicgb/PolyRing.cpp
index f8e0540..0c6ce37 100755
--- a/src/mathicgb/PolyRing.cpp
+++ b/src/mathicgb/PolyRing.cpp
@@ -222,20 +222,6 @@ void PolyRing::monomialGreatestCommonDivisor(ConstMonomial a,
setWeightsOnly(g);
}
-bool PolyRing::monomialIsLeastCommonMultipleNoWeights(
- ConstMonomial a,
- ConstMonomial b,
- ConstMonomial l) const
-{
- if (*l != 0)
- return false;
- for (size_t i = 1; i <= mNumVars; ++i)
- if (l[i] != std::max(a[i], b[i]))
- return false;
- return true;
-}
-
-
void PolyRing::mysteriousSPairMonomialRoutine(ConstMonomial newSig,
ConstMonomial newLead,
ConstMonomial baseDivSig,
diff --git a/src/mathicgb/PolyRing.hpp b/src/mathicgb/PolyRing.hpp
index e3cef31..a192659 100755
--- a/src/mathicgb/PolyRing.hpp
+++ b/src/mathicgb/PolyRing.hpp
@@ -18,6 +18,7 @@
#define EQ 0
#define GT 1
+
template<class T>
PrimeField<
typename std::make_unsigned<
@@ -91,6 +92,7 @@ typedef long coefficient;
typedef MonoMonoid<exponent> Monoid;
typedef PrimeField<unsigned long> Field;
+
typedef exponent* vecmonomial; // includes a component
typedef coefficient const_coefficient;
@@ -195,6 +197,9 @@ struct term {
class PolyRing {
public:
+ typedef MonoMonoid<exponent> Monoid;
+ typedef PrimeField<unsigned long> Field;
+
PolyRing(coefficient charac, int nvars, const std::vector<exponent>& weights);
PolyRing(coefficient charac, int nvars, int nweights);
~PolyRing() {}
@@ -452,18 +457,10 @@ public:
ConstMonomial b,
Monomial& l) const;
- inline void monomialLeastCommonMultipleNoWeights(ConstMonomial a,
- ConstMonomial b,
- Monomial& l) const;
-
bool monomialIsLeastCommonMultiple(ConstMonomial a,
ConstMonomial b,
ConstMonomial l) const;
- bool monomialIsLeastCommonMultipleNoWeights(ConstMonomial a,
- ConstMonomial b,
- ConstMonomial l) const;
-
// Returns true if there is a variable var such that hasLarger raises var to
// a strictly greater exponent than both smaller1 and smaller2 does.
inline bool monomialHasStrictlyLargerExponent(
@@ -498,6 +495,9 @@ public:
const coefficientStats & getCoefficientStats() const { return mStats; }
void resetCoefficientStats() const;
+ const Monoid& monoid() const {return mMonoid;}
+ const Field field() const {return mField;}
+
private:
inline HashValue computeHashValue(const_monomial a1) const;
@@ -518,10 +518,7 @@ private:
bool mTotalDegreeGradedOnly;
- const Monoid& monoid() const {return mMonoid;}
Monoid mMonoid;
-
- const Field field() const {return mField;}
Field mField;
};
@@ -651,14 +648,6 @@ inline void PolyRing::monomialLeastCommonMultiple(
monoid().lcm(a, b, l);
}
-inline void PolyRing::monomialLeastCommonMultipleNoWeights(
- ConstMonomial a,
- ConstMonomial b,
- Monomial& l) const
-{
- monoid().lcmRaw(a, b, l);
-}
-
inline bool PolyRing::monomialHasStrictlyLargerExponent(
ConstMonomial hasLarger,
ConstMonomial smaller1,
diff --git a/src/mathicgb/SPairs.cpp b/src/mathicgb/SPairs.cpp
index 5ca5f43..01aecbe 100755
--- a/src/mathicgb/SPairs.cpp
+++ b/src/mathicgb/SPairs.cpp
@@ -23,9 +23,12 @@ MATHICGB_DEFINE_LOG_ALIAS(
);
SPairs::SPairs(const PolyBasis& basis, bool preferSparseSPairs):
- mQueue(QueueConfiguration(basis, preferSparseSPairs)),
+ mQueue(QueueConfiguration(basis, basis.ring().monoid(), preferSparseSPairs)),
mBasis(basis),
- mRing(basis.ring()) {}
+ mRing(basis.ring()),
+ mMonoid(basis.ring().monoid()),
+ mOrderMonoid(mMonoid),
+ mBareMonoid(mMonoid) {}
std::pair<size_t, size_t> SPairs::pop() {
MATHICGB_LOG_TIME(SPairLate);
@@ -34,23 +37,24 @@ std::pair<size_t, size_t> SPairs::pop() {
MATHICGB_ASSERT(mEliminated.columnCount() == mBasis.size());
while (!mQueue.empty()) {
- std::pair<size_t, size_t> p;
- p = mQueue.topPair();
+ const auto p = mQueue.topPair();
if (mBasis.retired(p.first) || mBasis.retired(p.second)) {
mQueue.pop();
continue;
}
- const_monomial lcm = mQueue.topPairData();
- MATHICGB_ASSERT(mRing.monomialIsLeastCommonMultiple
- (mBasis.leadMonomial(p.first),
- mBasis.leadMonomial(p.second), lcm));
- // Can't pop before done with lcm as popping overwrites lcm.
+ auto lcm = bareMonoid().alloc(); // todo: just keep one around instead
+ bareMonoid().copy(monoid(), mQueue.topPairData(), lcm);
+ mQueue.pop();
+
+ MATHICGB_ASSERT(bareMonoid().isLcm(
+ monoid(), mBasis.leadMonomial(p.first),
+ monoid(), mBasis.leadMonomial(p.second),
+ lcm
+ ));
if (!advancedBuchbergerLcmCriterion(p.first, p.second, lcm)) {
- mQueue.pop();
mEliminated.setBit(p.first, p.second, true);
return p;
}
- mQueue.pop();
}
return std::make_pair(static_cast<size_t>(-1), static_cast<size_t>(-1));
}
@@ -61,25 +65,23 @@ std::pair<size_t, size_t> SPairs::pop(exponent& w) {
// Must call addPairs for new elements before popping.
MATHICGB_ASSERT(mEliminated.columnCount() == mBasis.size());
- while (!mQueue.empty()) {
- std::pair<size_t, size_t> p;
- p = mQueue.topPair();
- if (mBasis.retired(p.first) || mBasis.retired(p.second)) {
- mQueue.pop();
+ for (; !mQueue.empty(); mQueue.pop()) {
+ const auto p = mQueue.topPair();
+ if (mBasis.retired(p.first) || mBasis.retired(p.second))
continue;
- }
- const_monomial lcm = mQueue.topPairData();
- MATHICGB_ASSERT(mRing.monomialIsLeastCommonMultiple
- (mBasis.leadMonomial(p.first),
- mBasis.leadMonomial(p.second), lcm));
- // Can't pop before done with lcm as popping overwrites lcm.
- if (advancedBuchbergerLcmCriterion(p.first, p.second, lcm)) {
- mQueue.pop();
+ auto lcm = bareMonoid().alloc(); // todo: just keep one around instead
+ bareMonoid().copy(monoid(), mQueue.topPairData(), lcm);
+
+ MATHICGB_ASSERT(bareMonoid().isLcm(
+ monoid(), mBasis.leadMonomial(p.first),
+ monoid(), mBasis.leadMonomial(p.second),
+ lcm
+ ));
+ if (advancedBuchbergerLcmCriterion(p.first, p.second, lcm))
continue;
- }
if (w == 0)
- w = mRing.weight(lcm);
- else if (w != mRing.weight(lcm))
+ w = bareMonoid().degree(lcm);
+ else if (w != bareMonoid().degree(lcm))
break;
mQueue.pop();
mEliminated.setBit(p.first, p.second, true);
@@ -202,7 +204,7 @@ void SPairs::addPairs(size_t newGen) {
}
- typedef std::pair<monomial, Queue::Index> PrePair;
+ typedef std::pair<Mono, Queue::Index> PrePair;
std::vector<PrePair> prePairs;
MATHICGB_ASSERT(prePairs.empty());
@@ -210,29 +212,30 @@ void SPairs::addPairs(size_t newGen) {
throw std::overflow_error
("Too large basis element index in constructing S-pairs.");
- const_monomial const newLead = mBasis.leadMonomial(newGen);
- monomial lcm = mBasis.ring().allocMonomial();
+ Monoid::MonoPool pool(mMonoid);
+ ConstMonoRef newLead = mBasis.leadMonomial(newGen);
+ auto lcm = mBareMonoid.alloc();
for (size_t oldGen = 0; oldGen < newGen; ++oldGen) {
if (mBasis.retired(oldGen))
continue;
- const_monomial const oldLead = mBasis.leadMonomial(oldGen);
- if (mRing.monomialRelativelyPrime(newLead, oldLead)) {
+ ConstMonoRef oldLead = mBasis.leadMonomial(oldGen);
+ if (mMonoid.relativelyPrime(newLead, oldLead)) {
++mStats.relativelyPrimeHits;
mEliminated.setBit(newGen, oldGen, true);
continue;
}
- mRing.monomialLeastCommonMultipleNoWeights(newLead, oldLead, lcm);
+ mBareMonoid.lcm(mMonoid, newLead, mMonoid, oldLead, lcm);
if (simpleBuchbergerLcmCriterion(newGen, oldGen, lcm)) {
mEliminated.setBit(newGen, oldGen, true);
continue;
}
- mRing.setWeightsOnly(lcm);
- auto newLcm = mBasis.ring().allocMonomial();
- mBasis.ring().monomialCopy(lcm, newLcm);
- prePairs.emplace_back(newLcm, static_cast<Queue::Index>(oldGen));
+ auto orderLcm = monoid().alloc(); // todo: use orderMonoid
+ // todo: convert lcm instead of re-computing
+ monoid().lcm(monoid(), newLead, monoid(), oldLead, orderLcm);
+ prePairs.emplace_back
+ (std::move(orderLcm), static_cast<Queue::Index>(oldGen));
}
- mBasis.ring().freeMonomial(lcm);
std::sort(prePairs.begin(), prePairs.end(),
[&](const PrePair& a, const PrePair& b)
@@ -244,32 +247,41 @@ void SPairs::addPairs(size_t newGen) {
(makeSecondIterator(prePairs.begin()), makeSecondIterator(prePairs.end()));
for (auto it = prePairs.begin(); it != prePairs.end(); ++it)
- mRing.freeMonomial(it->first);
+ mMonoid.free(std::move(it->first));
}
size_t SPairs::getMemoryUse() const {
return mQueue.getMemoryUse();
}
-bool SPairs::simpleBuchbergerLcmCriterion
-(size_t a, size_t b, const_monomial lcmAB) const
-{
+bool SPairs::simpleBuchbergerLcmCriterion(
+ size_t a,
+ size_t b,
+ BareMonoid::ConstMonoRef lcmAB
+) const {
MATHICGB_ASSERT(a < mBasis.size());
MATHICGB_ASSERT(b < mBasis.size());
MATHICGB_ASSERT(a != b);
MATHICGB_ASSERT(!mBasis.retired(a));
MATHICGB_ASSERT(!mBasis.retired(b));
- MATHICGB_ASSERT(mRing.monomialIsLeastCommonMultipleNoWeights
- (mBasis.leadMonomial(a), mBasis.leadMonomial(b), lcmAB));
+ MATHICGB_ASSERT(bareMonoid().isLcm
+ (monoid(), mBasis.leadMonomial(a), monoid(), mBasis.leadMonomial(b), lcmAB)
+ );
MATHICGB_ASSERT(mEliminated.columnCount() == mBasis.size());
class Criterion : public DivisorLookup::EntryOutput {
public:
- Criterion(size_t a, size_t b, const_monomial lcmAB, const SPairs& sPairs):
+ Criterion(
+ const size_t a,
+ const size_t b,
+ BareMonoid::ConstMonoRef lcmAB,
+ const SPairs& sPairs
+ ):
mA(a), mB(b),
mLcmAB(lcmAB),
mSPairs(sPairs),
- mRing(sPairs.ring()),
+ mMonoid(sPairs.monoid()),
+ mBareMonoid(sPairs.bareMonoid()),
mBasis(sPairs.basis()),
mHit(static_cast<size_t>(-1)),
mAlmostApplies(false) {}
@@ -277,41 +289,45 @@ bool SPairs::simpleBuchbergerLcmCriterion
virtual bool proceed(size_t index) {
MATHICGB_ASSERT(index < mBasis.size());
MATHICGB_ASSERT(!applies()); // should have stopped search in this case
- MATHICGB_ASSERT(mRing.monomialIsDivisibleBy(mLcmAB, mBasis.leadMonomial(index)));
+ MATHICGB_ASSERT
+ (mBareMonoid.divides(mMonoid, mBasis.leadMonomial(index), mLcmAB));
if (index == mA || index == mB)
return true;
mAlmostApplies = true;
// check lcm(a,index) != lcm(a,b) <=>
// exists i such that max(a[i], c[i]) != max(a[i],b[i]) <=>
- // exists i such that b[i] > a[i] && b[i] > c[i]
+ // exists i such that b[i] > a[i] && b[i] > c[i] <=>
+ // exists i such that b[i] > max(a[i], c[i]) <=>
+ // b does not divide lcm(a[i], c[i])
const_monomial leadA = mBasis.leadMonomial(mA);
const_monomial leadB = mBasis.leadMonomial(mB);
const_monomial leadC = mBasis.leadMonomial(index);
if (!mSPairs.eliminated(index, mA) &&
- !mRing.monomialHasStrictlyLargerExponent(leadB, leadC, leadA))
- return true;
+ mMonoid.dividesLcm(leadB, leadC, leadA))
+ return true; // we had lcm(a,index) == lcm(a,b)
// check lcm(b,index) != lcm(a,b)
if (!mSPairs.eliminated(index, mB) &&
- !mRing.monomialHasStrictlyLargerExponent(leadA, leadC, leadB))
- return true;
+ mMonoid.dividesLcm(leadA, leadC, leadB))
+ return true; // we had lcm(b,index) == lcm(a,b)
mHit = index;
return false; // stop search
}
- const_monomial lcmAB() const {return mLcmAB;}
+ BareMonoid::ConstMonoRef lcmAB() const {return mLcmAB;}
bool almostApplies() const {return mAlmostApplies;}
bool applies() const {return mHit != static_cast<size_t>(-1);}
size_t hit() const {return mHit;}
private:
- size_t mA;
- size_t mB;
- const_monomial mLcmAB;
+ const size_t mA;
+ const size_t mB;
+ BareMonoid::ConstMonoRef mLcmAB;
const SPairs& mSPairs;
- const PolyRing& mRing;
+ const Monoid& mMonoid;
+ const BareMonoid& mBareMonoid;
const PolyBasis& mBasis;
size_t mHit; // the divisor that made the criterion apply
bool mAlmostApplies; // applies ignoring lcm(a,b)=lcm(a,c) complication
@@ -342,30 +358,35 @@ bool SPairs::simpleBuchbergerLcmCriterion
// That idea seems to be right since it worked better in the one
// test I did.
size_t cacheB = mBuchbergerLcmHitCache[b];
- if (!applies && !mBasis.retired(cacheB) &&
- mRing.monomialIsDivisibleBy
- (criterion.lcmAB(), mBasis.leadMonomial(cacheB)))
+ if (
+ !applies &&
+ !mBasis.retired(cacheB) &&
+ mBareMonoid.divides
+ (mMonoid, mBasis.leadMonomial(cacheB), criterion.lcmAB())
+ )
applies = !criterion.Criterion::proceed(cacheB);
size_t cacheA = mBuchbergerLcmHitCache[a];
- if (!applies && !mBasis.retired(cacheA) &&
- mRing.monomialIsDivisibleBy
- (criterion.lcmAB(), mBasis.leadMonomial(cacheA))) {
+ if (
+ !applies &&
+ !mBasis.retired(cacheA) &&
+ mBareMonoid.divides
+ (mMonoid, mBasis.leadMonomial(cacheA), criterion.lcmAB())
+ ) {
applies = !criterion.Criterion::proceed(cacheA);
if (applies)
mBuchbergerLcmHitCache[b] = cacheA;
}
}
- if (applies)
- {
- if (mStats.late)
- ++mStats.buchbergerLcmCacheHitsLate;
- else
- ++mStats.buchbergerLcmCacheHits;
- }
- else {
+ if (applies) {
+ if (mStats.late)
+ ++mStats.buchbergerLcmCacheHitsLate;
+ else
+ ++mStats.buchbergerLcmCacheHits;
+ } else {
MATHICGB_ASSERT(!criterion.applies());
- mBasis.divisorLookup().divisors(criterion.lcmAB(), criterion);
+ mBasis.divisorLookup().divisors
+ (BareMonoid::toOld(criterion.lcmAB()), criterion);
applies = criterion.applies();
if (mUseBuchbergerLcmHitCache && applies) {
@@ -408,7 +429,7 @@ bool SPairs::simpleBuchbergerLcmCriterionSlow(size_t a, size_t b) const {
for (; i < stop; ++i) {
if (mBasis.retired(i))
continue;
- if (!mRing.monomialIsDivisibleBy(lcmAB, mBasis.leadMonomial(i)))
+ if (!mMonoid.divides(mBasis.leadMonomial(i), lcmAB))
continue;
if (i == a || i == b)
continue;
@@ -432,13 +453,16 @@ bool SPairs::simpleBuchbergerLcmCriterionSlow(size_t a, size_t b) const {
return i != stop;
}
-bool SPairs::advancedBuchbergerLcmCriterion
- (size_t a, size_t b, const_monomial lcmAB) const
-{
+bool SPairs::advancedBuchbergerLcmCriterion(
+ size_t a,
+ size_t b,
+ BareMonoid::ConstMonoRef lcmAB
+) const {
MATHICGB_ASSERT(a != b);
MATHICGB_ASSERT(mEliminated.columnCount() == mBasis.size());
- MATHICGB_ASSERT(mRing.monomialIsLeastCommonMultipleNoWeights
- (mBasis.leadMonomial(a), mBasis.leadMonomial(b), lcmAB));
+ MATHICGB_ASSERT(bareMonoid().isLcm
+ (monoid(), mBasis.leadMonomial(a), monoid(), mBasis.leadMonomial(b), lcmAB)
+ );
mStats.late = true;
if (simpleBuchbergerLcmCriterion(a, b, lcmAB)) {
@@ -467,7 +491,7 @@ bool SPairs::advancedBuchbergerLcmCriterion
Graph& graph = mAdvancedBuchbergerLcmCriterionGraph;
graph.clear();
GraphBuilder builder(graph);
- mBasis.divisorLookup().divisors(lcmAB, builder);
+ mBasis.divisorLookup().divisors(BareMonoid::toOld(lcmAB), builder);
if (graph.size() <= 3) {
// For the graph approach to be better than the simpler approach of
@@ -513,12 +537,15 @@ bool SPairs::advancedBuchbergerLcmCriterion
// Note that
// lcm(c,d) != lcmAB <=>
// exists i such that max(c[i], d[i]) < lcmAB[i] <=>
- // exists i such that lcmAB[i] > c[i] && lcmAB[i] > d[i]
- if (!eliminated(currentIndex, otherIndex) &&
- !mRing.monomialHasStrictlyLargerExponent(lcmAB, currentLead, otherLead))
- {
+ // exists i such that lcmAB[i] > c[i] && lcmAB[i] > d[i] <=>
+ // exists i such that lcmAB[i] > max(c[i], d[i]) <=>
+ // lcmAB does not divide lcm(c[i], d[i])
+ if (
+ !eliminated(currentIndex, otherIndex) &&
+ monoid().dividesLcm
+ (bareMonoid(), lcmAB, monoid(), currentLead, otherLead)
+ )
continue; // not an edge in G
- }
if (otherConnect == NotConnected) {
other->second = currentConnect;
@@ -547,10 +574,9 @@ bool SPairs::advancedBuchbergerLcmCriterionSlow(size_t a, size_t b) const {
MATHICGB_ASSERT(a != b);
MATHICGB_ASSERT(mEliminated.columnCount() == mBasis.size());
- monomial lcmAB = mRing.allocMonomial();
- monomial lcm = mRing.allocMonomial();
- mRing.monomialLeastCommonMultiple
- (mBasis.leadMonomial(a), mBasis.leadMonomial(b), lcmAB);
+ auto lcmAB = mMonoid.alloc();
+ auto lcm = mMonoid.alloc();
+ mMonoid.lcm(mBasis.leadMonomial(a), mBasis.leadMonomial(b), lcmAB);
size_t stop = mBasis.size();
// *** Build the graph vertices
@@ -562,7 +588,7 @@ bool SPairs::advancedBuchbergerLcmCriterionSlow(size_t a, size_t b) const {
for (size_t i = 0; i != stop; ++i) {
if (mBasis.retired(i))
continue;
- if (!mRing.monomialIsDivisibleBy(lcmAB, mBasis.leadMonomial(i)))
+ if (!mMonoid.divides(mBasis.leadMonomial(i), lcmAB))
continue;
Connection con = NotConnected;
if (i == a) {
@@ -586,7 +612,7 @@ bool SPairs::advancedBuchbergerLcmCriterionSlow(size_t a, size_t b) const {
MATHICGB_ASSERT(node.second != NotConnected);
// loop through all potential edges (node.first, i)
- const_monomial leadNode = mBasis.leadMonomial(node.first);
+ ConstMonoRef leadNode = mBasis.leadMonomial(node.first);
for (size_t i = 0; i < graph.size(); ++i) {
if (node.second == graph[i].second)
continue;
@@ -594,8 +620,8 @@ bool SPairs::advancedBuchbergerLcmCriterionSlow(size_t a, size_t b) const {
size_t const other = graph[i].first;
const_monomial const leadOther = mBasis.leadMonomial(other);
- mRing.monomialLeastCommonMultiple(leadNode, leadOther, lcm);
- if (!eliminated(node.first, other) && mRing.monomialEQ(lcm, lcmAB))
+ mMonoid.lcm(leadNode, leadOther, lcm);
+ if (!eliminated(node.first, other) && mMonoid.equal(lcm, lcmAB))
continue; // not an edge in G
if (graph[i].second == NotConnected) {
@@ -610,8 +636,6 @@ bool SPairs::advancedBuchbergerLcmCriterionSlow(size_t a, size_t b) const {
}
}
}
- mRing.freeMonomial(lcmAB);
- mRing.freeMonomial(lcm);
MATHICGB_ASSERT(applies || !simpleBuchbergerLcmCriterionSlow(a, b));
return applies;
@@ -630,9 +654,8 @@ std::string SPairs::name() const {
void SPairs::QueueConfiguration::computePairData(
size_t a,
size_t b,
- monomial orderBy
+ MonoRef orderBy
) const {
- MATHICGB_ASSERT(!orderBy.isNull());
MATHICGB_ASSERT(a != b);
MATHICGB_ASSERT(a < mBasis.size());
MATHICGB_ASSERT(b < mBasis.size());
@@ -640,8 +663,8 @@ void SPairs::QueueConfiguration::computePairData(
// todo: do something special here?
return; //return false;
}
- const_monomial const leadA = mBasis.leadMonomial(a);
- const_monomial const leadB = mBasis.leadMonomial(b);
- mBasis.ring().monomialLeastCommonMultiple(leadA, leadB, orderBy);
+ ConstMonoRef leadA = mBasis.leadMonomial(a);
+ ConstMonoRef leadB = mBasis.leadMonomial(b);
+ mMonoid.lcm(leadA, leadB, orderBy);
return; //todo: return true;
}
diff --git a/src/mathicgb/SPairs.hpp b/src/mathicgb/SPairs.hpp
index a81f3f4..2cfa8b1 100755
--- a/src/mathicgb/SPairs.hpp
+++ b/src/mathicgb/SPairs.hpp
@@ -12,6 +12,28 @@ class PolyBasis;
// algorithm. Also eliminates useless S-pairs and orders the S-pairs.
class SPairs {
public:
+ typedef PolyRing::Monoid Monoid;
+ typedef Monoid::Exponent Exponent;
+ typedef Monoid::Mono Mono;
+ typedef Monoid::MonoPtr MonoPtr;
+ typedef Monoid::ConstMonoPtr ConstMonoPtr;
+ typedef Monoid::MonoRef MonoRef;
+ typedef Monoid::ConstMonoRef ConstMonoRef;
+
+ /// This monoid is used for computations to determine whether to eliminate
+ /// an S-pair. These computations do not require anything beyond just
+ /// considering the exponents. Since recomputing characteristics of
+ /// lcms such as hash and degree is expensive (unlike for a product), it is
+ /// worthwhile to disable those characteristics that we do not need.
+ typedef MonoMonoid<Exponent, true, false, false> BareMonoid;
+ //typedef Monoid BareMonoid;
+
+ /// This monoid is used to order S-pairs by their lcm. Here we need to
+ /// to store the ordering data for fast comparison, but we do not need
+ /// hashes.
+ typedef MonoMonoid<Exponent, true, false, true> OrderMonoid;
+ //typedef Monoid OrderMonoid;
+
SPairs(const PolyBasis& basis, bool preferSparseSPairs);
// Returns the number of S-pairs in the data structure.
@@ -58,6 +80,7 @@ public:
return mEliminated.bitUnordered(a, b);
}
+ const Monoid& monoid() const {return ring().monoid();}
const PolyRing& ring() const {return mRing;}
const PolyBasis& basis() const {return mBasis;}
@@ -89,6 +112,9 @@ public:
std::string name() const;
private:
+ const BareMonoid& bareMonoid() const {return mBareMonoid;}
+ const OrderMonoid& orderMonoid() const {return mOrderMonoid;}
+
// Returns true if Buchberger's second criterion for eliminating useless
// S-pairs applies to the pair (a,b). Define
// l(a,b) = lcm(lead(a), lead(b)).
@@ -107,8 +133,11 @@ private:
// cases to the advanced criterion, except if an S-pair has already been
// eliminated - in that case we do not check to see if the lcm's are the same
// as it is not necessary to do so.
- bool simpleBuchbergerLcmCriterion
- (size_t a, size_t b, const_monomial lcmAB) const;
+ bool simpleBuchbergerLcmCriterion(
+ size_t a,
+ size_t b,
+ BareMonoid::ConstMonoRef lcmAB
+ ) const;
// As the non-slow version, but uses simpler and slower code.
bool simpleBuchbergerLcmCriterionSlow(size_t a, size_t b) const;
@@ -123,24 +152,31 @@ private:
// been eliminated. It is a theorem that if there is a path from a to b
// in G then (a,b) is a useless S-pair that can be eliminated.
bool advancedBuchbergerLcmCriterion
- (size_t a, size_t b, const_monomial lcmAB) const;
+ (size_t a, size_t b, BareMonoid::ConstMonoRef lcmAB) const;
// As the non-slow version, but uses simpler and slower code.
bool advancedBuchbergerLcmCriterionSlow(size_t a, size_t b) const;
class QueueConfiguration {
public:
- QueueConfiguration(const PolyBasis& basis, const bool preferSparseSPairs):
- mBasis(basis), mPreferSparseSPairs(preferSparseSPairs) {}
-
- typedef monomial PairData;
- void computePairData(size_t col, size_t row, monomial m) const;
+ QueueConfiguration(
+ const PolyBasis& basis,
+ const Monoid& monoid,
+ const bool preferSparseSPairs
+ ):
+ mBasis(basis),
+ mMonoid(mBasis.ring().monoid()),
+ mPreferSparseSPairs(preferSparseSPairs) {}
+
+ typedef Mono PairData;
+ void computePairData(size_t col, size_t row, MonoRef m) const;
typedef bool CompareResult;
- bool compare(size_t colA, size_t rowA, const_monomial a,
- size_t colB, size_t rowB, const_monomial b) const
- {
- const auto cmp = mBasis.ring().monomialCompare(a, b);
+ bool compare(
+ size_t colA, size_t rowA, ConstMonoRef a,
+ size_t colB, size_t rowB, ConstMonoRef b
+ ) const {
+ const auto cmp = monoid().compare(a, b);
if (cmp == GT)
return true;
if (cmp == LT)
@@ -162,13 +198,17 @@ private:
}
bool cmpLessThan(bool v) const {return v;}
- // these are not required for a configuration but we will use
- // them from this code. TODO
- monomial allocPairData() {return mBasis.ring().allocMonomial();}
- void freePairData(monomial m) {return mBasis.ring().freeMonomial(m);}
+ // The following methods are not required of a configuration.
+ Mono allocPairData() {return monoid().alloc();}
+ void freePairData(Mono&& mono) {
+ return monoid().free(std::move(mono));
+ }
private:
+ const Monoid& monoid() const {return mMonoid;}
+
const PolyBasis& mBasis;
+ const Monoid& mMonoid;
const bool mPreferSparseSPairs;
};
typedef mathic::PairQueue<QueueConfiguration> Queue;
@@ -183,6 +223,10 @@ private:
const PolyRing& mRing;
mutable Stats mStats;
+ const Monoid& mMonoid;
+ OrderMonoid mOrderMonoid;
+ BareMonoid mBareMonoid;
+
static const bool mUseBuchbergerLcmHitCache = true;
mutable std::vector<size_t> mBuchbergerLcmHitCache;
@@ -198,26 +242,34 @@ private:
friend void mathic::PairQueueNamespace::constructPairData<QueueConfiguration>
(void*, Index, Index, QueueConfiguration&);
friend void mathic::PairQueueNamespace::destructPairData<QueueConfiguration>
- (monomial*, Index, Index, QueueConfiguration&);
+ (Mono*, Index, Index, QueueConfiguration&);
};
namespace mathic {
namespace PairQueueNamespace {
template<>
- inline void constructPairData<SPairs::QueueConfiguration>
- (void* memory, Index col, Index row, SPairs::QueueConfiguration& conf) {
+ inline void constructPairData<SPairs::QueueConfiguration>(
+ void* memory,
+ const Index col,
+ const Index row,
+ SPairs::QueueConfiguration& conf
+ ) {
MATHICGB_ASSERT(memory != 0);
MATHICGB_ASSERT(col > row);
- monomial* pd = new (memory) monomial(conf.allocPairData());
+ auto pd = new (memory) SPairs::Mono(conf.allocPairData());
conf.computePairData(col, row, *pd);
}
template<>
- inline void destructPairData
- (monomial* pd, Index col, Index row, SPairs::QueueConfiguration& conf) {
+ inline void destructPairData(
+ SPairs::Mono* pd,
+ const Index col,
+ const Index row,
+ SPairs::QueueConfiguration& conf
+ ) {
MATHICGB_ASSERT(pd != 0);
MATHICGB_ASSERT(col > row);
- conf.freePairData(*pd);
+ conf.freePairData(std::move(*pd));
}
}
}
diff --git a/src/test/MonoMonoid.cpp b/src/test/MonoMonoid.cpp
index 890c40f..0763d7e 100755
--- a/src/test/MonoMonoid.cpp
+++ b/src/test/MonoMonoid.cpp
@@ -385,17 +385,24 @@ TYPED_TEST(Monoid, MultiplyDivide) {
ASSERT_TRUE(m.compare(c, mono) == Monoid::EqualTo);
ASSERT_EQ(m.hash(c), m.hash(mono));
- // check properties that mono=a*b should have
+ // divides, check properties that mono=a*b should have
ASSERT_TRUE(m.divides(mono, c));
ASSERT_TRUE(m.divides(c, mono));
ASSERT_TRUE(m.divides(a, mono));
ASSERT_TRUE(m.divides(b, mono));
+ // divides, general
+ ASSERT_TRUE(m.divides(m, mono, c));
+ ASSERT_TRUE(m.divides(m, c, mono));
+ ASSERT_TRUE(m.divides(m, a, mono));
+ ASSERT_TRUE(m.divides(m, b, mono));
+
if (!m.isIdentity(a)) {
ASSERT_TRUE(m.lessThan(b, mono));
ASSERT_FALSE(m.lessThan(mono, b));
ASSERT_TRUE(m.compare(mono, b) == Monoid::GreaterThan);
ASSERT_FALSE(m.divides(mono, b));
+ ASSERT_FALSE(m.divides(m, mono, b));
ASSERT_FALSE(m.isProductOf(a, c, b));
ASSERT_FALSE(m.isProductOfHintTrue(a, c, b));
@@ -406,6 +413,7 @@ TYPED_TEST(Monoid, MultiplyDivide) {
ASSERT_TRUE(m.equal(b, mono));
ASSERT_TRUE(m.compare(b, mono) == Monoid::EqualTo);
ASSERT_TRUE(m.divides(mono, b));
+ ASSERT_TRUE(m.divides(m, mono, b));
}
if (!m.isIdentity(b)) {
@@ -413,6 +421,7 @@ TYPED_TEST(Monoid, MultiplyDivide) {
ASSERT_FALSE(m.lessThan(mono, a));
ASSERT_TRUE(m.compare(mono, a) == Monoid::GreaterThan);
ASSERT_FALSE(m.divides(mono, a));
+ ASSERT_FALSE(m.divides(m, mono, a));
ASSERT_FALSE(m.isProductOf(c, b, a));
ASSERT_FALSE(m.isProductOfHintTrue(b, c, a));
@@ -422,7 +431,7 @@ TYPED_TEST(Monoid, MultiplyDivide) {
} else {
ASSERT_TRUE(m.equal(a, mono));
ASSERT_TRUE(m.compare(a, mono) == Monoid::EqualTo);
- ASSERT_TRUE(m.divides(mono, a));
+ ASSERT_TRUE(m.divides(m, mono, a));
}
// Check that aliased parameters work.
@@ -453,7 +462,8 @@ TYPED_TEST(Monoid, MultiplyDivide) {
TYPED_TEST(Monoid, LcmColon) {
typedef TypeParam Monoid;
- Monoid m(49);
+ Monoid mNonConst(49);
+ auto& m = mNonConst;
typename Monoid::MonoPool pool(m);
auto mono = pool.alloc();
auto mono2 = pool.alloc();
@@ -466,11 +476,13 @@ TYPED_TEST(Monoid, LcmColon) {
const auto& b = *++v.begin();
const auto& lcm = v.back();
- // isLcm
+ // isLcm (+general)
ASSERT_TRUE(m.isLcm(a, b, lcm));
+ ASSERT_TRUE(m.isLcm(m, a, m, b, lcm));
m.copy(lcm, mono);
m.setExponent(1, m.exponent(mono, 1) + 1, mono);
ASSERT_FALSE(m.isLcm(a, b, mono));
+ ASSERT_FALSE(m.isLcm(m, a, m, b, mono));
// dividesLcm
ASSERT_TRUE(m.dividesLcm(lcm, a, b));
@@ -480,14 +492,22 @@ TYPED_TEST(Monoid, LcmColon) {
ASSERT_TRUE(m.dividesLcm(b, b, b));
ASSERT_TRUE(m.dividesLcm(b, b, a));
+ // dividesLcm, general
+ ASSERT_TRUE(m.dividesLcm(m, lcm, m, a, b));
+ ASSERT_FALSE(m.dividesLcm(m, mono, m, a, b));
+ ASSERT_TRUE(m.dividesLcm(m, a, m, a, a));
+ ASSERT_TRUE(m.dividesLcm(m, a, m, a, b));
+ ASSERT_TRUE(m.dividesLcm(m, b, m, b, b));
+ ASSERT_TRUE(m.dividesLcm(m, b, m, b, a));
+
// lcm(a, b)
m.lcm(a, b, mono);
ASSERT_TRUE(m.equal(mono, lcm));
ASSERT_TRUE(m.compare(mono, lcm) == Monoid::EqualTo);
ASSERT_EQ(m.hash(lcm), m.hash(mono));
- // lcm(b, a)
- m.lcm(b, a, mono);
+ // lcm(b, a), general
+ m.lcm(m, b, m, a, mono);
ASSERT_TRUE(m.equal(mono, lcm));
ASSERT_TRUE(m.compare(mono, lcm) == Monoid::EqualTo);
ASSERT_EQ(m.hash(lcm), m.hash(mono));
@@ -616,3 +636,101 @@ TYPED_TEST(Monoid, HasAmpleCapacityTotalDegree) {
}
}
}
+
+TYPED_TEST(Monoid, CopyEqualConversion) {
+ typedef TypeParam Monoid;
+ typedef typename Monoid::Exponent Exponent;
+ typedef typename Monoid::VarIndex VarIndex;
+ static const bool HasComponent = Monoid::HasComponent;
+ typedef MonoMonoid<Exponent, HasComponent, false, false> MonoidNone;
+ typedef MonoMonoid<Exponent, HasComponent, true, true> MonoidAll;
+ for (VarIndex varCount = 1; varCount < 33; ++varCount) {
+ MonoidNone none(varCount);
+ Monoid some(none);
+ MonoidAll all(some);
+
+ auto none1 = none.alloc();
+ auto none2 = none.alloc();
+ auto none3 = none.alloc();
+ auto some1 = some.alloc();
+ auto some2 = some.alloc();
+ auto some3 = some.alloc();
+ auto all1 = all.alloc();
+ auto all2 = all.alloc();
+ auto all3 = all.alloc();
+
+ none.setExponent(0, 1, none1);
+ none.setExponent(varCount / 2, 2, none1);
+ none.setExponent(varCount - 1, 3, none1);
+ none.copy(none1, none2);
+ none.setExponent(0, 4, none2);
+
+ some.setExponent(0, 1, some1);
+ some.setExponent(varCount / 2, 2, some1);
+ some.setExponent(varCount - 1, 3, some1);
+ some.copy(some1, some2);
+ some.setExponent(0, 4, some2);
+
+ all.setExponent(0, 1, all1);
+ all.setExponent(varCount / 2, 2, all1);
+ all.setExponent(varCount - 1, 3, all1);
+ all.copy(all1, all2);
+ all.setExponent(0, 4, all2);
+
+ // compare on none
+ ASSERT_TRUE(none.equal(none, none1, none1));
+ ASSERT_TRUE(none.equal(some, some1, none1));
+ ASSERT_TRUE(none.equal(all, all1, none1));
+ ASSERT_FALSE(none.equal(none, none1, none2));
+ ASSERT_FALSE(none.equal(some, some1, none2));
+ ASSERT_FALSE(none.equal(all, all1, none2));
+
+ // compare on some
+ ASSERT_TRUE(some.equal(none, none1, some1));
+ ASSERT_TRUE(some.equal(some, some1, some1));
+ ASSERT_TRUE(some.equal(all, all1, some1));
+ ASSERT_FALSE(some.equal(none, none1, some2));
+ ASSERT_FALSE(some.equal(some, some1, some2));
+ ASSERT_FALSE(some.equal(all, all1, some2));
+
+ // compare on all
+ ASSERT_TRUE(all.equal(none, none1, all1));
+ ASSERT_TRUE(all.equal(some, some1, all1));
+ ASSERT_TRUE(all.equal(all, all1, all1));
+ ASSERT_FALSE(all.equal(none, none1, all2));
+ ASSERT_FALSE(all.equal(some, some1, all2));
+ ASSERT_FALSE(all.equal(all, all1, all2));
+
+ // convert some->none
+ none.copy(some, some1, none3);
+ ASSERT_TRUE(none.equal(none1, none3));
+ ASSERT_FALSE(none.equal(none2, none3));
+ none.copy(some, some2, none3);
+ ASSERT_FALSE(none.equal(none1, none3));
+ ASSERT_TRUE(none.equal(none2, none3));
+
+ /// convert some->all
+ all.copy(some, some1, all3);
+ ASSERT_TRUE(all.equal(all1, all3));
+ ASSERT_FALSE(all.equal(all2, all3));
+ all.copy(some, some2, all3);
+ ASSERT_FALSE(all.equal(all1, all3));
+ ASSERT_TRUE(all.equal(all2, all3));
+
+ // convert none->some
+ some.copy(none, none1, some3);
+ ASSERT_TRUE(some.equal(some1, some3));
+ ASSERT_FALSE(some.equal(some2, some3));
+ some.copy(none, none2, some3);
+ ASSERT_FALSE(some.equal(some1, some3));
+ ASSERT_TRUE(some.equal(some2, some3));
+
+ // convert Y->some
+ some.copy(none, none1, some3);
+ ASSERT_TRUE(some.equal(some1, some3));
+ ASSERT_FALSE(some.equal(some2, some3));
+ some.copy(none, none2, some3);
+ ASSERT_FALSE(some.equal(some1, some3));
+ ASSERT_TRUE(some.equal(some2, some3));
+ }
+}
--
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