[mathicgb] 70/393: Now uses exponent and coefficient typedefs in more places where they should be used. Also added OMPIndex typedef for use in OpenMP loops where the loop variable is (perversely) not allowed to be unsigned by MSVC 2012. I tried to make it possible to have coefficient be unsigned, but there is some bug somewhere, so for now it remains signed.
Doug Torrance
dtorrance-guest at moszumanska.debian.org
Fri Apr 3 15:58:34 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 433cf6e5f62e3e85062d6fd83346693fd7f6c88d
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date: Tue Oct 23 16:15:19 2012 +0200
Now uses exponent and coefficient typedefs in more places where they should be used. Also added OMPIndex typedef for use in OpenMP loops where the loop variable is (perversely) not allowed to be unsigned by MSVC 2012. I tried to make it possible to have coefficient be unsigned, but there is some bug somewhere, so for now it remains signed.
---
src/mathicgb/F4MatrixReducer.cpp | 21 +++++--
src/mathicgb/Ideal.cpp | 16 +++---
src/mathicgb/MonomialMap.hpp | 1 -
src/mathicgb/Poly.cpp | 108 +++++++++++++++++------------------
src/mathicgb/Poly.hpp | 10 ++--
src/mathicgb/PolyRing.cpp | 76 ++++++++++---------------
src/mathicgb/PolyRing.hpp | 120 ++++++++++++++++++++++++---------------
src/mathicgb/stdinc.h | 6 ++
8 files changed, 188 insertions(+), 170 deletions(-)
diff --git a/src/mathicgb/F4MatrixReducer.cpp b/src/mathicgb/F4MatrixReducer.cpp
index 38fecab..3279dd1 100755
--- a/src/mathicgb/F4MatrixReducer.cpp
+++ b/src/mathicgb/F4MatrixReducer.cpp
@@ -329,7 +329,9 @@ void myReduce
std::vector<SparseMatrix::RowIndex> rowOrder(rowCount);
#pragma omp parallel for num_threads(threadCount) schedule(dynamic)
- for (long row = 0; row < rowCount; ++row) {
+ for (OMPIndex rowOMP = 0;
+ rowOMP < static_cast<OMPIndex>(rowCount); ++rowOMP) {
+ const size_t row = static_cast<size_t>(rowOMP);
#ifdef _OPENMP
auto& denseRow = denseRowPerThread[omp_get_thread_num()];
#endif
@@ -367,9 +369,9 @@ void myReduce
}
}
-
#pragma omp parallel for num_threads(threadCount) schedule(dynamic)
- for (long i = 0; i < rowCount; ++i) {
+ for (OMPIndex iOMP = 0; iOMP < static_cast<OMPIndex>(rowCount); ++iOMP) {
+ const size_t i = static_cast<size_t>(iOMP);
#ifdef _OPENMP
auto& denseRow = denseRowPerThread[omp_get_thread_num()];
#endif
@@ -435,7 +437,9 @@ void myReduceToEchelonForm5
// dense representation
std::vector<DenseRow<uint64> > dense(rowCount);
#pragma omp parallel for num_threads(threadCount) schedule(dynamic)
- for (long row = 0; row < rowCount; ++row) {
+ for (OMPIndex rowOMP = 0;
+ rowOMP < static_cast<OMPIndex>(rowCount); ++rowOMP) {
+ const size_t row = static_cast<size_t>(rowOMP);
MATHICGB_ASSERT(!toReduce.emptyRow(row));
dense[row].reset(colCount);
dense[row].addRow(toReduce, row);
@@ -466,7 +470,9 @@ void myReduceToEchelonForm5
//std::cout << "reducing " << reduced.rowCount() << " out of " << toReduce.rowCount() << std::endl;
#pragma omp parallel for num_threads(threadCount) schedule(dynamic)
- for (long row = 0; row < rowCount; ++row) {
+ for (OMPIndex rowOMP = 0;
+ rowOMP < static_cast<OMPIndex>(rowCount); ++rowOMP) {
+ const size_t row = static_cast<size_t>(rowOMP);
MATHICGB_ASSERT(leadCols[row] <= colCount);
DenseRow<uint64>& denseRow = dense[row];
if (denseRow.empty())
@@ -534,8 +540,11 @@ void myReduceToEchelonForm5
}
#pragma omp parallel for num_threads(threadCount) schedule(dynamic)
- for (long row = 0; row < rowCount; ++row)
+ for (OMPIndex rowOMP = 0;
+ rowOMP < static_cast<OMPIndex>(rowCount); ++rowOMP) {
+ const size_t row = static_cast<size_t>(rowOMP);
dense[row].takeModulus(modulus);
+ }
toReduce.clear(colCount);
for (size_t row = 0; row < rowCount; ++row)
diff --git a/src/mathicgb/Ideal.cpp b/src/mathicgb/Ideal.cpp
index 2f868da..09163ab 100755
--- a/src/mathicgb/Ideal.cpp
+++ b/src/mathicgb/Ideal.cpp
@@ -58,16 +58,14 @@ std::unique_ptr<Ideal> Ideal::parse(std::istream& in)
return result;
}
-void Ideal::display(std::ostream &o, bool print_comp) const
+void Ideal::display(std::ostream& out, bool printComponent) const
{
- mRing.write(o);
- o << std::endl;
- o << mGenerators.size() << std::endl;
- for (size_t i = 0; i<mGenerators.size(); i++)
- {
- mGenerators[i]->display(o, print_comp);
- o << std::endl;
- }
+ mRing.write(out);
+ out << '\n' << mGenerators.size() << '\n';
+ for (size_t i = 0; i < mGenerators.size(); ++i) {
+ mGenerators[i]->display(out, printComponent);
+ out << '\n';
+ }
}
// Local Variables:
diff --git a/src/mathicgb/MonomialMap.hpp b/src/mathicgb/MonomialMap.hpp
index cd43d16..6b22b2f 100755
--- a/src/mathicgb/MonomialMap.hpp
+++ b/src/mathicgb/MonomialMap.hpp
@@ -44,7 +44,6 @@ namespace MonomialMapInternal {
typedef MTT mapped_type;
private:
- typedef exponent HashValue;
struct Node {
Node* next;
mapped_type value;
diff --git a/src/mathicgb/Poly.cpp b/src/mathicgb/Poly.cpp
index ca6f86c..ca4f30e 100755
--- a/src/mathicgb/Poly.cpp
+++ b/src/mathicgb/Poly.cpp
@@ -243,47 +243,49 @@ void Poly::dump() const
void Poly::parseDoNotOrder(std::istream& i)
{
- char next = i.peek();
- if (next == '0') {
+ if (i.peek() == '0') {
i.get();
return;
}
- for (;;) {
- bool is_neg = false;
- next = i.peek();
+ while (true) {
+ bool preceededByMinus = false;
+ char next = i.peek();
if (next == '+') {
i.get();
next = i.peek();
} else if (next == '-') {
- is_neg = true;
- i.get();
- next = i.peek();
- } if (isdigit(next) || isalpha(next) || next == '<') {
- int a = 1;
- coefficient b;
- if (isdigit(next))
- {
- i >> a;
- next = i.peek();
- }
- if (is_neg) a = -a;
- size_t firstloc = monoms.size();
- monoms.resize(firstloc + R->maxMonomialSize());
- monomial m = &monoms[firstloc];
- R->coefficientFromInt(b,a);
- coeffs.push_back(b);
- if (isalpha(next) || next == '<')
- R->monomialParse(i, m);
- else
- R->monomialSetIdentity(m); // have to do this to set hash value
- MATHICGB_ASSERT(ring().hashValid(m));
+ preceededByMinus = true;
+ i.get();
+ next = i.peek();
+ }
+
+ if (!isdigit(next) && !isalpha(next) && next != '<')
+ break;
+
+ { // read coefficient
+ int64 bigCoef = 1;
+ if (isdigit(next)) {
+ i >> bigCoef;
next = i.peek();
- if (next == '>')
- i.get();
}
+ if (preceededByMinus)
+ bigCoef = -bigCoef;
+ coeffs.push_back(R->toCoefficient(bigCoef));
+ }
+
+ // read monic monomial
+ const size_t firstLocation = monoms.size();
+ monoms.resize(firstLocation + R->maxMonomialSize());
+ monomial m = &monoms[firstLocation];
+ if (isalpha(next) || next == '<')
+ R->monomialParse(i, m);
else
- break;
+ R->monomialSetIdentity(m); // have to do this to set hash value
+ MATHICGB_ASSERT(ring().hashValid(m));
+ next = i.peek();
+ if (next == '>')
+ i.get();
}
}
@@ -292,30 +294,26 @@ void Poly::parse(std::istream& in) {
sortTermsDescending();
}
-void Poly::display(std::ostream &o, bool print_comp) const
+void Poly::display(std::ostream& out, const bool printComponent) const
{
- long p = R->charac();
- int maxpos = (p+1)/2;
- if (nTerms() == 0)
- {
- o << "0";
- return;
- }
- int nterms = 0;
- for (const_iterator i = begin(); i != end(); ++i)
- {
- nterms++;
- coefficient a = i.getCoefficient();
- bool is_neg = (a > maxpos);
- if (is_neg) a = p-a;
- if (is_neg)
- o << "-";
- else if (nterms > 1)
- o << "+";
- bool print_one = (a == 1);
- if (a != 1) o << a;
- R->monomialDisplay(o, i.getMonomial(), print_comp, print_one);
- }
+ const coefficient p = R->charac();
+ const coefficient maxPositive = (p + 1) / 2; // half rounded up
+ if (isZero()) {
+ out << "0";
+ return;
+ }
+
+ for (const_iterator i = begin(); i != end(); ++i) {
+ coefficient coef = i.getCoefficient();
+ if (coef > maxPositive) {
+ out << "-";
+ R->coefficientNegateTo(coef);
+ } else if (i != begin())
+ out << '+';
+ if (coef != 1)
+ out << coef;
+ R->monomialDisplay(out, i.getMonomial(), printComponent, coef == 1);
+ }
}
void Poly::display(FILE* file, bool printComponent) const
@@ -326,11 +324,11 @@ void Poly::display(FILE* file, bool printComponent) const
}
const auto characteristic = R->charac();
- const exponent maxPositiveExponent = (characteristic + 1) / 2;
+ const coefficient maxPositiveCoefficient = (characteristic + 1) / 2;
bool firstTerm = true;
for (auto it = begin(); it != end(); ++it) {
coefficient coef = it.getCoefficient();
- if (coef > maxPositiveExponent) {
+ if (coef > maxPositiveCoefficient) {
coef = characteristic - coef;
fputc('-', file);
} else if (!firstTerm)
diff --git a/src/mathicgb/Poly.hpp b/src/mathicgb/Poly.hpp
index 81c1b7b..d3be160 100755
--- a/src/mathicgb/Poly.hpp
+++ b/src/mathicgb/Poly.hpp
@@ -16,14 +16,14 @@ public:
void parse(std::istream &i); // reads into this, sorts terms
void parseDoNotOrder(std::istream &i); // reads into this, does not sort terms
void display(FILE* file, bool printComponent = true) const;
- void display(std::ostream &o, bool print_comp = true) const;
+ void display(std::ostream& out, bool printComponent = true) const;
void see(bool print_comp) const;
class iterator {
// only for const objects...
size_t monsize;
std::vector<coefficient>::iterator ic;
- std::vector<int>::iterator im;
+ std::vector<exponent>::iterator im;
friend class Poly;
iterator(Poly& f) : monsize(f.getRing()->maxMonomialSize()), ic(f.coeffs.begin()), im(f.monoms.begin()) {}
@@ -45,7 +45,7 @@ public:
// only for const objects...
size_t monsize;
std::vector<coefficient>::const_iterator ic;
- std::vector<int>::const_iterator im;
+ std::vector<exponent>::const_iterator im;
friend class Poly;
const_iterator(const Poly& f) : monsize(f.getRing()->maxMonomialSize()), ic(f.coeffs.begin()), im(f.monoms.begin()) {}
@@ -109,7 +109,7 @@ public:
const_monomial getLeadMonomial() const { return &(monoms[0]); }
const_coefficient getLeadCoefficient() const { return coeffs[0]; }
- long getLeadComponent() const { return R->monomialGetComponent(&(monoms[0])); }
+ exponent getLeadComponent() const { return R->monomialGetComponent(&(monoms[0])); }
bool isZero() const { return coeffs.empty(); }
size_t nTerms() const { return coeffs.size(); }
@@ -133,7 +133,7 @@ public:
private:
const PolyRing *R;
std::vector<coefficient> coeffs;
- std::vector<int> monoms;
+ std::vector<exponent> monoms;
};
std::ostream& operator<<(std::ostream& out, const Poly& p);
diff --git a/src/mathicgb/PolyRing.cpp b/src/mathicgb/PolyRing.cpp
index d402d4c..b6eecda 100755
--- a/src/mathicgb/PolyRing.cpp
+++ b/src/mathicgb/PolyRing.cpp
@@ -15,7 +15,7 @@ bool PolyRing::hashValid(const_monomial m) const {
return monomialHashValue(m) == computeHashValue(m);
}
-PolyRing::PolyRing(long p0,
+PolyRing::PolyRing(coefficient p0,
int nvars,
int nweights)
: mCharac(p0),
@@ -29,9 +29,8 @@ PolyRing::PolyRing(long p0,
mTotalDegreeGradedOnly(false)
{
resetCoefficientStats();
- // rand();
for (size_t i=0; i<mNumVars; i++)
- mHashVals.push_back(rand());
+ mHashVals.push_back(static_cast<HashValue>(rand()));
}
void PolyRing::displayHashValues() const
@@ -64,10 +63,10 @@ bool PolyRing::weightsCorrect(ConstMonomial a1) const
{
exponent const *a = a1.mValue;
++a;
- const int *wts = &mWeights[0];
- for (int i=0; i < mNumWeights; ++i) {
- int result = 0;
- for (size_t j=0; j < mNumVars; ++j)
+ auto wts = &mWeights[0];
+ for (size_t i = 0; i < mNumWeights; ++i) {
+ exponent result = 0;
+ for (size_t j = 0; j < mNumVars; ++j)
result += *wts++ * a[j];
if (a[mNumVars + i] != result)
return false;
@@ -80,7 +79,7 @@ int PolyRing::monomialCompare(ConstMonomial sig, ConstMonomial m2, ConstMonomial
{
for (size_t i = mTopIndex; i != static_cast<size_t>(-1); --i)
{
- int cmp = sig[i] - m2[i] - sig2[i];
+ auto cmp = sig[i] - m2[i] - sig2[i];
if (cmp < 0) return GT;
if (cmp > 0) return LT;
}
@@ -91,15 +90,15 @@ void PolyRing::monomialSetIdentity(Monomial& result) const
{
for (size_t i = 0; i <= mTopIndex; ++i)
result[i] = 0;
- result[mHashIndex] = 0;
+ result[mHashIndex] = static_cast<HashValue>(0);
}
void PolyRing::monomialEi(size_t i, Monomial &result) const
{
for (size_t j=mTopIndex; j != static_cast<size_t>(-1); --j)
result[j] = 0;
- *result = static_cast<int>(i); // todo: handle overflow or change representation
- result[mHashIndex] = static_cast<int>(i); // todo: handle overflow or change representation
+ *result = static_cast<exponent>(i); // todo: handle overflow or change representation
+ result[mHashIndex] = static_cast<HashValue>(i); // todo: handle overflow or change representation
}
void PolyRing::monomialMultTo(Monomial &a, ConstMonomial b) const
@@ -126,7 +125,7 @@ void PolyRing::monomialQuotientAndMult(ConstMonomial a,
for (size_t i = 0; i <= mNumVars; ++i)
{
result[i] = c[i];
- int cmp = a[i] - b[i];
+ auto cmp = a[i] - b[i];
if (cmp > 0) result[i] += cmp;
}
setWeightsAndHash(result);
@@ -225,7 +224,8 @@ void PolyRing::monomialParse(std::istream &i,
// first initialize result:
for (size_t j=0; j<mMaxMonomialSize; j++) result[j] = 0;
- int v, e, x;
+ exponent e;
+ int v, x;
// now look at the next char
for (;;)
{
@@ -321,7 +321,7 @@ void PolyRing::printMonomialFrobbyM2Format(std::ostream& out, ConstMonomial m) c
out << " ";
bool isOne = true;
for (size_t i = 0; i < mNumVars; ++i) {
- int e = m[i + 1];
+ const auto e = m[i + 1];
if (e == 0)
continue;
if (!isOne)
@@ -359,9 +359,9 @@ bool PolyRing::monomialIsLeastCommonMultiple(
bool PolyRing::weightsCorrect(const_monomial a) const
{
++a;
- const int *wts = &mWeights[0];
+ auto wts = &mWeights[0];
for (int i=0; i < mNumWeights; ++i) {
- int result = 0;
+ exponent result = 0;
for (size_t j=0; j < mNumVars; ++j)
result += *wts++ * a[j];
if (a[mNumVars + i] != result)
@@ -388,7 +388,7 @@ int PolyRing::monomialCompare(const_monomial sig, const_monomial m2, const_monom
{
for (size_t i = mTopIndex; i != static_cast<size_t>(-1); --i)
{
- int cmp = sig[i] - m2[i] - sig2[i];
+ auto cmp = sig[i] - m2[i] - sig2[i];
if (cmp < 0) return GT;
if (cmp > 0) return LT;
}
@@ -411,8 +411,8 @@ void PolyRing::monomialSetIdentity(monomial& result) const {
void PolyRing::monomialEi(size_t i, monomial &result) const
{
for (size_t j=mTopIndex; j != static_cast<size_t>(-1); --j) result[j] = 0;
- *result = static_cast<int>(i); // todo: handle overflow or change representation
- result[mHashIndex] = static_cast<int>(i); // todo: handle overflow or change representation
+ *result = static_cast<exponent>(i); // todo: handle overflow or change representation
+ result[mHashIndex] = static_cast<exponent>(i); // todo: handle overflow or change representation
}
void PolyRing::monomialMult(const_monomial a, const_monomial b, monomial &result) const
@@ -443,7 +443,7 @@ void PolyRing::monomialQuotientAndMult(const_monomial a,
for (size_t i = 0; i <= mNumVars; ++i)
{
result[i] = c[i];
- int cmp = a[i] - b[i];
+ auto cmp = a[i] - b[i];
if (cmp > 0) result[i] += cmp;
}
setWeightsAndHash(result);
@@ -452,7 +452,7 @@ void PolyRing::monomialQuotientAndMult(const_monomial a,
void PolyRing::setWeightsAndHash(monomial a) const
{
setWeightsOnly(a);
- int hash = *a;
+ auto hash = *a;
a++;
for (size_t i = 0; i < mNumVars; ++i)
hash += a[i] * mHashVals[i];
@@ -542,7 +542,8 @@ void PolyRing::monomialRead(std::istream &I, monomial &result) const
result[i] = 0;
for (int i=len; i>0; --i)
{
- int v,e;
+ int v;
+ exponent e;
I >> v;
I >> e;
result[v+1] = e;
@@ -571,7 +572,8 @@ void PolyRing::monomialParse(std::istream &i, monomial &result) const
// first initialize result:
for (size_t j=0; j<mMaxMonomialSize; j++) result[j] = 0;
- int v, e, x;
+ exponent e;
+ int v, x;
// now look at the next char
for (;;)
{
@@ -636,7 +638,7 @@ void PolyRing::printMonomialFrobbyM2Format(std::ostream& out, const_monomial m)
out << " ";
bool isOne = true;
for (size_t i = 0; i < mNumVars; ++i) {
- int e = m[i + 1];
+ auto e = m[i + 1];
if (e == 0)
continue;
if (!isOne)
@@ -652,7 +654,7 @@ void PolyRing::printMonomialFrobbyM2Format(std::ostream& out, const_monomial m)
PolyRing *PolyRing::read(std::istream &i)
{
- long charac;
+ coefficient charac;
int mNumVars, mNumWeights;
i >> charac;
@@ -664,7 +666,7 @@ PolyRing *PolyRing::read(std::istream &i)
R->mTotalDegreeGradedOnly = (mNumWeights == 1);
for (int j=0; j <mNumVars * mNumWeights; j++)
{
- int a;
+ exponent a;
i >> a;
R->mWeights[j] = -a;
if (R->mWeights[j] != -1)
@@ -676,7 +678,7 @@ PolyRing *PolyRing::read(std::istream &i)
void PolyRing::write(std::ostream &o) const
{
o << mCharac << " " << mNumVars << " " << mNumWeights << std::endl;
- const int *wts = &mWeights[0];
+ auto wts = &mWeights[0];
for (int i=0; i<mNumWeights; i++)
{
for (size_t j=0; j<mNumVars; j++)
@@ -685,26 +687,6 @@ void PolyRing::write(std::ostream &o) const
}
}
-
-void gcd_extended(long a, long b, long &u, long &v, long &g)
-{
- long q ;
- long u1, v1, g1;
- long utemp, vtemp, gtemp;
-
- g1 = b; u1 = 0; v1 = 1;
- g = a; u = 1; v = 0;
- while (g1 != 0)
- {
- q = g / g1 ;
- gtemp= g - q * g1;
- utemp= u - q * u1;
- vtemp= v - q * v1;
- g = g1; u = u1; v = v1 ;
- g1 = gtemp; u1 = utemp; v1 = vtemp;
- }
-}
-
// Local Variables:
// compile-command: "make -C .. "
// indent-tabs-mode: nil
diff --git a/src/mathicgb/PolyRing.hpp b/src/mathicgb/PolyRing.hpp
index 808e3a3..d683988 100755
--- a/src/mathicgb/PolyRing.hpp
+++ b/src/mathicgb/PolyRing.hpp
@@ -28,7 +28,7 @@ T modularInverse(T a, T modulus) {
#ifdef MATHICGB_DEBUG
T origA = a;
#endif
- T b = modulus; // note that we actually only need modulus for asserts
+ T b = modulus;
T minusLastX = 0;
T x = 1;
while (true) {
@@ -97,9 +97,11 @@ T modularNegativeNonZero(T a, T modulus) {
return modulus - a;
}
-typedef int exponent ;
+typedef int32 exponent ;
+typedef uint32 HashValue;
typedef long coefficient;
-typedef exponent *vecmonomial; // includes a component
+
+typedef exponent* vecmonomial; // includes a component
typedef coefficient const_coefficient;
class Monomial;
@@ -200,7 +202,7 @@ struct term {
class PolyRing {
public:
- PolyRing(long charac,
+ PolyRing(coefficient charac,
int nvars,
int nweights
);
@@ -208,7 +210,7 @@ public:
memt::BufferPool &getMonomialPool() const { return mMonomialPool; }
- long charac() const { return mCharac; }
+ coefficient charac() const { return mCharac; }
size_t getNumVars() const { return mNumVars; }
// const std::vector<int> °s,
// const std::string &monorder);
@@ -288,6 +290,9 @@ public:
+ coefficient toCoefficient(int64 a) const;
+ coefficient coefficientNegate(coefficient result) const;
+
void coefficientFromInt(coefficient &result, int a) const;
void coefficientSetOne(coefficient &result) const { result = 1; }
void coefficientAddOneTo(coefficient &result) const;
@@ -322,7 +327,7 @@ public:
// Monomial Routines //////////////////////
///////////////////////////////////////////
- exponent monomialHashValue(ConstMonomial m) const {return m[mHashIndex];}
+ HashValue monomialHashValue(ConstMonomial m) const {return static_cast<HashValue>(m[mHashIndex]);}
exponent monomialExponent(ConstMonomial m, size_t var) const {
return m[var+1];
@@ -346,7 +351,7 @@ public:
inline void setHashOnly(Monomial& a) const;
- void setGivenHash(Monomial& a, exponent hash) const {a[mHashIndex] = hash;}
+ void setGivenHash(Monomial& a, HashValue hash) const {a[mHashIndex] = static_cast<exponent>(hash);}
bool hashValid(const_monomial m) const;
@@ -364,7 +369,7 @@ public:
bool monomialLT(ConstMonomial a, ConstMonomial b) const {
for (size_t i = mTopIndex; i != static_cast<size_t>(-1); --i)
{
- int cmp = a[i] - b[i];
+ exponent cmp = a[i] - b[i];
if (cmp == 0) continue;
if (cmp < 0) return false;
return true;
@@ -379,11 +384,11 @@ public:
size_t monomialSize(ConstMonomial) const { return mMaxMonomialSize; }
- int monomialGetComponent(ConstMonomial a) const { return *a.mValue; }
+ exponent monomialGetComponent(ConstMonomial a) const { return *a.mValue; }
void monomialChangeComponent(Monomial a, int x) const {
- a[mHashIndex] -= *a.mValue;
- a[mHashIndex] += x;
+ a[mHashIndex] -= static_cast<HashValue>(*a.mValue);
+ a[mHashIndex] += static_cast<HashValue>(x);
*a = x;
}
@@ -414,8 +419,9 @@ public:
(ConstMonomial a, ConstMonomial b, ConstMonomial ab) const;
/// Returns the hash of the product of a and b.
- exponent monomialHashOfProduct(ConstMonomial a, ConstMonomial b) const {
- return a[mHashIndex] + b[mHashIndex];
+ HashValue monomialHashOfProduct(ConstMonomial a, ConstMonomial b) const {
+ return static_cast<HashValue>(a[mHashIndex]) +
+ static_cast<HashValue>(b[mHashIndex]);
}
void monomialCopy(ConstMonomial a, Monomial &result) const;
@@ -488,7 +494,7 @@ public:
bool monomialLT(const_monomial a, const_monomial b) const {
for (size_t i = mTopIndex; i != static_cast<size_t>(-1); --i)
{
- int cmp = a[i] - b[i];
+ auto cmp = a[i] - b[i];
if (cmp == 0) continue;
if (cmp < 0) return false;
return true;
@@ -574,18 +580,18 @@ public:
void resetCoefficientStats() const;
private:
- inline exponent computeHashValue(const_monomial a1) const;
+ inline HashValue computeHashValue(const_monomial a1) const;
- long mCharac; // p=mCharac: ring is ZZ/p
+ coefficient mCharac; // p=mCharac: ring is ZZ/p
size_t mNumVars;
int mNumWeights; // stored as negative of weight vectors
size_t mTopIndex;
size_t mHashIndex; // 1 more than mTopIndex. Where the has value is stored.
size_t mMaxMonomialSize;
size_t mMaxMonomialByteSize;
- std::vector<int> mWeights; // 0..mNumWeights * mNumVars - 1.
+ std::vector<exponent> mWeights; // 0..mNumWeights * mNumVars - 1.
- std::vector<int> mHashVals; // one for each variable 0..mNumVars-1
+ std::vector<HashValue> mHashVals; // one for each variable 0..mNumVars-1
// stored as weightvec1 weightvec2 ...
mutable memt::BufferPool mMonomialPool;
@@ -679,22 +685,22 @@ inline void PolyRing::setWeightsOnly(Monomial& a1) const
{
exponent *a = a1.unsafeGetRepresentation();
a++;
- const int *wts = &mWeights[0];
- for (int i=0; i<mNumWeights; i++)
+ auto wts = &mWeights[0];
+ for (size_t i = 0; i < mNumWeights; ++i)
{
- int result = 0;
- for (size_t j=0; j<mNumVars; j++)
+ exponent result = 0;
+ for (size_t j = 0; j < mNumVars; ++j)
result += *wts++ * a[j];
a[mNumVars+i] = result;
}
}
-inline exponent PolyRing::computeHashValue(const_monomial a1) const {
+inline HashValue PolyRing::computeHashValue(const_monomial a1) const {
const exponent* a = a1.unsafeGetRepresentation();
- int hash = *a;
+ HashValue hash = static_cast<HashValue>(*a);
a++;
for (size_t i = 0; i < mNumVars; ++i)
- hash += a[i] * mHashVals[i];
+ hash += static_cast<HashValue>(a[i]) * mHashVals[i];
return hash;
}
@@ -709,7 +715,7 @@ inline int PolyRing::monomialCompare(ConstMonomial a, ConstMonomial b) const
{
for (size_t i = mTopIndex; i != static_cast<size_t>(-1); --i)
{
- int cmp = a[i] - b[i];
+ auto cmp = a[i] - b[i];
if (cmp < 0) return GT;
if (cmp > 0) return LT;
}
@@ -740,7 +746,7 @@ inline bool PolyRing::monomialDivide(ConstMonomial a,
size_t i;
for (i = 1; i <= mNumVars; i++)
{
- int c = a[i] - b[i];
+ exponent c = a[i] - b[i];
if (c >= 0)
result[i] = c;
else
@@ -759,7 +765,7 @@ inline void PolyRing::monomialDivideToNegative(ConstMonomial a,
{
for (size_t i = 0; i <= mHashIndex; ++i)
result[i] = a[i] - b[i];
- MATHICGB_ASSERT(result[mHashIndex] == a[mHashIndex] - b[mHashIndex]);
+ MATHICGB_ASSERT(monomialHashValue(result) == monomialHashValue(a) - monomialHashValue(b));
MATHICGB_ASSERT(!hashValid(a) || !hashValid(b) || hashValid(result));
MATHICGB_ASSERT(computeHashValue(result) ==
computeHashValue(a) - computeHashValue(b));
@@ -850,7 +856,7 @@ inline bool PolyRing::monomialDivide(const_monomial a, const_monomial b, monomia
size_t i;
for (i = 1; i <= mNumVars; i++)
{
- int c = a[i] - b[i];
+ auto c = a[i] - b[i];
if (c >= 0)
result[i] = c;
else
@@ -872,7 +878,7 @@ inline int PolyRing::monomialCompare(const_monomial a, const_monomial b) const
// if (a[i] < b[i]) return GT;
// return LT;
// if (a[i] > b[i]) return LT;
- int cmp = a[i] - b[i];
+ auto cmp = a[i] - b[i];
if (cmp < 0) return GT;
if (cmp > 0) return LT;
}
@@ -911,14 +917,13 @@ inline bool PolyRing::monomialHasStrictlyLargerExponent(
inline void PolyRing::setWeightsOnly(monomial a) const
{
a++;
- const int *wts = &mWeights[0];
- for (int i=0; i<mNumWeights; i++)
- {
- int result = 0;
- for (size_t j=0; j<mNumVars; j++)
- result += *wts++ * a[j];
- a[mNumVars+i] = result;
- }
+ auto wts = &mWeights[0];
+ for (size i = 0; i < mNumWeights; ++i) {
+ exponent result = 0;
+ for (size_t j=0; j<mNumVars; ++j)
+ result += *wts++ * a[j];
+ a[mNumVars+i] = result;
+ }
}
inline bool PolyRing::monomialRelativelyPrime(const_monomial a, const_monomial b) const
@@ -956,29 +961,49 @@ inline void PolyRing::coefficientDivide(coefficient a, coefficient b, coefficien
MATHICGB_ASSERT(result < mCharac);
}
+inline coefficient PolyRing::toCoefficient(int64 a) const {
+ auto modLong = a % mCharac;
+ if (modLong < 0)
+ modLong += mCharac;
+ const coefficient mod = static_cast<coefficient>(modLong);
+ MATHICGB_ASSERT(0 <= mod);
+ MATHICGB_ASSERT(mod < mCharac);
+ return mod;
+}
+
inline void PolyRing::coefficientFromInt(coefficient &result, int a) const
{
- result = a % mCharac;
- if (result < 0) result += mCharac;
+ result = toCoefficient(a);
}
inline void PolyRing::coefficientAddOneTo(coefficient &result) const
{
- result++;
- if (result == mCharac) result = 0;
+ ++result;
+ if (result == mCharac)
+ result = 0;
}
inline void PolyRing::coefficientNegateTo(coefficient &result) const
// result = -result
{
- result = mCharac - result;
+ if (result != 0)
+ result = mCharac - result;
+}
+
+inline coefficient PolyRing::coefficientNegate(coefficient coef) const
+ // result = -result
+{
+ if (coef == 0)
+ return 0;
+ else
+ return mCharac - coef;
}
inline void PolyRing::coefficientAddTo(coefficient &result, coefficient a, coefficient b) const
// result += a*b
{
mStats.n_addmult++;
- long c = a * b + result;
+ auto c = a * b + result;
result = c % mCharac;
}
@@ -987,21 +1012,22 @@ inline void PolyRing::coefficientAddTo(coefficient &result, coefficient a) const
{
mStats.n_add++;
result += a;
- if (result >= mCharac) result -= mCharac;
+ if (result >= mCharac)
+ result -= mCharac;
}
inline void PolyRing::coefficientMultTo(coefficient &result, coefficient a) const
// result *= a
{
mStats.n_mult++;
- long b = result * a;
+ coefficient b = result * a;
result = b % mCharac;
}
inline void PolyRing::coefficientMult(coefficient a, coefficient b, coefficient &result) const
{
mStats.n_mult++;
- long c = b * a;
+ coefficient c = b * a;
result = c % mCharac;
}
diff --git a/src/mathicgb/stdinc.h b/src/mathicgb/stdinc.h
index c9114be..39ce87f 100755
--- a/src/mathicgb/stdinc.h
+++ b/src/mathicgb/stdinc.h
@@ -172,6 +172,12 @@ typedef signed int int32;
typedef signed short int16;
typedef signed char int8;
+/// Bizarrely, OpenMP 2 only supports signed loop
+/// variables. This defect is fixed in OpenMP 3. MSVC 2012 only supports
+/// OpenMP 2. So signed loop indices are forced for loops that are
+/// parallelized using OpenMP.
+typedef signed long OMPIndex;
+
static const size_t BitsPerByte = 8;
static const size_t MemoryAlignment = sizeof(void*);
--
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