[mathicgb] 216/393: PolyRing can now forward to MonoMonoid and PrimeField for almost everything. Added lots of new methods on those two to make that possible.
Doug Torrance
dtorrance-guest at moszumanska.debian.org
Fri Apr 3 15:59:04 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 d241a72b61720190e65bf860e1e30a4a2e64ad35
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date: Sun Mar 31 21:43:45 2013 +0200
PolyRing can now forward to MonoMonoid and PrimeField for almost everything. Added lots of new methods on those two to make that possible.
---
src/mathicgb/GroebnerBasis.cpp | 17 ++--
src/mathicgb/GroebnerBasis.hpp | 10 ++-
src/mathicgb/MonoMonoid.hpp | 101 ++++++++++++++++++++++++
src/mathicgb/PolyRing.cpp | 6 ++
src/mathicgb/PolyRing.hpp | 174 +++++++++++++++++++++++++++++++++++------
src/mathicgb/PrimeField.hpp | 14 ++++
src/mathicgb/SignatureGB.cpp | 2 +-
src/mathicgb/io-util.cpp | 2 +-
src/test/poly-test.cpp | 56 +++++++------
9 files changed, 321 insertions(+), 61 deletions(-)
diff --git a/src/mathicgb/GroebnerBasis.cpp b/src/mathicgb/GroebnerBasis.cpp
index 9ccab34..98a1313 100755
--- a/src/mathicgb/GroebnerBasis.cpp
+++ b/src/mathicgb/GroebnerBasis.cpp
@@ -1,5 +1,6 @@
// Copyright 2011 Michael E. Stillman
+#define MATHICGB_SLOW_DEBUG
#include "stdinc.h"
#include <iostream>
#include <iomanip>
@@ -173,7 +174,8 @@ size_t GroebnerBasis::regularReducer(
MATHICGB_SLOW_ASSERT(debugValue != static_cast<size_t>(-1));
monomial m = ring().allocMonomial();
MATHICGB_SLOW_ASSERT
- (ring().monomialDivide(term, getLeadMonomial(reducer), m));
+ (ring().monomialIsDivisibleBy(term, getLeadMonomial(reducer)));
+ ring().monomialDivide(term, getLeadMonomial(reducer), m);
ring().monomialMultTo(m, getSignature(reducer));
MATHICGB_SLOW_ASSERT(order().signatureCompare(sig, m) == GT);
ring().freeMonomial(m);
@@ -189,8 +191,9 @@ size_t GroebnerBasis::regularReducerSlow(
monomial m = ring().allocMonomial();
const size_t stop = size();
for (size_t be = 0; be < stop; ++be) {
- if (!ring().monomialDivide(term, getLeadMonomial(be), m))
+ if (!ring().monomialIsDivisibleBy(term, getLeadMonomial(be)))
continue;
+ ring().monomialDivide(term, getLeadMonomial(be), m);
ring().monomialMultTo(m, getSignature(be));
if (order().signatureCompare(sig, m) == GT) {
ring().freeMonomial(m);
@@ -308,8 +311,9 @@ size_t GroebnerBasis::minimalLeadInSigSlow(const_monomial sig) const {
for (size_t gen = 0; gen < genCount; ++gen) {
if (ring().monomialGetComponent(getSignature(gen)) != sigComponent)
continue;
- if (!ring().monomialDivide(sig, getSignature(gen), multiplier))
+ if (!ring().monomialIsDivisibleBy(sig, getSignature(gen)))
continue;
+ ring().monomialDivide(sig, getSignature(gen), multiplier);
if (minLeadGen != static_cast<size_t>(-1)) {
const_monomial genLead = getLeadMonomial(gen);
int leadCmp = order().signatureCompare(minLead, multiplier, genLead);
@@ -345,8 +349,8 @@ size_t GroebnerBasis::minimalLeadInSigSlow(const_monomial sig) const {
return minLeadGen;
}
-bool GroebnerBasis::isSingularTopReducible
- (const Poly& poly, const_monomial sig) const {
+bool GroebnerBasis::isSingularTopReducibleSlow
+(const Poly& poly, const_monomial sig) const {
MATHICGB_ASSERT( ! sig.isNull() );
if (poly.isZero())
return false;
@@ -355,8 +359,9 @@ bool GroebnerBasis::isSingularTopReducible
const size_t genCount = size();
const_monomial polyLead = poly.getLeadMonomial();
for (size_t i = 0; i < genCount; ++i) {
- if (!ring().monomialDivide(polyLead, getLeadMonomial(i), multiplier))
+ if (!ring().monomialIsDivisibleBy(polyLead, getLeadMonomial(i)))
continue;
+ ring().monomialDivide(polyLead, getLeadMonomial(i), multiplier);
if (order().signatureCompare(sig, multiplier, getSignature(i)) == EQ)
return true;
}
diff --git a/src/mathicgb/GroebnerBasis.hpp b/src/mathicgb/GroebnerBasis.hpp
old mode 100644
new mode 100755
index 174fcab..4ec1686
--- a/src/mathicgb/GroebnerBasis.hpp
+++ b/src/mathicgb/GroebnerBasis.hpp
@@ -86,10 +86,12 @@ public:
size_t minimalLeadInSig(const_monomial sig) const;
// Returns true if poly can be singular reduced in signature sig.
- // In other words, returns true if there is a basis element with lead
- // term M and signature S such that M divides the lead term N of
- // poly and such that N/M*S == sig.
- bool isSingularTopReducible(const Poly& poly, const_monomial sig) const;
+ // In other words, returns true if there is a basis element with
+ // lead term M and signature S such that M divides the lead term N
+ // of poly and such that N/M*S == sig. This is slow because it is
+ // currently only used for asserts - implement a fast version if
+ // that changes.
+ bool isSingularTopReducibleSlow(const Poly& poly, const_monomial sig) const;
void display(std::ostream &o) const;
void displayBrief(std::ostream &o) const;
diff --git a/src/mathicgb/MonoMonoid.hpp b/src/mathicgb/MonoMonoid.hpp
index 9e7a070..2a6b2ad 100755
--- a/src/mathicgb/MonoMonoid.hpp
+++ b/src/mathicgb/MonoMonoid.hpp
@@ -310,6 +310,27 @@ public:
return true;
}
+ /// Returns true if div divides lcm(a, b).
+ bool dividesLcm(ConstMonoRef div, ConstMonoRef a, ConstMonoRef b) const {
+ for (auto i = exponentsIndexBegin(); i != exponentsIndexEnd(); ++i) {
+ const auto dive = access(div, i);
+ if (access(div, i) > access(a, i) && access(div, i) > access(b, i))
+ return false;
+ }
+ return true;
+ }
+
+ /// Returns true if lcm(a,b) == lcmAB.
+ bool isLcm(ConstMonoRef a, ConstMonoRef b, ConstMonoRef lcmAB) const {
+ MATHICGB_ASSERT(component(a) == component(b));
+
+ // Loop also checks component.
+ for (auto i = componentIndex(); i != exponentsIndexEnd(); ++i)
+ if (access(lcmAB, i) != std::max(access(a, i), access(b, i)))
+ return false;
+ return true;
+ }
+
// Graded reverse lexicographic order. The grading is total degree.
CompareResult compare(ConstMonoRef a, ConstMonoRef b) const {
MATHICGB_ASSERT(debugOrderValid(a));
@@ -328,6 +349,34 @@ public:
return compare(a, b) == LessThan;
}
+ /// Returns true if gcd(a, b) == 1.
+ bool relativelyPrime(ConstMonoRef a, ConstMonoRef b) const {
+ for (auto i = exponentsIndexBegin(); i != exponentsIndexEnd(); ++i)
+ if (access(a, i) > 0 && access(b, i) > 0)
+ return false;
+ return true;
+ }
+
+ // If this method returns true for monomials a and b then it is
+ // guaranteed that multiplying a and b together will not overflow
+ // the integers in the representation.
+ bool hasAmpleCapacity(ConstMonoRef mono) const {
+ const auto halfMin = std::numeric_limits<Exponent>::min() / 2;
+ const auto halfMax = std::numeric_limits<Exponent>::max() / 2;
+ for (VarIndex i = exponentsIndexBegin(); i != exponentsIndexEnd(); ++i) {
+ const auto value = access(mono, i);
+ if (!(halfMin <= value && value <= halfMax))
+ return false;
+ }
+ return true;
+ }
+
+ /// Returns the degree of mono using the grading on the monoid.
+ Exponent degree(ConstMonoRef mono) const {
+ MATHICGB_ASSERT(debugOrderValid(mono));
+ return access(mono, orderIndexBegin());
+ }
+
// *** Monomial mutating computations
@@ -354,6 +403,17 @@ public:
MATHICGB_ASSERT(debugValid(mono));
}
+ /// Sets all the exponents of mono. exponents must point to an array
+ /// of size varCount().
+ void setExponents(const Exponent* exponents, MonoRef mono) const {
+ MATHICGB_ASSERT(exponents != 0);
+ access(mono, componentIndex()) = 0;
+ std::copy_n(exponents, varCount(), ptr(mono, exponentsIndexBegin()));
+ setOrderData(mono);
+ setHash(mono);
+ MATHICGB_ASSERT(debugValid(mono));
+ }
+
/// Sets mono to 1, which is the identity for multiplication.
void setIdentity(MonoRef mono) const {
std::fill_n(rawPtr(mono), entryCount(), static_cast<Exponent>(0));
@@ -428,6 +488,47 @@ public:
MATHICGB_ASSERT(debugValid(quo));
}
+ /// Sets aColonB to a:b and bColonA to b:a.
+ void colons(
+ ConstMonoRef a,
+ ConstMonoRef b,
+ MonoRef aColonB,
+ MonoRef bColonA
+ ) const {
+ MATHICGB_ASSERT(debugValid(a));
+ MATHICGB_ASSERT(debugValid(b));
+ MATHICGB_ASSERT(component(a) == component(b));
+
+ access(aColonB, componentIndex()) = 0;
+ access(bColonA, componentIndex()) = 0;
+ for (auto i = exponentsIndexBegin(); i != exponentsIndexEnd(); ++i) {
+ const auto ae = access(a, i);
+ const auto be = access(b, i);
+ const auto max = std::max(ae, be);
+ access(aColonB, i) = max - be;
+ access(bColonA, i) = max - ae;
+ }
+ setOrderData(aColonB);
+ setHash(aColonB);
+ setOrderData(bColonA);
+ setHash(bColonA);
+
+ MATHICGB_ASSERT(debugValid(aColonB));
+ MATHICGB_ASSERT(debugValid(bColonA));
+ }
+
+ void lcm(ConstMonoRef a, ConstMonoRef b, MonoRef lcmAB) const {
+ MATHICGB_ASSERT(component(a) == component(b));
+
+ // Loop also sets component.
+ for (auto i = componentIndex(); i != exponentsIndexEnd(); ++i)
+ access(lcmAB, i) = std::max(access(a, i), access(b, i));
+ setOrderData(lcmAB);
+ setHash(lcmAB);
+
+ MATHICGB_ASSERT(debugValid(lcmAB));
+ }
+
/// 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
diff --git a/src/mathicgb/PolyRing.cpp b/src/mathicgb/PolyRing.cpp
index 3fa8ba7..10a9a7a 100755
--- a/src/mathicgb/PolyRing.cpp
+++ b/src/mathicgb/PolyRing.cpp
@@ -32,6 +32,9 @@ PolyRing::PolyRing(
#ifdef MATHICGB_USE_MONOID
, mMonoid(weights)
#endif
+#ifdef MATHICGB_USE_FIELD
+ , mField(p0)
+#endif
{
MATHICGB_ASSERT(weights.size() == nvars);
mTotalDegreeGradedOnly = true;
@@ -63,6 +66,9 @@ PolyRing::PolyRing(coefficient p0,
#ifdef MATHICGB_USE_MONOID
, mMonoid(nvars)
#endif
+#ifdef MATHICGB_USE_FIELD
+ , mField(p0)
+#endif
{
MATHICGB_ASSERT(nweights == 1);
diff --git a/src/mathicgb/PolyRing.hpp b/src/mathicgb/PolyRing.hpp
index 956686e..4eb5bba 100755
--- a/src/mathicgb/PolyRing.hpp
+++ b/src/mathicgb/PolyRing.hpp
@@ -174,6 +174,9 @@ typedef long coefficient;
#ifdef MATHICGB_USE_MONOID
typedef MonoMonoid<exponent> Monoid;
#endif
+#ifdef MATHICGB_USE_FIELD
+typedef PrimeField<unsigned long> Field;
+#endif
typedef exponent* vecmonomial; // includes a component
typedef coefficient const_coefficient;
@@ -351,18 +354,13 @@ public:
P.free(m.unsafeGetRepresentation());
}
-
-
- // Free monomials allocated here by calling freeMonomial().
- Monomial allocMonomial1() const {
- return static_cast<exponent *>(mMonomialPool.alloc());
- }
-
// Only call this method for monomials returned by allocMonomial().
- void freeMonomial(Monomial m) const {mMonomialPool.free(m.unsafeGetRepresentation());}
+ void freeMonomial(Monomial m) const {
+ mMonomialPool.free(m.unsafeGetRepresentation());
+ }
-// Free monomials allocated here by calling freeMonomial().
-monomial allocMonomial() const { return allocMonomial(mMonomialPool); }
+ // Free monomials allocated here by calling freeMonomial().
+ monomial allocMonomial() const { return allocMonomial(mMonomialPool); }
@@ -405,13 +403,26 @@ monomial allocMonomial() const { return allocMonomial(mMonomialPool); }
// Monomial Routines //////////////////////
///////////////////////////////////////////
- HashValue monomialHashValue(ConstMonomial m) const {return static_cast<HashValue>(m[mHashIndex]);}
+ HashValue monomialHashValue(ConstMonomial m) const {
+#ifdef MATHICGB_USE_MONOID
+ return monoid().hash(m);
+#else
+ return static_cast<HashValue>(m[mHashIndex]);
+#endif
+ }
void monomialSetExponent(Monomial m, size_t var, exponent c) const {
+#ifdef MATHICGB_USE_MONOID
+ monoid().setExponent(var, c, m);
+#else
m[var+1] = c;
+#endif
}
void monomialSetExponents(Monomial m, exponent* exponents) const {
+#ifdef MATHICGB_USE_MONOID
+ monoid().setExponents(exponents, m);
+#else
*m = 0;
std::memcpy(
m.unsafeGetRepresentation() + 1,
@@ -419,10 +430,15 @@ monomial allocMonomial() const { return allocMonomial(mMonomialPool); }
mNumVars * sizeof(exponent)
);
setWeightsAndHash(m);
+#endif
}
exponent monomialExponent(ConstMonomial m, size_t var) const {
+#ifdef MATHICGB_USE_MONOID
+ return monoid().exponent(m, var);
+#else
return m[var+1];
+#endif
}
// This function only sets component and the monomial itself. NOT weights, degree, or hash value
@@ -443,8 +459,6 @@ monomial allocMonomial() const { return allocMonomial(mMonomialPool); }
inline void setHashOnly(Monomial& a) const;
- void setGivenHash(Monomial& a, HashValue hash) const {a[mHashIndex] = static_cast<exponent>(hash);}
-
bool hashValid(const_monomial m) const;
bool weightsCorrect(ConstMonomial a) const;
@@ -463,6 +477,9 @@ monomial allocMonomial() const { return allocMonomial(mMonomialPool); }
bool monomialHasAmpleCapacity(ConstMonomial mono) const;
bool monomialLT(ConstMonomial a, ConstMonomial b) const {
+#ifdef MATHICGB_USE_MONOID
+ return monoid().lessThan(a, b);
+#else
for (size_t i = mTopIndex; i != static_cast<size_t>(-1); --i)
{
exponent cmp = a[i] - b[i];
@@ -471,6 +488,7 @@ monomial allocMonomial() const { return allocMonomial(mMonomialPool); }
return true;
}
return false;
+#endif
}
bool monomialEQ(ConstMonomial a, ConstMonomial b) const;
@@ -483,9 +501,13 @@ monomial allocMonomial() const { return allocMonomial(mMonomialPool); }
exponent monomialGetComponent(ConstMonomial a) const { return *a.mValue; }
void monomialChangeComponent(Monomial a, int x) const {
+#ifdef MATHICGB_USE_MONOID
+ monoid().setComponent(x, a);
+#else
a[mHashIndex] -= static_cast<HashValue>(*a.mValue);
a[mHashIndex] += static_cast<HashValue>(x);
*a = x;
+#endif
}
void monomialSetIdentity(Monomial& result) const;
@@ -496,8 +518,8 @@ monomial allocMonomial() const { return allocMonomial(mMonomialPool); }
void monomialMultTo(Monomial &a, ConstMonomial b) const; // a *= b
- /// returns true if b divides a, in this case, result is set to b//a.
- bool monomialDivide(ConstMonomial a, ConstMonomial b, Monomial &result) const;
+ /// Result is set to b/a. a must divide b.
+ void monomialDivide(ConstMonomial a, ConstMonomial b, Monomial &result) const;
/// sets result to a/b, even if that produces negative exponents.
void monomialDivideToNegative(ConstMonomial a, ConstMonomial b, Monomial &result) const;
@@ -630,11 +652,19 @@ private:
const Monoid& monoid() const {return mMonoid;}
Monoid mMonoid;
#endif
+#ifdef MATHICGB_USE_FIELD
+ const Field field() const {return mField;}
+ Field mField;
+#endif
};
inline exponent PolyRing::weight(ConstMonomial a) const {
+#ifdef MATHICGB_USE_MONOID
+ monoid().degree(a);
+#else
MATHICGB_ASSERT(weightsCorrect(a));
return a[mNumVars + 1];
+#endif
}
////////////////////////////////////////////////
@@ -857,25 +887,27 @@ inline bool PolyRing::monomialIsDivisibleBy(ConstMonomial a,
#endif
}
-inline bool PolyRing::monomialDivide(ConstMonomial a,
+inline void PolyRing::monomialDivide(ConstMonomial a,
ConstMonomial b,
Monomial& result) const
{
+#ifdef MATHICGB_USE_MONOID
+ return monoid().divide(b, a, result);
+#else
//// returns true if b divides a, in this case, result is set to b//a.
size_t i;
for (i = 1; i <= mNumVars; i++)
{
exponent c = a[i] - b[i];
- if (c >= 0)
- result[i] = c;
- else
- return false;
+ if (c < 0)
+ return;
+ result[i] = c;
}
// at this point we have divisibility, so need to fill in the rest of the monomial
*result = *a.mValue - *b.mValue; // component
for ( ; i<=mHashIndex; i++)
result[i] = a[i] - b[i];
- return true;
+#endif
}
inline void PolyRing::monomialColons(
@@ -884,6 +916,9 @@ inline void PolyRing::monomialColons(
monomial aColonB,
monomial bColonA
) const {
+#ifdef MATHICGB_USE_MONOID
+ monoid().colons(a, b, aColonB, bColonA);
+#else
*aColonB = *a;
*bColonA = *b;
for (size_t i = 1; i <= mNumVars; i++) {
@@ -893,12 +928,16 @@ inline void PolyRing::monomialColons(
}
setWeightsAndHash(aColonB);
setWeightsAndHash(bColonA);
+#endif
}
inline void PolyRing::monomialDivideToNegative(ConstMonomial a,
ConstMonomial b,
Monomial& result) const
{
+#ifdef MATHICGB_USE_MONOID
+ monoid().divideToNegative(b, a, result);
+#else
for (size_t i = 0; i <= mHashIndex; ++i)
result[i] = a[i] - b[i];
MATHICGB_ASSERT(monomialHashValue(result) ==
@@ -906,15 +945,20 @@ inline void PolyRing::monomialDivideToNegative(ConstMonomial a,
MATHICGB_ASSERT(!hashValid(a) || !hashValid(b) || hashValid(result));
MATHICGB_ASSERT(computeHashValue(result) == static_cast<exponent>
(computeHashValue(a) - computeHashValue(b)));
+#endif
}
inline bool PolyRing::monomialRelativelyPrime(ConstMonomial a,
ConstMonomial b) const
{
+#ifdef MATHICGB_USE_MONOID
+ return monoid().relativelyPrime(a, b);
+#else
for (size_t i = 1; i <= mNumVars; ++i)
if (a[i] > 0 && b[i] > 0)
return false;
return true;
+#endif
}
inline void PolyRing::monomialLeastCommonMultiple(
@@ -922,8 +966,12 @@ inline void PolyRing::monomialLeastCommonMultiple(
ConstMonomial b,
Monomial& l) const
{
+#ifdef MATHICGB_USE_MONOID
+ monoid().lcm(a, b, l);
+#else
monomialLeastCommonMultipleNoWeights(a, b, l);
setWeightsAndHash(l);
+#endif
}
inline void PolyRing::monomialLeastCommonMultipleNoWeights(
@@ -931,9 +979,13 @@ inline void PolyRing::monomialLeastCommonMultipleNoWeights(
ConstMonomial b,
Monomial& l) const
{
+#ifdef MATHICGB_USE_MONOID
+ monoid().lcm(a, b, l);
+#else
*l = 0;
for (size_t i = 1; i <= mNumVars; ++i)
l[i] = std::max(a[i], b[i]);
+#endif
}
inline bool PolyRing::monomialHasStrictlyLargerExponent(
@@ -941,10 +993,14 @@ inline bool PolyRing::monomialHasStrictlyLargerExponent(
ConstMonomial smaller1,
ConstMonomial smaller2) const
{
+#ifdef MATHICGB_USE_MONOID
+ return !monoid().dividesLcm(hasLarger, smaller1, smaller2);
+#else
for (size_t i = 1; i <= mNumVars; ++i)
if (hasLarger[i] > smaller1[i] && hasLarger[i] > smaller2[i])
return true;
return false;
+#endif
}
@@ -957,45 +1013,73 @@ inline bool PolyRing::monomialIsLeastCommonMultiple(
ConstMonomial b,
ConstMonomial l) const
{
+#ifdef MATHICGB_USE_MONOID
+ return monoid().isLcm(a, b, l);
+#else
return monomialIsLeastCommonMultipleNoWeights(a, b, l) && weightsCorrect(l);
+#endif
}
inline void PolyRing::coefficientReciprocalTo(coefficient& result) const
{
+#ifdef MATHICGB_USE_FIELD
+ result = field().inverse(field().toElementInRange(result)).value();
+#else
MATHICGB_ASSERT(result != 0);
mStats.n_recip++;
result = modularInverse(result, mCharac);
+#endif
}
inline void PolyRing::coefficientDivide(coefficient a, coefficient b, coefficient &result) const
// result = a/b
{
+#ifdef MATHICGB_USE_FIELD
+ result = field().quotient
+ (field().toElementInRange(a), field().toElementInRange(b)).value();
+#else
mStats.n_divide++;
result = (a * modularInverse(b, mCharac)) % mCharac;
MATHICGB_ASSERT((result * b) % mCharac == a);
MATHICGB_ASSERT(result >= 0);
MATHICGB_ASSERT(result < mCharac);
+#endif
}
inline void PolyRing::coefficientFromInt(coefficient &result, int a) const
{
+#ifdef MATHICGB_USE_FIELD
+ result = field().toElement(a).value();
+#else
result = toCoefficient(a);
+#endif
}
inline void PolyRing::coefficientAddOneTo(coefficient &result) const
{
+#ifdef MATHICGB_USE_FIELD
+ result = field().plusOne(field().toElementInRange(result)).value();
+#else
++result;
if (result == mCharac)
result = 0;
+#endif
}
inline void PolyRing::coefficientNegateTo(coefficient& result) const {
+#ifdef MATHICGB_USE_FIELD
+ result = field().negative(field().toElementInRange(result)).value();
+#else
MATHICGB_ASSERT(result < mCharac);
if (result != 0)
result = coefficientNegateNonZero(result);
+#endif
}
inline coefficient PolyRing::toCoefficient(const int64 value) const {
+#ifdef MATHICGB_USE_FIELD
+ return field().toElement(value).value();
+#else
auto modLong = value % mCharac;
if (modLong < 0)
modLong += mCharac;
@@ -1005,71 +1089,113 @@ inline coefficient PolyRing::toCoefficient(const int64 value) const {
MATHICGB_ASSERT(0 <= mod);
MATHICGB_ASSERT(mod < mCharac);
return mod;
+#endif
}
inline coefficient PolyRing::coefficientNegate(const coefficient coeff) const {
+#ifdef MATHICGB_USE_FIELD
+ return field().negative(field().toElementInRange(coeff)).value();
+#else
MATHICGB_ASSERT(coeff < mCharac);
return coeff == 0 ? 0 : coefficientNegateNonZero(coeff);
+#endif
}
inline coefficient PolyRing::coefficientNegateNonZero(
const coefficient coeff
) const {
+#ifdef MATHICGB_USE_FIELD
+ return field().negativeNonZero(field().toElementInRange(coeff)).value();
+#else
MATHICGB_ASSERT(coeff != 0);
MATHICGB_ASSERT(coeff < mCharac);
return mCharac - coeff;
+#endif
}
inline coefficient PolyRing::coefficientSubtract(
const coefficient a,
const coefficient b
) const {
+#ifdef MATHICGB_USE_FIELD
+ return field().difference
+ (field().toElementInRange(a), field().toElementInRange(b)).value();
+#else
MATHICGB_ASSERT(a < mCharac);
MATHICGB_ASSERT(b < mCharac);
const auto diff = a < b ? a + (mCharac - b) : a - b;
MATHICGB_ASSERT(diff < mCharac);
MATHICGB_ASSERT((diff + b) % mCharac == a);
return diff;
+#endif
}
-inline void PolyRing::coefficientAddTo(coefficient &result, coefficient a, coefficient b) const
+inline void PolyRing::coefficientAddTo
+(coefficient &result, coefficient a, coefficient b) const
// result += a*b
{
+#ifdef MATHICGB_USE_FIELD
+ const auto prod =
+ field().product(field().toElementInRange(a), field().toElementInRange(b));
+ result = field().sum(field().toElementInRange(result), prod).value();
+#else
mStats.n_addmult++;
auto c = a * b + result;
result = c % mCharac;
+#endif
}
inline void PolyRing::coefficientAddTo(coefficient &result, coefficient a) const
// result += a
{
+#ifdef MATHICGB_USE_FIELD
+ result = field().sum
+ (field().toElementInRange(result), field().toElementInRange(a)).value();
+#else
mStats.n_add++;
result += a;
if (result >= mCharac)
result -= mCharac;
+#endif
}
-inline void PolyRing::coefficientMultTo(coefficient &result, coefficient a) const
+inline void PolyRing::coefficientMultTo
+(coefficient &result, coefficient a) const
// result *= a
{
+#ifdef MATHICGB_USE_FIELD
+ result = field().product
+ (field().toElementInRange(result), field().toElementInRange(a)).value();
+#else
mStats.n_mult++;
coefficient b = result * a;
result = b % mCharac;
+#endif
}
-inline void PolyRing::coefficientMult(coefficient a, coefficient b, coefficient &result) const
+inline void PolyRing::coefficientMult
+(coefficient a, coefficient b, coefficient &result) const
{
+#ifdef MATHICGB_USE_FIELD
+ result = field().product
+ (field().toElementInRange(a), field().toElementInRange(b)).value();
+#else
mStats.n_mult++;
coefficient c = b * a;
result = c % mCharac;
+#endif
}
inline bool PolyRing::monomialHasAmpleCapacity(ConstMonomial mono) const {
+#ifdef MATHICGB_USE_MONOID
+ return monoid().hasAmpleCapacity(mono);
+#else
const auto halfMax = std::numeric_limits<exponent>::max() / 2;
for (size_t i = mTopIndex; i != 0; --i)
if (mono[i] > halfMax)
return false;
return true;
+#endif
}
#endif
diff --git a/src/mathicgb/PrimeField.hpp b/src/mathicgb/PrimeField.hpp
index a1f8cd9..a981b4b 100755
--- a/src/mathicgb/PrimeField.hpp
+++ b/src/mathicgb/PrimeField.hpp
@@ -95,6 +95,7 @@ public:
// 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.
+ // @todo: extend precision to rule out overflow without a branch.
if (a.value() > s || s >= charac()) {
MATHICGB_ASSERT(s - charac() < charac());
return Element(s - charac());
@@ -102,6 +103,14 @@ public:
return Element(s);
}
+ Element plusOne(const Element a) const {
+ const auto s = a.value() + 1;
+ if (s == charac())
+ return Element(0);
+ 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());
@@ -124,6 +133,11 @@ public:
Element product(const Element a, const Element b) const;
+ /// Returns a times the inverse of b.
+ Element quotient(const Element a, const Element b) const {
+ return product(a, inverse(b));
+ }
+
/// Returns the multiplicative inverse a^-1 mod charac(). a must not be zero.
Element inverse(const Element a) const;
diff --git a/src/mathicgb/SignatureGB.cpp b/src/mathicgb/SignatureGB.cpp
index 15147b0..31bfdc9 100755
--- a/src/mathicgb/SignatureGB.cpp
+++ b/src/mathicgb/SignatureGB.cpp
@@ -158,7 +158,7 @@ bool SignatureGB::processSPair
}
// new basis element
- MATHICGB_ASSERT(!GB->isSingularTopReducible(*f, sig));
+ MATHICGB_ASSERT(!GB->isSingularTopReducibleSlow(*f, sig));
{
std::unique_ptr<Poly> autoF(f);
GB->insert(sig, std::move(autoF));
diff --git a/src/mathicgb/io-util.cpp b/src/mathicgb/io-util.cpp
index 1e30596..514ed1c 100755
--- a/src/mathicgb/io-util.cpp
+++ b/src/mathicgb/io-util.cpp
@@ -47,7 +47,7 @@ std::unique_ptr<PolyRing> ringFromString(std::string ringinfo)
Monomial stringToMonomial(const PolyRing *R, std::string mon)
{
- Monomial result = R->allocMonomial1();
+ Monomial result = R->allocMonomial();
std::stringstream ifil(mon);
R->monomialParse(ifil, result);
return result;
diff --git a/src/test/poly-test.cpp b/src/test/poly-test.cpp
index f92fc3c..0795583 100755
--- a/src/test/poly-test.cpp
+++ b/src/test/poly-test.cpp
@@ -155,7 +155,7 @@ TEST(Monomial,mult) {
Monomial m2 = stringToMonomial(R.get(), "a2b<0>");
Monomial m3ans = stringToMonomial(R.get(), "a3b3<0>");
- Monomial m3 = R->allocMonomial1();
+ Monomial m3 = R->allocMonomial();
R->monomialMult(m1,m2,m3);
EXPECT_TRUE(R->monomialEQ(m3ans,m3));
@@ -186,8 +186,8 @@ TEST(Monomial, divide) {
Monomial m1 = stringToMonomial(R.get(), "ab2<0>");
Monomial m2 = stringToMonomial(R.get(), "a2b<0>");
Monomial m3ans = stringToMonomial(R.get(), "a3b3<0>");
- Monomial m3 = R->allocMonomial1();
- Monomial m1a = R->allocMonomial1();
+ Monomial m3 = R->allocMonomial();
+ Monomial m1a = R->allocMonomial();
R->monomialMult(m1,m2,m3);
EXPECT_TRUE(R->monomialIsDivisibleBy(m3,m2));
@@ -210,9 +210,9 @@ TEST(Monomial, monomialQuotientAndMult) {
Monomial m2 = stringToMonomial(R.get(), "af<0>");
Monomial m3 = stringToMonomial(R.get(), "<2>");
- Monomial n = R->allocMonomial1();
- Monomial n1 = R->allocMonomial1();
- Monomial na = R->allocMonomial1();
+ Monomial n = R->allocMonomial();
+ Monomial n1 = R->allocMonomial();
+ Monomial na = R->allocMonomial();
R->monomialQuotientAndMult(m1,m2,m3,n);
R->monomialDivide(m1,m2,n1); // m1//m2
@@ -236,15 +236,15 @@ void testMonomialOps(const PolyRing* R, std::string s1, std::string s2)
// m1 * m2 == lcm(m1,m2) * gcd(m1,m2)
- Monomial m4 = R->allocMonomial1();
- Monomial lcm = R->allocMonomial1();
- Monomial gcd = R->allocMonomial1();
- Monomial m7 = R->allocMonomial1();
- Monomial m8 = R->allocMonomial1();
- Monomial m1a = R->allocMonomial1();
- Monomial m2a = R->allocMonomial1();
- Monomial m1b = R->allocMonomial1();
- Monomial m2b = R->allocMonomial1();
+ Monomial m4 = R->allocMonomial();
+ Monomial lcm = R->allocMonomial();
+ Monomial gcd = R->allocMonomial();
+ Monomial m7 = R->allocMonomial();
+ Monomial m8 = R->allocMonomial();
+ Monomial m1a = R->allocMonomial();
+ Monomial m2a = R->allocMonomial();
+ Monomial m1b = R->allocMonomial();
+ Monomial m2b = R->allocMonomial();
R->monomialMult(m1,m2,m4);
R->monomialLeastCommonMultiple(m1,m2,lcm);
@@ -254,15 +254,21 @@ void testMonomialOps(const PolyRing* R, std::string s1, std::string s2)
EXPECT_TRUE(R->monomialEQ(m4, m7));
// lcm(m1,m2)/m1, lcm(m1,m2)/m2: relatively prime
- EXPECT_TRUE(R->monomialDivide(lcm, m1, m1a));
- EXPECT_TRUE(R->monomialDivide(lcm, m2, m2a));
+ EXPECT_TRUE(R->monomialIsDivisibleBy(lcm, m1));
+ EXPECT_TRUE(R->monomialIsDivisibleBy(lcm, m2));
+ R->monomialDivide(lcm, m1, m1a);
+ R->monomialDivide(lcm, m2, m2a);
EXPECT_TRUE(R->monomialRelativelyPrime(m1a,m2a));
- EXPECT_TRUE(R->monomialDivide(lcm, m1a, m1b));
- EXPECT_TRUE(R->monomialDivide(lcm, m2a, m2b));
+ EXPECT_TRUE(R->monomialIsDivisibleBy(lcm, m1a));
+ EXPECT_TRUE(R->monomialIsDivisibleBy(lcm, m2a));
+ R->monomialDivide(lcm, m1a, m1b);
+ R->monomialDivide(lcm, m2a, m2b);
EXPECT_TRUE(R->monomialEQ(m1, m1b));
EXPECT_TRUE(R->monomialEQ(m2, m2b));
- EXPECT_TRUE(R->monomialDivide(lcm,gcd,m8));
+
+ EXPECT_TRUE(R->monomialIsDivisibleBy(lcm,gcd));
+ R->monomialDivide(lcm,gcd,m8);
size_t supp1 = R->monomialSizeOfSupport(m1a);
size_t supp2 = R->monomialSizeOfSupport(m2a);
@@ -290,7 +296,7 @@ TEST(Monomial, ei)
std::unique_ptr<PolyRing> R(ringFromString("32003 6 1\n1 1 1 1 1 1"));
Monomial m1 = stringToMonomial(R.get(), "<1>");
- Monomial m1a = R->allocMonomial1();
+ Monomial m1a = R->allocMonomial();
R->monomialEi(1, m1a);
EXPECT_TRUE(R->monomialEQ(m1,m1a));
@@ -327,9 +333,9 @@ TEST(Monomial, divideToNegative)
Monomial m1 = stringToMonomial(R.get(), "ab100<0>");
Monomial m2 = stringToMonomial(R.get(), "ab2c3d4<0>");
- Monomial m3 = R->allocMonomial1();
- Monomial m4 = R->allocMonomial1();
- Monomial m5 = R->allocMonomial1();
+ Monomial m3 = R->allocMonomial();
+ Monomial m4 = R->allocMonomial();
+ Monomial m5 = R->allocMonomial();
Monomial mone = stringToMonomial(R.get(), "<0>");
R->monomialDivideToNegative(m1,m2,m3);
@@ -352,7 +358,7 @@ TEST(Monomial, findSignature)
Monomial v1 = stringToMonomial(R.get(), "abef");
Monomial v2 = stringToMonomial(R.get(), "acdf2");
Monomial u1 = stringToMonomial(R.get(), "f5<13>");
- Monomial t1 = R->allocMonomial1();
+ Monomial t1 = R->allocMonomial();
Monomial t1ans = stringToMonomial(R.get(), "cdf6<13>");
R->monomialFindSignature(v1,v2,u1,t1);
--
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