diff --git a/.github/workflows/test_build.yml b/.github/workflows/test_build.yml index d679a854847..d4dc91192bb 100644 --- a/.github/workflows/test_build.yml +++ b/.github/workflows/test_build.yml @@ -10,12 +10,21 @@ name: Build and Test Macaulay2 on: workflow_dispatch: - pull_request: + push: branches: - stable + - Dmodules + - Truncations + - Varieties - development - schedule: - - cron: '0 6 * * SUN' # runs tests on the main branch every Sunday at 06:00 UTC + - feature/functors + - feature/monoids + - feature/torsion + - feature/varieties + - refactor/Core + - refactor/engine + - refactor/monoids + - research concurrency: # Cancel in-progress runs when a new workflow with the same group name is triggered @@ -28,7 +37,6 @@ defaults: jobs: build: - if: github.repository == 'Macaulay2/M2' || contains(github.ref, 'global') name: ${{ matrix.build-system }}-${{ matrix.os }}-${{ matrix.compiler }} runs-on: ${{ matrix.os }} strategy: diff --git a/M2/.clang-format b/M2/.clang-format index 06b09d42b57..5cc5aad9346 100644 --- a/M2/.clang-format +++ b/M2/.clang-format @@ -80,7 +80,7 @@ KeepEmptyLinesAtTheStartOfBlocks: false MacroBlockBegin: '' MacroBlockEnd: '' MaxEmptyLinesToKeep: 1 -NamespaceIndentation: Inner +NamespaceIndentation: All ObjCBlockIndentWidth: 2 ObjCSpaceAfterProperty: false ObjCSpaceBeforeProtocolList: false diff --git a/M2/Macaulay2/CMakeLists.txt b/M2/Macaulay2/CMakeLists.txt index bf199311246..02fadb39390 100644 --- a/M2/Macaulay2/CMakeLists.txt +++ b/M2/Macaulay2/CMakeLists.txt @@ -46,6 +46,10 @@ add_subdirectory(m2) # rename to Core? add_subdirectory(packages) add_subdirectory(editors) +# uncomment when developing +#set_target_properties(M2-engine PROPERTIES C_INCLUDE_WHAT_YOU_USE "${IWYU};--transitive_includes_only") +#set_target_properties(M2-engine PROPERTIES CXX_INCLUDE_WHAT_YOU_USE "${IWYU};--transitive_includes_only") + if(BUILD_TESTING) add_subdirectory(tests) endif() diff --git a/M2/Macaulay2/d/Makefile.files.in b/M2/Macaulay2/d/Makefile.files.in index d0d07b16630..10d785969a8 100644 --- a/M2/Macaulay2/d/Makefile.files.in +++ b/M2/Macaulay2/d/Makefile.files.in @@ -153,7 +153,6 @@ interface.o : \ @srcdir@/../e/buffer.hpp \ @srcdir@/../e/mem.hpp \ @srcdir@/../e/../system/mutex.h \ - @srcdir@/../e/imonorder.hpp \ @srcdir@/../e/aring.hpp \ @srcdir@/../e/exceptions.hpp \ @srcdir@/../e/ZZ.hpp \ @@ -246,7 +245,6 @@ engine.o : \ @srcdir@/../e/ZZ.hpp \ @srcdir@/../e/matrix.hpp \ @srcdir@/../e/monoid.hpp \ - @srcdir@/../e/imonorder.hpp \ @srcdir@/../e/interface/monomial-ordering.h \ @srcdir@/../e/freemod.hpp \ @srcdir@/../e/schorder.hpp \ diff --git a/M2/Macaulay2/d/interface.dd b/M2/Macaulay2/d/interface.dd index b25a030e876..8d2a26cf795 100644 --- a/M2/Macaulay2/d/interface.dd +++ b/M2/Macaulay2/d/interface.dd @@ -957,12 +957,13 @@ export rawInverse(e:Expr):Expr := ( setupfun("rawInverse",rawInverse); export rawMultiDegree(e:Expr):Expr := ( - when e - is x:RawRingElementCell do toExpr( Ccode(RawArrayIntOrNull, "IM2_RingElement_multidegree(", x.p, ")" ) ) - is x:RawFreeModuleCell do toExpr( Ccode(RawArrayIntOrNull, "IM2_FreeModule_get_degrees(", x.p, ")" ) ) - is x:RawMatrixCell do toExpr( Ccode(RawArrayIntOrNull, "IM2_Matrix_get_degree(", x.p, ")" ) ) - else WrongArg("a raw ring element, matrix, or free module")); -setupfun("rawMultiDegree",rawMultiDegree); + when e + -- TODO: unify these function names + is x:RawRingElementCell do toExpr( Ccode(RawArrayIntOrNull, "IM2_RingElement_multidegree(", x.p, ")" ) ) + is x:RawFreeModuleCell do toExpr( Ccode(RawArrayIntOrNull, "IM2_FreeModule_get_degrees(", x.p, ")" ) ) + is x:RawMatrixCell do toExpr( Ccode(RawArrayIntOrNull, "IM2_Matrix_get_degree(", x.p, ")" ) ) + else WrongArg("a raw ring element, matrix, or free module")); +setupfun("rawMultiDegree", rawMultiDegree); export rawDiscreteLog(e:Expr):Expr := ( when e diff --git a/M2/Macaulay2/e/CMakeLists.txt b/M2/Macaulay2/e/CMakeLists.txt index 72fd6ae7a45..9631355a6f7 100644 --- a/M2/Macaulay2/e/CMakeLists.txt +++ b/M2/Macaulay2/e/CMakeLists.txt @@ -155,6 +155,8 @@ set(CPPONLYFILES # these files all have .hpp and .cpp files # MES TODO: put these NCAlgebra files (and friends) into alpha order in list set(SRCLIST + MonomialOrder + MonomialOrderOld BasicPoly BasicPolyList BasicPolyListParser @@ -231,7 +233,6 @@ set(SRCLIST gbweight hermite hilb - imonorder int-bag interreduce interrupted @@ -463,6 +464,8 @@ if(BUILD_TESTING) unit-tests/basics-test.cpp unit-tests/MatrixIOTest.cpp + unit-tests/MonomialOrderTest.cpp + unit-tests/fromStream.cpp unit-tests/testMain.cpp # not needed, except for GC_INIT ) diff --git a/M2/Macaulay2/e/ExponentList.hpp b/M2/Macaulay2/e/ExponentList.hpp index 3b0f21bfa3a..f8d7b92410f 100644 --- a/M2/Macaulay2/e/ExponentList.hpp +++ b/M2/Macaulay2/e/ExponentList.hpp @@ -150,6 +150,7 @@ class ExponentList static void mult(ConstExponents a, ConstExponents b, Vector& result); // compute the quotient a:b static void quotient(ConstExponents a, ConstExponents b, Vector& result); + static void quotient(ConstExponents a, ConstExponents b, Exponents result); static void power(ConstExponents a, Exponent n, Vector& result); static void monsyz(ConstExponents a, diff --git a/M2/Macaulay2/e/Makefile.files b/M2/Macaulay2/e/Makefile.files index f75b786285a..2fb2313eefd 100644 --- a/M2/Macaulay2/e/Makefile.files +++ b/M2/Macaulay2/e/Makefile.files @@ -9,6 +9,8 @@ INTERFACE2 = interrupted # NCAlgebras/SuffixTree \ INTERFACE = \ + MonomialOrder \ + MonomialOrderOld \ myalloc \ BasicPoly \ BasicPolyList \ @@ -104,7 +106,6 @@ INTERFACE = \ SLP \ overflow \ memory-status \ - imonorder \ sagbi \ monideal-minprimes \ lapack \ diff --git a/M2/Macaulay2/e/MonomialOrder.cpp b/M2/Macaulay2/e/MonomialOrder.cpp new file mode 100644 index 00000000000..5e5a63fb322 --- /dev/null +++ b/M2/Macaulay2/e/MonomialOrder.cpp @@ -0,0 +1,606 @@ +#include +#include + +#include "MonomialOrder.hpp" +#include "exceptions.hpp" + +namespace MO { + MonomialBlock::MonomialBlock(order_type t, + const std::vector& data, + int packing, + int& next_var + ) + : mType(t), + mPacking(packing), + mNumVars(0), + mStartVariable(next_var) // will change for MO::Weights. + { + if (t == MO::PositionUp or t == MO::PositionDown) + { + mPacking = 1; + mStartVariable = 0; + return; + } + + if (data.size() == 0) + { + throw exc::engine_error("expected each monomial block to have positive number of variables"); + } + + if (t == MO::Weights) + { + mNumVars = 0; + mStartVariable = data[0]; + std::copy(data.cbegin()+1, data.cend(), std::back_inserter(mWeights)); + return; + } + + mNumVars = (data.size() >= 2 ? data.size() : data[0]); + if (mNumVars <= 0) + { + throw exc::engine_error("expected each monomial block to have positive number of variables"); + } + + switch (t) { + case MO::GroupLex: + case MO::GroupRevLex: + mPacking = 1; + break; + case MO::GRevLexWeights: + //TODO: check that these are all positive. + std::copy(data.cbegin(), data.cend(), std::back_inserter(mWeights)); + break; + case MO::Lex: + case MO::RevLex: + case MO::GRevLex: + // The following should not occur here (they are handled above) + case MO::Weights: + case MO::PositionUp: + case MO::PositionDown: + break; + // The following should not occur here (they are never placed into this type). + case MO::Packing: + throw exc::engine_error("internal error: received incorrect monomial block"); + }; + + next_var += mNumVars; + } + + std::vector prepend(int a, const std::vector& b) + { + std::vector result; + result.push_back(a); + for (auto c : b) result.push_back(c); + return result; + } + + std::pair> MonomialBlock::toMonomialType() const + { + std::pair> result; + std::vector newwts; // I don't know how to easily get around this? + switch (mType) + { + case MO::Lex: return {MO::Lex, {mNumVars}}; + case MO::RevLex: return {MO::RevLex, {mNumVars}}; + case MO::GroupLex: return {MO::GroupLex, {mNumVars}}; + case MO::GroupRevLex: return {MO::GroupRevLex, {mNumVars}}; + case MO::GRevLex: return {MO::GRevLex, {mNumVars}}; + case MO::GRevLexWeights: return {MO::GRevLexWeights, mWeights}; + case MO::Weights: return {MO::Weights, prepend(mStartVariable, mWeights)}; + case MO::PositionUp: return {MO::PositionUp, {}}; + case MO::PositionDown: return {MO::PositionUp, {}}; + case MO::Packing: return {MO::Packing, {100}}; // This should never occur. TODO: error here + } + } + + std::string toString(MO::order_type t) + { + switch (t) { + case MO::Lex: return "Lex"; + case MO::RevLex: return "RevLex"; + case MO::GroupLex: return "GroupLex"; + case MO::GroupRevLex: return "GroupRevLex"; + case MO::GRevLex: return "GRevLex"; + case MO::GRevLexWeights: return "GRevLexWeights"; + case MO::Weights: return "Weights"; + case MO::PositionUp: return "PositionUp"; + case MO::PositionDown: return "PositionDown"; + case MO::Packing: return "Packing"; + } + } + + void MonomialBlock::debugDisplay(std::ostream& o) const + { + o << " " << mType; + o << "[#vars=" << mNumVars << " start=" << mStartVariable; + if (mWeights.size() > 0) + o << " weights=" << mWeights; + // if (mPacking != 1) + o << " packing " << mPacking; + o << "]"; + } + + EncodingBlock::EncodingBlock(MonomialBlock b, + int& next_slot + ) + : mBlock(b), + mStartEncoded(next_slot) + { + switch (b.type()) { + case MO::Lex: + case MO::RevLex: + mNumEncoded = b.numVars() / b.packing(); + break; + case MO::GroupLex: // Lex + case MO::GroupRevLex: + mNumEncoded = b.numVars(); + break; + case MO::GRevLex: + mNumEncoded = 1 + b.numVars() / b.packing(); + break; + case MO::GRevLexWeights: + mNumEncoded = 1 + b.numVars() / b.packing(); + break; + case MO::Weights: + mNumEncoded = 1; + break; + case MO::PositionUp: + case MO::PositionDown: + case MO::Packing: + // this is handled higher up the food chain + break; + }; + next_slot += mNumEncoded; + } + + void EncodingBlock::debugDisplay(std::ostream& o) const + { + mBlock.debugDisplay(o); + o << " encoded [start,#] = [" << mStartEncoded << ", " << mNumEncoded << "]"; + } + +}; // end of MO namespace + +std::ostream& operator<<(std::ostream& o, MO::order_type t) +{ + return (o << MO::toString(t)); +} + +std::ostream& operator<<(std::ostream& o, MO::MonomialBlock b) +{ + o << b.type() << " => "; + switch (b.type()) { + case MO::Lex: + case MO::RevLex: + case MO::GroupLex: + case MO::GroupRevLex: + case MO::GRevLex: + o << b.numVars(); + break; + case MO::GRevLexWeights: + case MO::Weights: + o << b.weights(); + break; + case MO::PositionUp: + case MO::PositionDown: + o << "{}"; + break; + case MO::Packing: + o << b.packing(); + break; + } + return o; + // return (o << MO::toString(t)); +} + +std::ostream& operator<<(std::ostream& o, const std::pair>& a) +{ + o << "{" << a.first; //<< " => " << a.second << "}"; + return o; +} + +std::ostream& operator<<(std::ostream& o, const NewAgeMonomialOrder::MOInput& a) + { + for (const auto& a1 : a) + { + o << "{" << a1.first << ", " << a1.second << "}" << std::endl; + } + return o; + } + +std::ostream& operator<<(std::ostream& o, const std::vector& a) +{ + o << "{"; + for (int i = 0; i < a.size(); i++) + { + if (i > 0) o << ","; + o << a[i]; + } + o << "}"; + return o; +} + +NewAgeMonomialOrder::NewAgeMonomialOrder(int numvars, const MOInput& input) + : mNumVars(numvars) +{ + int next_var = 0; + int packing = 1; + + // Create each block in the monomial order in turn. + for (const auto &i : input) + { + // We do not store directly the Packing elements + if (i.first == MO::Packing) + { + if (i.second.size() != 1) + throw exc::engine_error("packing size not given correctly"); + packing = i.second[0]; + if (not (packing == 1 or packing == 2 or packing == 4)) + throw exc::engine_error("packing size must be 1, 2 or 4"); + } + else + mParts.push_back(MO::MonomialBlock(i.first, i.second, packing, next_var)); + } + + if (numvars != next_var) + throw exc::engine_error("wrong number of variables"); +} + +auto NewAgeMonomialOrder::toMonomialType() const -> MOInput +{ + MOInput result; + int packing = 1; + for (const auto &m : mParts) + { + if (m.packing() != packing + and (m.type() == MO::Lex or m.type() == MO::GRevLex or m.type() == MO::GRevLexWeights)) + { + result.push_back({MO::Packing, {m.packing()}}); + packing = m.packing(); + } + result.push_back(m.toMonomialType()); + } + return result; +} + +std::ostream& operator<<(std::ostream& o, const NewAgeMonomialOrder& mo) { + int packing = 1; + o << "MonomialOrder => {"; + for (int i = 0; i < mo.mParts.size(); ++i) + { + const auto& b = mo.mParts[i]; + if (i == 0) + o << "\n "; + else + o << ",\n "; + if (b.packing() != packing + and (b.type() == MO::Lex or b.type() == MO::GRevLex or b.type() == MO::GRevLexWeights)) + { + o << "Packing => " << b.packing() << ",\n "; + } + o << b; + } + o << "\n }"; + return o; +} + +void NewAgeMonomialOrder::debugDisplay(std::ostream& o) const +{ + o << "MonomialOrder:" << std::endl; + o << " #vars = " << mNumVars << std::endl; + for (const auto& m : mParts) + { + m.debugDisplay(o); + o << std::endl; + } +} + +std::vector NewAgeMonomialOrder::invertibleVariables() const +{ + auto isInvertible = invertibleVariablesPredicate(); + std::vector result; + for (int i=0; i < isInvertible.size(); ++i) + { + if (isInvertible[i]) result.push_back(i); + } + return result; +} + +std::vector NewAgeMonomialOrder::invertibleVariablesPredicate() const +{ + std::vector result; + for (int i=0; i < numVars(); ++i) + result.push_back(false); + for (const auto& b : mParts) + { + if ((b.type() == MO::GroupRevLex) or (b.type() == MO::GroupLex)) + { + for (int i = 0; i < b.numVars(); ++i) + result[b.firstVar() + i] = true; + } + } + return result; +} + +std::vector NewAgeMonomialOrder::localVariables() const +{ + auto isLocal = localVariablesPredicate(); + std::vector result; + for (int i=0; i < isLocal.size(); ++i) + { + if (isLocal[i]) result.push_back(i); + } + return result; +} + +std::vector NewAgeMonomialOrder::localVariablesPredicate() const +{ + std::vector status(numVars(), 0); // initialized to 0's. meaning: 0: don't know yet, 1: is global, -1: is local. + for (const auto& b : mParts) + { + switch (b.type()) + { + case MO::Lex: + case MO::GroupLex: + case MO::GRevLex: + case MO::GRevLexWeights: + // these are global variables + for (int i=0; i 0) status[v] = 1; + else if (b.weights()[i] < 0) status[v] = -1; + // don't change value if weight is 0. + } + } + break; + case MO::PositionUp: + case MO::PositionDown: + case MO::Packing: + // these are ignored + break; + } + } + // At this point status[v] should be +1 or -1 for each v. + // now we create the result + std::vector result; + for (int i=0; i &grading, + int nvars, + int which, + int value) +{ + for (int i = 0; i < nvars; i++) + if (i == which) + grading.push_back(value); + else + grading.push_back(0); +} +static void write_weights(std::vector &grading, + int nvars, + int firstvar, + const std::vector& wts, + int nwts) +// place nvars ints into grading: 0 ... 0 wts[0] wts[1] ... wts[nwts-1] 0 .... +// 0 +// where wts[0] is in the 'firstvar' location. If wts is NULL, treat it as the +// vector with nwts '1's. +{ + for (int i = 0; i < firstvar; i++) grading.push_back(0); + if (wts.size() == 0) + for (int i = 0; i < nwts; i++) grading.push_back(1); + else + for (int i = 0; i < nwts; i++) grading.push_back(wts[i]); + for (int i = firstvar + nwts; i < nvars; i++) grading.push_back(0); +} + +bool NewAgeMonomialOrder::monomialOrderingToMatrix( + std::vector& mat, + bool& base_is_revlex, + int& component_direction, // -1 is Down, +1 is Up, 0 is not present + int& component_is_before_row // -1 means: at the end. 0 means before the + // order. + // and r means considered before row 'r' of the matrix. + ) const +{ + // a false return value means: this order cannot be used in GB computations (i.e. there are invertible variables). + int nvars = numVars(); + base_is_revlex = true; + enum LastBlock { LEX, REVLEX, WEIGHTS, NONE }; + LastBlock last = NONE; + int nrows = 0; + int firstvar = 0; + component_direction = 0; + component_is_before_row = -2; // what should the default value be? Probably: -1. + size_t last_element = 0; // The vector 'mat' will be resized back to this + // value if the last part of the order is lex or + // revlex. + for (const auto& b : mParts) + { + switch (b.type()) + { + case MO::Lex: + // printf("lex %d\n", p->nvars); + last_element = mat.size(); + for (int j = 0; j < b.numVars(); ++j) + { + write_row(mat, nvars, firstvar + j, 1); + } + last = LEX; + firstvar += b.numVars(); + nrows += b.numVars(); + break; + case MO::GRevLex: + // printf("grevlex %d %ld\n", p->nvars, p->wts); + // TODO: write the 1's here directly, as b.weights() is not set... + write_weights(mat, nvars, firstvar, b.weights(), b.numVars()); + last_element = mat.size(); + for (int j = b.numVars() - 1; j >= 1; --j) + { + write_row(mat, nvars, firstvar + j, -1); + } + last = REVLEX; + firstvar += b.numVars(); + nrows += b.numVars(); + break; + case MO::GRevLexWeights: + // printf("grevlex_wts %d %ld\n", p->nvars, p->wts); + write_weights(mat, nvars, firstvar, b.weights(), b.numVars()); + last_element = mat.size(); + for (int j = b.numVars() - 1; j >= 1; --j) + { + write_row(mat, nvars, firstvar + j, -1); + } + last = REVLEX; + firstvar += b.numVars(); + nrows += b.numVars(); + break; + case MO::RevLex: + // printf("revlex %d\n", p->nvars); + last_element = mat.size(); + for (int j = b.numVars() - 1; j >= 0; --j) + { + write_row(mat, nvars, firstvar + j, -1); + } + last = REVLEX; + firstvar += b.numVars(); + nrows += b.numVars(); + break; + case MO::Weights: + // printf("matsize= %d weights %d p->wts=%lu\n", mat.size(), + // p->nvars, p->wts); + write_weights(mat, nvars, b.firstVar(), b.weights(), b.weights().size()); + nrows++; + last_element = mat.size(); + last = WEIGHTS; + break; + case MO::GroupLex: + case MO::GroupRevLex: + return false; + case MO::PositionUp: + component_direction = 1; + component_is_before_row = nrows; + break; + case MO::PositionDown: + component_direction = -1; + component_is_before_row = nrows; + break; + case MO::Packing: + // DO nothing + break; + } + } + if (last == LEX) + { + // last block was lex, so use lex tie-breaker + mat.resize(last_element); + if (nrows == component_is_before_row) component_is_before_row = -1; + base_is_revlex = false; + } + else if (last == REVLEX) + { + // last block was revlex, so use revlex tie-breaker + if (nrows == component_is_before_row) component_is_before_row = -1; + mat.resize(last_element); + } + else + { + // last block is a weight vector, so use revlex as the tie-breaker. + // nothing to change here. + } + return true; +} + +//////////////////////////// Encoder code ///////////////////////////////////// +MonomialEncoder::MonomialEncoder(const NewAgeMonomialOrder& mo) +{ + int next_var = 0; + int next_slot = 0; + int packing = 1; + + mComponentOccursBeforeThisBlock = -1; +#if 0 + // TODO: reinstate the code in this block. + // Create each block in the monomial order in turn. + for (const auto &i : input) + { + if (i.first == MO::PositionUp) + { + mComponentOccursBeforeThisBlock = mParts.size(); + mEncodedComponentLocation = next_slot; + mPositionUp = true; + next_slot++; + continue; + } + else if (i.first == MO::PositionDown) + { + mComponentOccursBeforeThisBlock = mParts.size(); + mEncodedComponentLocation = next_slot; + mPositionUp = false; + next_slot++; + continue; + } + mParts.push_back(MO::Block(i.first, i.second, packing, next_var, next_slot)); + } + + if (mComponentOccursBeforeThisBlock == -1) + { + // In this case, Position wasn't given. We add in Position => Up here. + mComponentOccursBeforeThisBlock = mParts.size(); + mEncodedComponentLocation = next_slot; + mPositionUp = true; + next_slot++; + } + if (numvars != next_var) + throw exc::engine_error("wrong number of variables"); +#endif + mEncodedSize = next_slot; +} + +void MonomialEncoder::debugDisplay(std::ostream& o) const +{ + o << "MonomialEncoder:" << std::endl; + o << " #vars = " << mNumVars << std::endl; + o << " #encoded = " << mEncodedSize << std::endl; + o << " #compslot = " << mEncodedComponentLocation << std::endl; + o << " #compblock = " << mComponentOccursBeforeThisBlock << std::endl; + o << " #comp up = " << (mPositionUp ? "true" : "false") << std::endl; + for (const auto& m : mParts) + { + m.debugDisplay(o); + o << std::endl; + } +} + +// Local Variables: +// compile-command: "make -C $M2BUILDDIR/Macaulay2/e " +// indent-tabs-mode: nil +// End: + diff --git a/M2/Macaulay2/e/MonomialOrder.hpp b/M2/Macaulay2/e/MonomialOrder.hpp new file mode 100644 index 00000000000..2bef84fc2ec --- /dev/null +++ b/M2/Macaulay2/e/MonomialOrder.hpp @@ -0,0 +1,250 @@ +#pragma once + +#include +#include +#include + +using ExponentVector = int*; +using EncodedMonomial = int*; + +///TODO: maybe we need 2 classes here. +/// 1. MonomialOrder +/// 2. MonomialEncoder (and decoder). +/// this is like the current intermal monomial ordering code. +/// it has a list of parts (monomial order parts, and (firstslot, numslot, packing, overflow) info for each. +/// operations +/// - encodeExponentVector +/// - decodeToExponentVector +/// - mult (with overflow checks) (this needs the encoder for the over flow fields) +/// - comparison (this might not need the Encoder at all). + +/// TODO: this is in some sense a todo list for my current modifications +/// In Macaulay2, a MonomialOrder is a data structure that allows teh following options: +/// - Create a monomial order from its parts +/// - Encode a monomial into a list of ints for which comparison and multipliciation is easy. +/// - Decode a monomial back to an exponent vector, also a "varpower" product (sparse monomial). +/// - Comparison function +/// - Multiplication of two monomials, with overflow checking. +/// - to/from a "matrix" description of the monomial order. This is lossy, in that +/// the matrix representation doesn't know about the packing information in the order. +/// Additionally, there are methods for querying aspects of the monomial order. +/// - How many slots (i.e. ints in the encoded monomial) correspond to the first d parts of the monomial order?) +/// - Is this a global ordering? +/// - Which variables are local? +/// - How many (and which) variables allow negative exponents? +/// - Is this graded reverse lex order? Is it Lex order? +/// - What else is needed? +/// Questions remaining: +/// - where best to put torsion in monoids? In the mult function? + +/// MO = monomial order types, internal functions and types for monomial orders */ +namespace MO { + enum order_type { + Packing, /**< # exponents per (currently 32 bit) word. Only affects Lex, + GRevLex, GRevLexWeights (but not the weight part). Does not + affect the monomial order. */ + Lex, /**< lexicographic order on a specified number of variables */ + GRevLex, /**< graded reverse lexicographic order + on a specified number of variables, + with weights 1 for each variable */ + GRevLexWeights, /**< graded reverse lexicographic order + with specified positive weights, on a number of variables + */ + RevLex, /**< reverse lexicographic order + on a specified number of variables. + Not a global order! */ + Weights, + GroupLex, // previously MO_LAURENT + GroupRevLex, // previously MO_LAURENT_REVLEX + PositionUp, + PositionDown, + }; + // Previous names and values + // MO_LEX = 1, + // MO_GREVLEX = 4, + // MO_GREVLEX_WTS = 7, + // MO_REVLEX = 10, + // MO_WEIGHTS = 11, + // MO_LAURENT = 12, /* Lex order here */ + // MO_LAURENT_REVLEX = 13, /* Rev lex order here */ + // MO_POSITION_UP = 15, + // MO_POSITION_DOWN = 16 + + std::string toString(MO::order_type t); + + // Information about a piece of the monomial order + class MonomialBlock + { + private: + enum order_type mType; + + /** This is the number of exponents to pack in to a word. This is + the one piece of data that isn't part of the monomial order. + But we keep it, because it is easier to collect the + information once from the front end + */ + int mPacking; + + /** number of variables for this block + * for Lex, GRevLex, RevLex, GroupLex, GroupRevLex: number of variables + * for GRevLexWeights: number of variables + * PositionUp, PositionDown: 0 + * Weights: 0. + */ + int mNumVars; + + /// index of the first variable for this block + int mStartVariable; + + // Only used for MO::Weights, MO::GRevLex + std::vector mWeights; + + public: + /** constructor */ + MonomialBlock(enum order_type t, + const std::vector& data, + int packing, + int& next_var + ); + + enum order_type type() const { return mType; } + int numVars() const { return mNumVars; } + int firstVar() const { return mStartVariable; } + auto weights() const -> const std::vector& { return mWeights; } + int packing() const { return mPacking; } + void debugDisplay(std::ostream& o) const; + std::pair> toMonomialType() const; + }; + + // Information for encoding and decoding monomials + class EncodingBlock + { + private: + MonomialBlock mBlock; + + /// encoded locations (start, length). + int mStartEncoded; + int mNumEncoded; + public: + /** constructor */ + EncodingBlock(MonomialBlock b, + int& next_slot + ); + + //TODO: reinstate the ones needed here. + // enum order_type type() const { return mType; } + // int numVars() const { return mNumVars; } + // auto weights() const -> const std::vector& { return mWeights; } + void debugDisplay(std::ostream& o) const; + }; +}; // end namespace MO + + +/////////////////////////////////////////////////////////////////////// + +class NewAgeMonomialOrder +{ + friend std::ostream& operator<<(std::ostream&, const NewAgeMonomialOrder& mo); +private: + /// Number of variables + int mNumVars; + + // The sequence of blocks for the monomial order + std::vector mParts; + + // Informational about the order + std::vector mIsLocalVariable; + std::vector mIsGroupVariable; + +public: + using MOInput = std::vector>>; + NewAgeMonomialOrder(int nvars, const MOInput& input); + + MOInput toMonomialType() const; + + /// TODO + bool isWellDefined(int verbose_level) const; + + void debugDisplay(std::ostream& o) const; + + int numVars() const { return mNumVars; } +public: + ////// informational functions /////// + + // Is the monomial order Lex or GRevLex? + + + // The list of indices of variables which are < 1. + std::vector localVariables() const; + std::vector invertibleVariables() const; + + std::vector invertibleVariablesPredicate() const; // aka laurentVariables + std::vector localVariablesPredicate() const; // aka non term order variables. + + bool hasLocalVariables() const { return localVariables().size() > 0; } + bool hasInvertibleVariables() const { return invertibleVariables().size() > 0; } +public: + //TODO unencoded comparison. + int compareExponentVectors(const ExponentVector &a, const ExponentVector* b) const; + + //TODO: is an unencoded monomial in the subring where the first n parts are 0? + + //TODO: to and from std::vector data for mgb and the front end. + + bool monomialOrderingToMatrix( + std::vector& mat, + bool& base_is_revlex, + int& component_direction, // -1 is Down, +1 is Up, 0 is not present + int& component_is_before_row // -1 means: at the end. 0 means before the + // order. + // and r means considered before row 'r' of the matrix. + ) const; +}; + +std::ostream& operator<<(std::ostream& o, const std::vector& a); +std::ostream& operator<<(std::ostream& o, MO::order_type t); +std::ostream& operator<<(std::ostream& o, const std::pair>& a); +std::ostream& operator<<(std::ostream& o, const NewAgeMonomialOrder& mo); +std::ostream& operator<<(std::ostream& o, const NewAgeMonomialOrder::MOInput& a); + + +/////////////// Encoder code ///////////////////////////// +class MonomialEncoder +{ +private: + std::vector mParts; + + /// Number of variables + int mNumVars; + + /// Size of encoded monomial (includes component, if present). + int mEncodedSize; + + /// the component comes after the following block #. + int mComponentOccursBeforeThisBlock; + + /// the component value in the encoded monomial is at this location + int mEncodedComponentLocation; + + /// if true, components are encoded via negation, so comparison will be faster + bool mPositionUp; + +public: + MonomialEncoder(const NewAgeMonomialOrder& mo); + + //TODO. + int compareEncodedMonomials(); + void multEncodedMonomials(); // can overflow + void encode(const ExponentVector& a, EncodedMonomial& result) const; + void decode(const EncodedMonomial& a, ExponentVector& b) const; + + void debugDisplay(std::ostream& o) const; +}; + + +/* +// Local Variables: +// compile-command: "make -C $M2BUILDDIR/Macaulay2/e " +// indent-tabs-mode: nil +// End: +*/ diff --git a/M2/Macaulay2/e/MonomialOrderOld.cpp b/M2/Macaulay2/e/MonomialOrderOld.cpp new file mode 100644 index 00000000000..333040cf503 --- /dev/null +++ b/M2/Macaulay2/e/MonomialOrderOld.cpp @@ -0,0 +1,1680 @@ +#include "MonomialOrderOld.hpp" +#include "exceptions.hpp" +#include "error.h" +#include +#include + +static struct mon_part_rec_ *mo_make(enum MonomialOrdering_type type, + int nvars, + const int *wts) +{ + mon_part result; + result = getmemstructtype(mon_part); + result->type = type; + result->nvars = nvars; + if (wts != nullptr) + { + int i; + result->wts = getmematomicvectortype(int, nvars); + for (i = 0; i < nvars; i++) result->wts[i] = wts[i]; + } + else + result->wts = nullptr; + return result; +} + +static MonomialOrdering *make_mon_order(int n) +{ + static unsigned int next_hash = 23023421; + MonomialOrdering *z = getmemarraytype(MonomialOrdering *, n); + z->len = n; + z->_hash = next_hash++; + int i; + for (i = 0; i < n; i++) z->array[i] = nullptr; + return z; +} + +static MonomialOrdering *M2_mo_offset(const MonomialOrdering *mo, int offset) +{ + int i, j; + MonomialOrdering *result = make_mon_order(mo->len); + for (i = 0; i < mo->len; i++) + { + mon_part p = mo->array[i]; + if (p->type != MO_WEIGHTS) + result->array[i] = mo_make(p->type, p->nvars, p->wts); + else + { + mon_part q = mo_make(MO_WEIGHTS, offset + p->nvars, nullptr); + q->wts = getmemvectortype(int, q->nvars); + for (j = 0; j < offset; j++) q->wts[j] = 0; + for (; j < q->nvars; j++) q->wts[j] = p->wts[j - offset]; + } + } + return result; +} + +static bool is_good(mon_part p) +{ + switch (p->type) + { + case MO_LEX: + case MO_LEX2: + case MO_LEX4: + case MO_GREVLEX: + case MO_GREVLEX2: + case MO_GREVLEX4: + case MO_GREVLEX_WTS: + case MO_GREVLEX2_WTS: + case MO_GREVLEX4_WTS: + case MO_LAURENT: + case MO_NC_LEX: + case MO_LAURENT_REVLEX: + case MO_REVLEX: + case MO_WEIGHTS: + return (p->nvars > 0); + case MO_POSITION_UP: + case MO_POSITION_DOWN: + return true; + } + return false; +} + + +std::ostringstream &toString(std::ostringstream &o, int len, int *p) +{ + o << "{"; + for (int i = 0; i < len; i++) + { + if (i > 0) o << ","; + o << p[i]; + } + o << "}"; + return o; +} + +std::ostringstream &ones(std::ostringstream &o, int len) +{ + o << "{"; + for (int i = 0; i < len; i++) + { + if (i > 0) o << ","; + o << 1; + } + o << "}"; + return o; +} + +std::string MonomialOrderings::toString(const MonomialOrdering *mo) +{ + std::ostringstream o; + o << "MonomialOrder => {"; + for (int i = 0; i < mo->len; i++) + { + mon_part p = mo->array[i]; + bool p_ones = false; + if (i == 0) + o << "\n "; + else + o << ",\n "; + switch (p->type) + { + case MO_LEX: + o << "Lex => " << p->nvars; + break; + case MO_LEX2: + o << "LexSmall => " << p->nvars; + break; + case MO_LEX4: + o << "LexTiny => " << p->nvars; + break; + case MO_GREVLEX: + o << "GRevLex => "; + p_ones = true; + break; + case MO_GREVLEX2: + o << "GRevLexSmall => "; + p_ones = true; + break; + case MO_GREVLEX4: + o << "GRevLexTiny => "; + p_ones = true; + break; + case MO_GREVLEX_WTS: + o << "GRevLex => "; + break; + case MO_GREVLEX2_WTS: + o << "GRevLexSmall => "; + break; + case MO_GREVLEX4_WTS: + o << "GRevLexTiny => "; + break; + case MO_REVLEX: + o << "RevLex => " << p->nvars; + break; + case MO_WEIGHTS: + o << "Weights => "; + break; + case MO_LAURENT: + o << "GroupLex => " << p->nvars; + break; + case MO_LAURENT_REVLEX: + o << "GroupRevLex => " << p->nvars; + break; + case MO_NC_LEX: + o << "NCLex => " << p->nvars; + break; + case MO_POSITION_UP: + o << "Position => Up"; + break; + case MO_POSITION_DOWN: + o << "Position => Down"; + break; + default: + o << "UNKNOWN"; + break; + } + if (p->wts != nullptr) { ::toString(o, p->nvars, p->wts); } + else if (p_ones) + { + ::ones(o, p->nvars); + } + } + o << "\n }"; + return o.str(); +} + +// Many monomial ordering routines are in monordering.c +// Here are some that use c++ features, so cannot be there. +// @todo Make monordering.{h,c} into a c++ class. + +static void write_row(std::vector &grading, + int nvars, + int which, + int value) +{ + for (int i = 0; i < nvars; i++) + if (i == which) + grading.push_back(value); + else + grading.push_back(0); +} +static void write_weights(std::vector &grading, + int nvars, + int firstvar, + int *wts, + int nwts) +// place nvars ints into grading: 0 ... 0 wts[0] wts[1] ... wts[nwts-1] 0 .... +// 0 +// where wts[0] is in the 'firstvar' location. If wts is NULL, treat it as the +// vector with nwts '1's. +{ + for (int i = 0; i < firstvar; i++) grading.push_back(0); + if (wts == nullptr) + for (int i = 0; i < nwts; i++) grading.push_back(1); + else + for (int i = 0; i < nwts; i++) grading.push_back(wts[i]); + for (int i = firstvar + nwts; i < nvars; i++) grading.push_back(0); +} + +bool monomialOrderingToMatrix( + const struct MonomialOrdering& mo, + std::vector& mat, + bool& base_is_revlex, + int& component_direction, // -1 is Down, +1 is Up, 0 is not present + int& component_is_before_row // -1 means: at the end. 0 means before the + // order. + // and r means considered before row 'r' of the matrix. + ) +{ + // a false return value means an error has occurred. + int nvars = rawNumberOfVariables(&mo); + base_is_revlex = true; + enum LastBlock { LEX, REVLEX, WEIGHTS, NONE }; + LastBlock last = NONE; + int nwts = 0; // local var used in MO_WEIGHTS section + int nrows = 0; + int firstvar = 0; + component_direction = 0; + component_is_before_row = + -2; // what should the default value be? Probably: -1. + size_t last_element = 0; // The vector 'mat' will be resized back to this + // value if the last part of the order is lex or + // revlex. + for (int i = 0; i < mo.len; i++) + { + mon_part p = mo.array[i]; + switch (p->type) + { + case MO_LEX: + case MO_LEX2: + case MO_LEX4: + // printf("lex %d\n", p->nvars); + last_element = mat.size(); + for (int j = 0; j < p->nvars; j++) + { + write_row(mat, nvars, firstvar + j, 1); + } + last = LEX; + firstvar += p->nvars; + nrows += p->nvars; + break; + case MO_GREVLEX: + case MO_GREVLEX2: + case MO_GREVLEX4: + // printf("grevlex %d %ld\n", p->nvars, p->wts); + write_weights(mat, nvars, firstvar, p->wts, p->nvars); + last_element = mat.size(); + for (int j = p->nvars - 1; j >= 1; --j) + { + write_row(mat, nvars, firstvar + j, -1); + } + last = REVLEX; + firstvar += p->nvars; + nrows += p->nvars; + break; + case MO_GREVLEX_WTS: + case MO_GREVLEX2_WTS: + case MO_GREVLEX4_WTS: + // printf("grevlex_wts %d %ld\n", p->nvars, p->wts); + write_weights(mat, nvars, firstvar, p->wts, p->nvars); + last_element = mat.size(); + for (int j = p->nvars - 1; j >= 1; --j) + { + write_row(mat, nvars, firstvar + j, -1); + } + last = REVLEX; + firstvar += p->nvars; + nrows += p->nvars; + break; + case MO_REVLEX: + // printf("revlex %d\n", p->nvars); + last_element = mat.size(); + for (int j = p->nvars - 1; j >= 0; --j) + { + write_row(mat, nvars, firstvar + j, -1); + } + last = REVLEX; + firstvar += p->nvars; + nrows += p->nvars; + break; + case MO_WEIGHTS: + // printf("matsize= %d weights %d p->wts=%lu\n", mat.size(), + // p->nvars, p->wts); + nwts = (p->nvars > nvars ? nvars : p->nvars); + write_weights(mat, nvars, 0, p->wts, nwts); + nrows++; + last_element = mat.size(); + last = WEIGHTS; + break; + case MO_LAURENT: + case MO_LAURENT_REVLEX: + case MO_NC_LEX: + return false; + break; + case MO_POSITION_UP: + component_direction = 1; + component_is_before_row = nrows; + break; + case MO_POSITION_DOWN: + component_direction = -1; + component_is_before_row = nrows; + break; + default: + // DO nothing + break; + } + } + if (last == LEX) + { + // last block was lex, so use lex tie-breaker + mat.resize(last_element); + if (nrows == component_is_before_row) component_is_before_row = -1; + base_is_revlex = false; + } + else if (last == REVLEX) + { + // last block was revlex, so use revlex tie-breaker + if (nrows == component_is_before_row) component_is_before_row = -1; + mat.resize(last_element); + } + else + { + // last block is a weight vector, so use revlex as the tie-breaker. + // nothing to change here. + } + return true; +} + +int moIsLex(const MonomialOrdering *mo) +{ + // The monomial order is lex if what? + // one lex block, no grevlex blocks, no weightvector blocks. + // only: lex block and position blocks are allowed. + int nlex = 0; + int i; + for (i = 0; i < mo->len; i++) + { + enum MonomialOrdering_type typ = mo->array[i]->type; + switch (typ) + { + case MO_LEX: + case MO_LEX2: + case MO_LEX4: + nlex++; + break; + case MO_POSITION_UP: + case MO_POSITION_DOWN: + break; + default: + return 0; + } + } + return (nlex == 1); +} + +int moIsGRevLex(const MonomialOrdering *mo) +{ + // The monomial order is lex if what? + // one lex block, no grevlex blocks, no weightvector blocks. + // only: lex block and position blocks are allowed. + int ngrevlex = 0; + int i; + for (i = 0; i < mo->len; i++) + { + enum MonomialOrdering_type typ = mo->array[i]->type; + switch (typ) + { + case MO_GREVLEX: + case MO_GREVLEX2: + case MO_GREVLEX4: + case MO_GREVLEX_WTS: + case MO_GREVLEX2_WTS: + case MO_GREVLEX4_WTS: + ngrevlex++; + break; + case MO_POSITION_UP: + case MO_POSITION_DOWN: + break; + default: + return 0; + } + } + return (ngrevlex == 1); +} + +int rawNumberOfVariables(const MonomialOrdering *mo) +{ + int i, sum = 0; + for (i = 0; i < mo->len; i++) + if (mo->array[i]->type != MO_WEIGHTS) sum += mo->array[i]->nvars; + return sum; +} + +M2_arrayint moGetWeightValues(const MonomialOrdering *mo) +{ + int nvars = rawNumberOfVariables(mo); + // grab the first weight vector + if (mo->len == 0) return nullptr; + if (mo->array[0]->type == MO_WEIGHTS) + { + int i; + M2_arrayint result = M2_makearrayint(nvars); + int *wts = mo->array[0]->wts; + for (i = 0; i < mo->array[0]->nvars; i++) result->array[i] = wts[i]; + for (; i < nvars; i++) result->array[i] = 0; + return result; + } + return nullptr; +} + +int rawNumberOfInvertibleVariables(const MonomialOrdering *mo) +{ + int i, sum = 0; + for (i = 0; i < mo->len; i++) + if (mo->array[i]->type == MO_LAURENT || + mo->array[i]->type == MO_LAURENT_REVLEX) + sum += mo->array[i]->nvars; + return sum; +} + +M2_arrayint rawNonTermOrderVariables(const MonomialOrdering *mo) +// returns a list of the indices of those variables which are less than 1 in +// the given monomial order. +{ + int i, j, sum, nextvar; + int nvars = rawNumberOfVariables(mo); + int *gt = getmematomicvectortype(int, nvars); + for (i = 0; i < nvars; i++) + gt[i] = + 0; // 0 means undecided, -1 means non term order, 1 means term order + // Now we loop through the parts of the monomial order + nextvar = 0; + for (i = 0; i < mo->len; i++) + { + mon_part p = mo->array[i]; + switch (p->type) + { + case MO_LEX: + case MO_LEX2: + case MO_LEX4: + case MO_GREVLEX: + case MO_GREVLEX2: + case MO_GREVLEX4: + case MO_GREVLEX_WTS: + case MO_GREVLEX2_WTS: + case MO_GREVLEX4_WTS: + case MO_LAURENT: + case MO_NC_LEX: + for (j = 0; j < p->nvars; j++, nextvar++) + if (gt[nextvar] == 0) gt[nextvar] = 1; + break; + case MO_LAURENT_REVLEX: + case MO_REVLEX: + for (j = 0; j < p->nvars; j++, nextvar++) + if (gt[nextvar] == 0) gt[nextvar] = -1; + break; + case MO_WEIGHTS: + for (j = nextvar; j < p->nvars; j++) + if (gt[j] == 0) + { + if (p->wts[j] > 0) + gt[j] = 1; + else if (p->wts[j] < 0) + gt[j] = -1; + } + break; + case MO_POSITION_UP: + case MO_POSITION_DOWN: + break; + } + } + // At this point every variables' gt should be 1 or -1. + sum = 0; + for (i = 0; i < nvars; i++) + { + if (gt[i] == 0) INTERNAL_ERROR("gt[i] should not be 0"); + if (gt[i] < 0) sum++; + } + + // Make an array of this length. + M2_arrayint result = M2_makearrayint(sum); + nextvar = 0; + for (i = 0; i < nvars; i++) + if (gt[i] < 0) result->array[nextvar++] = i; + return result; +} + +MonomialOrdering *rawLexMonomialOrdering(int nvars, int packing) +{ + MonomialOrdering *result; + mon_part p; + enum MonomialOrdering_type typ; + + if (packing == 2) + typ = MO_LEX2; + else if (packing == 4) + typ = MO_LEX4; + else + typ = MO_LEX; + + p = mo_make(typ, nvars, nullptr); + result = make_mon_order(1); + result->array[0] = p; + return result; +} + +MonomialOrdering /* or null */ *rawGRevLexMonomialOrdering(M2_arrayint degs, + int packing) +{ + MonomialOrdering *result; + mon_part p; + enum MonomialOrdering_type typ; + int *wts; + int all_one = 1; + unsigned int i; + for (i = 0; i < degs->len; i++) + if (degs->array[i] <= 0) + { + ERROR("grevlex: expected all degrees to be positive"); + return nullptr; + } + else if (degs->array[i] > 1) + all_one = 0; + + if (all_one) + { + if (packing == 2) + typ = MO_GREVLEX2; + else if (packing == 4) + typ = MO_GREVLEX4; + else + typ = MO_GREVLEX; + wts = nullptr; + } + else + { + if (packing == 2) + typ = MO_GREVLEX2_WTS; + else if (packing == 4) + typ = MO_GREVLEX4_WTS; + else + typ = MO_GREVLEX_WTS; + wts = degs->array; + } + + p = mo_make(typ, degs->len, wts); + result = make_mon_order(1); + result->array[0] = p; + return result; +} + +MonomialOrdering *rawRevLexMonomialOrdering(int nvars) +{ + mon_part p = mo_make(MO_REVLEX, nvars, nullptr); + MonomialOrdering *result = make_mon_order(1); + result->array[0] = p; + return result; +} + +MonomialOrdering *rawWeightsMonomialOrdering(M2_arrayint wts) +{ + mon_part p = mo_make(MO_WEIGHTS, wts->len, wts->array); + MonomialOrdering *result = make_mon_order(1); + result->array[0] = p; + return result; +} +MonomialOrdering *rawGroupLexMonomialOrdering(int nvars) +{ + mon_part p = mo_make(MO_LAURENT, nvars, nullptr); + MonomialOrdering *result = make_mon_order(1); + result->array[0] = p; + return result; +} +MonomialOrdering *rawGroupRevLexMonomialOrdering(int nvars) +{ + mon_part p = mo_make(MO_LAURENT_REVLEX, nvars, nullptr); + MonomialOrdering *result = make_mon_order(1); + result->array[0] = p; + return result; +} +MonomialOrdering *rawNClexMonomialOrdering(int nvars) +{ + mon_part p = mo_make(MO_NC_LEX, nvars, nullptr); + MonomialOrdering *result = make_mon_order(1); + result->array[0] = p; + return result; +} +MonomialOrdering *rawPositionMonomialOrdering(M2_bool up_or_down) +{ + mon_part p = + mo_make((up_or_down ? MO_POSITION_UP : MO_POSITION_DOWN), 0, nullptr); + MonomialOrdering *result = make_mon_order(1); + result->array[0] = p; + return result; +} + +MonomialOrdering *rawJoinMonomialOrdering(engine_RawMonomialOrderingArray M) +{ + MonomialOrdering *result; + const MonomialOrdering *mo; + int i, j, sum, next; + int nvars_so_far = 0; + + // sum = 0; + // for (i=0; ilen; i++) + // sum += M->array[i]->len; + + sum = 0; + for (i = 0; i < M->len; i++) + { + mo = M->array[i]; + for (j = 0; j < mo->len; j++) + if (is_good(mo->array[j])) sum++; + } + + result = make_mon_order(sum); + next = 0; + for (i = 0; i < M->len; i++) + { + mo = M->array[i]; + for (j = 0; j < mo->len; j++) + { + mon_part p = mo->array[j]; + if (!is_good(p)) continue; + if (p->type != MO_WEIGHTS) + nvars_so_far += p->nvars; + else + { + /* Shift the weights over by nvars_so_far */ + mon_part q = mo_make(MO_WEIGHTS, nvars_so_far + p->nvars, nullptr); + q->wts = getmemvectortype(int, q->nvars); + for (j = 0; j < nvars_so_far; j++) q->wts[j] = 0; + for (; j < q->nvars; j++) q->wts[j] = p->wts[j - nvars_so_far]; + p = q; + } + result->array[next++] = p; + } + } + return result; +} + +MonomialOrdering *rawProductMonomialOrdering(engine_RawMonomialOrderingArray M) +{ + MonomialOrdering *result; + int i, j, sum, next, offset; + sum = 0; + for (i = 0; i < M->len; i++) sum += M->array[i]->len; + + result = make_mon_order(sum); + next = 0; + offset = 0; + for (i = 0; i < M->len; i++) + { + int nvars = rawNumberOfVariables(M->array[i]); + MonomialOrdering *mo = M2_mo_offset(M->array[i], offset); + for (j = 0; j < mo->len; j++) result->array[next++] = mo->array[j]; + offset += nvars; + } + return result; +} + +M2_string intarray_to_string(int len, int *p) +{ + int i; + char s[200]; + M2_string result = M2_tostring("{"); + for (i = 0; i < len; i++) + { + if (i > 0) result = M2_join(result, M2_tostring(",")); + sprintf(s, "%d", p[i]); + result = M2_join(result, M2_tostring(s)); + } + result = M2_join(result, M2_tostring("}")); + return result; +} + +M2_string ones_to_string(int len) +{ + int i; + char s[200]; + M2_string one; + M2_string result = M2_tostring("{"); + sprintf(s, "1"); + one = M2_tostring(s); + for (i = 0; i < len; i++) + { + if (i > 0) result = M2_join(result, M2_tostring(",")); + result = M2_join(result, one); + } + result = M2_join(result, M2_tostring("}")); + return result; +} + +unsigned int rawMonomialOrderingHash(const MonomialOrdering *mo) +{ + return mo->_hash; +} + +M2_string IM2_MonomialOrdering_to_string(const MonomialOrdering *mo) +{ + int i; + char s[200]; + int p_ones = 0; + M2_string result = M2_tostring("MonomialOrder => {"); + for (i = 0; i < mo->len; i++) + { + mon_part p = mo->array[i]; + p_ones = 0; + if (i == 0) + result = M2_join(result, M2_tostring("\n ")); + else + result = M2_join(result, M2_tostring(",\n ")); + switch (p->type) + { + case MO_LEX: + sprintf(s, "Lex => %d", p->nvars); + break; + case MO_LEX2: + sprintf(s, "LexSmall => %d", p->nvars); + break; + case MO_LEX4: + sprintf(s, "LexTiny => %d", p->nvars); + break; + case MO_GREVLEX: + sprintf(s, "GRevLex => "); + p_ones = 1; + break; + case MO_GREVLEX2: + sprintf(s, "GRevLexSmall => "); + p_ones = 1; + break; + case MO_GREVLEX4: + sprintf(s, "GRevLexTiny => "); + p_ones = 1; + break; + case MO_GREVLEX_WTS: + sprintf(s, "GRevLex => "); + break; + case MO_GREVLEX2_WTS: + sprintf(s, "GRevLexSmall => "); + break; + case MO_GREVLEX4_WTS: + sprintf(s, "GRevLexTiny => "); + break; + case MO_REVLEX: + sprintf(s, "RevLex => %d", p->nvars); + break; + case MO_WEIGHTS: + sprintf(s, "Weights => "); + break; + case MO_LAURENT: + sprintf(s, "GroupLex => %d", p->nvars); + break; + case MO_LAURENT_REVLEX: + sprintf(s, "GroupRevLex => %d", p->nvars); + break; + case MO_NC_LEX: + sprintf(s, "NCLex => %d", p->nvars); + break; + case MO_POSITION_UP: + sprintf(s, "Position => Up"); + break; + case MO_POSITION_DOWN: + sprintf(s, "Position => Down"); + break; + default: + sprintf(s, "UNKNOWN"); + break; + } + result = M2_join(result, M2_tostring(s)); + if (p->wts != nullptr) + result = M2_join(result, intarray_to_string(p->nvars, p->wts)); + else if (p_ones) + result = M2_join(result, ones_to_string(p->nvars)); + } + result = M2_join(result, M2_tostring("\n }")); + return result; +} + +#include "overflow.hpp" + +std::vector laurentVariables(const MonomialOrder* mo) +{ + std::vector result; + for (auto i = 0; i < mo->nvars; ++i) + result.push_back(mo->is_laurent[i] == 1); + return result; +} +/* TODO: + -- negative exponent versions need to be included (at least for MO_LEX) + -- non-commutative blocks should be added in +*/ + +static void mo_block_revlex(struct mo_block *b, int nvars) +{ + b->typ = MO_REVLEX; + b->nvars = nvars; + b->nslots = nvars; + b->first_exp = 0; /* will be set later */ + b->first_slot = 0; /* will be set later */ + b->nweights = 0; + b->weights = nullptr; +} + +static void mo_block_grevlex(struct mo_block *b, int nvars) +{ + b->typ = MO_GREVLEX; + b->nvars = nvars; + b->nslots = nvars; + b->first_exp = 0; /* will be set later */ + b->first_slot = 0; /* will be set later */ + b->nweights = 0; + b->weights = nullptr; +} + +static void mo_block_grevlex2(struct mo_block *b, int nvars) +{ + b->typ = MO_GREVLEX2; + b->nvars = nvars; + b->nslots = (nvars + 1) / 2; /* 2 per word */ + b->first_exp = 0; /* will be set later */ + b->first_slot = 0; /* will be set later */ + b->nweights = 0; + b->weights = nullptr; +} + +static void mo_block_grevlex4(struct mo_block *b, int nvars) +{ + b->typ = MO_GREVLEX4; + b->nvars = nvars; + b->nslots = (nvars + 3) / 4; /* 4 per word */ + b->first_exp = 0; /* will be set later */ + b->first_slot = 0; /* will be set later */ + b->nweights = 0; + b->weights = nullptr; +} + +static void mo_block_grevlex_wts(struct mo_block *b, int nvars) +{ + b->typ = MO_GREVLEX_WTS; + b->nvars = nvars; + b->nslots = nvars; + b->first_exp = 0; /* will be set later */ + b->first_slot = 0; /* will be set later */ + b->nweights = nvars; + b->weights = nullptr; /* will be set later */ +} + +static void mo_block_grevlex2_wts(struct mo_block *b, int nvars) +{ + b->typ = MO_GREVLEX2_WTS; + b->nvars = nvars; + b->nslots = (nvars + 1) / 2; /* 2 per word */ + b->first_exp = 0; /* will be set later */ + b->first_slot = 0; /* will be set later */ + b->nweights = nvars; + b->weights = nullptr; /* will be set later */ +} + +static void mo_block_grevlex4_wts(struct mo_block *b, int nvars) +{ + b->typ = MO_GREVLEX4_WTS; + b->nvars = nvars; + b->nslots = (nvars + 3) / 4; /* 4 per word */ + b->first_exp = 0; /* will be set later */ + b->first_slot = 0; /* will be set later */ + b->nweights = nvars; + b->weights = nullptr; /* will be set later */ +} + +static void mo_block_lex(struct mo_block *b, int nvars) +{ + b->typ = MO_LEX; + b->nvars = nvars; + b->nslots = nvars; + b->first_exp = 0; /* will be set later */ + b->first_slot = 0; /* will be set later */ + b->nweights = 0; + b->weights = nullptr; +} + +static void mo_block_lex2(struct mo_block *b, int nvars) +{ + b->typ = MO_LEX2; + b->nvars = nvars; + b->nslots = (nvars + 1) / 2; /* 2 per word */ + b->first_exp = 0; /* will be set later */ + b->first_slot = 0; /* will be set later */ + b->nweights = 0; + b->weights = nullptr; +} + +static void mo_block_lex4(struct mo_block *b, int nvars) +{ + b->typ = MO_LEX4; + b->nvars = nvars; + b->nslots = (nvars + 3) / 4; /* 4 per word */ + b->first_exp = 0; /* will be set later */ + b->first_slot = 0; /* will be set later */ + b->nweights = 0; + b->weights = nullptr; +} + +static void mo_block_group_lex(struct mo_block *b, int nvars) +{ + b->typ = MO_LAURENT; + b->nvars = nvars; + b->nslots = nvars; + b->first_exp = 0; /* will be set later */ + b->first_slot = 0; /* will be set later */ + b->nweights = 0; + b->weights = nullptr; +} + +static void mo_block_group_revlex(struct mo_block *b, int nvars) +{ + b->typ = MO_LAURENT_REVLEX; + b->nvars = nvars; + b->nslots = nvars; + b->first_exp = 0; /* will be set later */ + b->first_slot = 0; /* will be set later */ + b->nweights = 0; + b->weights = nullptr; +} + +static void mo_block_wt_function(struct mo_block *b, int nvars, deg_t *wts) +{ + b->typ = MO_WEIGHTS; + b->nvars = 0; + b->nslots = 1; + b->first_exp = 0; + b->first_slot = 0; /* will be set later */ + b->nweights = nvars; + b->weights = wts; +} + +MonomialOrder *monomialOrderMake(const MonomialOrdering *mo) +{ + MonomialOrder *result; + int i, j, nv, this_block; + deg_t *wts = nullptr; + /* Determine the number of variables, the number of blocks, and the location + of the component */ + int nblocks = 0; + int nvars = 0; + int hascomponent = 0; + for (i = 0; i < mo->len; i++) + { + struct mon_part_rec_ *t = mo->array[i]; + nblocks++; + if (t->type == MO_POSITION_DOWN || t->type == MO_POSITION_UP) + hascomponent++; + else if (t->type == MO_NC_LEX) + { + // Currently, do nothing. + } + if (t->type != MO_WEIGHTS) nvars += t->nvars; + } + nblocks -= hascomponent; + + /* Now create the blocks, and fill them in. Also fill in the deg vector */ + result = getmemstructtype(MonomialOrder *); + result->nvars = nvars; + result->nslots = 0; + result->nblocks = nblocks; + result->blocks = + (struct mo_block *)getmem(nblocks * sizeof(result->blocks[0])); + result->degs = (deg_t *)getmem_atomic(nvars * sizeof(result->degs[0])); + if (hascomponent == 0) result->nblocks_before_component = nblocks; + + this_block = 0; + nvars = 0; + for (i = 0; i < mo->len; i++) + { + struct mon_part_rec_ *t = mo->array[i]; + if (t->type != MO_WEIGHTS) + { + if (t->wts == nullptr) + for (j = 0; j < t->nvars; j++) result->degs[nvars++] = 1; + else + for (j = 0; j < t->nvars; j++) result->degs[nvars++] = t->wts[j]; + } + else + { + wts = (deg_t *)getmem_atomic(t->nvars * sizeof(wts[0])); + for (j = 0; j < t->nvars; j++) wts[j] = t->wts[j]; + } + switch (t->type) + { + case MO_REVLEX: + mo_block_revlex(result->blocks + this_block++, t->nvars); + break; + case MO_GREVLEX: + mo_block_grevlex(result->blocks + this_block++, t->nvars); + break; + case MO_GREVLEX2: + mo_block_grevlex2(result->blocks + this_block++, t->nvars); + break; + case MO_GREVLEX4: + mo_block_grevlex4(result->blocks + this_block++, t->nvars); + break; + case MO_GREVLEX_WTS: + mo_block_grevlex_wts(result->blocks + this_block++, t->nvars); + break; + case MO_GREVLEX2_WTS: + mo_block_grevlex2_wts(result->blocks + this_block++, t->nvars); + break; + case MO_GREVLEX4_WTS: + mo_block_grevlex4_wts(result->blocks + this_block++, t->nvars); + break; + case MO_LEX: + mo_block_lex(result->blocks + this_block++, t->nvars); + break; + case MO_LEX2: + mo_block_lex2(result->blocks + this_block++, t->nvars); + break; + case MO_LEX4: + mo_block_lex4(result->blocks + this_block++, t->nvars); + break; + case MO_WEIGHTS: + // if extra weight values are given (more than "nvars", ignore the + // rest. + mo_block_wt_function( + result->blocks + this_block++, + (t->nvars <= result->nvars ? t->nvars : result->nvars), + wts); + break; + case MO_LAURENT: + mo_block_group_lex(result->blocks + this_block++, t->nvars); + break; + case MO_LAURENT_REVLEX: + mo_block_group_revlex(result->blocks + this_block++, t->nvars); + break; + case MO_NC_LEX: + /* MES */ + break; + case MO_POSITION_UP: + if (--hascomponent == 0) + { + // Set the information about the component + result->component_up = 1; + result->nblocks_before_component = this_block; + } + // mo_block_position_up(result->blocks + this_block); + break; + case MO_POSITION_DOWN: + if (--hascomponent == 0) + { + // Set the information about the component + result->component_up = 0; + result->nblocks_before_component = this_block; + } + // mo_block_position_down(result->blocks + this_block); + break; + } + } + + /* Go back and fill in the 'slots' information */ + /* Now fix the first_exp, first_slot values, and also result->{nslots,nvars}; + */ + nv = 0; + result->nslots = 0; + result->nslots_before_component = 0; + for (i = 0; i < nblocks; i++) + { + enum MonomialOrdering_type typ = result->blocks[i].typ; + + result->blocks[i].first_exp = nv; + result->blocks[i].first_slot = result->nslots; + nv += result->blocks[i].nvars; + result->nslots += result->blocks[i].nslots; + + if (typ == MO_WEIGHTS) + { + result->blocks[i].first_exp = 0; + + /* divide the wt vector by the degree vector */ + for (j = 0; j < result->blocks[i].nvars; j++) + safe::div_by(result->blocks[i].weights[j], result->degs[j]); + ; + } + else if (typ == MO_GREVLEX_WTS || typ == MO_GREVLEX2_WTS || + typ == MO_GREVLEX4_WTS) + { + result->blocks[i].weights = + result->degs + result->blocks[i].first_exp; + } + + if (i == result->nblocks_before_component - 1) + { + result->nslots_before_component = result->nslots; + } + } + + /* Set is_laurent */ + result->is_laurent = (int *)getmem_atomic(result->nvars * sizeof(int)); + for (i = 0; i < result->nvars; i++) result->is_laurent[i] = 0; + + for (i = 0; i < result->nblocks; i++) + if (result->blocks[i].typ == MO_LAURENT || + result->blocks[i].typ == MO_LAURENT_REVLEX) + { + for (j = 0; j < result->blocks[i].nvars; j++) + result->is_laurent[result->blocks[i].first_exp + j] = 1; + } + + return result; +} + +extern void monomialOrderFree(MonomialOrder *mo) {} +static void MO_pack4(int nvars, const int *expon, int *slots) +{ + int32_t i; + if (nvars == 0) return; + while (1) + { + i = safe::fits_7(*expon++) << 24; + if (--nvars == 0) break; + i |= safe::fits_7(*expon++) << 16; + if (--nvars == 0) break; + i |= safe::fits_7(*expon++) << 8; + if (--nvars == 0) break; + i |= safe::fits_7(*expon++); + if (--nvars == 0) break; + *slots++ = i; + } + *slots++ = i; +} + +static void MO_pack2(int nvars, const int *expon, int *slots) +{ + int32_t i; + if (nvars == 0) return; + while (1) + { + i = safe::fits_15(*expon++) << 16; + if (--nvars == 0) break; + i |= safe::fits_15(*expon++); + if (--nvars == 0) break; + *slots++ = i; + } + *slots++ = i; +} + +static void MO_unpack4(int nvars, const int *slots, int *expon) +{ + int32_t i; + if (nvars == 0) return; + while (1) + { + i = *slots++; + *expon++ = (i >> 24); + if (--nvars == 0) break; + *expon++ = (i >> 16) & 0x7f; + if (--nvars == 0) break; + *expon++ = (i >> 8) & 0x7f; + if (--nvars == 0) break; + *expon++ = i & 0x7f; + if (--nvars == 0) break; + } +} + +static void MO_unpack2(int nvars, const int *slots, int *expon) +{ + int32_t i; + if (nvars == 0) return; + while (1) + { + i = *slots++; + *expon++ = i >> 16; + if (--nvars == 0) break; + *expon++ = i & 0x7fff; + if (--nvars == 0) break; + } +} + +std::vector MonomialOrder::overflow_flags() const +{ + std::vector overflow; + enum overflow_type flag; + int i = 0; + int k = 0; + for (; i < nblocks; i++) + { + mo_block *b = &blocks[i]; + switch (b->typ) + { + case MO_REVLEX: + case MO_WEIGHTS: + case MO_LAURENT: + case MO_LAURENT_REVLEX: + case MO_NC_LEX: + flag = OVER; + goto fillin; + case MO_POSITION_UP: + case MO_POSITION_DOWN: + ERROR( + "internal error - MO_POSITION_DOWN or MO_POSITION_UP " + "encountered"); + assert(0); + break; + case MO_LEX: + case MO_GREVLEX: + case MO_GREVLEX_WTS: + flag = OVER1; + goto fillin; + case MO_LEX2: + case MO_GREVLEX2: + case MO_GREVLEX2_WTS: + flag = OVER2; + goto fillin; + case MO_LEX4: + case MO_GREVLEX4: + case MO_GREVLEX4_WTS: + flag = OVER4; + goto fillin; + fillin: + assert(b->first_slot == k); + for (int p = b->nslots; p > 0; p--) + { + assert(k < nslots); + overflow.push_back(flag); + k++; + } + break; + default: + ERROR("internal error - missing case"); + assert(0); + break; + } + } + assert(k == nslots); + return overflow; +} + + + +void monomialOrderEncodeFromActualExponents(const MonomialOrder *mo, + const_exponents expon, + monomial result_psums) +/* Given 'expon', compute the encoded partial sums value */ +{ + if (mo == nullptr) return; + int *tmpexp = static_cast(alloca((mo->nvars + 1) * sizeof(int))); + int i, j, nvars, s; + int *p1; + deg_t *degs; + struct mo_block *b = mo->blocks; + int nblocks = mo->nblocks; + const int *e = expon; + int *p = result_psums; + for (i = nblocks; i > 0; --i, b++) switch (b->typ) + { + case MO_LEX: + case MO_LAURENT: + nvars = b->nvars; + for (j = 0; j < nvars; j++) *p++ = *e++; + break; + case MO_REVLEX: + case MO_LAURENT_REVLEX: + nvars = b->nvars; + for (j = 0; j < nvars; j++) *p++ = safe::minus(*e++); + break; + case MO_GREVLEX: + nvars = b->nvars; + p += b->nslots; + p1 = p; + *--p1 = *e++; + for (j = 1; j < nvars; j++) + { + --p1; + *p1 = safe::add(*e++, p1[1]); + } + break; + case MO_GREVLEX_WTS: + nvars = b->nvars; + degs = mo->degs + b->first_exp; + p += b->nslots; + p1 = p; + *--p1 = safe::mult(*e++, *degs++); + for (j = 1; j < nvars; j++) + { + --p1; + int tmp = safe::mult(*e++, *degs++); + *p1 = safe::add(tmp, p1[1]); + } + break; + case MO_GREVLEX4: + nvars = b->nvars; + p1 = tmpexp + b->nvars; + *--p1 = *e++; + for (j = 1; j < nvars; j++) + { + --p1; + *p1 = safe::add(*e++, p1[1]); + } + MO_pack4(nvars, p1, p); + p += b->nslots; + break; + case MO_GREVLEX4_WTS: + nvars = b->nvars; + degs = mo->degs + b->first_exp; + p1 = tmpexp + b->nvars; + *--p1 = safe::mult(*e++, *degs++); + for (j = 1; j < nvars; j++) + { + --p1; + int tmp = safe::mult(*e++, *degs++); + *p1 = safe::add(tmp, p1[1]); + } + MO_pack4(nvars, p1, p); + p += b->nslots; + break; + case MO_GREVLEX2: + nvars = b->nvars; + p1 = tmpexp + b->nvars; + *--p1 = *e++; + for (j = 1; j < nvars; j++) + { + --p1; + *p1 = safe::add(*e++, p1[1]); + } + MO_pack2(nvars, p1, p); + p += b->nslots; + break; + case MO_GREVLEX2_WTS: + nvars = b->nvars; + degs = mo->degs + b->first_exp; + p1 = tmpexp + b->nvars; + *--p1 = safe::mult(*e++, *degs++); + for (j = 1; j < nvars; j++) + { + --p1; + int tmp = safe::mult(*e++, *degs++); + *p1 = safe::add(tmp, p1[1]); + } + MO_pack2(nvars, p1, p); + p += b->nslots; + break; + case MO_LEX4: + nvars = b->nvars; + MO_pack4(nvars, e, p); + p += b->nslots; + e += nvars; + break; + case MO_LEX2: + nvars = b->nvars; + MO_pack2(nvars, e, p); + p += b->nslots; + e += nvars; + break; + case MO_WEIGHTS: + if (b->nweights == 0) + { + s = 0; + } + else + { + s = safe::mult(b->weights[0], expon[0]); + for (j = 1; j < b->nweights; j++) + s = safe::add(s, safe::mult(b->weights[j], expon[j])); + } + *p++ = s; + break; + case MO_POSITION_UP: + case MO_POSITION_DOWN: + /* nothing to do here */ + break; + case MO_NC_LEX: + /* nothing to do here */ + break; + } +} + +void monomialOrderDecodeToActualExponents(const MonomialOrder *mo, + const_monomial psums, + exponents_t expon) +{ + if (mo == nullptr) return; + int *tmpexp = static_cast(alloca((mo->nvars + 1) * sizeof(int))); + int i, j, nvars; + deg_t *degs = mo->degs; + deg_t *d; + struct mo_block *b = mo->blocks; + int nblocks = mo->nblocks; + int *e = expon; + const int *p = psums; + for (i = nblocks; i > 0; --i, b++) switch (b->typ) + { + case MO_LEX: + case MO_LAURENT: + nvars = b->nvars; + p = psums + b->first_slot; + e = expon + b->first_exp; + for (j = 0; j < nvars; j++) *e++ = *p++; + break; + case MO_REVLEX: + case MO_LAURENT_REVLEX: + nvars = b->nvars; + p = psums + b->first_slot; + e = expon + b->first_exp; + for (j = 0; j < nvars; j++) *e++ = safe::minus(*p++); + break; + case MO_GREVLEX: + nvars = b->nvars; + p = psums + b->first_slot + nvars - 1; + e = expon + b->first_exp; + *e++ = *p--; + for (j = nvars - 1; j >= 1; --j, --p) *e++ = safe::sub(*p, p[1]); + break; + case MO_GREVLEX_WTS: + nvars = b->nvars; + d = degs + b->first_exp; + p = psums + b->first_slot + nvars - 1; + e = expon + b->first_exp; + *e++ = *p-- / *d++; + for (j = nvars - 1; j >= 1; --j, --p) + *e++ = safe::sub(*p, p[1]) / *d++; + break; + case MO_GREVLEX4: + nvars = b->nvars; + MO_unpack4(nvars, psums + b->first_slot, tmpexp); + p = tmpexp + nvars - 1; + e = expon + b->first_exp; + *e++ = *p--; + for (j = nvars - 1; j >= 1; --j, --p) *e++ = safe::sub(*p, p[1]); + break; + case MO_GREVLEX4_WTS: + nvars = b->nvars; + d = degs + b->first_exp; + MO_unpack4(nvars, psums + b->first_slot, tmpexp); + p = tmpexp + nvars - 1; + e = expon + b->first_exp; + *e++ = *p-- / *d++; + for (j = nvars - 1; j >= 1; --j, --p) + *e++ = safe::sub(*p, p[1]) / *d++; + break; + case MO_GREVLEX2: + nvars = b->nvars; + MO_unpack2(nvars, psums + b->first_slot, tmpexp); + p = tmpexp + nvars - 1; + e = expon + b->first_exp; + *e++ = *p--; + for (j = nvars - 1; j >= 1; --j, --p) *e++ = safe::sub(*p, p[1]); + break; + case MO_GREVLEX2_WTS: + nvars = b->nvars; + d = degs + b->first_exp; + MO_unpack2(nvars, psums + b->first_slot, tmpexp); + p = tmpexp + nvars - 1; + e = expon + b->first_exp; + *e++ = *p-- / *d++; + for (j = nvars - 1; j >= 1; --j, --p) + *e++ = safe::sub(*p, p[1]) / *d++; + break; + case MO_LEX4: + nvars = b->nvars; + e = expon + b->first_exp; + MO_unpack4(nvars, psums + b->first_slot, e); + break; + case MO_LEX2: + nvars = b->nvars; + e = expon + b->first_exp; + MO_unpack2(nvars, psums + b->first_slot, e); + break; + case MO_WEIGHTS: + break; + case MO_POSITION_UP: + case MO_POSITION_DOWN: + /* should not occur, but do nothing in any case */ + break; + case MO_NC_LEX: + /* nothing to do here */ + break; + } +} + +MonomialOrdering *MonomialOrderings::Lex(int nvars) +{ + return rawLexMonomialOrdering(nvars, 1); +} +MonomialOrdering *MonomialOrderings::Lex2(int nvars) +{ + return rawLexMonomialOrdering(nvars, 2); +} +MonomialOrdering *MonomialOrderings::Lex4(int nvars) +{ + return rawLexMonomialOrdering(nvars, 4); +} +MonomialOrdering *MonomialOrderings::GRevLex(int nvars) +{ + std::vector w; + for (int i = 0; i < nvars; i++) w.push_back(1); + return GRevLex(w); +} +MonomialOrdering *MonomialOrderings::GRevLex2(int nvars) +{ + std::vector w; + for (int i = 0; i < nvars; i++) w.push_back(1); + return GRevLex2(w); +} +MonomialOrdering *MonomialOrderings::GRevLex4(int nvars) +{ + std::vector w; + for (int i = 0; i < nvars; i++) w.push_back(1); + return GRevLex4(w); +} +MonomialOrdering *MonomialOrderings::GRevLex(const std::vector &wts) +{ + return GRevLex(wts, 1); +} +MonomialOrdering *MonomialOrderings::GRevLex2(const std::vector &wts) +{ + return GRevLex(wts, 2); +} +MonomialOrdering *MonomialOrderings::GRevLex4(const std::vector &wts) +{ + return GRevLex(wts, 4); +} +MonomialOrdering *MonomialOrderings::RevLex(int nvars) +{ + return rawRevLexMonomialOrdering(nvars); +} +MonomialOrdering *MonomialOrderings::Weights(const std::vector &wts) +{ + mon_part p = mo_make(MO_WEIGHTS, wts.size(), wts.data()); + MonomialOrdering *result = make_mon_order(1); + result->array[0] = p; + return result; +} +MonomialOrdering *MonomialOrderings::GroupLex(int nvars) +{ + return rawGroupLexMonomialOrdering(nvars); +} +MonomialOrdering *MonomialOrderings::GroupRevLex(int nvars) +{ + return rawGroupRevLexMonomialOrdering(nvars); +} +MonomialOrdering *MonomialOrderings::PositionUp() +{ + return rawPositionMonomialOrdering(true); +} +MonomialOrdering *MonomialOrderings::PositionDown() +{ + return rawPositionMonomialOrdering(false); +} + +MonomialOrdering *MonomialOrderings::GRevLex(const std::vector °s, + int packing) +{ + MonomialOrdering *result; + mon_part p; + enum MonomialOrdering_type typ; + const int *wts; + bool all_one = true; + for (int i = 0; i < degs.size(); i++) + if (degs[i] <= 0) + { + ERROR("grevlex: expected all degrees to be positive"); + return nullptr; + } + else if (degs[i] > 1) + all_one = false; + + if (all_one) + { + if (packing == 2) + typ = MO_GREVLEX2; + else if (packing == 4) + typ = MO_GREVLEX4; + else + typ = MO_GREVLEX; + wts = nullptr; + } + else + { + if (packing == 2) + typ = MO_GREVLEX2_WTS; + else if (packing == 4) + typ = MO_GREVLEX4_WTS; + else + typ = MO_GREVLEX_WTS; + wts = degs.data(); + } + + p = mo_make(typ, degs.size(), wts); + result = make_mon_order(1); + result->array[0] = p; + return result; +} + +MonomialOrdering *MonomialOrderings::join( + const std::vector &M) +{ + MonomialOrdering *result; + const MonomialOrdering *mo; + int i, j, sum, next; + int nvars_so_far = 0; + + sum = 0; + for (i = 0; i < M.size(); i++) + { + mo = M[i]; + for (j = 0; j < mo->len; j++) + if (is_good(mo->array[j])) sum++; + } + + result = make_mon_order(sum); + next = 0; + for (i = 0; i < M.size(); i++) + { + mo = M[i]; + for (j = 0; j < mo->len; j++) + { + mon_part p = mo->array[j]; + if (!is_good(p)) continue; + if (p->type != MO_WEIGHTS) + nvars_so_far += p->nvars; + else + { + /* Shift the weights over by nvars_so_far */ + mon_part q = mo_make(MO_WEIGHTS, nvars_so_far + p->nvars, nullptr); + q->wts = getmemvectortype(int, q->nvars); + for (j = 0; j < nvars_so_far; j++) q->wts[j] = 0; + for (; j < q->nvars; j++) q->wts[j] = p->wts[j - nvars_so_far]; + p = q; + } + result->array[next++] = p; + } + } + return result; +} + +MonomialOrdering *MonomialOrderings::product( + const std::vector &M) +{ + MonomialOrdering *result; + int i, j, sum, next, offset; + sum = 0; + for (i = 0; i < M.size(); i++) sum += M[i]->len; + + result = make_mon_order(sum); + next = 0; + offset = 0; + for (i = 0; i < M.size(); i++) + { + int nvars = rawNumberOfVariables(M[i]); + MonomialOrdering *mo = M2_mo_offset(M[i], offset); + for (j = 0; j < mo->len; j++) result->array[next++] = mo->array[j]; + offset += nvars; + } + return result; +} + + +// Local Variables: +// compile-command: "make -C $M2BUILDDIR/Macaulay2/e " +// indent-tabs-mode: nil +// End: + diff --git a/M2/Macaulay2/e/MonomialOrderOld.hpp b/M2/Macaulay2/e/MonomialOrderOld.hpp new file mode 100644 index 00000000000..b24ee2e2c58 --- /dev/null +++ b/M2/Macaulay2/e/MonomialOrderOld.hpp @@ -0,0 +1,123 @@ +#pragma once + +#include +#include +#include + +#include "interface/monomial-ordering.h" + +#include "ExponentVector.hpp" + +typedef struct mon_part_rec_ +{ + enum MonomialOrdering_type type; + int nvars; + int *wts; +} * mon_part; + +struct MonomialOrdering +{ + unsigned int _hash; + unsigned int len; + mon_part array[1]; +}; + +class MonomialOrderings { +public: + static std::string toString(const MonomialOrdering *mo); + static MonomialOrdering* join(const std::vector& M); + static MonomialOrdering* product(const std::vector& M); + static MonomialOrdering* Lex(int nvars); + static MonomialOrdering* Lex2(int nvars); + static MonomialOrdering* Lex4(int nvars); + static MonomialOrdering* GRevLex(int nvars); + static MonomialOrdering* GRevLex2(int nvars); + static MonomialOrdering* GRevLex4(int nvars); + static MonomialOrdering* GRevLex(const std::vector& wts); + static MonomialOrdering* GRevLex2(const std::vector& wts); + static MonomialOrdering* GRevLex4(const std::vector& wts); + static MonomialOrdering* RevLex(int nvars); + static MonomialOrdering* Weights(const std::vector& wts); + static MonomialOrdering* GroupLex(int nvars); + static MonomialOrdering* GroupRevLex(int nvars); + static MonomialOrdering* PositionUp(); + static MonomialOrdering* PositionDown(); + + static MonomialOrdering* GRevLex(const std::vector& wts, int packing); +}; + +// This is currently located in interface/monomial-ordering.cpp +bool monomialOrderingToMatrix( + const struct MonomialOrdering &mo, + std::vector &mat, + bool &base_is_revlex, + int &component_direction, // -1 is Down, +1 is Up, 0 is not present + int &component_is_before_row); // -1 means at the end, 0 means before the + // order, and r means considered before row + // 'r' of the matrix. + +typedef int *monomial; + +typedef int32_t + deg_t; // this is the integer type to use for degrees and weights + +typedef const int *const_monomial; + +struct mo_block +{ + enum MonomialOrdering_type typ; + int nvars; + int nslots; + int first_exp; + int first_slot; + int nweights; + deg_t *weights; +}; + +struct MonomialOrder +{ + int nvars; + int nslots; + int nblocks; + int nblocks_before_component; + int nslots_before_component; + int component_up; /* bool */ + deg_t + *degs; /* 0..nvars: heuristic degree of each variable. degs[nvars] = 1. + Assumption: degs[i] >= 1, for all i, and should be an integer. + Any graded rev lex block assumes graded wrt these degrees. */ + struct mo_block *blocks; /* 0..nblocks-1 with each entry a struct mo_block */ + int *is_laurent; /* 0..nvars-1: 0 or 1: 1 means negative exponents allowed */ + +public: + enum overflow_type { OVER, OVER1, OVER2, OVER4 } * overflow; + std::vector overflow_flags() const; +}; + +MonomialOrder *monomialOrderMake(const MonomialOrdering *mo); +void monomialOrderFree(MonomialOrder *mo); +void monomialOrderEncodeFromActualExponents(const MonomialOrder *mo, + const_exponents a, + monomial b); +void monomialOrderDecodeToActualExponents(const MonomialOrder *mo, + const_monomial a, + exponents_t b); + +std::vector laurentVariables(const MonomialOrder* mo); + +bool monomialOrderingToMatrix( + const struct MonomialOrdering& mo, + std::vector& mat, + bool& base_is_revlex, + int& component_direction, // -1 is Down, +1 is Up, 0 is not present + int& component_is_before_row // -1 means: at the end. 0 means before the + // order. + // and r means considered before row 'r' of the matrix. + ); + +/* +// Local Variables: +// compile-command: "make -C $M2BUILDDIR/Macaulay2/e " +// indent-tabs-mode: nil +// End: +*/ diff --git a/M2/Macaulay2/e/engine-includes.hpp b/M2/Macaulay2/e/engine-includes.hpp index 7d6a0bb7e37..b65293035f9 100644 --- a/M2/Macaulay2/e/engine-includes.hpp +++ b/M2/Macaulay2/e/engine-includes.hpp @@ -9,6 +9,8 @@ #include +// TODO: eventually the engine should not depend on +// types from the interpreter, making this unnecessary. #if !defined(SAFEC_EXPORTS) #include #endif diff --git a/M2/Macaulay2/e/exptable.h b/M2/Macaulay2/e/exptable.h index 803049a5107..d7ee306dfbb 100644 --- a/M2/Macaulay2/e/exptable.h +++ b/M2/Macaulay2/e/exptable.h @@ -3,6 +3,8 @@ /* Implementation of a hashtable [exponent vectors, of fixed length] --> unsigned long int. */ /* The implementation uses table.{h,c}, which was written by David R. Hanson */ +// TODO: replace with std::unordered_map + typedef int * exponent; typedef struct exponent_table exponent_table; diff --git a/M2/Macaulay2/e/f4/moninfo.cpp b/M2/Macaulay2/e/f4/moninfo.cpp index dea242c225b..9e8f2d6efea 100644 --- a/M2/Macaulay2/e/f4/moninfo.cpp +++ b/M2/Macaulay2/e/f4/moninfo.cpp @@ -1,7 +1,7 @@ // Copyright 2005-2021 Michael E. Stillman #include "f4/moninfo.hpp" -#include "monordering.hpp" // for monomialOrderingToMatrix +#include "MonomialOrder.hpp" // for monomialOrderingToMatrix #include "newdelete.hpp" // for freemem, newarray_atomic #include // for fprintf, stderr, stdout diff --git a/M2/Macaulay2/e/freemod.cpp b/M2/Macaulay2/e/freemod.cpp index 6094ee99d1f..89441c3cd43 100644 --- a/M2/Macaulay2/e/freemod.cpp +++ b/M2/Macaulay2/e/freemod.cpp @@ -388,6 +388,7 @@ FreeModule *FreeModule::symm(int n) const return result; } +// TODO: adjust this for arbitrary cone boxes static bool degree_in_box(int len, const_exponents deg, M2_arrayint lo, M2_arrayint hi) { if (lo->len != 0) diff --git a/M2/Macaulay2/e/imonorder.cpp b/M2/Macaulay2/e/imonorder.cpp deleted file mode 100644 index cbc82ec67cd..00000000000 --- a/M2/Macaulay2/e/imonorder.cpp +++ /dev/null @@ -1,656 +0,0 @@ -// Copyright 2009 Michael E. Stillman - -#include "imonorder.hpp" - -#include "engine-includes.hpp" - -#ifdef HAVE_ALLOCA_H -# include -#else -# include -#endif - -#include "ExponentVector.hpp" -#include "overflow.hpp" - -std::vector laurentVariables(const MonomialOrder* mo) -{ - std::vector result; - for (auto i = 0; i < mo->nvars; ++i) - result.push_back(mo->is_laurent[i] == 1); - return result; -} -/* TODO: - -- negative exponent versions need to be included (at least for MO_LEX) - -- non-commutative blocks should be added in -*/ - -static void mo_block_revlex(struct mo_block *b, int nvars) -{ - b->typ = MO_REVLEX; - b->nvars = nvars; - b->nslots = nvars; - b->first_exp = 0; /* will be set later */ - b->first_slot = 0; /* will be set later */ - b->nweights = 0; - b->weights = nullptr; -} - -static void mo_block_grevlex(struct mo_block *b, int nvars) -{ - b->typ = MO_GREVLEX; - b->nvars = nvars; - b->nslots = nvars; - b->first_exp = 0; /* will be set later */ - b->first_slot = 0; /* will be set later */ - b->nweights = 0; - b->weights = nullptr; -} - -static void mo_block_grevlex2(struct mo_block *b, int nvars) -{ - b->typ = MO_GREVLEX2; - b->nvars = nvars; - b->nslots = (nvars + 1) / 2; /* 2 per word */ - b->first_exp = 0; /* will be set later */ - b->first_slot = 0; /* will be set later */ - b->nweights = 0; - b->weights = nullptr; -} - -static void mo_block_grevlex4(struct mo_block *b, int nvars) -{ - b->typ = MO_GREVLEX4; - b->nvars = nvars; - b->nslots = (nvars + 3) / 4; /* 4 per word */ - b->first_exp = 0; /* will be set later */ - b->first_slot = 0; /* will be set later */ - b->nweights = 0; - b->weights = nullptr; -} - -static void mo_block_grevlex_wts(struct mo_block *b, int nvars) -{ - b->typ = MO_GREVLEX_WTS; - b->nvars = nvars; - b->nslots = nvars; - b->first_exp = 0; /* will be set later */ - b->first_slot = 0; /* will be set later */ - b->nweights = nvars; - b->weights = nullptr; /* will be set later */ -} - -static void mo_block_grevlex2_wts(struct mo_block *b, int nvars) -{ - b->typ = MO_GREVLEX2_WTS; - b->nvars = nvars; - b->nslots = (nvars + 1) / 2; /* 2 per word */ - b->first_exp = 0; /* will be set later */ - b->first_slot = 0; /* will be set later */ - b->nweights = nvars; - b->weights = nullptr; /* will be set later */ -} - -static void mo_block_grevlex4_wts(struct mo_block *b, int nvars) -{ - b->typ = MO_GREVLEX4_WTS; - b->nvars = nvars; - b->nslots = (nvars + 3) / 4; /* 4 per word */ - b->first_exp = 0; /* will be set later */ - b->first_slot = 0; /* will be set later */ - b->nweights = nvars; - b->weights = nullptr; /* will be set later */ -} - -static void mo_block_lex(struct mo_block *b, int nvars) -{ - b->typ = MO_LEX; - b->nvars = nvars; - b->nslots = nvars; - b->first_exp = 0; /* will be set later */ - b->first_slot = 0; /* will be set later */ - b->nweights = 0; - b->weights = nullptr; -} - -static void mo_block_lex2(struct mo_block *b, int nvars) -{ - b->typ = MO_LEX2; - b->nvars = nvars; - b->nslots = (nvars + 1) / 2; /* 2 per word */ - b->first_exp = 0; /* will be set later */ - b->first_slot = 0; /* will be set later */ - b->nweights = 0; - b->weights = nullptr; -} - -static void mo_block_lex4(struct mo_block *b, int nvars) -{ - b->typ = MO_LEX4; - b->nvars = nvars; - b->nslots = (nvars + 3) / 4; /* 4 per word */ - b->first_exp = 0; /* will be set later */ - b->first_slot = 0; /* will be set later */ - b->nweights = 0; - b->weights = nullptr; -} - -static void mo_block_group_lex(struct mo_block *b, int nvars) -{ - b->typ = MO_LAURENT; - b->nvars = nvars; - b->nslots = nvars; - b->first_exp = 0; /* will be set later */ - b->first_slot = 0; /* will be set later */ - b->nweights = 0; - b->weights = nullptr; -} - -static void mo_block_group_revlex(struct mo_block *b, int nvars) -{ - b->typ = MO_LAURENT_REVLEX; - b->nvars = nvars; - b->nslots = nvars; - b->first_exp = 0; /* will be set later */ - b->first_slot = 0; /* will be set later */ - b->nweights = 0; - b->weights = nullptr; -} - -static void mo_block_wt_function(struct mo_block *b, int nvars, deg_t *wts) -{ - b->typ = MO_WEIGHTS; - b->nvars = 0; - b->nslots = 1; - b->first_exp = 0; - b->first_slot = 0; /* will be set later */ - b->nweights = nvars; - b->weights = wts; -} - -MonomialOrder *monomialOrderMake(const MonomialOrdering *mo) -{ - MonomialOrder *result; - int i, j, nv, this_block; - deg_t *wts = nullptr; - /* Determine the number of variables, the number of blocks, and the location - of the component */ - int nblocks = 0; - int nvars = 0; - int hascomponent = 0; - for (i = 0; i < mo->len; i++) - { - struct mon_part_rec_ *t = mo->array[i]; - nblocks++; - if (t->type == MO_POSITION_DOWN || t->type == MO_POSITION_UP) - hascomponent++; - else if (t->type == MO_NC_LEX) - { - // Currently, do nothing. - } - if (t->type != MO_WEIGHTS) nvars += t->nvars; - } - nblocks -= hascomponent; - - /* Now create the blocks, and fill them in. Also fill in the deg vector */ - result = getmemstructtype(MonomialOrder *); - result->nvars = nvars; - result->nslots = 0; - result->nblocks = nblocks; - result->blocks = - (struct mo_block *)getmem(nblocks * sizeof(result->blocks[0])); - result->degs = (deg_t *)getmem_atomic(nvars * sizeof(result->degs[0])); - if (hascomponent == 0) result->nblocks_before_component = nblocks; - - this_block = 0; - nvars = 0; - for (i = 0; i < mo->len; i++) - { - struct mon_part_rec_ *t = mo->array[i]; - if (t->type != MO_WEIGHTS) - { - if (t->wts == nullptr) - for (j = 0; j < t->nvars; j++) result->degs[nvars++] = 1; - else - for (j = 0; j < t->nvars; j++) result->degs[nvars++] = t->wts[j]; - } - else - { - wts = (deg_t *)getmem_atomic(t->nvars * sizeof(wts[0])); - for (j = 0; j < t->nvars; j++) wts[j] = t->wts[j]; - } - switch (t->type) - { - case MO_REVLEX: - mo_block_revlex(result->blocks + this_block++, t->nvars); - break; - case MO_GREVLEX: - mo_block_grevlex(result->blocks + this_block++, t->nvars); - break; - case MO_GREVLEX2: - mo_block_grevlex2(result->blocks + this_block++, t->nvars); - break; - case MO_GREVLEX4: - mo_block_grevlex4(result->blocks + this_block++, t->nvars); - break; - case MO_GREVLEX_WTS: - mo_block_grevlex_wts(result->blocks + this_block++, t->nvars); - break; - case MO_GREVLEX2_WTS: - mo_block_grevlex2_wts(result->blocks + this_block++, t->nvars); - break; - case MO_GREVLEX4_WTS: - mo_block_grevlex4_wts(result->blocks + this_block++, t->nvars); - break; - case MO_LEX: - mo_block_lex(result->blocks + this_block++, t->nvars); - break; - case MO_LEX2: - mo_block_lex2(result->blocks + this_block++, t->nvars); - break; - case MO_LEX4: - mo_block_lex4(result->blocks + this_block++, t->nvars); - break; - case MO_WEIGHTS: - // if extra weight values are given (more than "nvars", ignore the - // rest. - mo_block_wt_function( - result->blocks + this_block++, - (t->nvars <= result->nvars ? t->nvars : result->nvars), - wts); - break; - case MO_LAURENT: - mo_block_group_lex(result->blocks + this_block++, t->nvars); - break; - case MO_LAURENT_REVLEX: - mo_block_group_revlex(result->blocks + this_block++, t->nvars); - break; - case MO_NC_LEX: - /* MES */ - break; - case MO_POSITION_UP: - if (--hascomponent == 0) - { - // Set the information about the component - result->component_up = 1; - result->nblocks_before_component = this_block; - } - // mo_block_position_up(result->blocks + this_block); - break; - case MO_POSITION_DOWN: - if (--hascomponent == 0) - { - // Set the information about the component - result->component_up = 0; - result->nblocks_before_component = this_block; - } - // mo_block_position_down(result->blocks + this_block); - break; - } - } - - /* Go back and fill in the 'slots' information */ - /* Now fix the first_exp, first_slot values, and also result->{nslots,nvars}; - */ - nv = 0; - result->nslots = 0; - result->nslots_before_component = 0; - for (i = 0; i < nblocks; i++) - { - enum MonomialOrdering_type typ = result->blocks[i].typ; - - result->blocks[i].first_exp = nv; - result->blocks[i].first_slot = result->nslots; - nv += result->blocks[i].nvars; - result->nslots += result->blocks[i].nslots; - - if (typ == MO_WEIGHTS) - { - result->blocks[i].first_exp = 0; - - /* divide the wt vector by the degree vector */ - for (j = 0; j < result->blocks[i].nvars; j++) - safe::div_by(result->blocks[i].weights[j], result->degs[j]); - ; - } - else if (typ == MO_GREVLEX_WTS || typ == MO_GREVLEX2_WTS || - typ == MO_GREVLEX4_WTS) - { - result->blocks[i].weights = - result->degs + result->blocks[i].first_exp; - } - - if (i == result->nblocks_before_component - 1) - { - result->nslots_before_component = result->nslots; - } - } - - /* Set is_laurent */ - result->is_laurent = (int *)getmem_atomic(result->nvars * sizeof(int)); - for (i = 0; i < result->nvars; i++) result->is_laurent[i] = 0; - - for (i = 0; i < result->nblocks; i++) - if (result->blocks[i].typ == MO_LAURENT || - result->blocks[i].typ == MO_LAURENT_REVLEX) - { - for (j = 0; j < result->blocks[i].nvars; j++) - result->is_laurent[result->blocks[i].first_exp + j] = 1; - } - - return result; -} - -extern void monomialOrderFree(MonomialOrder *mo) {} -static void MO_pack4(int nvars, const int *expon, int *slots) -{ - int32_t i; - if (nvars == 0) return; - while (1) - { - i = safe::fits_7(*expon++) << 24; - if (--nvars == 0) break; - i |= safe::fits_7(*expon++) << 16; - if (--nvars == 0) break; - i |= safe::fits_7(*expon++) << 8; - if (--nvars == 0) break; - i |= safe::fits_7(*expon++); - if (--nvars == 0) break; - *slots++ = i; - } - *slots++ = i; -} - -static void MO_pack2(int nvars, const int *expon, int *slots) -{ - int32_t i; - if (nvars == 0) return; - while (1) - { - i = safe::fits_15(*expon++) << 16; - if (--nvars == 0) break; - i |= safe::fits_15(*expon++); - if (--nvars == 0) break; - *slots++ = i; - } - *slots++ = i; -} - -static void MO_unpack4(int nvars, const int *slots, int *expon) -{ - int32_t i; - if (nvars == 0) return; - while (1) - { - i = *slots++; - *expon++ = (i >> 24); - if (--nvars == 0) break; - *expon++ = (i >> 16) & 0x7f; - if (--nvars == 0) break; - *expon++ = (i >> 8) & 0x7f; - if (--nvars == 0) break; - *expon++ = i & 0x7f; - if (--nvars == 0) break; - } -} - -static void MO_unpack2(int nvars, const int *slots, int *expon) -{ - int32_t i; - if (nvars == 0) return; - while (1) - { - i = *slots++; - *expon++ = i >> 16; - if (--nvars == 0) break; - *expon++ = i & 0x7fff; - if (--nvars == 0) break; - } -} - -void monomialOrderEncodeFromActualExponents(const MonomialOrder *mo, - const_exponents expon, - monomial result_psums) -/* Given 'expon', compute the encoded partial sums value */ -{ - if (mo == nullptr) return; - int *tmpexp = static_cast(alloca((mo->nvars + 1) * sizeof(int))); - int i, j, nvars, s; - int *p1; - deg_t *degs; - struct mo_block *b = mo->blocks; - int nblocks = mo->nblocks; - const int *e = expon; - int *p = result_psums; - for (i = nblocks; i > 0; --i, b++) switch (b->typ) - { - case MO_LEX: - case MO_LAURENT: - nvars = b->nvars; - for (j = 0; j < nvars; j++) *p++ = *e++; - break; - case MO_REVLEX: - case MO_LAURENT_REVLEX: - nvars = b->nvars; - for (j = 0; j < nvars; j++) *p++ = safe::minus(*e++); - break; - case MO_GREVLEX: - nvars = b->nvars; - p += b->nslots; - p1 = p; - *--p1 = *e++; - for (j = 1; j < nvars; j++) - { - --p1; - *p1 = safe::add(*e++, p1[1]); - } - break; - case MO_GREVLEX_WTS: - nvars = b->nvars; - degs = mo->degs + b->first_exp; - p += b->nslots; - p1 = p; - *--p1 = safe::mult(*e++, *degs++); - for (j = 1; j < nvars; j++) - { - --p1; - int tmp = safe::mult(*e++, *degs++); - *p1 = safe::add(tmp, p1[1]); - } - break; - case MO_GREVLEX4: - nvars = b->nvars; - p1 = tmpexp + b->nvars; - *--p1 = *e++; - for (j = 1; j < nvars; j++) - { - --p1; - *p1 = safe::add(*e++, p1[1]); - } - MO_pack4(nvars, p1, p); - p += b->nslots; - break; - case MO_GREVLEX4_WTS: - nvars = b->nvars; - degs = mo->degs + b->first_exp; - p1 = tmpexp + b->nvars; - *--p1 = safe::mult(*e++, *degs++); - for (j = 1; j < nvars; j++) - { - --p1; - int tmp = safe::mult(*e++, *degs++); - *p1 = safe::add(tmp, p1[1]); - } - MO_pack4(nvars, p1, p); - p += b->nslots; - break; - case MO_GREVLEX2: - nvars = b->nvars; - p1 = tmpexp + b->nvars; - *--p1 = *e++; - for (j = 1; j < nvars; j++) - { - --p1; - *p1 = safe::add(*e++, p1[1]); - } - MO_pack2(nvars, p1, p); - p += b->nslots; - break; - case MO_GREVLEX2_WTS: - nvars = b->nvars; - degs = mo->degs + b->first_exp; - p1 = tmpexp + b->nvars; - *--p1 = safe::mult(*e++, *degs++); - for (j = 1; j < nvars; j++) - { - --p1; - int tmp = safe::mult(*e++, *degs++); - *p1 = safe::add(tmp, p1[1]); - } - MO_pack2(nvars, p1, p); - p += b->nslots; - break; - case MO_LEX4: - nvars = b->nvars; - MO_pack4(nvars, e, p); - p += b->nslots; - e += nvars; - break; - case MO_LEX2: - nvars = b->nvars; - MO_pack2(nvars, e, p); - p += b->nslots; - e += nvars; - break; - case MO_WEIGHTS: - if (b->nweights == 0) - { - s = 0; - } - else - { - s = safe::mult(b->weights[0], expon[0]); - for (j = 1; j < b->nweights; j++) - s = safe::add(s, safe::mult(b->weights[j], expon[j])); - } - *p++ = s; - break; - case MO_POSITION_UP: - case MO_POSITION_DOWN: - /* nothing to do here */ - break; - case MO_NC_LEX: - /* nothing to do here */ - break; - } -} - -void monomialOrderDecodeToActualExponents(const MonomialOrder *mo, - const_monomial psums, - exponents_t expon) -{ - if (mo == nullptr) return; - int *tmpexp = static_cast(alloca((mo->nvars + 1) * sizeof(int))); - int i, j, nvars; - deg_t *degs = mo->degs; - deg_t *d; - struct mo_block *b = mo->blocks; - int nblocks = mo->nblocks; - int *e = expon; - const int *p = psums; - for (i = nblocks; i > 0; --i, b++) switch (b->typ) - { - case MO_LEX: - case MO_LAURENT: - nvars = b->nvars; - p = psums + b->first_slot; - e = expon + b->first_exp; - for (j = 0; j < nvars; j++) *e++ = *p++; - break; - case MO_REVLEX: - case MO_LAURENT_REVLEX: - nvars = b->nvars; - p = psums + b->first_slot; - e = expon + b->first_exp; - for (j = 0; j < nvars; j++) *e++ = safe::minus(*p++); - break; - case MO_GREVLEX: - nvars = b->nvars; - p = psums + b->first_slot + nvars - 1; - e = expon + b->first_exp; - *e++ = *p--; - for (j = nvars - 1; j >= 1; --j, --p) *e++ = safe::sub(*p, p[1]); - break; - case MO_GREVLEX_WTS: - nvars = b->nvars; - d = degs + b->first_exp; - p = psums + b->first_slot + nvars - 1; - e = expon + b->first_exp; - *e++ = *p-- / *d++; - for (j = nvars - 1; j >= 1; --j, --p) - *e++ = safe::sub(*p, p[1]) / *d++; - break; - case MO_GREVLEX4: - nvars = b->nvars; - MO_unpack4(nvars, psums + b->first_slot, tmpexp); - p = tmpexp + nvars - 1; - e = expon + b->first_exp; - *e++ = *p--; - for (j = nvars - 1; j >= 1; --j, --p) *e++ = safe::sub(*p, p[1]); - break; - case MO_GREVLEX4_WTS: - nvars = b->nvars; - d = degs + b->first_exp; - MO_unpack4(nvars, psums + b->first_slot, tmpexp); - p = tmpexp + nvars - 1; - e = expon + b->first_exp; - *e++ = *p-- / *d++; - for (j = nvars - 1; j >= 1; --j, --p) - *e++ = safe::sub(*p, p[1]) / *d++; - break; - case MO_GREVLEX2: - nvars = b->nvars; - MO_unpack2(nvars, psums + b->first_slot, tmpexp); - p = tmpexp + nvars - 1; - e = expon + b->first_exp; - *e++ = *p--; - for (j = nvars - 1; j >= 1; --j, --p) *e++ = safe::sub(*p, p[1]); - break; - case MO_GREVLEX2_WTS: - nvars = b->nvars; - d = degs + b->first_exp; - MO_unpack2(nvars, psums + b->first_slot, tmpexp); - p = tmpexp + nvars - 1; - e = expon + b->first_exp; - *e++ = *p-- / *d++; - for (j = nvars - 1; j >= 1; --j, --p) - *e++ = safe::sub(*p, p[1]) / *d++; - break; - case MO_LEX4: - nvars = b->nvars; - e = expon + b->first_exp; - MO_unpack4(nvars, psums + b->first_slot, e); - break; - case MO_LEX2: - nvars = b->nvars; - e = expon + b->first_exp; - MO_unpack2(nvars, psums + b->first_slot, e); - break; - case MO_WEIGHTS: - break; - case MO_POSITION_UP: - case MO_POSITION_DOWN: - /* should not occur, but do nothing in any case */ - break; - case MO_NC_LEX: - /* nothing to do here */ - break; - } -} - -/* -// Local Variables: -// compile-command: "make -C $M2BUILDDIR/Macaulay2/e " -// indent-tabs-mode: nil -// End: -*/ diff --git a/M2/Macaulay2/e/imonorder.hpp b/M2/Macaulay2/e/imonorder.hpp deleted file mode 100644 index 76e7907ef5f..00000000000 --- a/M2/Macaulay2/e/imonorder.hpp +++ /dev/null @@ -1,67 +0,0 @@ -// Copyright 2009 Michael E. Stillman - -#ifndef __imonorder_h_ -#define __imonorder_h_ - -/* This is the internal form of the monomial ordering */ -/* Used in monomial encoding/decoding/comparison */ - -#include - -#include "interface/monomial-ordering.h" - -#include "ExponentVector.hpp" - -typedef int *monomial; - -typedef int32_t - deg_t; // this is the integer type to use for degrees and weights - -typedef const int *const_monomial; - -struct mo_block -{ - enum MonomialOrdering_type typ; - int nvars; - int nslots; - int first_exp; - int first_slot; - int nweights; - deg_t *weights; -}; - -struct MonomialOrder_rec -{ - int nvars; - int nslots; - int nblocks; - int nblocks_before_component; - int nslots_before_component; - int component_up; /* bool */ - deg_t - *degs; /* 0..nvars: heuristic degree of each variable. degs[nvars] = 1. - Assumption: degs[i] >= 1, for all i, and should be an integer. - Any graded rev lex block assumes graded wrt these degrees. */ - struct mo_block *blocks; /* 0..nblocks-1 with each entry a struct mo_block */ - int *is_laurent; /* 0..nvars-1: 0 or 1: 1 means negative exponents allowed */ -}; - -typedef struct MonomialOrder_rec MonomialOrder; -MonomialOrder *monomialOrderMake(const MonomialOrdering *mo); -void monomialOrderFree(MonomialOrder *mo); -void monomialOrderEncodeFromActualExponents(const MonomialOrder *mo, - const_exponents a, - monomial b); -void monomialOrderDecodeToActualExponents(const MonomialOrder *mo, - const_monomial a, - exponents_t b); - -std::vector laurentVariables(const MonomialOrder* mo); -#endif - -/* -// Local Variables: -// compile-command: "make -C $M2BUILDDIR/Macaulay2/e " -// indent-tabs-mode: nil -// End: -*/ diff --git a/M2/Macaulay2/e/int-bag.hpp b/M2/Macaulay2/e/int-bag.hpp index fdd664b7ffa..333c65e5b09 100644 --- a/M2/Macaulay2/e/int-bag.hpp +++ b/M2/Macaulay2/e/int-bag.hpp @@ -4,6 +4,7 @@ #include "newdelete.hpp" +// TODO: replace with a templated tagged varpower struct class int_bag : public our_new_delete { union diff --git a/M2/Macaulay2/e/interface/monomial-ordering.cpp b/M2/Macaulay2/e/interface/monomial-ordering.cpp index 117ed6085cb..7b7a4ce08d3 100644 --- a/M2/Macaulay2/e/interface/monomial-ordering.cpp +++ b/M2/Macaulay2/e/interface/monomial-ordering.cpp @@ -1,990 +1,6 @@ -// TODO: there are many duplicate methods in C and C++ style in this file - #include "interface/monomial-ordering.h" -// TODO: remove this when to_string methods are moved together -#include - -#include -#include -#include - -#include "M2mem.h" // for getmemvectortype, getmematomicvectortype -#include "error.h" -#include "monordering.hpp" // TODO: where can this go? it only defines one class - -static struct mon_part_rec_ *mo_make(enum MonomialOrdering_type type, - int nvars, - const int *wts) -{ - mon_part result; - result = getmemstructtype(mon_part); - result->type = type; - result->nvars = nvars; - if (wts != nullptr) - { - int i; - result->wts = getmematomicvectortype(int, nvars); - for (i = 0; i < nvars; i++) result->wts[i] = wts[i]; - } - else - result->wts = nullptr; - return result; -} - -static MonomialOrdering *make_mon_order(int n) -{ - static unsigned int next_hash = 23023421; - MonomialOrdering *z = getmemarraytype(MonomialOrdering *, n); - z->len = n; - z->_hash = next_hash++; - int i; - for (i = 0; i < n; i++) z->array[i] = nullptr; - return z; -} - -static MonomialOrdering *M2_mo_offset(const MonomialOrdering *mo, int offset) -{ - int i, j; - MonomialOrdering *result = make_mon_order(mo->len); - for (i = 0; i < mo->len; i++) - { - mon_part p = mo->array[i]; - if (p->type != MO_WEIGHTS) - result->array[i] = mo_make(p->type, p->nvars, p->wts); - else - { - mon_part q = mo_make(MO_WEIGHTS, offset + p->nvars, nullptr); - q->wts = getmemvectortype(int, q->nvars); - for (j = 0; j < offset; j++) q->wts[j] = 0; - for (; j < q->nvars; j++) q->wts[j] = p->wts[j - offset]; - } - } - return result; -} - -static bool is_good(mon_part p) -{ - switch (p->type) - { - case MO_LEX: - case MO_LEX2: - case MO_LEX4: - case MO_GREVLEX: - case MO_GREVLEX2: - case MO_GREVLEX4: - case MO_GREVLEX_WTS: - case MO_GREVLEX2_WTS: - case MO_GREVLEX4_WTS: - case MO_LAURENT: - case MO_NC_LEX: - case MO_LAURENT_REVLEX: - case MO_REVLEX: - case MO_WEIGHTS: - return (p->nvars > 0); - case MO_POSITION_UP: - case MO_POSITION_DOWN: - return true; - } - return false; -} - -MonomialOrdering *MonomialOrderings::Lex(int nvars) -{ - return rawLexMonomialOrdering(nvars, 1); -} -MonomialOrdering *MonomialOrderings::Lex2(int nvars) -{ - return rawLexMonomialOrdering(nvars, 2); -} -MonomialOrdering *MonomialOrderings::Lex4(int nvars) -{ - return rawLexMonomialOrdering(nvars, 4); -} -MonomialOrdering *MonomialOrderings::GRevLex(int nvars) -{ - std::vector w; - for (int i = 0; i < nvars; i++) w.push_back(1); - return GRevLex(w); -} -MonomialOrdering *MonomialOrderings::GRevLex2(int nvars) -{ - std::vector w; - for (int i = 0; i < nvars; i++) w.push_back(1); - return GRevLex2(w); -} -MonomialOrdering *MonomialOrderings::GRevLex4(int nvars) -{ - std::vector w; - for (int i = 0; i < nvars; i++) w.push_back(1); - return GRevLex4(w); -} -MonomialOrdering *MonomialOrderings::GRevLex(const std::vector &wts) -{ - return GRevLex(wts, 1); -} -MonomialOrdering *MonomialOrderings::GRevLex2(const std::vector &wts) -{ - return GRevLex(wts, 2); -} -MonomialOrdering *MonomialOrderings::GRevLex4(const std::vector &wts) -{ - return GRevLex(wts, 4); -} -MonomialOrdering *MonomialOrderings::RevLex(int nvars) -{ - return rawRevLexMonomialOrdering(nvars); -} -MonomialOrdering *MonomialOrderings::Weights(const std::vector &wts) -{ - mon_part p = mo_make(MO_WEIGHTS, wts.size(), wts.data()); - MonomialOrdering *result = make_mon_order(1); - result->array[0] = p; - return result; -} -MonomialOrdering *MonomialOrderings::GroupLex(int nvars) -{ - return rawGroupLexMonomialOrdering(nvars); -} -MonomialOrdering *MonomialOrderings::GroupRevLex(int nvars) -{ - return rawGroupRevLexMonomialOrdering(nvars); -} -MonomialOrdering *MonomialOrderings::PositionUp() -{ - return rawPositionMonomialOrdering(true); -} -MonomialOrdering *MonomialOrderings::PositionDown() -{ - return rawPositionMonomialOrdering(false); -} - -MonomialOrdering *MonomialOrderings::GRevLex(const std::vector °s, - int packing) -{ - MonomialOrdering *result; - mon_part p; - enum MonomialOrdering_type typ; - const int *wts; - bool all_one = true; - for (int i = 0; i < degs.size(); i++) - if (degs[i] <= 0) - { - ERROR("grevlex: expected all degrees to be positive"); - return nullptr; - } - else if (degs[i] > 1) - all_one = false; - - if (all_one) - { - if (packing == 2) - typ = MO_GREVLEX2; - else if (packing == 4) - typ = MO_GREVLEX4; - else - typ = MO_GREVLEX; - wts = nullptr; - } - else - { - if (packing == 2) - typ = MO_GREVLEX2_WTS; - else if (packing == 4) - typ = MO_GREVLEX4_WTS; - else - typ = MO_GREVLEX_WTS; - wts = degs.data(); - } - - p = mo_make(typ, degs.size(), wts); - result = make_mon_order(1); - result->array[0] = p; - return result; -} - -MonomialOrdering *MonomialOrderings::join( - const std::vector &M) -{ - MonomialOrdering *result; - const MonomialOrdering *mo; - int i, j, sum, next; - int nvars_so_far = 0; - - sum = 0; - for (i = 0; i < M.size(); i++) - { - mo = M[i]; - for (j = 0; j < mo->len; j++) - if (is_good(mo->array[j])) sum++; - } - - result = make_mon_order(sum); - next = 0; - for (i = 0; i < M.size(); i++) - { - mo = M[i]; - for (j = 0; j < mo->len; j++) - { - mon_part p = mo->array[j]; - if (!is_good(p)) continue; - if (p->type != MO_WEIGHTS) - nvars_so_far += p->nvars; - else - { - /* Shift the weights over by nvars_so_far */ - mon_part q = mo_make(MO_WEIGHTS, nvars_so_far + p->nvars, nullptr); - q->wts = getmemvectortype(int, q->nvars); - for (j = 0; j < nvars_so_far; j++) q->wts[j] = 0; - for (; j < q->nvars; j++) q->wts[j] = p->wts[j - nvars_so_far]; - p = q; - } - result->array[next++] = p; - } - } - return result; -} - -MonomialOrdering *MonomialOrderings::product( - const std::vector &M) -{ - MonomialOrdering *result; - int i, j, sum, next, offset; - sum = 0; - for (i = 0; i < M.size(); i++) sum += M[i]->len; - - result = make_mon_order(sum); - next = 0; - offset = 0; - for (i = 0; i < M.size(); i++) - { - int nvars = rawNumberOfVariables(M[i]); - MonomialOrdering *mo = M2_mo_offset(M[i], offset); - for (j = 0; j < mo->len; j++) result->array[next++] = mo->array[j]; - offset += nvars; - } - return result; -} - -std::ostringstream &toString(std::ostringstream &o, int len, int *p) -{ - o << "{"; - for (int i = 0; i < len; i++) - { - if (i > 0) o << ","; - o << p[i]; - } - o << "}"; - return o; -} - -std::ostringstream &ones(std::ostringstream &o, int len) -{ - o << "{"; - for (int i = 0; i < len; i++) - { - if (i > 0) o << ","; - o << 1; - } - o << "}"; - return o; -} - -std::string MonomialOrderings::toString(const MonomialOrdering *mo) -{ - std::ostringstream o; - o << "MonomialOrder => {"; - for (int i = 0; i < mo->len; i++) - { - mon_part p = mo->array[i]; - bool p_ones = false; - if (i == 0) - o << "\n "; - else - o << ",\n "; - switch (p->type) - { - case MO_LEX: - o << "Lex => " << p->nvars; - break; - case MO_LEX2: - o << "LexSmall => " << p->nvars; - break; - case MO_LEX4: - o << "LexTiny => " << p->nvars; - break; - case MO_GREVLEX: - o << "GRevLex => "; - p_ones = true; - break; - case MO_GREVLEX2: - o << "GRevLexSmall => "; - p_ones = true; - break; - case MO_GREVLEX4: - o << "GRevLexTiny => "; - p_ones = true; - break; - case MO_GREVLEX_WTS: - o << "GRevLex => "; - break; - case MO_GREVLEX2_WTS: - o << "GRevLexSmall => "; - break; - case MO_GREVLEX4_WTS: - o << "GRevLexTiny => "; - break; - case MO_REVLEX: - o << "RevLex => " << p->nvars; - break; - case MO_WEIGHTS: - o << "Weights => "; - break; - case MO_LAURENT: - o << "GroupLex => " << p->nvars; - break; - case MO_LAURENT_REVLEX: - o << "GroupRevLex => " << p->nvars; - break; - case MO_NC_LEX: - o << "NCLex => " << p->nvars; - break; - case MO_POSITION_UP: - o << "Position => Up"; - break; - case MO_POSITION_DOWN: - o << "Position => Down"; - break; - default: - o << "UNKNOWN"; - break; - } - if (p->wts != nullptr) { ::toString(o, p->nvars, p->wts); } - else if (p_ones) - { - ::ones(o, p->nvars); - } - } - o << "\n }"; - return o.str(); -} - -///////// Below this is from monordering.h //////////////// - -int moIsLex(const MonomialOrdering *mo) -{ - // The monomial order is lex if what? - // one lex block, no grevlex blocks, no weightvector blocks. - // only: lex block and position blocks are allowed. - int nlex = 0; - int i; - for (i = 0; i < mo->len; i++) - { - enum MonomialOrdering_type typ = mo->array[i]->type; - switch (typ) - { - case MO_LEX: - case MO_LEX2: - case MO_LEX4: - nlex++; - break; - case MO_POSITION_UP: - case MO_POSITION_DOWN: - break; - default: - return 0; - } - } - return (nlex == 1); -} - -int moIsGRevLex(const MonomialOrdering *mo) -{ - // The monomial order is lex if what? - // one lex block, no grevlex blocks, no weightvector blocks. - // only: lex block and position blocks are allowed. - int ngrevlex = 0; - int i; - for (i = 0; i < mo->len; i++) - { - enum MonomialOrdering_type typ = mo->array[i]->type; - switch (typ) - { - case MO_GREVLEX: - case MO_GREVLEX2: - case MO_GREVLEX4: - case MO_GREVLEX_WTS: - case MO_GREVLEX2_WTS: - case MO_GREVLEX4_WTS: - ngrevlex++; - break; - case MO_POSITION_UP: - case MO_POSITION_DOWN: - break; - default: - return 0; - } - } - return (ngrevlex == 1); -} - -int rawNumberOfVariables(const MonomialOrdering *mo) -{ - int i, sum = 0; - for (i = 0; i < mo->len; i++) - if (mo->array[i]->type != MO_WEIGHTS) sum += mo->array[i]->nvars; - return sum; -} - -M2_arrayint moGetWeightValues(const MonomialOrdering *mo) -{ - int nvars = rawNumberOfVariables(mo); - // grab the first weight vector - if (mo->len == 0) return nullptr; - if (mo->array[0]->type == MO_WEIGHTS) - { - int i; - M2_arrayint result = M2_makearrayint(nvars); - int *wts = mo->array[0]->wts; - for (i = 0; i < mo->array[0]->nvars; i++) result->array[i] = wts[i]; - for (; i < nvars; i++) result->array[i] = 0; - return result; - } - return nullptr; -} - -int rawNumberOfInvertibleVariables(const MonomialOrdering *mo) -{ - int i, sum = 0; - for (i = 0; i < mo->len; i++) - if (mo->array[i]->type == MO_LAURENT || - mo->array[i]->type == MO_LAURENT_REVLEX) - sum += mo->array[i]->nvars; - return sum; -} - -M2_arrayint rawNonTermOrderVariables(const MonomialOrdering *mo) -// returns a list of the indices of those variables which are less than 1 in -// the given monomial order. -{ - int i, j, sum, nextvar; - int nvars = rawNumberOfVariables(mo); - int *gt = getmematomicvectortype(int, nvars); - for (i = 0; i < nvars; i++) - gt[i] = - 0; // 0 means undecided, -1 means non term order, 1 means term order - // Now we loop through the parts of the monomial order - nextvar = 0; - for (i = 0; i < mo->len; i++) - { - mon_part p = mo->array[i]; - switch (p->type) - { - case MO_LEX: - case MO_LEX2: - case MO_LEX4: - case MO_GREVLEX: - case MO_GREVLEX2: - case MO_GREVLEX4: - case MO_GREVLEX_WTS: - case MO_GREVLEX2_WTS: - case MO_GREVLEX4_WTS: - case MO_LAURENT: - case MO_NC_LEX: - for (j = 0; j < p->nvars; j++, nextvar++) - if (gt[nextvar] == 0) gt[nextvar] = 1; - break; - case MO_LAURENT_REVLEX: - case MO_REVLEX: - for (j = 0; j < p->nvars; j++, nextvar++) - if (gt[nextvar] == 0) gt[nextvar] = -1; - break; - case MO_WEIGHTS: - for (j = nextvar; j < p->nvars; j++) - if (gt[j] == 0) - { - if (p->wts[j] > 0) - gt[j] = 1; - else if (p->wts[j] < 0) - gt[j] = -1; - } - break; - case MO_POSITION_UP: - case MO_POSITION_DOWN: - break; - } - } - // At this point every variables' gt should be 1 or -1. - sum = 0; - for (i = 0; i < nvars; i++) - { - if (gt[i] == 0) INTERNAL_ERROR("gt[i] should not be 0"); - if (gt[i] < 0) sum++; - } - - // Make an array of this length. - M2_arrayint result = M2_makearrayint(sum); - nextvar = 0; - for (i = 0; i < nvars; i++) - if (gt[i] < 0) result->array[nextvar++] = i; - return result; -} - -MonomialOrdering *rawLexMonomialOrdering(int nvars, int packing) -{ - MonomialOrdering *result; - mon_part p; - enum MonomialOrdering_type typ; - - if (packing == 2) - typ = MO_LEX2; - else if (packing == 4) - typ = MO_LEX4; - else - typ = MO_LEX; - - p = mo_make(typ, nvars, nullptr); - result = make_mon_order(1); - result->array[0] = p; - return result; -} - -MonomialOrdering /* or null */ *rawGRevLexMonomialOrdering(M2_arrayint degs, - int packing) -{ - MonomialOrdering *result; - mon_part p; - enum MonomialOrdering_type typ; - int *wts; - int all_one = 1; - unsigned int i; - for (i = 0; i < degs->len; i++) - if (degs->array[i] <= 0) - { - ERROR("grevlex: expected all degrees to be positive"); - return nullptr; - } - else if (degs->array[i] > 1) - all_one = 0; - - if (all_one) - { - if (packing == 2) - typ = MO_GREVLEX2; - else if (packing == 4) - typ = MO_GREVLEX4; - else - typ = MO_GREVLEX; - wts = nullptr; - } - else - { - if (packing == 2) - typ = MO_GREVLEX2_WTS; - else if (packing == 4) - typ = MO_GREVLEX4_WTS; - else - typ = MO_GREVLEX_WTS; - wts = degs->array; - } - - p = mo_make(typ, degs->len, wts); - result = make_mon_order(1); - result->array[0] = p; - return result; -} - -MonomialOrdering *rawRevLexMonomialOrdering(int nvars) -{ - mon_part p = mo_make(MO_REVLEX, nvars, nullptr); - MonomialOrdering *result = make_mon_order(1); - result->array[0] = p; - return result; -} - -MonomialOrdering *rawWeightsMonomialOrdering(M2_arrayint wts) -{ - mon_part p = mo_make(MO_WEIGHTS, wts->len, wts->array); - MonomialOrdering *result = make_mon_order(1); - result->array[0] = p; - return result; -} -MonomialOrdering *rawGroupLexMonomialOrdering(int nvars) -{ - mon_part p = mo_make(MO_LAURENT, nvars, nullptr); - MonomialOrdering *result = make_mon_order(1); - result->array[0] = p; - return result; -} -MonomialOrdering *rawGroupRevLexMonomialOrdering(int nvars) -{ - mon_part p = mo_make(MO_LAURENT_REVLEX, nvars, nullptr); - MonomialOrdering *result = make_mon_order(1); - result->array[0] = p; - return result; -} -MonomialOrdering *rawNClexMonomialOrdering(int nvars) -{ - mon_part p = mo_make(MO_NC_LEX, nvars, nullptr); - MonomialOrdering *result = make_mon_order(1); - result->array[0] = p; - return result; -} -MonomialOrdering *rawPositionMonomialOrdering(M2_bool up_or_down) -{ - mon_part p = - mo_make((up_or_down ? MO_POSITION_UP : MO_POSITION_DOWN), 0, nullptr); - MonomialOrdering *result = make_mon_order(1); - result->array[0] = p; - return result; -} - -MonomialOrdering *rawJoinMonomialOrdering(engine_RawMonomialOrderingArray M) -{ - MonomialOrdering *result; - const MonomialOrdering *mo; - int i, j, sum, next; - int nvars_so_far = 0; - - // sum = 0; - // for (i=0; ilen; i++) - // sum += M->array[i]->len; - - sum = 0; - for (i = 0; i < M->len; i++) - { - mo = M->array[i]; - for (j = 0; j < mo->len; j++) - if (is_good(mo->array[j])) sum++; - } - - result = make_mon_order(sum); - next = 0; - for (i = 0; i < M->len; i++) - { - mo = M->array[i]; - for (j = 0; j < mo->len; j++) - { - mon_part p = mo->array[j]; - if (!is_good(p)) continue; - if (p->type != MO_WEIGHTS) - nvars_so_far += p->nvars; - else - { - /* Shift the weights over by nvars_so_far */ - mon_part q = mo_make(MO_WEIGHTS, nvars_so_far + p->nvars, nullptr); - q->wts = getmemvectortype(int, q->nvars); - for (j = 0; j < nvars_so_far; j++) q->wts[j] = 0; - for (; j < q->nvars; j++) q->wts[j] = p->wts[j - nvars_so_far]; - p = q; - } - result->array[next++] = p; - } - } - return result; -} - -MonomialOrdering *rawProductMonomialOrdering(engine_RawMonomialOrderingArray M) -{ - MonomialOrdering *result; - int i, j, sum, next, offset; - sum = 0; - for (i = 0; i < M->len; i++) sum += M->array[i]->len; - - result = make_mon_order(sum); - next = 0; - offset = 0; - for (i = 0; i < M->len; i++) - { - int nvars = rawNumberOfVariables(M->array[i]); - MonomialOrdering *mo = M2_mo_offset(M->array[i], offset); - for (j = 0; j < mo->len; j++) result->array[next++] = mo->array[j]; - offset += nvars; - } - return result; -} - -M2_string intarray_to_string(int len, int *p) -{ - int i; - const int N = 200; - char s[N]; - M2_string result = M2_tostring("{"); - for (i = 0; i < len; i++) - { - if (i > 0) result = M2_join(result, M2_tostring(",")); - snprintf(s, N, "%d", p[i]); - result = M2_join(result, M2_tostring(s)); - } - result = M2_join(result, M2_tostring("}")); - return result; -} - -M2_string ones_to_string(int len) -{ - int i; - const int N = 200; - char s[N]; - M2_string one; - M2_string result = M2_tostring("{"); - snprintf(s, N, "1"); - one = M2_tostring(s); - for (i = 0; i < len; i++) - { - if (i > 0) result = M2_join(result, M2_tostring(",")); - result = M2_join(result, one); - } - result = M2_join(result, M2_tostring("}")); - return result; -} - -unsigned int rawMonomialOrderingHash(const MonomialOrdering *mo) -{ - return mo->_hash; -} - -M2_string IM2_MonomialOrdering_to_string(const MonomialOrdering *mo) -{ - int i; - const int N = 200; - char s[N]; - int p_ones = 0; - M2_string result = M2_tostring("MonomialOrder => {"); - for (i = 0; i < mo->len; i++) - { - mon_part p = mo->array[i]; - p_ones = 0; - if (i == 0) - result = M2_join(result, M2_tostring("\n ")); - else - result = M2_join(result, M2_tostring(",\n ")); - switch (p->type) - { - case MO_LEX: - snprintf(s, N, "Lex => %d", p->nvars); - break; - case MO_LEX2: - snprintf(s, N, "LexSmall => %d", p->nvars); - break; - case MO_LEX4: - snprintf(s, N, "LexTiny => %d", p->nvars); - break; - case MO_GREVLEX: - snprintf(s, N, "GRevLex => "); - p_ones = 1; - break; - case MO_GREVLEX2: - snprintf(s, N, "GRevLexSmall => "); - p_ones = 1; - break; - case MO_GREVLEX4: - snprintf(s, N, "GRevLexTiny => "); - p_ones = 1; - break; - case MO_GREVLEX_WTS: - snprintf(s, N, "GRevLex => "); - break; - case MO_GREVLEX2_WTS: - snprintf(s, N, "GRevLexSmall => "); - break; - case MO_GREVLEX4_WTS: - snprintf(s, N, "GRevLexTiny => "); - break; - case MO_REVLEX: - snprintf(s, N, "RevLex => %d", p->nvars); - break; - case MO_WEIGHTS: - snprintf(s, N, "Weights => "); - break; - case MO_LAURENT: - snprintf(s, N, "GroupLex => %d", p->nvars); - break; - case MO_LAURENT_REVLEX: - snprintf(s, N, "GroupRevLex => %d", p->nvars); - break; - case MO_NC_LEX: - snprintf(s, N, "NCLex => %d", p->nvars); - break; - case MO_POSITION_UP: - snprintf(s, N, "Position => Up"); - break; - case MO_POSITION_DOWN: - snprintf(s, N, "Position => Down"); - break; - default: - snprintf(s, N, "UNKNOWN"); - break; - } - result = M2_join(result, M2_tostring(s)); - if (p->wts != nullptr) - result = M2_join(result, intarray_to_string(p->nvars, p->wts)); - else if (p_ones) - result = M2_join(result, ones_to_string(p->nvars)); - } - result = M2_join(result, M2_tostring("\n }")); - return result; -} - -// Many monomial ordering routines are in monordering.c -// Here are some that use c++ features, so cannot be there. -// @todo Make monordering.{h,c} into a c++ class. - -static void write_row(std::vector &grading, - int nvars, - int which, - int value) -{ - for (int i = 0; i < nvars; i++) - if (i == which) - grading.push_back(value); - else - grading.push_back(0); -} -static void write_weights(std::vector &grading, - int nvars, - int firstvar, - int *wts, - int nwts) -// place nvars ints into grading: 0 ... 0 wts[0] wts[1] ... wts[nwts-1] 0 .... -// 0 -// where wts[0] is in the 'firstvar' location. If wts is NULL, treat it as the -// vector with nwts '1's. -{ - for (int i = 0; i < firstvar; i++) grading.push_back(0); - if (wts == nullptr) - for (int i = 0; i < nwts; i++) grading.push_back(1); - else - for (int i = 0; i < nwts; i++) grading.push_back(wts[i]); - for (int i = firstvar + nwts; i < nvars; i++) grading.push_back(0); -} - -bool monomialOrderingToMatrix( - const struct MonomialOrdering &mo, - std::vector &mat, - bool &base_is_revlex, - int &component_direction, // -1 is Down, +1 is Up, 0 is not present - int &component_is_before_row) // -1 means at the end, 0 means before the - // order, and r means considered before row - // 'r' of the matrix. -{ - // a false return value means an error has occurred. - int nvars = rawNumberOfVariables(&mo); - base_is_revlex = true; - enum LastBlock { LEX, REVLEX, WEIGHTS, NONE }; - LastBlock last = NONE; - int nwts = 0; // local var used in MO_WEIGHTS section - int nrows = 0; - int firstvar = 0; - component_direction = 0; - component_is_before_row = - -2; // what should the default value be? Probably: -1. - size_t last_element = 0; // The vector 'mat' will be resized back to this - // value if the last part of the order is lex or - // revlex. - for (int i = 0; i < mo.len; i++) - { - mon_part p = mo.array[i]; - switch (p->type) - { - case MO_LEX: - case MO_LEX2: - case MO_LEX4: - // printf("lex %d\n", p->nvars); - last_element = mat.size(); - for (int j = 0; j < p->nvars; j++) - { - write_row(mat, nvars, firstvar + j, 1); - } - last = LEX; - firstvar += p->nvars; - nrows += p->nvars; - break; - case MO_GREVLEX: - case MO_GREVLEX2: - case MO_GREVLEX4: - // printf("grevlex %d %ld\n", p->nvars, p->wts); - write_weights(mat, nvars, firstvar, p->wts, p->nvars); - last_element = mat.size(); - for (int j = p->nvars - 1; j >= 1; --j) - { - write_row(mat, nvars, firstvar + j, -1); - } - last = REVLEX; - firstvar += p->nvars; - nrows += p->nvars; - break; - case MO_GREVLEX_WTS: - case MO_GREVLEX2_WTS: - case MO_GREVLEX4_WTS: - // printf("grevlex_wts %d %ld\n", p->nvars, p->wts); - write_weights(mat, nvars, firstvar, p->wts, p->nvars); - last_element = mat.size(); - for (int j = p->nvars - 1; j >= 1; --j) - { - write_row(mat, nvars, firstvar + j, -1); - } - last = REVLEX; - firstvar += p->nvars; - nrows += p->nvars; - break; - case MO_REVLEX: - // printf("revlex %d\n", p->nvars); - last_element = mat.size(); - for (int j = p->nvars - 1; j >= 0; --j) - { - write_row(mat, nvars, firstvar + j, -1); - } - last = REVLEX; - firstvar += p->nvars; - nrows += p->nvars; - break; - case MO_WEIGHTS: - // printf("matsize= %d weights %d p->wts=%lu\n", mat.size(), - // p->nvars, p->wts); - nwts = (p->nvars > nvars ? nvars : p->nvars); - write_weights(mat, nvars, 0, p->wts, nwts); - nrows++; - last_element = mat.size(); - last = WEIGHTS; - break; - case MO_LAURENT: - case MO_LAURENT_REVLEX: - case MO_NC_LEX: - return false; - break; - case MO_POSITION_UP: - component_direction = 1; - component_is_before_row = nrows; - break; - case MO_POSITION_DOWN: - component_direction = -1; - component_is_before_row = nrows; - break; - default: - // DO nothing - break; - } - } - if (last == LEX) - { - // last block was lex, so use lex tie-breaker - mat.resize(last_element); - if (nrows == component_is_before_row) component_is_before_row = -1; - base_is_revlex = false; - } - else if (last == REVLEX) - { - // last block was revlex, so use revlex tie-breaker - if (nrows == component_is_before_row) component_is_before_row = -1; - mat.resize(last_element); - } - else - { - // last block is a weight vector, so use revlex as the tie-breaker. - // nothing to change here. - } - return true; -} +#include "MonomialOrderOld.hpp" M2_arrayint rawMonomialOrderingToMatrix(const struct MonomialOrdering *mo) { @@ -1005,3 +21,7 @@ M2_arrayint rawMonomialOrderingToMatrix(const struct MonomialOrdering *mo) } return result; } + +// Local Variables: +// indent-tabs-mode: nil +// End: diff --git a/M2/Macaulay2/e/interface/monomial-ordering.h b/M2/Macaulay2/e/interface/monomial-ordering.h index 647c3ee0a15..8977b7adebe 100644 --- a/M2/Macaulay2/e/interface/monomial-ordering.h +++ b/M2/Macaulay2/e/interface/monomial-ordering.h @@ -3,12 +3,19 @@ # include "engine-includes.hpp" +# if defined(__cplusplus) +class MonomialOrdering; +# else typedef struct MonomialOrdering MonomialOrdering; +# endif + +# if defined(__cplusplus) +extern "C" { +# endif /** MonomialOrdering interface routines */ - enum MonomialOrdering_type { MO_LEX = 1, MO_LEX2 = 2, @@ -28,24 +35,6 @@ enum MonomialOrdering_type { MO_POSITION_DOWN = 16 }; -typedef struct mon_part_rec_ -{ - enum MonomialOrdering_type type; - int nvars; - int *wts; -} * mon_part; - -struct MonomialOrdering -{ - unsigned int _hash; - unsigned int len; - mon_part array[1]; -}; - -# if defined(__cplusplus) -extern "C" { -# endif - MonomialOrdering *rawLexMonomialOrdering(int nvars, int packing); /* drg: connected rawMonomialOrdering*/ /* Lex, LexSmall, LexTiny */ @@ -89,6 +78,8 @@ MonomialOrdering *rawJoinMonomialOrdering(engine_RawMonomialOrderingArray mo); /* drg: connected rawMonomialOrdering*/ /* default, when making monoids and polynomial rings */ +//////// Informational //////////////////////////////// + int rawNumberOfVariables(const MonomialOrdering *mo); /* drg: connected rawNumberOfVariables*/ @@ -134,7 +125,7 @@ M2_arrayintOrNull rawMonomialOrderingToMatrix( } # endif -#endif /* _monomial-ordering_h_ */ +#endif /* _monomial_ordering_h_ */ // Local Variables: // indent-tabs-mode: nil diff --git a/M2/Macaulay2/e/interface/ringelement.cpp b/M2/Macaulay2/e/interface/ringelement.cpp index 6aed3d0b018..46ce4cdb8f1 100644 --- a/M2/Macaulay2/e/interface/ringelement.cpp +++ b/M2/Macaulay2/e/interface/ringelement.cpp @@ -350,6 +350,7 @@ M2_bool IM2_RingElement_is_graded(const RingElement *a) return a->is_homogeneous(); } +// TODO: should this be M2_arrayint or M2_arrayintOrNull? M2_arrayint IM2_RingElement_multidegree(const RingElement *a) { try diff --git a/M2/Macaulay2/e/interface/ringelement.h b/M2/Macaulay2/e/interface/ringelement.h index 593cb9ea631..3600fc92f42 100644 --- a/M2/Macaulay2/e/interface/ringelement.h +++ b/M2/Macaulay2/e/interface/ringelement.h @@ -91,6 +91,7 @@ const RingElement /* or null */ *IM2_RingElement_lift(int *success_return, M2_bool IM2_RingElement_is_graded(const RingElement *a); +// FIXME: why are there duplicates? which return type is correct? M2_arrayint IM2_RingElement_multidegree(const RingElement *a); const RingElement * /* or null */ rawRingElementAntipode(const RingElement *f); diff --git a/M2/Macaulay2/e/monideal.cpp b/M2/Macaulay2/e/monideal.cpp index 59c6e6edfbe..b4242d1ff69 100644 --- a/M2/Macaulay2/e/monideal.cpp +++ b/M2/Macaulay2/e/monideal.cpp @@ -277,7 +277,8 @@ void MonomialIdeal::find_all_divisors(const_exponents exp, VECTOR(Bag *)& b) con int MonomialIdeal::search(const_varpower m, Bag *&b) const { - exponents_t exp = ARRAY_ON_STACK(int, get_ring()->n_vars()); + // TODO: unify best practices; why alloca instead of getmem_atomic? + exponents_t exp = ALLOCATE_EXPONENTS(EXPONENT_BYTE_SIZE(get_ring()->n_vars())); varpower::to_expvector(get_ring()->n_vars(), m, exp); return search_expvector(exp, b); } diff --git a/M2/Macaulay2/e/monoid.hpp b/M2/Macaulay2/e/monoid.hpp index 418adcf136d..4060e99ebda 100644 --- a/M2/Macaulay2/e/monoid.hpp +++ b/M2/Macaulay2/e/monoid.hpp @@ -11,7 +11,7 @@ #include "ExponentList.hpp" #include "ExponentVector.hpp" #include "hash.hpp" -#include "imonorder.hpp" +#include "MonomialOrderOld.hpp" #include "newdelete.hpp" #include "style.hpp" @@ -36,6 +36,7 @@ typedef const int *const_monomial; // TODO: make sure monoid is deconstructed and not garbage collected class Monoid : public MutableEngineObject { + // TODO: perhaps mDegreeMonoid should just be a ZZ-module (potentially with torsion?) const Monoid *mDegreeMonoid; const PolynomialRing *mDegreeRing; diff --git a/M2/Macaulay2/e/monomial.hpp b/M2/Macaulay2/e/monomial.hpp index cd4f2553596..f995136106b 100644 --- a/M2/Macaulay2/e/monomial.hpp +++ b/M2/Macaulay2/e/monomial.hpp @@ -35,6 +35,7 @@ class EngineMonomial : public EngineObject // [2n+1, v1, e1, v2, e2, ..., vn, en] // with each ei != 0. + // TODO: change to const_varpower? int * ints() { return val.data(); } const int * ints() const { return val.data(); } diff --git a/M2/Macaulay2/e/monordering.hpp b/M2/Macaulay2/e/monordering.hpp deleted file mode 100644 index f533aec51ab..00000000000 --- a/M2/Macaulay2/e/monordering.hpp +++ /dev/null @@ -1,50 +0,0 @@ -// TODO: move this file: these are interface routines - -#ifndef __monordering_hpp_ -#define __monordering_hpp_ - -#include -#include - -class MonomialOrderings { -public: - static std::string toString(const MonomialOrdering *mo); - static MonomialOrdering* join(const std::vector& M); - static MonomialOrdering* product(const std::vector& M); - static MonomialOrdering* Lex(int nvars); - static MonomialOrdering* Lex2(int nvars); - static MonomialOrdering* Lex4(int nvars); - static MonomialOrdering* GRevLex(int nvars); - static MonomialOrdering* GRevLex2(int nvars); - static MonomialOrdering* GRevLex4(int nvars); - static MonomialOrdering* GRevLex(const std::vector& wts); - static MonomialOrdering* GRevLex2(const std::vector& wts); - static MonomialOrdering* GRevLex4(const std::vector& wts); - static MonomialOrdering* RevLex(int nvars); - static MonomialOrdering* Weights(const std::vector& wts); - static MonomialOrdering* GroupLex(int nvars); - static MonomialOrdering* GroupRevLex(int nvars); - static MonomialOrdering* PositionUp(); - static MonomialOrdering* PositionDown(); - - static MonomialOrdering* GRevLex(const std::vector& wts, int packing); -}; - -// This is currently located in interface/monomial-ordering.cpp -bool monomialOrderingToMatrix( - const struct MonomialOrdering &mo, - std::vector &mat, - bool &base_is_revlex, - int &component_direction, // -1 is Down, +1 is Up, 0 is not present - int &component_is_before_row); // -1 means at the end, 0 means before the - // order, and r means considered before row - // 'r' of the matrix. - -#endif - -/* -// Local Variables: -// compile-command: "make -C $M2BUILDDIR/Macaulay2/e " -// indent-tabs-mode: nil -// End: -*/ diff --git a/M2/Macaulay2/e/ring.cpp b/M2/Macaulay2/e/ring.cpp index 4b9b94b483e..6145a2a6e2a 100644 --- a/M2/Macaulay2/e/ring.cpp +++ b/M2/Macaulay2/e/ring.cpp @@ -26,6 +26,7 @@ const CoefficientRingR *Ring::getCoefficientRingR() const return cR; } +// TODO: get rid of this function and directly use constructors void Ring::initialize_ring(long P0, const PolynomialRing *DR, const std::vector &heft_vector) diff --git a/M2/Macaulay2/e/schreyer-resolution/res-f4-m2-interface.cpp b/M2/Macaulay2/e/schreyer-resolution/res-f4-m2-interface.cpp index 72809f84565..35ee28147bc 100644 --- a/M2/Macaulay2/e/schreyer-resolution/res-f4-m2-interface.cpp +++ b/M2/Macaulay2/e/schreyer-resolution/res-f4-m2-interface.cpp @@ -270,6 +270,7 @@ MutableMatrix* ResF4toM2Interface::to_M2_MutableMatrix(SchreyerFrame& C, Nterm** last = newarray(Nterm*, nrows); monomial m1 = M->make_one(); + // TODO: should be exponents // FIXME: is exp a monomial or exponent vector? int* exp = new int[M->n_vars() + 1]; diff --git a/M2/Macaulay2/e/schreyer-resolution/res-moninfo-dense.hpp b/M2/Macaulay2/e/schreyer-resolution/res-moninfo-dense.hpp index 48be26f5e25..578f74e1445 100644 --- a/M2/Macaulay2/e/schreyer-resolution/res-moninfo-dense.hpp +++ b/M2/Macaulay2/e/schreyer-resolution/res-moninfo-dense.hpp @@ -397,18 +397,22 @@ class ResMonoidDense int (ResMonoidDense::*compare)(res_const_packed_monomial m, res_const_packed_monomial n) const; #endif + // TODO: remove this void variable_as_vp(int v, res_varpower_monomial result) const { - result[0] = 1; + result[0] = 1; // ATTN result[1] = v; result[2] = 1; } + // TODO: remove this int degree_of_vp(res_const_varpower_monomial a) const { return static_cast(res_varpower_monomials::weight(a, mVarDegrees)); } + // TODO: should this be merged elsewhere? + // TODO: return the consumed size void quotient_as_vp(res_const_packed_monomial a, res_const_packed_monomial b, res_varpower_monomial result) const diff --git a/M2/Macaulay2/e/unit-tests/MonomialOrderTest.cpp b/M2/Macaulay2/e/unit-tests/MonomialOrderTest.cpp new file mode 100644 index 00000000000..d22950e91aa --- /dev/null +++ b/M2/Macaulay2/e/unit-tests/MonomialOrderTest.cpp @@ -0,0 +1,172 @@ +#include +#include + +#include "MonomialOrder.hpp" + +TEST(MonomialOrder, createFromArray) +{ + std::vector>> array1 { + {MO::Lex, {3}}, + {MO::Weights, {0, 1, 2, 3, 4}}, + {MO::GRevLex, {4}}, + {MO::GRevLexWeights, {1, 1, 1, 1, 2}}, + {MO::PositionUp, {}} + }; + + /* + //day 1: std::vector array1 coming from the top level; + //day 2: MonomialOrder {BLOCK{...}, BLOCK{...}, ..., list of beginning points, is it local order?, etc.}; + //day 3: M->to_monomial(MO, exp, mon); + + // hpp + template + int encode(m, int* vec, int* iter, ...); + + // cpp + encoder(exp, monorder) + { + exp; + int *iter = monomial; + for (auto& mo : array1 + ...) + { + // e.g. mo = {MO::Weights, {0, 1, 2, 3, 4}} + switch mo.type: + ... + //iter += encode(&exp, mo.content, iter) + } + } + + inline int encode_Lex(exp, int* vec, int* iter, ...) + { + e++... + } + */ + + NewAgeMonomialOrder mo1 {12, array1}; + auto array2 = mo1.toMonomialType(); + EXPECT_EQ(array1, array2); + std::cout << "--- display monomialType---" << std::endl; + std::cout << mo1.toMonomialType() << std::endl; + std::cout << "--- debug display ---" << std::endl; + mo1.debugDisplay(std::cout); + + NewAgeMonomialOrder mo2 {12, { + {MO::Lex, {3}}, + {MO::Weights, {0, 1, 2, 3, 4}}, + {MO::GRevLex, {4}}, + {MO::GRevLexWeights, {1, 1, 1, 1, 2}}, + {MO::PositionUp, {}} + }}; + + std::cout << "mo2 = " << mo2 << std::endl; + mo2.debugDisplay(std::cout); + + + NewAgeMonomialOrder mo3 {9, { + {MO::Lex, {3}}, + {MO::RevLex, {2}}, + {MO::GroupRevLex, {4}}, + {MO::PositionUp, {}} + }}; + + std::cout << "mo3 = " << mo3 << std::endl; + mo3.debugDisplay(std::cout); + +} + +TEST(MonomialOrder, packing) +{ + std::vector>> array1 { + {MO::Packing, {4}}, + {MO::Lex, {3}}, + {MO::Weights, {0, 1, 2, 3, 4}}, + {MO::Packing, {2}}, + {MO::GRevLex, {4}}, + {MO::GRevLexWeights, {1, 1, 1, 1, 2}}, + {MO::PositionUp, {}} + }; + NewAgeMonomialOrder mo1 {12, array1}; + auto array2 = mo1.toMonomialType(); + EXPECT_EQ(array1, array2); + std::cout << mo1.toMonomialType() << std::endl; + std::cout << "--- debug display ---" << std::endl; + mo1.debugDisplay(std::cout); + + std::cout << "mo1 = " << mo1 << std::endl; + +} + +TEST(MonomialOrder, localVariables) +{ + NewAgeMonomialOrder mo1 {10, { + {MO::Packing, {4}}, + {MO::Lex, {3}}, + {MO::RevLex, {2}}, + {MO::GroupRevLex, {2}}, + {MO::PositionUp, {}}, + {MO::Weights, {0, 0,0,0,0,0,0,0,1, -1, 1}}, + {MO::Lex, {3}} + }}; + + std::cout << mo1 << std::endl; + + auto locals = mo1.localVariablesPredicate(); + std::vector locals_ans { + false, false, false, + true, true, + true, true, + false, true, false + }; + EXPECT_EQ(locals, locals_ans); + + auto inverts = mo1.invertibleVariablesPredicate(); + std::vector inverts_ans { + false, false, false, + false, false, + true, true, + false, false, false + }; + EXPECT_EQ(inverts, inverts_ans); +} + +TEST(MonomialOrder, matrixFromMonomialOrder) +{ + NewAgeMonomialOrder mo1 {6, { + {MO::Packing, {4}}, + {MO::Lex, {3}}, + {MO::PositionUp, {}}, + {MO::Weights, {0, 0,0,0,1, -1, 1}}, + {MO::Lex, {3}} + }}; + + std::vector mat; + bool base_is_revlex; + int component_dir; + int comp_is_before_row; + EXPECT_TRUE(mo1.monomialOrderingToMatrix(mat, base_is_revlex, component_dir, comp_is_before_row)); + std::cout << mat << std::endl; + EXPECT_FALSE(base_is_revlex); + EXPECT_EQ(component_dir, 1); + EXPECT_EQ(comp_is_before_row, 3); +} + +TEST(MonomialOrder, matrixFromMonomialOrder2) +{ + NewAgeMonomialOrder mo1 {6, { + {MO::Packing, {4}}, + {MO::Lex, {3}}, + {MO::PositionUp, {}}, + {MO::GRevLex, {3}} + }}; + + std::vector mat; + bool base_is_revlex; + int component_dir; + int comp_is_before_row; + EXPECT_TRUE(mo1.monomialOrderingToMatrix(mat, base_is_revlex, component_dir, comp_is_before_row)); + std::cout << mat << std::endl; + EXPECT_TRUE(base_is_revlex); + EXPECT_EQ(component_dir, 1); + EXPECT_EQ(comp_is_before_row, 3); +} + diff --git a/M2/Macaulay2/e/unit-tests/NCGroebnerTest.cpp b/M2/Macaulay2/e/unit-tests/NCGroebnerTest.cpp index 34228d8c868..0e5abe4c46a 100644 --- a/M2/Macaulay2/e/unit-tests/NCGroebnerTest.cpp +++ b/M2/Macaulay2/e/unit-tests/NCGroebnerTest.cpp @@ -14,7 +14,7 @@ #include "NCAlgebras/OverlapTable.hpp" #include "NCAlgebras/SuffixTree.hpp" #include "NCAlgebras/NCReduction.hpp" -#include "monordering.hpp" +#include "MonomialOrderOld.hpp" #include "monoid.hpp" #include "util-polyring-creation.hpp" diff --git a/M2/Macaulay2/e/unit-tests/util-polyring-creation.hpp b/M2/Macaulay2/e/unit-tests/util-polyring-creation.hpp index 95f55f69c3b..61a298847e9 100644 --- a/M2/Macaulay2/e/unit-tests/util-polyring-creation.hpp +++ b/M2/Macaulay2/e/unit-tests/util-polyring-creation.hpp @@ -6,7 +6,6 @@ #include #include "interface/monomial-ordering.h" -#include "monordering.hpp" #include "interface/ring.h" #include "interface/aring.h" #include "monoid.hpp" diff --git a/M2/Macaulay2/m2/engine.m2 b/M2/Macaulay2/m2/engine.m2 index 842007a53fb..ac31de138ea 100644 --- a/M2/Macaulay2/m2/engine.m2 +++ b/M2/Macaulay2/m2/engine.m2 @@ -41,181 +41,7 @@ gcd(RawMonomial,RawMonomial) := (x,y) -> rawGCD(x,y) setAttribute(RawMonomialOrdering,ReverseDictionary,symbol RawMonomialOrdering) RawMonomialOrdering.synonym = "raw monomial ordering" -Eliminate = new SelfInitializingType of BasicList -new Eliminate from ZZ := (Eliminate,n) -> Eliminate {n} -expression Eliminate := v -> ( - if #v === 1 - then new FunctionApplication from {Eliminate, v#0} - else new FunctionApplication from {Eliminate, toList v}) -ProductOrder = new SelfInitializingType of BasicList - -isSmall := i -> class i === ZZ and i < 2^15 and i > -2^15 -isCount := i -> class i === ZZ and i >= 0 and i < 2^15 -isListOfIntegers = x -> instance(x, List) and all(x,i -> class i === ZZ) -isListOfListsOfIntegers = x -> instance(x, List) and all(x,isListOfIntegers) -listZ = listZZ = v -> if isListOfIntegers(v = toList splice v) then v else error "expected a list of integers" -checkCount := i -> if not isCount i then error "expected a small positive integer" - -fixup1 := method(Dispatch => Thing) -- stage 1, everything except Tiny and Small -fixup2 := method(Dispatch => Thing) -- stage 2, Tiny and Small -err := o -> error ("unrecognized ordering item " | toString o) - -deglist := {} -- filled in below each time -getdegs := (m,n) -> ( - -- get m-th through n-th degree, using 1 when there aren't enough - join ( take(deglist, {m,n}), (n - (#deglist - 1)):1) - ) -numvars := 0 -- re-initialized below each time -varcount := 0 -- re-initialized below each time -MonSize := 0 -- re-initialized below each time -invert := false -- re-initialized below each time -bump := n -> varcount = varcount + n - -processMonSize := monsize -> ( - oldMonSize := MonSize; - monsize === null or ( - if class monsize =!= ZZ then error "expected an integer as MonomialSize option"; - if monsize <= 8 then MonSize = 8 - else if monsize <= 16 then MonSize = 16 - else MonSize = 32; - MonSize != oldMonSize ) ) - -optionSizes := hashTable { - -- this table omits the default of 32 on purpose, so that MonSize=0 will work out. - (Lex,8) => LexTiny, - (Lex,16) => LexSmall, - (GRevLex,8) => GRevLexTiny, - (GRevLex,16) => GRevLexSmall - } -fixSize := key -> if optionSizes#?(key,MonSize) then optionSizes#(key,MonSize) else key -optionInverses := hashTable { - Lex => key -> GroupLex, - RevLex => key -> GroupRevLex, - GroupRevLex => key -> GroupRevLex, - GroupLex => key -> GroupLex --- GRevLex => key -> ( --- newkey := GroupRevLex; --- stderr << "--warning: for rawMonomialOrdering, replacing monomial ordering " << key << " by " << newkey << endl; --- newkey) - } -fix1 := key -> ( - if invert then ( - if not optionInverses#?key then error ("Inverses => true not compatible with ",toString key); - optionInverses#key key) - else key) -intOption := (key,n) -> ( - key = fix1 key; - checkCount n; - bump n; - key => n) -fixList := v -> if instance(v,List) then spliceInside v else v -grevOption := (key,v) -> ( - key = fix1 key; - if instance(v,ZZ) - then v = getdegs ( varcount, varcount + v - 1 ) - else ( - v = fixList v; - if not isListOfIntegers v then error "expected an integer or a list of integers"; - ); - bump(#v); - key => v) -optionFixes := hashTable { - Weights => (key,val) -> key => fixList val, - Position => (key,val) -> key => val, - MonomialSize => (key,val) -> if processMonSize val then key => val, - Lex => intOption, - RevLex => intOption, - GroupLex => intOption, - GroupRevLex => intOption, - NCLex => intOption, - GRevLex => grevOption - } - -ordOption := o -> if numvars > varcount then fixup1 ( o => numvars - varcount ) -symbolFixes := hashTable { - GLex => o -> (Weights => toList (numvars:1), fixup1 Lex), - RevLex => ordOption, - GRevLex => ordOption, - Lex => ordOption - } - -fixup1 Thing := err -fixup1 List := v -> splice (toSequence v / fixup1) -fixup1 Sequence := v -> v / fixup1 -fixup1 Eliminate := e -> Weights => toList (e#0 : 1) -fixup1 ProductOrder := o -> fixup1(toSequence o / (i -> GRevLex => i)) -fixup1 Symbol := o -> if symbolFixes#?o then symbolFixes#o o else err o -fixup1 ZZ := n -> fixup1( GRevLex => n ) -fixup1 Option := o -> ( - key := o#0; - val := o#1; - if optionFixes#?key then optionFixes#key (key,val) - else error ("unrecognized ordering item keyword : " | toString key) - ) -intOption2 := (key,n) -> ( - bump n; - fixSize key => n - ) -grevOption2 := (key,v) -> fixSize key => v -optionFixes2 := hashTable { - MonomialSize => (key,val) -> (processMonSize val; null), - Lex => intOption2, - GRevLex => grevOption2 - } - -fixup2 Thing := identity -fixup2 List := v -> toSequence v / fixup2 -fixup2 Sequence := v -> v / fixup2 -fixup2 Option := o -> ( - key := o#0; - val := o#1; - if optionFixes2#?key then optionFixes2#key (key,val) else o - ) - -oops = () -> error "expected a list of integers or a list of lists of small integers" - -processWeights = (nvars,weights) -> ( - weights = spliceInside toList weights; - if #weights > 0 and isListOfIntegers weights then weights = {weights} - else ( - weights = apply(weights, w -> ( if not instance(w,List) then oops(); deepSplice w)); - if not isListOfListsOfIntegers weights then oops(); - ); - scan(weights, - wt -> ( - if # wt != nvars - then error("Weights: expected weight vector of length ",toString nvars," but got ",toString (#wt)))); - weights); - -makeMonomialOrdering = (monsize,inverses,nvars,degs,weights,ordering) -> ( - -- 'monsize' is the old MonomialSize option, usually 8 or 16, or 'null' if unset - -- 'inverses' is true or false, and tells whether the old "Inverses => true" option was used. - -- 'nvars' tells the total number of variables. Any extra variables will be ordered with GRevLex or GroupLex. - -- 'degs' is a list of integers, the first components of the multi-degrees of the variables - -- if it's too short, additional degrees are taken to be 1. Could be an empty list. - -- 'weights' is the list of integers or the list of lists of integers provided by the user under - -- the *separate* Weights option. Could be an empty list. - -- 'ordering' is a list of ordering options, e.g., { Lex => 4, GRevLex => 4 } - -- If it's not a list, we'll make a list of one element from it. - if monsize === null then monsize = null; - ordering = {MonomialSize => monsize, ordering}; - invert = inverses; - if not isListOfIntegers degs then error "expected a list of integers"; - deglist = degs; - varcount = 0; - MonSize = 0; - numvars = nvars; - weights = processWeights(nvars,weights); - ordering = join(weights / (i -> Weights => i), ordering); - t':= toList nonnull splice fixup1 ordering; - if varcount < nvars then t' = append(t',fixup1(GRevLex => nvars - varcount)); - if not any(t', x -> class x === Option and x#0 === Position) then t' = append(t', Position => Up); - MonSize = 0; - varcount = 0; - t := toList nonnull fixup2 t'; - logmo := new FunctionApplication from {rawMonomialOrdering,t}; - (t,t',value logmo, logmo)) - +-- TODO: also add rawJoinMonomialOrdering RawMonomialOrdering ** RawMonomialOrdering := RawMonomialOrdering => rawProductMonomialOrdering -- used for debugging mgb interface, moved from ofcm.m2 diff --git a/M2/Macaulay2/m2/enginering.m2 b/M2/Macaulay2/m2/enginering.m2 index 350c123ce85..f2398f71d94 100644 --- a/M2/Macaulay2/m2/enginering.m2 +++ b/M2/Macaulay2/m2/enginering.m2 @@ -485,7 +485,14 @@ leadCoefficient RingElement := RingElement => (f) -> ( k := coefficientRing ring f; promote(rawLeadCoefficient(raw k, raw f),k)) -degree RingElement := f -> if f == 0 then -infinity else rawMultiDegree raw f +-- TODO: would be better in the engine +degree RingElement := f -> if f == 0 then -infinity else ( + deg := rawMultiDegree raw f; + -- TODO: perhaps there should be a flag in the ring, + -- determining whether reduction should happen? + if not (M := monoid ring f).?degreeGroup + or isFreeModule(G := M.degreeGroup) then deg + else reduceDegree(G, deg)) leadTerm RingElement := RingElement => (f) -> someTerms(f,0,1) diff --git a/M2/Macaulay2/m2/modules.m2 b/M2/Macaulay2/m2/modules.m2 index f915b630ab6..a94d64110d0 100644 --- a/M2/Macaulay2/m2/modules.m2 +++ b/M2/Macaulay2/m2/modules.m2 @@ -263,6 +263,7 @@ toExternalString Module := M -> toString describe M -- TODO: where is it set before being cached? degrees Module := -*(cacheValue symbol degrees) (*-N -> ( r := degreeLength(R := ring N); + -- TODO: this check should be unnecessary once issue #3007 is fixed if r == 0 then toList(numgens N : {}) else ( degs := pack(r, rawMultiDegree raw cover N); if not (M := monoid R).?degreeGroup @@ -288,6 +289,11 @@ Ring ^ List := Module => (R, degs) -> ( else if isListOfListsOfIntegers degs then ( if any(degs, deg -> degrk != #deg) then error("expected each multidegree to be of length ", degrk)) else error "expected a list of integers or a list of lists of integers"; + -- TODO: perhaps there should be a flag in the ring, + -- determining whether reduction should happen? + degs = if not (M := monoid R).?degreeGroup + or isFreeModule(G := M.degreeGroup) then degs + else apply(degs, reduceDegree_G); -- then flatten the args fdegs := toSequence flatten degs; new Module from (R, rawFreeModule(R.RawRing, if #fdegs === 0 then #degs else fdegs))) diff --git a/M2/Macaulay2/m2/monoids.m2 b/M2/Macaulay2/m2/monoids.m2 index 64b3a0a7adc..98cd2665f8a 100644 --- a/M2/Macaulay2/m2/monoids.m2 +++ b/M2/Macaulay2/m2/monoids.m2 @@ -37,6 +37,12 @@ indices := (M, variables) -> apply(variables, x -> ( sameMonoid := (x, y) -> if (M := class x) === class y then M else error "expected elements in the same monoid" +isSmall := i -> class i === ZZ and i < 2^15 and i > -2^15 +isCount := i -> class i === ZZ and i < 2^15 and i >= 0 +isListOfIntegers = L -> instance(L, List) and all(L, i -> instance(i, ZZ)) +isListOfListsOfIntegers = L -> instance(L, List) and all(L, isListOfIntegers) +listZ = listZZ = v -> if isListOfIntegers(v = toList splice v) then v else error "expected a list of integers" + ----------------------------------------------------------------------------- -- MonoidElement type declarations and basic methods ----------------------------------------------------------------------------- @@ -376,7 +382,9 @@ degreesMonoid List := memoize lookup(degreesMonoid, List) -- findHeft ----------------------------------------------------------------------------- +-- TODO: this should check something stronger now that the grading group is arbitrary checkHeft = (degs, heftvec) -> all(degs, d -> sum(d, heftvec, times) > 0) + -- vector that zeros the torsion part of the degrees -- TODO: what should happen when there are zero-divisors that are not torsion? -- compare with freeComponents and torsionComponents in basis.m2 @@ -456,6 +464,156 @@ processDegrees = (degs, degrk, group, nvars) -> ( (degs, degrk, group)) else error "expected Degrees option to be list of degrees") +----------------------------------------------------------------------------- +-- Monomial Orderings +----------------------------------------------------------------------------- + +Eliminate = new SelfInitializingType of BasicList +new Eliminate from ZZ := (Eliminate, n) -> Eliminate {n} +expression Eliminate := v -> ( if #v === 1 + then new FunctionApplication from {Eliminate, v#0} + else new FunctionApplication from {Eliminate, toList v}) +ProductOrder = new SelfInitializingType of BasicList + +-- + +-- TODO: remove this +processMOkey := (opts, key) -> ( + -- omit the default of 32 on purpose, so that opts#"monsize"=0 will work out. + if key === Lex then ( + if opts#"monsize" === 8 then LexTiny else if opts#"monsize" === 16 then LexSmall else Lex ) else + if key === GRevLex then ( + if opts#"monsize" === 8 then GRevLexTiny else if opts#"monsize" === 16 then GRevLexSmall else GRevLex ) else key) + +bump := (opts, n) -> (opts#"offset" = opts#"offset" + n; n) + +intOption := (opts, k, n) -> ( + if not isCount n then error "expected small positive integers in monomial orders"; + k => bump(opts, n)) + +processMonSize := (opts, monsize) -> opts#"monsize" = ( + if not isCount monsize then error "expected a small positive integer for MonomialSize"; + if monsize <= 8 then 8 else if monsize <= 16 then 16 else 32) + +processMOfunctions := hashTable { + Weights => (opts, k, v) -> k => if instance(v, List) then spliceInside v else v, + Position => (opts, k, v) -> k => (opts.Position = v), + -- FIXME: MonomialSize => 16, Lex => 2 has two MonomialSizes + MonomialSize => (opts, k, v) -> k => processMonSize(opts, v), + NCLex => intOption, -- TODO: remove? + Lex => intOption, + RevLex => intOption, + GroupLex => intOption, + GroupRevLex => intOption, + GRevLex => (opts, k, v) -> ( + if instance(v, ZZ) + -- get m-th through n-th degree + then v = take(opts#"heftdegs", {opts#"offset", opts#"offset" + v - 1}) + else ( + v = if instance(v, List) then spliceInside v else v; + if not isListOfIntegers v then error "expected an integer or a list of integers"; + ); + bump(opts, #v); + k => v) + } + +processMO = method() -- stage 1, everything except Tiny and Small +processMO(MutableHashTable, Thing) := (opts, mo) -> error("unrecognized monomial ordering block: ", toString mo) + +fillMO := (opts, mo) -> if opts#"nvars" > opts#"offset" then processMO(opts, mo => opts#"nvars" - opts#"offset") +symbolFixes := hashTable { + Lex => fillMO, + GLex => (opts, mo) -> (Weights => toList(opts#"nvars" : 1), processMO(opts, Lex)), + RevLex => fillMO, + GRevLex => fillMO, + } +processMO(MutableHashTable, Symbol) := (opts, s) -> if symbolFixes#?s then symbolFixes#s(opts, s) else + error("unrecognized monomial ordering keyword: ", toString s) +processMO(MutableHashTable, Sequence) := (opts, v) -> v / (x -> processMO(opts, x)) +processMO(MutableHashTable, List) := (opts, v) -> splice(toSequence v / (x -> processMO(opts, x))) +processMO(MutableHashTable, ZZ) := (opts, n) -> processMO(opts, GRevLex => n) -- default setting +processMO(MutableHashTable, ProductOrder) := (opts, v) -> processMO(opts, toSequence v / (x -> processMO(opts, x))) +processMO(MutableHashTable, Eliminate) := (opts, w) -> Weights => toList(w#0 : 1) + +-- check that Inverses => true is compatible with the monomial order +checkInverses := key -> ( + if key === Lex then GroupLex else + if key === RevLex then GroupRevLex else + if isMember(key, {GroupLex, GroupRevLex}) then key else + if isMember(key, {Weights, Position, MonomialSize}) then key else + error("Inverses => true not compatible with ", toString key)) + +processMO(MutableHashTable, Option) := (opts, mo) -> ( + (key, val) := toSequence mo; + if opts.Inverses then key = checkInverses key; + if processMOfunctions#?key then processMOfunctions#key(opts, key, val) + else error("unrecognized monomial ordering block: ", toString key)) + +processWeights = (nvars, vecs) -> ( + vecs = spliceInside toList vecs; + -- e.g. allows {}, {1,1,1,1}, {{1,1,1,1}}, {{{0,0},1,1},{1,1,{0,0}}} + vecs = if isListOfIntegers vecs and #vecs > 0 then {vecs} else apply(vecs, deepSplice); + if not isListOfListsOfIntegers vecs then error "expected a list of integers or a list of lists of integers"; + apply(vecs, wt -> if #wt < nvars then join(wt, nvars - #wt : 0) else if #wt == nvars then wt + else error("encountered weight vector of length ", #wt))) + +-- TODO: remove these +stage2Options := hashTable { + MonomialSize => (opts, k, v) -> ( processMonSize(opts, v); null ), + Lex => (opts, k, n) -> processMOkey(opts, k) => bump(opts, n), + GRevLex => (opts, k, v) -> processMOkey(opts, k) => v + } + +stage2 := method() -- stage 2, Tiny and Small +stage2(MutableHashTable, Thing) := (opts, mo) -> mo +stage2(MutableHashTable, List) := (opts, v) -> toSequence v / (x -> stage2(opts, x)) +stage2(MutableHashTable, Sequence) := (opts, v) -> v / (x -> stage2(opts, x)) +stage2(MutableHashTable, Option) := (opts, mo) -> ( + (key, val) := toSequence mo; + if stage2Options#?key then stage2Options#key(opts, key, val) else mo) + +newMonomialOrder = (ordering, monsize, inverses, weights, heftdegs, nvars) -> ( + -- 'ordering' is a list of ordering options, e.g., { Lex => 4, GRevLex => 4 } + -- If it's not a list, we'll make a list of one element from it. + -- 'monsize' is the old MonomialSize option, usually 8 or 16, or 'null' if unset + -- 'inverses' is true or false, and tells whether the old "Inverses => true" option was used. + -- 'nvars' tells the total number of variables. Any extra variables will be ordered with GRevLex or GroupLex. + -- 'degs' is a list of integers, the first components of the multi-degrees of the variables + -- if it's too short, additional degrees are taken to be 1. Could be an empty list. + -- 'weights' is the list of integers or the list of lists of integers provided by the user under + -- the *separate* Weights option. Could be an empty list. + if not isListOfIntegers heftdegs then error "expected a list of integers"; + if not instance(ordering, List) then ordering = {ordering}; + if instance(monsize, ZZ) then ordering = prepend(MonomialSize => monsize, ordering); + -- combine these two + weights = processWeights(nvars, weights); + ordering = join(weights / (v -> Weights => v), ordering); + opts := new MutableHashTable from { + Inverses => inverses, + Position => null, + "heftdegs" => heftdegs, + "monsize" => 0, + "offset" => 0, + "nvars" => nvars}; + ordering = toList nonnull splice processMO(opts, ordering); + -- use default GRevLex for remaining variables and Position => Up if not specified + if opts#"offset" < nvars then ordering = append(ordering, processMO(opts, GRevLex)); + if opts#Position === null then ordering = append(ordering, processMO(opts, Position => Up)); + -- second pass + opts#"monsize" = opts#"offset" = 0; + MOopts := toList nonnull stage2(opts, ordering); + rawMO := rawMonomialOrdering MOopts; + (MOopts, ordering, rawMO)) + +-* +restart +errorDepth=1 +R = ZZ/101[a,b,c,d]/(a^4+b^4+c^4+d^4) + X = Proj R + errorDepth=1 + result = table(3,3,(p,q) -> timing ((p,q) => rank HH^q(cotangentSheaf(p,X)))) +*- + ----------------------------------------------------------------------------- makeVars = (n, var) -> toList( @@ -560,15 +718,14 @@ newMonoid = opts -> ( degs := M.degrees = opts.Degrees; nvars := M.numgens = #varlist; M.vars = M.generators = apply(nvars, i -> new M from rawVarMonomial(i, 1)); + -- TODO: why is this set to degrk:0? heftvec := M.heft = if opts.Heft === null then toList(degrk:0) else opts.Heft; + heftdegrees := if degrk == 0 or opts.Heft === null then toList(nvars:1) + else apply(degs, d -> if (w := dotprod(d, heftvec)) > 0 then w else 1); -- see engine.m2 - (MOopts, MOoptsint, rawMO, logMO) := makeMonomialOrdering( - opts.MonomialSize, opts.Inverses, nvars, - if degreeLength M == 0 or opts.Heft === null then toList(nvars:1) - else apply(degs, d -> if (w := dotprod(d, heftvec)) > 0 then w else 1), - opts.Weights, opts.MonomialOrder); + (MOopts, MOoptsint, rawMO) := newMonomialOrder( + opts.MonomialOrder, opts.MonomialSize, opts.Inverses, opts.Weights, heftdegrees, nvars); M.RawMonomialOrdering = rawMO; - M#"raw creation log" = new Bag from {logMO}; M#"options" = MOopts; -- these are exactly the arguments given to rawMonomialOrdering -- these are the options given to rawMonomialOrdering minus MonomialSize, -- except Tiny and Small aren't there yet and GRevLex=>nvars hasn't been expanded @@ -584,9 +741,13 @@ newMonoid = opts -> ( madeTrivialMonoid = true; rawMonoid()) else ( + -- TODO: when degree group has torsion, this isn't correct yet + -- TODO: pass a ZZ-module, potentially with torsion, instead? M.degreesRing = if opts.Heft =!= null then degreesRing heftvec else degreesRing degrk; -* shouldn't really be needed *- M.degreesMonoid = monoid M.degreesRing; rawMonoid( + -- TODO: fix the order + -- TODO: also pass heftdegrees M.RawMonomialOrdering, raw M.degreesRing, toSequence M.generators / toString, diff --git a/M2/Macaulay2/tests/engine/raw.m2 b/M2/Macaulay2/tests/engine/raw.m2 index 176279143c5..96dc68d1705 100644 --- a/M2/Macaulay2/tests/engine/raw.m2 +++ b/M2/Macaulay2/tests/engine/raw.m2 @@ -18,6 +18,10 @@ rawMonomialOrdering {Lex => 1, LexSmall => 2, LexTiny => 4} rawMonomialOrdering {GRevLex => {1,1,1}, GRevLexSmall => {2,2}, GRevLexTiny => {4,4,4,4}} rawMonomialOrdering {GroupLex => 3, NCLex => 4, Weights => {3,4,5}} +-- FIXME: fix this file +makeMonomialOrdering = (monsize,inverses,nvars,heftdegs,weights,ordering) -> + processMO(ordering,monsize,inverses,weights,heftdegs,nvars) + makeMonomialOrdering( null, false, 0, {1,2,3,4,5,6,7,8}, {}, {} ) makeMonomialOrdering( null, false, 4, {1,2,3,4,5,6,7,8}, {}, {} ) makeMonomialOrdering( null, false, 10, {1,2,3,4,5,6,7,8}, {}, {} ) diff --git a/M2/Macaulay2/tests/normal/monoids.m2 b/M2/Macaulay2/tests/normal/monoids.m2 index 9c29f77ef94..3b0612f7956 100644 --- a/M2/Macaulay2/tests/normal/monoids.m2 +++ b/M2/Macaulay2/tests/normal/monoids.m2 @@ -9,7 +9,7 @@ monoid[(a,b,c),(x,y,z)] monoid[{a,b,c},{x,y,z}] monoid[vars(0..2,23..25)] monoid[vars(0..2),vars(23..25)] -monoid[a_0..a_3,b_0..b_3] -- FIXME: this would fail if QQ[a,b,c] was defined before +monoid[a_0..a_3,b_0..b_3] --- monoid[a,b, WeylAlgebra => a=>b] monoid[a,b, WeylAlgebra => 0=>1] @@ -106,10 +106,28 @@ assert(heft B == {0,1,0,1}) S = (ZZ/3) A T = (ZZ/3) B +--- degrees monoids or rings with torsion in degree group +-- FIXME: these two should have torsion, and probably not be ordered either +degreesMonoid A +degreesRing A + -- module over ring with torsion degree group assert(degrees S^{{3,3}} == {{0,-3}}) assert(S^{{3,3}} == S^{{0,3}}) --- assert(S^{{3,3}} === S^{{0,3}}) -- TODO: should this be true? +assert(S^{{3,3}} === S^{{0,3}}) -- TODO: do in the engine + +-- elements of ring with torsion degree group +--restart +--debug Core +--G = coker matrix{{3}, {0}} +--S = QQ[x,y,z, +-- DegreeGroup => G, +-- Degrees => {{1,1}, {-1,1}, {0,1}}] +assert(degree S_0^3 == {0,3}) -- TODO: do in the engine +--prepend(-sum(d = degree S_2, heft S, times), d) +--code(degree, RingElement) +--code(degrees, Module) +--rawMultiDegree raw S_0^3 --- testing graphIdeal f = map(T, S, {T_0, T_1, T_2}, DegreeMap => d -> join(d, {0,0})) @@ -124,10 +142,12 @@ J = graphIdeal g assert isHomogeneous J assert(degreeGroup monoid S == degreeGroup monoid ring J) ---- TODO +--- test degreeGroup of rings use S R = S/(x-y) -monoid R +assert(degreeGroup S == G) +assert(degreeGroup R == G) +-- degreeGroup frac R -- TODO: what should this be? -- used to be in debugme.m2 A = ZZ[a] @@ -174,6 +194,9 @@ assert(degrees M == splice {{1, 0}, 2:{0, 1}}) --- test passing degrees from generators of the degree group H = subquotient(ambient A, relations G) M = monoid[a,b,c, DegreeGroup => H] +-- TODO: should we allow degrees to coefficients of the generators of H? +-- in that case, freeComponents and torsionComponents need to be fixed, because numgens H changes +-- M = monoid[a,b,c, DegreeGroup => H, Degrees => entries id_(ZZ^3)] assert(degreeGroup M == G) assert(degrees M == degrees S) diff --git a/M2/cmake/check-libraries.cmake b/M2/cmake/check-libraries.cmake index c06b88aaa90..3a7d5aa704c 100644 --- a/M2/cmake/check-libraries.cmake +++ b/M2/cmake/check-libraries.cmake @@ -327,6 +327,7 @@ if(FACTORY_FOUND) set(CMAKE_REQUIRED_INCLUDES "${FACTORY_INCLUDE_DIR}") # whether factory was built with --enable-streamio check_cxx_source_compiles([[#include + #include int main(){Variable x; x = Variable(); std::cout << x;return 0;}]] FACTORY_STREAMIO) else() unset(FACTORY_STREAMIO CACHE)