diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 00000000000..1fd7aa42b90 --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,326 @@ +# Macaulay2 Engine Development Guide + +## Project Overview + +Macaulay2 is a software system for algebraic geometry research. This guide focuses on the +C++ engine in `M2/Macaulay2/e/`, which is the computational core. + +## Repository Structure + +``` +M2/ Top-level build directory + CMakeLists.txt CMake entry point (requires 3.24+, C++17) + configure.ac Autotools configuration + Macaulay2/ + e/ Engine (C++ computational core, ~500 source files) + CMakeLists.txt Builds M2-engine static library and M2-unit-tests + Makefile.in Autotools build for engine + unit-tests/ GoogleTest-based unit tests + Makefile.in Autotools build for unit tests + Makefile.files List of test source files (shared by autotools) + interface/ Public C/C++ API (ring, matrix, groebner, etc.) + f4/ Faugere F4 algorithm + gb-f4/ GB F4 computation interface + NCAlgebras/ Non-commutative algebra support + schreyer-resolution/ Schreyer resolution (F4-style linear algebra) + bibasis/ Boolean Involutive Groebner Bases + d/ Interpreter + m2/ M2 language libraries + packages/ User packages + submodules/ Git submodules (googletest, flint, bdwgc, etc.) + cmake/ CMake modules (check-libraries, build-libraries, etc.) + libraries/ External library sources +``` + +## Engine Unit Tests + +### Location and Framework + +Tests are in `M2/Macaulay2/e/unit-tests/` using **GoogleTest** (gtest). + +### Key Files + +- `testMain.cpp` - Test entry point; calls `IM2_initialize()` then `RUN_ALL_TESTS()` +- `M2-cpp-replacement.cpp` - Stub for `system_interrupted()` (avoids linking interpreter) +- `ARingTest.hpp` - Templated test helpers for arithmetic ring operations (negate, add, + subtract, multiply, divide, reciprocal, power, axioms, coercions) +- `RingTest.hpp` - Test helpers for Ring interface +- `DMatTest.hpp` - Test helpers for dense matrices +- `RingElem.hpp/cpp` - Lightweight value-semantics wrapper for ring elements (see below) +- `util-polyring-creation.hpp/cpp` - Helpers for creating rings in tests: + - `simplePolynomialRing(p, names)` - polynomial ring over ZZ/p (or QQ if p=0) + - `simpleWeylAlgebra(p, names, comms, derivs)` - Weyl algebra + - `degreeRing(n)` - degree ring with n variables + +### Existing Test Files + +| Test File | What It Tests | +|-----------|---------------| +| `ARingZZTest.cpp` | Integers (flint) | +| `ARingZZpTest.cpp` | Z/p (multiple implementations) | +| `ARingQQFlintTest.cpp` | Rationals (flint) | +| `ARingQQGmpTest.cpp` | Rationals (GMP) | +| `ARingRRTest.cpp` | Machine reals | +| `ARingCCTest.cpp` | Machine complex | +| `ARingRRRTest.cpp` | Arbitrary-precision reals (MPFR) | +| `ARingRRiTest.cpp` | Arbitrary-precision real intervals (MPFI) | +| `ARingCCCTest.cpp` | Arbitrary-precision complex | +| `RingZZTest.cpp` | ZZ via Ring interface | +| `RingZZpTest.cpp` | ZZp via Ring interface | +| `RingQQTest.cpp` | QQ via Ring interface | +| `RingRRRTest.cpp` | RRR via Ring interface | +| `RingCCCTest.cpp` | CCC via Ring interface | +| `RingTowerTest.cpp` | Tower of polynomial rings | +| `DMatZZpTest.cpp` | Dense matrices over ZZp | +| `MonoidTest.cpp` | Monoid operations | +| `PolyRingTest.cpp` | Polynomial ring operations | +| `NCGroebnerTest.cpp` | Non-commutative Groebner bases | +| `WeylAlgebraTest.cpp` | Weyl algebra creation, commutators, binomial, multinomial, fromString | +| `NewF4Test.cpp` | New F4 algorithm | +| `ResTest.cpp` | Resolutions | +| `MatrixIOTest.cpp` | Matrix I/O | +| `SubsetTest.cpp` | Subset operations | +| `PointArray.cpp` | Point array operations | +| `basics-test.cpp` | Buffer, utility functions | +| `fromStream.cpp` | Stream parsing | + +### Excluded Test Files + +| Test File | Why Excluded | +|-----------|-------------| +| `ARingGFTest.cpp` | ARingGFFlint API changed: constructor now requires `PolynomialRing` + primitive element (was `(int p, int n)`), and `cardinality()` method is missing | + +### Linking Requirements + +The unit tests link against the M2 engine library. The only additional file needed is +`M2-cpp-replacement.cpp`, which stubs out `system_interrupted()` — a function normally +supplied by the Macaulay2 executable. `testMain.cpp` calls `IM2_initialize()` to +handle engine initialization (including GC). + +### Writing New Tests + +Pattern for arithmetic ring tests (standalone `TEST` macros): +```cpp +#include +#include "aring-zz-flint.hpp" // or whichever ring +#include "ARingTest.hpp" // templated test helpers + +TEST(MyRing, create) { + M2::ARingZZ R; + // ... test ring properties +} + +TEST(MyRing, arithmetic) { + M2::ARingZZ R; + testCoercions(R); + testNegate(R, ntrials); + testAdd(R, ntrials); + // ... +} +``` + +Pattern for tests using polynomial rings (use `util-polyring-creation.hpp`): +```cpp +#include "util-polyring-creation.hpp" +const PolynomialRing* R = simplePolynomialRing(101, {"a", "b", "c"}); +const WeylAlgebra* W = simpleWeylAlgebra(0, {"x","y","Dx","Dy"}, {0,1}, {2,3}); +``` + +Pattern for test fixtures with shared setup (use `TEST_F`): +```cpp +class MyTest : public ::testing::Test { + protected: + SomeRing* R = nullptr; + void SetUp() override { + R = /* create ring */; + } +}; + +TEST_F(MyTest, someTest) { + // R is available here, freshly created for each test +} +``` + +### Testing Private/Protected Members + +Use the friend class pattern. Add `friend class FooTestAccessor;` to the class header, +then define a test accessor in the test file with static forwarding methods: +```cpp +// In the .hpp file: +class Foo { + friend class FooTestAccessor; + // ... +}; + +// In the test .cpp file: +class FooTestAccessor { + public: + static int privateMethod(const Foo* f, int arg) { + return f->privateMethod(arg); + } +}; +``` +See `WeylAlgebraTest.cpp` and `weylalg.hpp` for a concrete example. + +### RingElem — Lightweight Value Wrapper for Tests + +`RingElem` (in `unit-tests/RingElem.hpp`) wraps `const Ring*` + `ring_elem` with +value semantics. Operators return values (not pointers), making test code concise: +```cpp +#include "RingElem.hpp" +auto x = RingElem::var(R, 0); +auto y = RingElem::var(R, 1); +auto one = RingElem::fromInt(R, 1); +EXPECT_EQ(x * y - y * x, one); // value comparison, prints elements on failure + +// Parse from string (requires explicit ^ and * in polynomial syntax): +auto f = RingElem::fromString(R, "x^2+3*x*y-1"); + +// Scalar multiplication: +auto g = x * 3; // RingElem * long +auto h = 3 * x; // long * RingElem +``` + +Factories: `RingElem::var(R, i)`, `RingElem::fromInt(R, n)`, `RingElem::fromString(R, s)`. +Arithmetic: `+`, `-`, `*`, `/`, unary `-`, `.power(n)`. +Output: `to_string()`, `operator<<` for gtest diagnostics. + +**fromString format**: Uses `parseBasicPoly` from `BasicPoly.hpp`. Requires `*` between +factors and `^` for exponents (e.g. `"3*x^2*y-1"`). Note: `to_string()` outputs a +different format (`x2y-1` without `*` or `^`), so round-tripping is not yet supported. + +**Naming convention**: Member fields use `m` prefix (e.g., `mRing`, `mValue`). + +### RingElement (Legacy Interface) + +`RingElement` operators (`*`, `+`, `-`, `/`) return `RingElement*` (pointers), not values. +Prefer `RingElem` for new test code. `RingElement` is still used in the interpreter interface. +```cpp +RingElement *x = new RingElement(W, W->var(0)); +RingElement *product = (*x) * (*y); // dereference, then multiply +EXPECT_TRUE(product->is_equal(*expected)); +``` +`is_equal` checks ring pointer equality first — both elements must come from the same +ring object. Use `IM2_Ring_QQ()` (QQGMP) consistently, not `rawARingQQFlint()`, when +creating rings whose elements will be compared. + +### M2_arrayint in Tests + +Many engine functions take `M2_arrayint` (a GC-allocated array). Convert from +`std::vector` using: +```cpp +#include "util.hpp" +M2_arrayint arr = stdvector_to_M2_arrayint(std::vector{2, 3, 5}); +``` + +## Building and Running Tests + +Build configurations are defined in `M2/BUILD/mike/Makefile`. Existing builds live +under `M2/BUILD/mike/builds.tmp/`. + +### CMake Build + +The active cmake build directory is: +`M2/BUILD/mike/builds.tmp/cmake-appleclang` (RelWithDebInfo, Ninja) + +```sh +# Build unit tests: +ninja -C M2/BUILD/mike/builds.tmp/cmake-appleclang M2-unit-tests + +# Run all engine unit tests: +ctest --test-dir M2/BUILD/mike/builds.tmp/cmake-appleclang -R unit-tests + +# Run a specific test by name: +ctest --test-dir M2/BUILD/mike/builds.tmp/cmake-appleclang -R "WeylAlgebra" --output-on-failure + +# Or run the executable directly: +M2/BUILD/mike/builds.tmp/cmake-appleclang/Macaulay2/e/M2-unit-tests +# With a gtest filter: +M2/BUILD/mike/builds.tmp/cmake-appleclang/Macaulay2/e/M2-unit-tests --gtest_filter="WeylAlgebra*" +``` + +The CMake build uses `gtest_discover_tests()` with prefix `unit-tests:`. + +### Autotools Build + +The active autotools build directory is: +`M2/BUILD/mike/builds.tmp/arm64-appleclang` + +```sh +# From the autotools build directory, in Macaulay2/e/unit-tests/: +gmake -k check +# Or run the executable directly: +./testMain +``` + +If `e/unit-tests/Makefile.in` changes, regenerate the build Makefile: +```sh +# From the autotools build root: +./config.status Macaulay2/e/unit-tests/Makefile +``` + +### Adding a New Test File + +1. Add the `.cpp` file to `e/unit-tests/` +2. Add it to `e/CMakeLists.txt` in the `add_executable(M2-unit-tests ...)` section +3. Add it to `e/unit-tests/Makefile.files` in the `UNITTEST_CCFILES` list +4. Build and run tests in **both** cmake and autotools to verify + +### Keeping Builds in Sync + +The CMake and autotools builds maintain separate lists of test files: +- CMake: `e/CMakeLists.txt` — the `add_executable(M2-unit-tests ...)` block +- Autotools: `e/unit-tests/Makefile.files` — the `UNITTEST_CCFILES` variable + +Both should include the same set of test files. Currently 199 tests pass in both builds. + +## Engine Dependencies + +All of the following are linked to the `M2-engine` target (see `e/CMakeLists.txt`). + +### Submodule Dependencies (built as part of engine) +- **memtailor** - Special-purpose memory allocators +- **mathic** - Symbolic algebra data structures +- **mathicgb** - Signature Groebner bases library +- **googletest** - Unit testing framework (for M2-unit-tests only) + +### Header-Only Libraries +- **Eigen3** (3.4.0+) - Linear algebra templates + +### Libraries (found via pkg-config) +- **FFLAS_FFPACK** (2.4.3+) - Finite field linear algebra routines (needs LAPACK, GIVARO) +- **GIVARO** (4.1.1+) - Prime field and algebraic computations + +### Libraries (found via find_package / FindXxx.cmake) +- **GMP** (6.0.0+) - GNU multiprecision arithmetic +- **MPFR** (4.0.1+) - Multiprecision floating-point (needs GMP) +- **MPFI** (1.5.1+) - Multiprecision floating-point intervals (needs GMP, MPFR) +- **FLINT** (3.0.0+) - Fast library for number theory (needs GMP, MPFR) +- **NTL** (10.5.0+) - Number theory library (needs GMP) +- **FACTORY** (4.4.0+) - Polynomial factorization (needs GMP, FLINT, NTL) +- **BDWGC** (7.6.4+) - Boehm-Demers-Weiser garbage collector +- **LAPACK** - Linear algebra (BLAS/LAPACK) +- **MPSOLVE** (3.2.0+) - Multiprecision polynomial solver +- **FROBBY** (0.9.0+) - Computations with monomial ideals +- **NORMALIZ** (3.8.0+) - Discrete convex geometry (needs GMP, Nauty) +Note: READLINE, HISTORY, and GDBM are in the top-level `LIBRARY_LIST` but are only +used by the interpreter (`d/`), not the engine. + +### Optional +- **OpenMP** - Parallel processing + +## Code Conventions + +- C++17 standard (`-std=gnu++17`) +- Arithmetic ring types live in `M2` namespace (e.g., `M2::ARingZZ`, `M2::ARingZZp`) +- Ring elements use init/clear pattern: `R.init(a); ... R.clear(a);` +- `buffer` class (from `buffer.hpp`) used for string building and text output +- Interface functions prefixed with `IM2_` (e.g., `IM2_initialize()`) +- Source files generally come in `.cpp`/`.hpp` pairs +- Use `IM2_Ring_QQ()` for the rationals (QQGMP) in tests, not `rawARingQQFlint()` +- New classes should use `m` prefix for member fields (e.g., `mRing`, `mValue`) + +## Branch Info + +- Current branch: `unit-testing` (improving engine unit test infrastructure) +- PR target: `stable` diff --git a/M2/Macaulay2/d/actors4.d b/M2/Macaulay2/d/actors4.d index d4e663e3497..f94f98f3d57 100644 --- a/M2/Macaulay2/d/actors4.d +++ b/M2/Macaulay2/d/actors4.d @@ -9,10 +9,10 @@ use pthread; use regex; header "// required for toString routines -#include // for IM2_GB_to_string, rawMuta... // TODO: remove this one +#include // for rawGBToString, rawMuta... // TODO: remove this one #include // for rawHomotopyToString, rawP... -#include // for IM2_FreeModule_to_string -#include // for IM2_Matrix_to_string +#include // for rawFreeModuleToString +#include // for rawMatrixToString #include // for rawMonoidToString #include // for IM2_MonomialOrdering_to_s... #include // for IM2_MutableMatrix_to_string @@ -998,8 +998,8 @@ tostringfun(e:Expr):Expr := ( else toExpr("<>")) is s:SpecialExpr do tostringfun(s.e) is x:RawMonomialCell do toExpr(tostring(x.p)) - is x:RawFreeModuleCell do toExpr(Ccode(string, "IM2_FreeModule_to_string(",x.p,")" )) - is x:RawMatrixCell do toExpr(Ccode(string, "IM2_Matrix_to_string(",x.p,")" )) + is x:RawFreeModuleCell do toExpr(Ccode(string, "rawFreeModuleToString(",x.p,")" )) + is x:RawMatrixCell do toExpr(Ccode(string, "rawMatrixToString(",x.p,")" )) is x:RawMutableMatrixCell do toExpr(Ccode(string, "IM2_MutableMatrix_to_string(",x.p,")" )) is x:RawMutableComplexCell do toExpr(Ccode(string, "rawMutableComplexToString(",x.p,")" )) -- NAG stuff begin @@ -1016,7 +1016,7 @@ tostringfun(e:Expr):Expr := ( is x:RawRingCell do toExpr(Ccode(string, "IM2_Ring_to_string(",x.p,")" )) is x:RawRingElementCell do toExpr( Ccode(string, "IM2_RingElement_to_string(",x.p,")" ) ) is x:RawMonomialIdealCell do toExpr( Ccode(string, "IM2_MonomialIdeal_to_string(",x.p,")" ) ) - is c:RawComputationCell do toExpr(Ccode(string, "IM2_GB_to_string(",c.p,")" )) + is c:RawComputationCell do toExpr(Ccode(string, "rawGBToString(",c.p,")" )) is pythonObjectCell do toExpr("<>") is x:xmlNodeCell do toExpr(toString(x.v)) is xmlAttrCell do toExpr("<>") diff --git a/M2/Macaulay2/d/equality.dd b/M2/Macaulay2/d/equality.dd index bf13bccfa93..b808ee24cb6 100644 --- a/M2/Macaulay2/d/equality.dd +++ b/M2/Macaulay2/d/equality.dd @@ -3,8 +3,8 @@ use util; use tokens; header "// required for equality checks -#include // for IM2_FreeModule_is_equal -#include // for IM2_Matrix_is_equal +#include // for rawFreeModuleIsEqual +#include // for rawMatrixIsEqual #include // for IM2_MonomialIdeal_is_equal #include // for IM2_RingElement_is_equal #include // for IM2_RingMap_is_equal @@ -208,14 +208,14 @@ export equal(lhs:Expr,rhs:Expr):Expr := ( is x:RawFreeModuleCell do ( when rhs is y:RawFreeModuleCell do ( - if Ccode(bool, "IM2_FreeModule_is_equal(",x.p,",",y.p,")") + if Ccode(bool, "rawFreeModuleIsEqual(",x.p,",",y.p,")") then True else False ) else False ) is x:RawMatrixCell do ( when rhs - is y:RawMatrixCell do toExpr(Ccode(bool, "1 == IM2_Matrix_is_equal(",x.p,",",y.p,")")) + is y:RawMatrixCell do toExpr(Ccode(bool, "1 == rawMatrixIsEqual(",x.p,",",y.p,")")) else False ) is x:MysqlConnectionWrapper do ( diff --git a/M2/Macaulay2/d/interface.dd b/M2/Macaulay2/d/interface.dd index 356ac209310..e94b2aebab7 100644 --- a/M2/Macaulay2/d/interface.dd +++ b/M2/Macaulay2/d/interface.dd @@ -24,7 +24,7 @@ header "// TODO: break apart this file so each piece includes only one: #include "; header "// TODO: remove the following headers -#include // for IM2_Computation_set_stop +#include // for rawComputationSetStop #include // for ERROR #include // for Z_mod #include // for engine_error, CATCH, TRY @@ -988,8 +988,8 @@ 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, ")" ) ) + is x:RawFreeModuleCell do toExpr( Ccode(RawArrayIntOrNull, "rawFreeModuleGetDegrees(", x.p, ")" ) ) + is x:RawMatrixCell do toExpr( Ccode(RawArrayIntOrNull, "rawMatrixDegree(", x.p, ")" ) ) else WrongArg("a raw ring element, matrix, or free module")); setupfun("rawMultiDegree",rawMultiDegree); @@ -1028,21 +1028,21 @@ setupfun("rawTermCount",rawTermCount); export rawIsHomogeneous(e:Expr):Expr := ( when e is x:RawRingElementCell do toExpr( Ccode( bool, "IM2_RingElement_is_graded(", x.p, ")" )) - is x:RawMatrixCell do toExpr( Ccode( bool, "IM2_Matrix_is_graded(", x.p, ")" )) + is x:RawMatrixCell do toExpr( Ccode( bool, "rawMatrixIsHomogeneous(", x.p, ")" )) else WrongArg("a raw ring element") ); setupfun("rawIsHomogeneous",rawIsHomogeneous); export rawIsDense(e:Expr):Expr := ( when e is x:RawMatrixCell do - toExpr(Ccode( bool, "IM2_Matrix_is_implemented_as_dense(", x.p, ")" )) + toExpr(Ccode( bool, "rawMatrixIsDense(", x.p, ")" )) else WrongArgMatrix()); setupfun("rawIsDense",rawIsDense); export rawIsZero(e:Expr):Expr := ( when e is x:RawRingElementCell do toExpr( Ccode( bool, "IM2_RingElement_is_zero(", x.p, ")" )) - is x:RawMatrixCell do toExpr( Ccode( bool, "IM2_Matrix_is_zero(", x.p, ")" )) + is x:RawMatrixCell do toExpr( Ccode( bool, "rawMatrixIsZero(", x.p, ")" )) is x:RawMutableMatrixCell do toExpr( Ccode( bool, "IM2_MutableMatrix_is_zero(", x.p, ")" )) else WrongArg("a raw ring element or matrix or mutable matrix") ); @@ -1262,7 +1262,7 @@ export rawPromote(e:Expr):Expr := ( is M:RawFreeModuleCell do ( -- M is the new target free module when a.1 is f:RawMatrixCell do toExpr( Ccode(RawMatrixOrNull, - "IM2_Matrix_promote(", + "rawMatrixPromote(", M.p, ",", f.p, ")" )) else WrongArg("a raw matrix")) @@ -1294,7 +1294,7 @@ export rawLift(e:Expr):Expr := ( is M:RawFreeModuleCell do ( -- M is the new target free module when a.1 is f:RawMatrixCell do ( m := Ccode(RawMatrixOrNull, - "IM2_Matrix_lift(", + "rawMatrixLift(", "&", success, ",", M.p, ",", f.p, ")" ); @@ -1310,7 +1310,7 @@ export rawRing(e:Expr):Expr := ( is x:RawRingElementCell do toExpr( Ccode(RawRing, "IM2_RingElement_ring(", x.p, ")" )) is x:RawFreeModuleCell do toExpr( - Ccode(RawRing, "IM2_FreeModule_ring(", x.p, ")" )) + Ccode(RawRing, "rawFreeModuleRing(", x.p, ")" )) else WrongArg("a raw ring element or free module") ); setupfun("rawRing", rawRing); @@ -1335,7 +1335,7 @@ export rawHomogenize(e:Expr):Expr := ( if !isSmallInt(s.1) then WrongArgSmallInteger(2) else if !isSequenceOfSmallIntegers(s.2) then WrongArg(3,"a sequence of small integers") else toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_homogenize(", + "rawMatrixHomogenize(", M.p, ",", getSmallInt(s.1), ",", -- var number v getSequenceOfSmallIntegers(s.2), -- weights @@ -1719,7 +1719,7 @@ setupfun("rawTowerTranslatePoly",rawTowerTranslatePoly); export rawRank(e:Expr):Expr := ( when e - is x:RawFreeModuleCell do toExpr( Ccode( int, "IM2_FreeModule_rank(", x.p, ")" )) + is x:RawFreeModuleCell do toExpr( Ccode( int, "rawFreeModuleRank(", x.p, ")" )) else WrongArg("a raw free module") ); setupfun("rawRank",rawRank); @@ -1727,7 +1727,7 @@ setupfun("rawRank",rawRank); export rawSchreyerSource(e:Expr):Expr := ( when e is m:RawMatrixCell do toExpr( - Ccode(RawFreeModuleOrNull, "IM2_FreeModule_make_schreyer(", + Ccode(RawFreeModuleOrNull, "rawFreeModuleMakeSchreyer(", m.p, ")" ) ) else WrongArgMatrix()); setupfun("rawSchreyerSource",rawSchreyerSource); @@ -1742,13 +1742,13 @@ export rawFreeModule(e:Expr):Expr := ( is rank:ZZcell do ( if isInt(rank) then toExpr( Ccode( RawFreeModuleOrNull, - "IM2_FreeModule_make(", + "rawFreeModuleMake(", R.p, ",", toInt(rank), ")" )) else WrongArg(2,"a small integer or a sequence of small integers")) is degs:Sequence do ( if isSequenceOfSmallIntegers(degs) then toExpr( Ccode( RawFreeModuleOrNull, - "IM2_FreeModule_make_degs(", + "rawFreeModuleMakeDegs(", R.p, ",", getSequenceOfSmallIntegers(degs), ")" )) @@ -1762,7 +1762,7 @@ setupfun("rawFreeModule",rawFreeModule); export rawGetSchreyer(e:Expr):Expr := ( when e is F:RawFreeModuleCell do toExpr( - Ccode(RawMatrix, "IM2_FreeModule_get_schreyer(", F.p, ")" ) ) + Ccode(RawMatrix, "rawFreeModuleGetSchreyer(", F.p, ")" ) ) else WrongArg("a raw free module")); setupfun("rawGetSchreyer",rawGetSchreyer); @@ -1772,7 +1772,7 @@ export rawZero(e:Expr):Expr := ( when s.1 is G:RawFreeModuleCell do when s.2 is preference:ZZcell do if !isInt(preference) then WrongArgSmallInteger(3) else toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_zero(", + "rawMatrixZero(", F.p, ",", G.p, ",", toInt(preference), @@ -1789,7 +1789,7 @@ export rawExteriorPower(e:Expr):Expr := ( when s.0 is n:ZZcell do if !isInt(n) then WrongArgSmallInteger(1) else when s.1 is F:RawFreeModuleCell do toExpr(Ccode(RawFreeModule, - "IM2_FreeModule_exterior(", + "rawFreeModuleExterior(", toInt(n), ",", F.p, ")" )) @@ -1802,7 +1802,7 @@ export rawExteriorPower(e:Expr):Expr := ( when s.2 is strategy:ZZcell do if !isInt(strategy) then WrongArgSmallInteger(3) else toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_exterior(", + "rawMatrixExteriorPower(", toInt(n), ",", M.p, ",", toInt(strategy), @@ -1820,12 +1820,12 @@ export rawSymmetricPower(e:Expr):Expr := ( if !isInt(n) then WrongArgSmallInteger(1) else when s.1 is F:RawFreeModuleCell do toExpr(Ccode(RawFreeModule, - "IM2_FreeModule_symm(", + "rawFreeModuleSymm(", toInt(n), ",", F.p, ")" )) is M:RawMatrixCell do toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_symm(", + "rawMatrixSymmetricPower(", toInt(n), ",", M.p, ")" )) @@ -1836,22 +1836,22 @@ setupfun("rawSymmetricPower",rawSymmetricPower); export rawDual(e:Expr):Expr := ( when e - is F:RawFreeModuleCell do toExpr(Ccode(RawFreeModuleOrNull, "IM2_FreeModule_dual(", F.p, ")" )) - is M:RawMatrixCell do toExpr(Ccode(RawMatrixOrNull, "IM2_Matrix_transpose(", M.p, ")" )) + is F:RawFreeModuleCell do toExpr(Ccode(RawFreeModuleOrNull, "rawFreeModuleDual(", F.p, ")" )) + is M:RawMatrixCell do toExpr(Ccode(RawMatrixOrNull, "rawMatrixDual(", M.p, ")" )) is M:RawMutableMatrixCell do toExpr(Ccode(RawMutableMatrixOrNull, "rawMutableMatrixTranspose(", M.p, ")" )) else WrongArg("a raw free module, matrix, or mutable matrix")); setupfun("rawDual",rawDual); export rawDirectSum(e:Expr):Expr := ( if isSequenceOfMatrices(e) then toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_direct_sum(", + "rawMatrixDirectSum(", getSequenceOfMatrices(e), ")")) else when e is s:Sequence do if length(s) == 2 then when s.0 is F:RawFreeModuleCell do when s.1 is G:RawFreeModuleCell do toExpr(Ccode(RawFreeModuleOrNull, - "IM2_FreeModule_sum(", F.p, ",", G.p, ")" )) + "rawFreeModuleSum(", F.p, ",", G.p, ")" )) else WrongArg(2,"a raw free module") else WrongArg(1,"a raw matrix or free module") else WrongArg("a sequence of raw matrices or a pair of raw free modules") @@ -1866,7 +1866,7 @@ export rawSubmodule(e:Expr):Expr := ( if !isSequenceOfSmallIntegers(s.1) then WrongArg(2,"a sequence of small integers") else ( selection := getSequenceOfSmallIntegers(s.1); toExpr(Ccode(RawFreeModuleOrNull, - "IM2_FreeModule_submodule(", + "rawFreeModuleSubmodule(", M.p, ",", selection, ")" ) ) ) @@ -1901,7 +1901,7 @@ export rawIsEqual(e:Expr):Expr := ( when e is s:Sequence do if length(s) != 2 then WrongNumArgs(2) else when s.0 is x:RawMatrixCell do when s.1 is y:RawMatrixCell do ( - r := Ccode(int, "IM2_Matrix_is_equal(",x.p,",",y.p,")"); + r := Ccode(int, "rawMatrixIsEqual(",x.p,",",y.p,")"); if r == -1 then engineErrorMessage() else toExpr(r == 1)) else WrongArgMatrix(2) else @@ -1915,14 +1915,14 @@ setupfun("rawIsEqual",rawIsEqual); export rawSource(e:Expr):Expr := ( when e - is M:RawMatrixCell do toExpr( Ccode( RawFreeModule, "IM2_Matrix_get_source(", M.p, ")" )) + is M:RawMatrixCell do toExpr( Ccode( RawFreeModule, "rawMatrixSource(", M.p, ")" )) else WrongArgMatrix() ); setupfun("rawSource",rawSource); export rawTarget(e:Expr):Expr := ( when e - is M:RawMatrixCell do toExpr( Ccode( RawFreeModule, "IM2_Matrix_get_target(", M.p, ")" )) + is M:RawMatrixCell do toExpr( Ccode( RawFreeModule, "rawMatrixTarget(", M.p, ")" )) is F:RawRingMapCell do toExpr( Ccode( RawRing, "IM2_RingMap_target(", F.p, ")" )) else WrongArgMatrix() ); @@ -1936,7 +1936,7 @@ export rawMatrix1(e:Expr):Expr := ( if isSequenceOfRingElements(s.2) then when s.3 is preference:ZZcell do if !isInt(preference) then WrongArgSmallInteger(4) else toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_make1(", + "rawMatrix1(", target.p, ",", toInt(ncols), ",", getSequenceOfRingElements(s.2), ",", -- entries @@ -1958,7 +1958,7 @@ export rawMatrix2(e:Expr):Expr := ( if !isSequenceOfRingElements(s.3) then WrongArg(4,"a sequence of raw ring elements") else when s.4 is preference:ZZcell do if !isInt(preference) then WrongArgSmallInteger(5) else toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_make2(", + "rawMatrix2(", target.p, ",", source.p, ",", getSequenceOfSmallIntegers(s.2), ",", -- deg @@ -1977,7 +1977,7 @@ export rawMatrixRemake1(e:Expr):Expr := ( when s.0 is target:RawFreeModuleCell do when s.1 is M:RawMatrixCell do when s.2 is preference:ZZcell do if !isInt(preference) then WrongArgSmallInteger(3) else - toExpr(Ccode(RawMatrixOrNull, "IM2_Matrix_remake1(", + toExpr(Ccode(RawMatrixOrNull, "rawMatrixRemake1(", target.p, ",", M.p, ",", toInt(preference), @@ -1997,7 +1997,7 @@ export rawMatrixRemake2(e:Expr):Expr := ( when s.3 is M:RawMatrixCell do when s.4 is preference:ZZcell do if !isInt(preference) then WrongArgSmallInteger(5) else toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_remake2(", + "rawMatrixRemake2(", target.p, ",", source.p, ",", getSequenceOfSmallIntegers(s.2), ",", -- deg @@ -2021,7 +2021,7 @@ export rawSparseMatrix1(e:Expr):Expr := ( if isSequenceOfRingElements(s.4) then when s.5 is preference:ZZcell do if !isInt(preference) then WrongArgSmallInteger(6) else toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_make_sparse1(", + "rawSparseMatrix1(", target.p, ",", toInt(ncols), ",", getSequenceOfSmallIntegers(s.2), ",", -- rows @@ -2047,7 +2047,7 @@ export rawRandomConstantMatrix(e:Expr):Expr := ( when s.4 is specialType:ZZcell do if !isInt(specialType) then WrongArgSmallInteger(5) else when s.5 is preference:ZZcell do if !isInt(preference) then WrongArgSmallInteger(6) else toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_random(", + "rawRandomConstantMatrix(", R.p, ",", toInt(r), ",", toInt(c), ",", @@ -2075,7 +2075,7 @@ export rawSparseMatrix2(e:Expr):Expr := ( if isSequenceOfRingElements(s.5) then when s.6 is preference:ZZcell do if !isInt(preference) then WrongArgSmallInteger(6) else toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_make_sparse2(", + "rawSparseMatrix2(", target.p, ",", source.p, ",", getSequenceOfSmallIntegers(s.2), ",", -- deg @@ -2131,7 +2131,7 @@ setupfun("rawMatrixReadMsolveFile",rawMatrixReadMsolveFile); export rawConcat(e:Expr):Expr := ( if isSequenceOfMatrices(e) then toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_concat(", + "rawMatrixConcat(", getSequenceOfMatrices(e), ")")) else WrongArg("a raw matrix or a sequence of raw matrices") @@ -2155,7 +2155,7 @@ export rawMatrixEntry(e:Expr):Expr := ( if !isInt(r) then WrongArgSmallInteger(2) else when s.2 is c:ZZcell do if !isInt(c) then WrongArgSmallInteger(3) else ( - toExpr(Ccode(RawRingElementOrNull, "IM2_Matrix_get_entry(", M.p, ",", toInt(r), ",", toInt(c), ")" ) ) ) + toExpr(Ccode(RawRingElementOrNull, "rawMatrixEntry(", M.p, ",", toInt(r), ",", toInt(c), ")" ) ) ) else WrongArgZZ(3) else WrongArgZZ(2) else WrongArg(1,"a raw matrix or mutable matrix") @@ -2177,7 +2177,7 @@ export rawMatrixEntries(e:Expr):Expr := ( is M:RawMutableMatrixCell do ( entriesArrayToExpr(Ccode(RawRingElementArrayArrayOrNull, "IM2_MutableMatrix_get_entries(", M.p, ")"))) is M:RawMatrixCell do ( - entriesArrayToExpr(Ccode(RawRingElementArrayArrayOrNull, "IM2_Matrix_get_entries(", M.p, ")"))) + entriesArrayToExpr(Ccode(RawRingElementArrayArrayOrNull, "rawMatrixEntries(", M.p, ")"))) else WrongArg("a raw matrix or mutable matrix") ); setupfun("rawMatrixEntries",rawMatrixEntries); @@ -2191,7 +2191,7 @@ export rawSortColumns(e:Expr):Expr := ( when s.2 is mon_order:ZZcell do if !isInt(mon_order) then WrongArgSmallInteger(3) else ( toExprSeq(Ccode(RawArrayInt, - "IM2_Matrix_sort_columns(", + "rawMatrixSortColumns(", M.p, ",", toInt(deg_order), ",", toInt(mon_order), @@ -2212,7 +2212,7 @@ export rawEliminateVariables(e:Expr):Expr := ( when s.0 is nparts:ZZcell do if !isInt(nparts) then WrongArgSmallInteger(1) else when s.1 is M:RawMatrixCell do ( toExpr(Ccode(RawArrayInt, - "IM2_Matrix_elim_vars(", + "rawMatrixEliminateVariables(", toInt(nparts), ",", M.p, ")" ) ) ) @@ -2228,7 +2228,7 @@ export rawKeepVariables(e:Expr):Expr := ( when s.0 is nparts:ZZcell do if !isInt(nparts) then WrongArgSmallInteger(1) else when s.1 is M:RawMatrixCell do ( toExpr(Ccode(RawArrayInt, - "IM2_Matrix_keep_vars(", + "rawMatrixKeepVariables(", toInt(nparts), ",", M.p, ")" ) ) ) @@ -2247,7 +2247,7 @@ export rawDivideByVariable(e:Expr):Expr := ( when s.2 is maxdegree:ZZcell do if !isInt(maxdegree) then WrongArgSmallInteger(3) else ( toExpr(Ccode(RawMatrixAndInt, - "IM2_Matrix_divide_by_var(", + "rawMatrixDivideByVariable(", M.p, ",", toInt(var), ",", toInt(maxdegree), @@ -2291,7 +2291,7 @@ export rawInitial(e:Expr):Expr := ( when s.0 is p:ZZcell do if !isInt(p) then WrongArgSmallInteger(1) else when s.1 is M:RawMatrixCell do ( toExpr(Ccode(RawMatrix, - "IM2_Matrix_initial(", + "rawMatrixInitial(", toInt(p), ",", M.p, ")" @@ -2310,7 +2310,7 @@ export rawPfaffians(e:Expr):Expr := ( when s.0 is p:ZZcell do if !isInt(p) then WrongArgSmallInteger(1) else when s.1 is M:RawMatrixCell do ( toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_pfaffians(", toInt(p), ",", M.p, ")" ) ) ) + "rawMatrixPfaffians(", toInt(p), ",", M.p, ")" ) ) ) else WrongArgMatrix(2) else WrongArgZZ(1) else WrongNumArgs(2) @@ -2320,7 +2320,7 @@ setupfun("rawPfaffians",rawPfaffians); export rawPfaffian(e:Expr):Expr := ( when e is M:RawMatrixCell do toExpr(Ccode(RawRingElementOrNull, - "IM2_Matrix_pfaffian(", M.p, ")")) + "rawMatrixPfaffian(", M.p, ")")) else WrongArgMatrix()); setupfun("rawPfaffian", rawPfaffian); @@ -2330,7 +2330,7 @@ export rawTensor(e:Expr):Expr := ( is f:RawMatrixCell do when s.1 is g:RawMatrixCell do ( toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_tensor(", + "rawMatrixTensor(", f.p, ",", g.p, ")" @@ -2341,7 +2341,7 @@ export rawTensor(e:Expr):Expr := ( is M:RawFreeModuleCell do when s.1 is N:RawFreeModuleCell do ( toExpr(Ccode(RawFreeModuleOrNull, - "IM2_FreeModule_tensor(", + "rawFreeModuleTensor(", M.p, ",", N.p, ")" @@ -2395,7 +2395,7 @@ export rawMatrixDiff(e:Expr):Expr := ( when s.0 is f:RawMatrixCell do when s.1 is g:RawMatrixCell do ( toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_diff(", + "rawMatrixDiff(", f.p, ",", g.p, ")" @@ -2413,7 +2413,7 @@ export rawMatrixContract(e:Expr):Expr := ( when s.0 is f:RawMatrixCell do when s.1 is g:RawMatrixCell do ( toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_contract(", + "rawMatrixContract(", f.p, ",", g.p, ")" @@ -2431,7 +2431,7 @@ export rawIdentity(e:Expr):Expr := ( when s.0 is F:RawFreeModuleCell do when s.1 is preference:ZZcell do if !isInt(preference) then WrongArgSmallInteger(2) else toExpr(Ccode(RawMatrix, - "IM2_Matrix_identity(", + "rawMatrixIdentity(", F.p, ",", toInt(preference), ")" )) @@ -2514,7 +2514,7 @@ export rawMonomials(e:Expr):Expr := ( when s.1 is M:RawMatrixCell do ( vars := getSequenceOfSmallIntegers(s.0); toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_monomials(", + "rawMatrixMonomials(", vars, ",", M.p, ")" )) @@ -2549,7 +2549,7 @@ export rawSubmatrix(e:Expr):Expr := ( rows := getSequenceOfSmallIntegers(s.1); cols := getSequenceOfSmallIntegers(s.2); toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_submatrix(", + "rawMatrixSubmatrix(", M.p, ",", rows, ",", cols, @@ -2573,7 +2573,7 @@ export rawSubmatrix(e:Expr):Expr := ( if !isSequenceOfSmallIntegers(s.1) then WrongArg(2,"a sequence of small integers") else ( cols := getSequenceOfSmallIntegers(s.1); toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_submatrix1(", + "rawMatrixSubmatrixColumns(", M.p, ",", cols, ")" ) ) ) @@ -2599,7 +2599,7 @@ export rawReshape(e:Expr):Expr := ( when s.0 is M:RawMatrixCell do when s.1 is F:RawFreeModuleCell do when s.2 is G:RawFreeModuleCell do toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_reshape(", + "rawMatrixReshape(", M.p, ",", F.p, ",", G.p, ")" )) else WrongArg(3,"a raw free module") else WrongArg(2,"a raw free module") @@ -2613,7 +2613,7 @@ export rawFlip(e:Expr):Expr := ( if length(s) == 2 then when s.0 is F:RawFreeModuleCell do when s.1 is G:RawFreeModuleCell do toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_flip(", + "rawMatrixFlip(", F.p, ",", G.p, ")" )) else WrongArg(2,"a raw free module") else WrongArg(1,"a raw free module") @@ -2627,7 +2627,7 @@ export rawKoszul(e:Expr):Expr := ( if !isInt(n) then WrongArgSmallInteger(1) else when s.1 is F:RawMatrixCell do toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_koszul(", toInt(n), ",", F.p, ")" )) + "rawMatrixKoszul(", toInt(n), ",", F.p, ")" )) else WrongArgMatrix(2) else WrongArgZZ(1) else WrongNumArgs(2)); @@ -2651,7 +2651,7 @@ export rawHilbert(e:Expr):Expr := ( when e is M:RawMatrixCell do ( toExpr(Ccode(RawRingElementOrNull, - "IM2_Matrix_Hilbert(", + "rawMatrixHilbert(", M.p, ")" ) ) ) is M:RawMonomialIdealCell do ( toExpr(Ccode(RawRingElementOrNull, @@ -3276,14 +3276,14 @@ setupfun("rawRingMapEval",rawRingMapEval); export rawNumberOfRows(e:Expr):Expr := ( when e - is M:RawMatrixCell do toExpr(Ccode( int, "IM2_Matrix_n_rows(", M.p, ")" )) + is M:RawMatrixCell do toExpr(Ccode( int, "rawMatrixNumRows(", M.p, ")" )) is M:RawMutableMatrixCell do toExpr(Ccode( int, "IM2_MutableMatrix_n_rows(", M.p, ")" )) else WrongArg("a raw matrix")); setupfun("rawNumberOfRows",rawNumberOfRows); export rawNumberOfColumns(e:Expr):Expr := ( when e - is M:RawMatrixCell do toExpr(Ccode( int, "IM2_Matrix_n_cols(", M.p, ")" )) + is M:RawMatrixCell do toExpr(Ccode( int, "rawMatrixNumColumns(", M.p, ")" )) is M:RawMutableMatrixCell do toExpr(Ccode( int, "IM2_MutableMatrix_n_cols(", M.p, ")" )) else WrongArg("a raw matrix")); setupfun("rawNumberOfColumns",rawNumberOfColumns); @@ -3438,7 +3438,7 @@ export rawStatusResolution(e:Expr):Expr := ( when e is G:RawComputationCell do ( completionDegree := 0; completionLevel := 0; - ret := Ccode(int,"IM2_Resolution_status(", + ret := Ccode(int,"rawResolutionStatus(", G.p,",", "&",completionDegree,",", "&",completionLevel, @@ -3464,7 +3464,7 @@ export rawGB(e:Expr):Expr := ( when s.8 is maxReductionCount:ZZcell do if !isInt(maxReductionCount) then WrongArgSmallInteger(9) else toExpr( Ccode(RawComputationOrNull, - "IM2_GB_make(", + "rawGBMake(", m.p,",", isTrue(collectSyz),",", toInt(nRowsToKeep),",", @@ -3504,7 +3504,7 @@ export rawResolution(e:Expr):Expr := ( if !isInt(algorithm) then WrongArgSmallInteger(8) else toExpr( Ccode(RawComputationOrNull, - "IM2_res_make(", + "rawResMake(", m.p,",", isTrue(resolveCokernel),",", toInt(maxLevel),",", @@ -3534,7 +3534,7 @@ export rawGBSetHilbertFunction(e:Expr):Expr := ( when a.0 is G:RawComputationCell do when a.1 is h:RawRingElementCell do toExpr( Ccode(RawComputationOrNull, - "IM2_GB_set_hilbert_function(", + "rawGBSetHilbertFunction(", G.p, ",", h.p, ")" )) @@ -3552,7 +3552,7 @@ export rawGBForce(e:Expr):Expr := ( when a.2 is change:RawMatrixCell do when a.3 is syz:RawMatrixCell do toExpr( Ccode(RawComputationOrNull, - "IM2_GB_force(", + "rawGBForce(", m.p, ",", gb.p, ",", change.p, @@ -3631,7 +3631,7 @@ export rawGBSetStop(e:Expr):Expr := ( when s.8 is just_min_gens:Boolean do if !isSequenceOfSmallIntegers(s.9) then WrongArg(10,"a sequence of small integers") else -- length_limit toExpr( - Ccode(RawComputationOrNull,"IM2_Computation_set_stop(", + Ccode(RawComputationOrNull,"rawComputationSetStop(", G.p, ",", True == always_stop, ",", getSequenceOfSmallIntegers(s.2), ",", @@ -3775,7 +3775,7 @@ export rawResolutionStatusLevel(e:Expr):Expr := ( if !isInt(level) then WrongArgSmallInteger(2) else when s.2 is minimize:Boolean do ( completionDegree := 0; - ret := Ccode(int,"IM2_Resolution_status_level(", + ret := Ccode(int,"rawResolutionStatusLevel(", G.p,",", toInt(level), ",", True == minimize, ",", "&",completionDegree, ")" ); if ret == -1 then engineErrorMessage() else Expr(Sequence(toExpr(ret),toExpr(completionDegree)))) @@ -3850,7 +3850,7 @@ export rawGBMatrixLift(e:Expr):Expr := ( when a.1 is m:RawMatrixCell do ( resultRemainder := RawMatrixOrNull(NULL); resultQuotient := RawMatrixOrNull(NULL); - cplt := Ccode(bool, "IM2_GB_matrix_lift(", + cplt := Ccode(bool, "rawGBMatrixLift(", G.p, ",", m.p, ",", "&", resultRemainder, ",", @@ -3869,7 +3869,7 @@ export rawGBContains(e:Expr):Expr := ( when a.0 is G:RawComputationCell do when a.1 is m:RawMatrixCell do toExpr( Ccode(int, - "IM2_GB_contains(", + "rawGBContains(", G.p, ",", m.p, ")" )) diff --git a/M2/Macaulay2/docs/Doxyfile.in b/M2/Macaulay2/docs/Doxyfile.in index 232315a477f..600bd50ca6c 100644 --- a/M2/Macaulay2/docs/Doxyfile.in +++ b/M2/Macaulay2/docs/Doxyfile.in @@ -195,7 +195,7 @@ SHORT_NAMES = NO # description.) # The default value is: NO. -JAVADOC_AUTOBRIEF = NO +JAVADOC_AUTOBRIEF = YES # If the QT_AUTOBRIEF tag is set to YES then doxygen will interpret the first # line (until the first dot) of a Qt-style comment as the brief description. If diff --git a/M2/Macaulay2/e/BasicPoly.hpp b/M2/Macaulay2/e/BasicPoly.hpp index e2714d29a13..ce820bed4f4 100644 --- a/M2/Macaulay2/e/BasicPoly.hpp +++ b/M2/Macaulay2/e/BasicPoly.hpp @@ -1,6 +1,8 @@ -// Current restriction: the coefficients must be an integral type. -// TODO: allow infinite precision integers too. -// TODO: how should we handle coefficients which are: GF(p^n), QQ, fraction fields? or even polynomials? +/** + * Current restriction: the coefficients must be an integral type. +*/ +// TODO: allow infinite precision integers too. +// TODO: how should we handle coefficients which are: GF(p^n), QQ, fraction fields? or even polynomials? #pragma once #include "exceptions.hpp" @@ -14,10 +16,10 @@ class BasicPoly { public: std::vector mCoefficients; - std::vector mComponents; // if zero length: all components are 0. - std::vector mMonomials; // a concatenated list of varpower monomials. Each first entry is its length. + std::vector mComponents; ///< if zero length: all components are 0. + std::vector mMonomials; ///< a concatenated list of varpower monomials. Each first entry is its length. - void clear(); // resets all data to represent the zero polynomial + void clear(); ///< resets all data to represent the zero polynomial ~BasicPoly() { clear(); } size_t termCount() const { return mCoefficients.size(); } @@ -81,16 +83,23 @@ class IdentifierHash } }; -// These will throw a parsing_error if there is a parsing error. The -// plan is that that will include the location in the string of the -// error. +// TODO: include location of the parse_error in the string +/** + * Parse a polynomial with integer coefficients from a string with variables `varnames`. + * Variables a valid idenifiers made up of alphanumeric characters and underscore. + * Terms are of the form integer * variable ^ exponent + ... + * Exponents are assumed to be positive. + * \throws parsing_error +*/ BasicPoly parseBasicPoly(std::string poly, std::vector varnames); -/// This version is a potentially faster alternative when reading many polynomials +/** + * This version is a potentially faster alternative when reading many polynomials + * \throws parsing_error +*/ void parseBasicPoly(const std::string_view& str, const IdentifierHash& idenHash, BasicPoly& result); - // TODO: we want an iterator type here. // TODO: // diff --git a/M2/Macaulay2/e/BasicPolyList.hpp b/M2/Macaulay2/e/BasicPolyList.hpp index 71a6921686c..e515fa04062 100644 --- a/M2/Macaulay2/e/BasicPolyList.hpp +++ b/M2/Macaulay2/e/BasicPolyList.hpp @@ -1,8 +1,11 @@ -// BasicPolyList is a vector of polynomials (with components) -// which we can easily translate to and from other polynomial and matrix types. -// This class really doesn't require any ring. -// Current restriction: the coefficients must be an integral type. TODO: allow infinite precision integers too. -// TODO: how should we handle coefficients which are: GF(p^n), QQ, fraction fields? or even polynomials? +/** + * BasicPolyList is a vector of BasicPoly, polynomials (with components), + * which we can easily translate to and from other polynomial and matrix types. + * This class really doesn't require any ring. + * Current restriction: the coefficients must be an integral type. + * TODO: allow infinite precision integers too. +*/ +// TODO: how should we handle coefficients which are: GF(p^n), QQ, fraction fields? or even polynomials? #pragma once #include "exceptions.hpp" @@ -98,8 +101,13 @@ void toStream(const BasicPolyList& Fs, S &str) const Matrix* toMatrix(const FreeModule *target, const BasicPolyList& Fs); -// The following can certainly throw an error. You need to check that! +/** + * \throws parsing_error +*/ auto basicPolyListFromString(std::vector varNames, std::string polyPerLine) -> BasicPolyList; +/** + * \throws parsing_error +*/ auto basicPolyListFromFile(std::vector varNames, std::string fileName) -> BasicPolyList; // Local Variables: diff --git a/M2/Macaulay2/e/BasicPolyListParser.cpp b/M2/Macaulay2/e/BasicPolyListParser.cpp index 2d59ab51bf7..979d140423a 100644 --- a/M2/Macaulay2/e/BasicPolyListParser.cpp +++ b/M2/Macaulay2/e/BasicPolyListParser.cpp @@ -1,3 +1,7 @@ +/** + * Parsing utilities for polynomials and list of polynomials in both the engine + * and MSolve's format + */ #include "BasicPolyListParser.hpp" #include diff --git a/M2/Macaulay2/e/BasicPolyListParser.hpp b/M2/Macaulay2/e/BasicPolyListParser.hpp index 17a19d11712..55867b569f3 100644 --- a/M2/Macaulay2/e/BasicPolyListParser.hpp +++ b/M2/Macaulay2/e/BasicPolyListParser.hpp @@ -1,5 +1,7 @@ -// This class implements parsing of polynomials from a string or file -// as well as Msolve format. +/** + * This class implements parsing of polynomials from a string or file + * as well as Msolve format. +*/ #pragma once diff --git a/M2/Macaulay2/e/CMakeLists.txt b/M2/Macaulay2/e/CMakeLists.txt index f9f8907d91e..803db970d42 100644 --- a/M2/Macaulay2/e/CMakeLists.txt +++ b/M2/Macaulay2/e/CMakeLists.txt @@ -446,6 +446,8 @@ if(BUILD_TESTING) unit-tests/util-polyring-creation.hpp unit-tests/util-polyring-creation.cpp + unit-tests/RingElem.hpp + unit-tests/RingElem.cpp unit-tests/NewF4Test.cpp unit-tests/MonoidTest.cpp unit-tests/PolyRingTest.cpp @@ -462,7 +464,8 @@ if(BUILD_TESTING) unit-tests/ARingRRiTest.cpp unit-tests/ARingCCCTest.cpp unit-tests/NCGroebnerTest.cpp -# unit-tests/WeylAlgebraTest.cpp # TODO: add this file + unit-tests/WeylAlgebraTest.cpp + unit-tests/QuotientRingTest.cpp unit-tests/RingTest.hpp unit-tests/RingZZTest.cpp diff --git a/M2/Macaulay2/e/Eschreyer.hpp b/M2/Macaulay2/e/Eschreyer.hpp index fb515cff859..8738156367f 100644 --- a/M2/Macaulay2/e/Eschreyer.hpp +++ b/M2/Macaulay2/e/Eschreyer.hpp @@ -1,4 +1,7 @@ -// Copyright 1999 Michael E. Stillman +/** + * \author Micahel E. Stillman + * \copyright Copyright 1999 Michael E. Stillman +*/ #ifndef _Eschreyer_hpp_ #define _Eschreyer_hpp_ @@ -31,19 +34,18 @@ class GBKernelComputation : public Computation const Ring *K; GBRing *GR; const Monoid *M; - const SchreyerOrder *SF; // order for F. - const SchreyerOrder *SG; // order for G. - const FreeModule *F; // This is where the action is... - const FreeModule *G; // This is where the resulting syzygies live. + const SchreyerOrder *SF; ///< order for F. + const SchreyerOrder *SG; ///< order for G. + const FreeModule *F; ///< This is where the action is... + const FreeModule *G; ///< This is where the resulting syzygies live. // This MUST be a Schreyer free module compatible with the input! - gc_vector mi; // Used in reduction. - gc_vector gb; // This is the "stripped" GB. - gc_vector syzygies; // This is basically the result. + gc_vector mi; ///< Used in reduction. + gc_vector gb; ///< This is the "stripped" GB. + gc_vector syzygies; ///< This is basically the result. - // byte sizes for allocating temp exp vectors and monomials on the stack - size_t exp_size; - size_t monom_size; + size_t exp_size; ///< byte size for allocating temp exp vectors on the stack + size_t monom_size; ///< byte size for allocating monomials on the stack int n_ones; int n_unique; @@ -54,21 +56,35 @@ class GBKernelComputation : public Computation void strip_gb(const gc_vector &m); void strip_gb(const GBMatrix *m); + /** + * This routine grabs 'c', and 'monom' should be the total monomial. + */ gbvector *make_syz_term(ring_elem c, const_monomial monom, int comp) const; - // This routine grabs 'c', and 'monom' should be the total monomial. bool find_ring_divisor(const_exponents exp, const gbvector *&result); + + /** + * Returns the index of the least element in the monomial order which divides. + */ int find_divisor(const MonomialIdeal *mi, const_exponents exp, int &result); - // Returns the index of the least element in the monomial order which divides. + /** + * removes every term of f which is not a lead term of some element of gb. + */ void wipe_unneeded_terms(gbvector *&f); - // removes every term of f which is not a lead term of some element of gb. gbvector *s_pair(gbvector *syz); - void reduce(gbvector *&g, - gbvector *&gsyz); // Reduces g to zero. gsyz is real result. + + /** + * Reduces g to zero. gsyz is real result. + */ + void reduce(gbvector *&g, gbvector *&gsyz); + + /** + * Reduces g to zero. gsyz is real result. + */ void geo_reduce(gbvector *&g, - gbvector *&gsyz); // Reduces g to zero. gsyz is real result. + gbvector *&gsyz); public: GBKernelComputation(const GBMatrix *m); diff --git a/M2/Macaulay2/e/ExponentVector.hpp b/M2/Macaulay2/e/ExponentVector.hpp index bb02f1fbb38..d2150a15390 100644 --- a/M2/Macaulay2/e/ExponentVector.hpp +++ b/M2/Macaulay2/e/ExponentVector.hpp @@ -35,19 +35,19 @@ class ExponentVector typedef const Exponent *ConstExponents; typedef typename std::make_unsigned::type HashExponent; - // result = a + /// result = a static inline void copy(int nvars, ConstExponents a, Exponents result) { memcpy(result, a, nvars * sizeof(Exponent)); } - // result = (0, 0, ..., 0) + /// result = (0, 0, ..., 0) static inline void one(int nvars, Exponents result) { for (int i = 0; i < nvars; i++) *result++ = 0; } - // returns whether a_i == 0, all i + /// returns whether a_i == 0, all i static inline bool is_one(int nvars, ConstExponents a) { for (int i = 0; i < nvars; i++) @@ -55,13 +55,13 @@ class ExponentVector return true; } - // returns whether a_i == b_i, all i + /// returns whether a_i == b_i, all i static inline bool equal(int nvars, ConstExponents a, ConstExponents b) { return std::equal(a, a + nvars, b); } - // result = a + b + /// result = a + b static inline void mult(int nvars, ConstExponents a, ConstExponents b, @@ -73,7 +73,7 @@ class ExponentVector for (int i = nvars; i > 0; i--) *result++ = *a++ + *b++; } - // result = n * a + /// result = n * a static inline void power(int nvars, ConstExponents a, const Exponent n, @@ -85,7 +85,7 @@ class ExponentVector for (int i = nvars; i > 0; i--) *result++ = *a++ * n; } - // result = a + n * b + /// result = a + n * b static inline void multpower(int nvars, ConstExponents a, ConstExponents b, @@ -99,7 +99,7 @@ class ExponentVector for (int i = nvars; i > 0; i--) *result++ = *a++ + *b++ * n; } - // returns whether b_i >= a_i, all i + /// returns whether b_i >= a_i, all i static inline bool divides(int nvars, ConstExponents a, ConstExponents b) { // we go upward, because some rings have unused variables at the end @@ -108,7 +108,7 @@ class ExponentVector return true; } - // result = a - b + /// result = a - b static inline void divide(int nvars, ConstExponents a, ConstExponents b, @@ -120,7 +120,7 @@ class ExponentVector for (int i = 0; i < nvars; i++) *result++ = *a++ - *b++; } - // result = max(a - b, 0) + /// result = max(a - b, 0) static inline void quotient(int nvars, ConstExponents a, ConstExponents b, @@ -142,7 +142,7 @@ class ExponentVector } } - // result = min(a_i, b_i), all i + /// result = min(a_i, b_i), all i static inline void gcd(int nvars, ConstExponents a, ConstExponents b, @@ -156,7 +156,7 @@ class ExponentVector } } - // result = max(a_i, b_i), all i + /// result = max(a_i, b_i), all i static inline void lcm(int nvars, ConstExponents a, ConstExponents b, @@ -170,7 +170,7 @@ class ExponentVector } } - // returns GT, LT, or EQ + /// returns GT, LT, or EQ static inline int lex_compare(int nvars, ConstExponents a, ConstExponents b) { for (int i = 0; i < nvars; i++) @@ -181,7 +181,7 @@ class ExponentVector return EQ; } - // returns sum_i a_i + /// returns sum_i a_i static inline Exponent simple_degree(int nvars, ConstExponents a) { // TODO: use std::accumulate? @@ -195,7 +195,7 @@ class ExponentVector return sum; } - // returns sum_i a_i * wt_i + /// returns sum_i a_i * wt_i static inline Exponent weight(int nvars, ConstExponents a, const std::vector &wts) @@ -214,11 +214,11 @@ class ExponentVector return weight(nvars, a, M2_arrayint_to_stdvector(wts)); } - // i-th bit of output is 1 if 0 < a[i + k * size_of(int)] for some k + /// i-th bit of output is 1 if 0 < a[i + k * size_of(int)] for some k static HashExponent mask(int nvars, ConstExponents a); // FIXME: merge diverging specializations - // assigns c and d such that a*c = d*b = lcm(a,b) + /// assigns c and d such that a*c = d*b = lcm(a,b) static inline void syz(int nvars, ConstExponents a, ConstExponents b, diff --git a/M2/Macaulay2/e/NAG.hpp b/M2/Macaulay2/e/NAG.hpp index dec8341327f..d96cc2e235b 100644 --- a/M2/Macaulay2/e/NAG.hpp +++ b/M2/Macaulay2/e/NAG.hpp @@ -116,8 +116,8 @@ class PointArray private: std::map mMap; std::vector mPoints; - Weight mEpsilon; // tolerance - RealVector mWeights; // random positive numbers summing up to 1 + Weight mEpsilon; /**< \brief tolerance */ + RealVector mWeights; /**< \brief random positive numbers summing up to 1 */ decltype(mMap)::const_iterator left(Weight key) const { @@ -132,11 +132,20 @@ class PointArray // patching defs and functions: ///////////////////////////////////////// // switching from CCC to ConcreteRing ///////////////////////// #define CCC M2::ConcreteRing +/** + * \brief Cast a ring to a concrete, complex ring + * \param R the ring to cast + */ inline const CCC* cast_to_CCC(const Ring* R) { return dynamic_cast(R); } +/** + * \brief Construct a ring element of ARingCC from a real and imaginary part + * \param re the real part + * \param im the imaginary part + */ inline ring_elem from_doubles(const CCC* C, double re, double im) { M2::ARingCC::Element a(C->ring()); @@ -146,6 +155,11 @@ inline ring_elem from_doubles(const CCC* C, double re, double im) return result; } +/** + * \brief Convert a ring element of ARingCC to GMP's complex type + * \param C the concrete ring `a` belongs to + * \param a the element to convert + */ inline gmp_CC toBigComplex(const CCC* C, ring_elem a) { M2::ARingCC::Element b(C->ring()); @@ -155,12 +169,14 @@ inline gmp_CC toBigComplex(const CCC* C, ring_elem a) } /////////////////////////////////////////////////////////////////////////// -// Simple complex number class +/** + * \brief Simple complex number class + */ class complex { private: - double real; // Real Part - double imag; // Imaginary Part + double real; /**< \brief real part */ + double imag; /**< \brief imaginary part */ public: complex(); complex(double); @@ -185,19 +201,30 @@ class complex // CONSTRUCTOR inline complex::complex() {} +/** + * \param r real part + */ inline complex::complex(double r) { real = r; imag = 0; } +/** + * + * \param r real part + * \param im imaginary part + */ inline complex::complex(double r, double im) { real = r; imag = im; } -// COPY CONSTRUCTOR +/** + * \brief Copy constructor + * \param c the complex number to copy from + */ inline complex::complex(const complex& c) { this->real = c.real; @@ -422,17 +449,19 @@ class SLP : public MutableEngineObject static SLP* catalog[MAX_NUM_SLPs]; // get rid of... !!! static int num_slps; - bool is_relative_position; // can use relative or absolute addressing - M2_arrayint program; // std::vector??? - element_type* nodes; // array of CCs - gc_vector node_index; // points to position in program (rel. to start) - // of operation corresponding to a node + bool is_relative_position; /**< \brief can use relative or absolute addressing */ + M2_arrayint program; /**< \brief std::vector??? */ + element_type* nodes; /**< \brief array of CCs */ + /** + * \brief points to position in program (relative to start) of operation corresponding to a node + */ + gc_vector node_index; int num_consts, num_inputs, num_operations, rows_out, cols_out; - void* handle; // dynamic library handle + void* handle; /**< \brief dynamic library handle */ void (*compiled_fn)(element_type*, element_type*); - clock_t eval_time; // accumulates time spent in evaluation - int n_calls; // number of times called + clock_t eval_time; /**< \brief accumulates time spent in evaluation */ + int n_calls; /**< \brief number of times called */ SLP(); @@ -528,13 +557,13 @@ class StraightLineProgram : public SLP // enum SolutionStatus { ... defined in SLP-imp.hpp ... }; struct Solution { - int n; // number of coordinates - complex* x; // array of n coordinates - double t; // last value of parameter t used - complex* start_x; // start of the path that produced x - double cond; // reverse condition number of Hx - SolutionStatus status; - int num_steps; // number of steps taken along the path + int n; /**< \brief number of coordinates */ + complex* x; /**< \brief array of n coordinates */ + double t; /**< \brief last value of parameter t used */ + complex* start_x; /**< \brief start of the path that produced x */ + double cond; /**< \brief reverse condition number of Hx */ + SolutionStatus status; /**< \brief current status of the solution */ + int num_steps; /**< \brief number of steps taken along the path */ Solution() { status = UNDETERMINED; } void make(int m, const complex* s_s); @@ -551,18 +580,22 @@ class PathTracker : public MutableEngineObject static PathTracker* catalog[MAX_NUM_PATH_TRACKERS]; static int num_path_trackers; - int number; // trackers are enumerated + int number; /**< \brief trackers are enumerated */ Matrix* target; - const Matrix *H, *S, *T; // homotopy, start, target - StraightLineProgram *slpH, *slpHxt, *slpHxtH, - *slpHxH, // slps for evaluating H, H_{x,t}, H_{x,t}|H, H_{x}|H - *slpS, *slpSx, *slpSxS, *slpT, *slpTx, - *slpTxT; // slps for S and T, needed if is_projective - double productST, // real part of the Bombieri-Weyl (hermitian) product - bigT; // length of arc between S and T - double* DMforPN; // multipliers used in ProjectiveNewton - double maxDegreeTo3halves; // max(degree of equation)^{3/2} + const Matrix *H, *S, *T; /**< \brief homotopy, start, target */ + /** + * \brief slps for evaluating H, H_{x,t}, H_{x,t}|H, H_{x}|H + */ + StraightLineProgram *slpH, *slpHxt, *slpHxtH, *slpHxH; + /** + * \brief slps for S and T, needed if is_projective + */ + StraightLineProgram *slpS, *slpSx, *slpSxS, *slpT, *slpTx, *slpTxT; + double productST; /**< \brief real part of the Bombieri-Weyl (hermitian) product */ + double bigT; /**< \brief length of the arc between S and T */ + double* DMforPN; /**< \brief multipliers used in ProjectiveNewton */ + double maxDegreeTo3halves; /**< \brief max(degree of equation)^{3/2} */ // inline functions needed by track void evaluate_slpHxt(int n, const complex* x0t0, complex* Hxt) { @@ -635,13 +668,15 @@ class PathTracker : public MutableEngineObject slpHxH->evaluate(n + 1, x0t0, HxH); } - const CCC* C; // coefficient field (complex numbers) - const PolyRing* homotopy_R; // polynomial ring where homotopy lives (does not - // include t if is_projective) + const CCC* C; /**< \brief coefficient field (complex numbers) */ + /** + * \brief polynomial rign where homotopy lives (does no include t if is_projective) + */ + const PolyRing* homotopy_R; int n_coords; int n_sols; - Solution* raw_solutions; // solutions + stats - Matrix* solutions; // Matrix of solutions passed to top level + Solution* raw_solutions; /**< \brief solutions + stats */ + Matrix* solutions; /**< \brief Matrix of solutions passed to top level */ // parameters M2_bool is_projective; @@ -654,7 +689,7 @@ class PathTracker : public MutableEngineObject gmp_RR infinity_threshold; int pred_type; - void make_slps(); // creates slpHxt and alpHxH + void make_slps(); /**< \brief creates slpHxt and alpHxH */ PathTracker(); @@ -715,6 +750,11 @@ class PathTracker : public MutableEngineObject // ------------ service functions // -------------------------------------------------- int degree_ring_elem(const PolyRing* R, ring_elem re); +/** + * \brief Prints a square complex matrix to stdout in human-readable form. + * \param size number of rows/columns + * \param A pointer to the matrix data stored in row-major order + */ void print_complex_matrix(int size, const double* A); #endif diff --git a/M2/Macaulay2/e/complex.h b/M2/Macaulay2/e/complex.h index aae494c8595..2d849523d92 100644 --- a/M2/Macaulay2/e/complex.h +++ b/M2/Macaulay2/e/complex.h @@ -1,13 +1,20 @@ // Copyright 2008 Michael E. Stillman +/** + * \file complex.h + * \author Michael E. Stillman + * \date 2008 + * \brief declares the engine's gmp_CC complex number primitives + * + * Initialisation, arithmetic, and special functions on arbitrary-precision complex values built from two MPFR reals. + * The interface is similar to mpfr: + * - Every gmp_CC struct needs to be initialized with init or init_set. + * - All rounding is MPFR_RNDN. + * - Resulting values are the first argument +*/ #ifndef _complex_h_ #define _complex_h_ -/* The interface is similar to mpfr: - Every gmp_CC struct needs to be initialized with init or init_set. - All rounding is MPFR_RNDN. - Resulting values are the first argument -*/ #if !defined(SAFEC_EXPORTS) //#include @@ -31,8 +38,8 @@ extern "C" { void mpfc_mul(gmp_CCmutable result, gmp_CCmutable a, gmp_CCmutable b); void mpfc_invert(gmp_CCmutable result, gmp_CCmutable v); + /// result -= a * b void mpfc_sub_mult(gmp_CCmutable result, gmp_CCmutable a, gmp_CCmutable b); - /* result -= a*b */ void mpfc_div(gmp_CCmutable result, gmp_CCmutable a, gmp_CCmutable b); void mpfc_abs(gmp_RRmutable result, gmp_CCmutable a); diff --git a/M2/Macaulay2/e/engine.h b/M2/Macaulay2/e/engine.h index e6f157944dc..e77d21d8dcb 100644 --- a/M2/Macaulay2/e/engine.h +++ b/M2/Macaulay2/e/engine.h @@ -87,7 +87,7 @@ extern "C" { /**** Groebner basis and resolution routines ******/ /**************************************************/ - Computation /* or null */* IM2_Computation_set_stop(Computation *G, + Computation /* or null */* rawComputationSetStop(Computation *G, M2_bool always_stop, /* 1 */ M2_arrayint degree_limit, /* 2*/ int basis_element_limit, /* 3 */ @@ -110,7 +110,7 @@ extern "C" { /* The computation is complete up to and including this degree. The exact meaning of 'degree' is computation specific */ - M2_string IM2_GB_to_string(Computation *C); /* drg: connected, in actors4.d */ + M2_string rawGBToString(Computation *C); /* drg: connected, in actors4.d */ unsigned int rawComputationHash(const Computation *C); /* drg: connected, in basic.d */ @@ -123,7 +123,7 @@ extern "C" { /* LongPolynomial, Sort, Primary, Inhomogeneous, Homogeneous */ /* Res: SortStrategy, 0, 1, 2, 3 ?? */ - Computation /* or null */ *IM2_res_make(const Matrix *m, + Computation /* or null */ *rawResMake(const Matrix *m, M2_bool resolve_cokernel, int max_level, M2_bool use_max_slanted_degree, @@ -170,7 +170,7 @@ extern "C" { */ /* I don't know what this is supposed to do (mike) */ - int IM2_Resolution_status(Computation *G, + int rawResolutionStatus(Computation *G, int * complete_up_through_this_degree, int * complete_up_through_this_level); /* drg: TODO */ /* -1: error condition, and the error message is set. @@ -181,7 +181,7 @@ extern "C" { 4: finished the computation completely */ - enum ComputationStatusCode IM2_Resolution_status_level(Computation *G, + enum ComputationStatusCode rawResolutionStatusLevel(Computation *G, int level, M2_bool minimize, int * complete_up_through_this_degree); diff --git a/M2/Macaulay2/e/interface/factory.cpp b/M2/Macaulay2/e/interface/factory.cpp index 7a827a7ee0a..504274a50e3 100644 --- a/M2/Macaulay2/e/interface/factory.cpp +++ b/M2/Macaulay2/e/interface/factory.cpp @@ -862,7 +862,7 @@ engine_RawMatrixArrayOrNull rawCharSeries(const Matrix *M) if (error()) return nullptr; } result->array[next++] = - IM2_Matrix_make1(M->rows(), u.length(), result1, false); + rawMatrix1(M->rows(), u.length(), result1, false); } return result; diff --git a/M2/Macaulay2/e/interface/freemodule.cpp b/M2/Macaulay2/e/interface/freemodule.cpp index 25b0264b601..840f9cbe08a 100644 --- a/M2/Macaulay2/e/interface/freemodule.cpp +++ b/M2/Macaulay2/e/interface/freemodule.cpp @@ -12,9 +12,9 @@ class Matrix; -const Ring *IM2_FreeModule_ring(const FreeModule *F) { return F->get_ring(); } -int IM2_FreeModule_rank(const FreeModule *F) { return F->rank(); } -M2_string IM2_FreeModule_to_string(const FreeModule *F) +const Ring *rawFreeModuleRing(const FreeModule *F) { return F->get_ring(); } +int rawFreeModuleRank(const FreeModule *F) { return F->rank(); } +M2_string rawFreeModuleToString(const FreeModule *F) { buffer o; F->text_out(o); @@ -22,7 +22,7 @@ M2_string IM2_FreeModule_to_string(const FreeModule *F) } unsigned int rawFreeModuleHash(const FreeModule *F) { return F->hash(); } -const FreeModule /* or null */ *IM2_FreeModule_make(const Ring *R, int rank) +const FreeModule /* or null */ *rawFreeModuleMake(const Ring *R, int rank) { try { @@ -39,32 +39,12 @@ const FreeModule /* or null */ *IM2_FreeModule_make(const Ring *R, int rank) } } -const FreeModule /* or null */ *IM2_FreeModule_make_degs(const Ring *R, +const FreeModule /* or null */ *rawFreeModuleMakeDegs(const Ring *R, M2_arrayint degs) { try { - auto D = R->degree_monoid(); - unsigned int eachdeg = D->n_vars(); - if (eachdeg == 0) - { - ERROR("rawFreeModule: degree rank 0, but sequence of degrees given"); - return nullptr; - } - unsigned int rank = degs->len / eachdeg; - if (rank * eachdeg != degs->len) - { - ERROR("inappropriate number of degrees"); - return nullptr; - } - monomial deg = D->make_one(); - FreeModule *F = R->make_FreeModule(); - for (unsigned int i = 0; i < rank; i++) - { - D->from_expvector(degs->array + i * eachdeg, deg); - F->append(deg); - } - return F; + return R->make_FreeModule(degs->len, degs->array); } catch (const exc::engine_error& e) { ERROR(e.what()); @@ -72,7 +52,7 @@ const FreeModule /* or null */ *IM2_FreeModule_make_degs(const Ring *R, } } -const FreeModule /* or null */ *IM2_FreeModule_make_schreyer(const Matrix *m) +const FreeModule /* or null */ *rawFreeModuleMakeSchreyer(const Matrix *m) { try { @@ -84,7 +64,7 @@ const FreeModule /* or null */ *IM2_FreeModule_make_schreyer(const Matrix *m) } } -M2_arrayint IM2_FreeModule_get_degrees(const FreeModule *F) +M2_arrayint rawFreeModuleGetDegrees(const FreeModule *F) { auto D = F->get_ring()->degree_monoid(); auto n = D->n_vars(); @@ -95,12 +75,12 @@ M2_arrayint IM2_FreeModule_get_degrees(const FreeModule *F) return result; } -const Matrix *IM2_FreeModule_get_schreyer(const FreeModule *F) +const Matrix *rawFreeModuleGetSchreyer(const FreeModule *F) { return F->get_induced_order(); } -M2_bool IM2_FreeModule_is_equal(const FreeModule *F, const FreeModule *G) +M2_bool rawFreeModuleIsEqual(const FreeModule *F, const FreeModule *G) /* Determines if F and G are the same graded module. If one has a Schreyer order and one does not, but their ranks and degrees are the same, then they are considered equal by this routine. */ @@ -108,13 +88,13 @@ M2_bool IM2_FreeModule_is_equal(const FreeModule *F, const FreeModule *G) return F->is_equal(G); } -const FreeModule /* or null */ *IM2_FreeModule_sum(const FreeModule *F, +const FreeModule /* or null */ *rawFreeModuleSum(const FreeModule *F, const FreeModule *G) { return F->direct_sum(G); } -const FreeModule /* or null */ *IM2_FreeModule_tensor(const FreeModule *F, +const FreeModule /* or null */ *rawFreeModuleTensor(const FreeModule *F, const FreeModule *G) { try @@ -127,7 +107,7 @@ const FreeModule /* or null */ *IM2_FreeModule_tensor(const FreeModule *F, } } -const FreeModule /* or null */ *IM2_FreeModule_dual(const FreeModule *F) +const FreeModule /* or null */ *rawFreeModuleDual(const FreeModule *F) { try { @@ -139,7 +119,7 @@ const FreeModule /* or null */ *IM2_FreeModule_dual(const FreeModule *F) } } -const FreeModule *IM2_FreeModule_symm(int n, const FreeModule *F) +const FreeModule *rawFreeModuleSymm(int n, const FreeModule *F) { try { @@ -151,7 +131,7 @@ const FreeModule *IM2_FreeModule_symm(int n, const FreeModule *F) } } -const FreeModule *IM2_FreeModule_exterior(int n, const FreeModule *F) +const FreeModule *rawFreeModuleExterior(int n, const FreeModule *F) { try { @@ -163,7 +143,7 @@ const FreeModule *IM2_FreeModule_exterior(int n, const FreeModule *F) } } -const FreeModule *IM2_FreeModule_submodule(const FreeModule *F, +const FreeModule *rawFreeModuleSubmodule(const FreeModule *F, M2_arrayint selection) { try diff --git a/M2/Macaulay2/e/interface/freemodule.h b/M2/Macaulay2/e/interface/freemodule.h index 27d23bb093f..73fc9a80956 100644 --- a/M2/Macaulay2/e/interface/freemodule.h +++ b/M2/Macaulay2/e/interface/freemodule.h @@ -38,22 +38,22 @@ typedef struct Ring Ring; extern "C" { # endif -const Ring *IM2_FreeModule_ring(const FreeModule *F); +const Ring *rawFreeModuleRing(const FreeModule *F); /* drg: connected rawRing*/ -int IM2_FreeModule_rank(const FreeModule *F); +int rawFreeModuleRank(const FreeModule *F); /* drg: connected rawRank*/ -M2_string IM2_FreeModule_to_string(const FreeModule *F); +M2_string rawFreeModuleToString(const FreeModule *F); /* drg: connected */ unsigned int rawFreeModuleHash(const FreeModule *F); /* not quite connected */ -const FreeModule /* or null */ *IM2_FreeModule_make(const Ring *R, int rank); +const FreeModule /* or null */ *rawFreeModuleMake(const Ring *R, int rank); /* drg: connected rawFreeModule*/ -const FreeModule /* or null */ *IM2_FreeModule_make_degs(const Ring *R, +const FreeModule /* or null */ *rawFreeModuleMakeDegs(const Ring *R, M2_arrayint degs); /* drg: connected rawFreeModule*/ /* Make a graded free module over R. 'degs' should be of length @@ -62,7 +62,7 @@ const FreeModule /* or null */ *IM2_FreeModule_make_degs(const Ring *R, * i = 0. */ -const FreeModule /* or null */ *IM2_FreeModule_make_schreyer(const Matrix *m); +const FreeModule /* or null */ *rawFreeModuleMakeSchreyer(const Matrix *m); /* drg: connected rawSchreyerSource */ /* Returns G, (a copy of) the source free module of 'm', modified to * use the induced order via m: compare two monomials of G via @@ -73,27 +73,27 @@ const FreeModule /* or null */ *IM2_FreeModule_make_schreyer(const Matrix *m); * handled efficiently. */ -M2_arrayint IM2_FreeModule_get_degrees(const FreeModule *F); +M2_arrayint rawFreeModuleGetDegrees(const FreeModule *F); /* drg: connected rawMultiDegree*/ -const Matrix *IM2_FreeModule_get_schreyer(const FreeModule *F); +const Matrix *rawFreeModuleGetSchreyer(const FreeModule *F); /* drg: connected rawGetSchreyer*/ -M2_bool IM2_FreeModule_is_equal(const FreeModule *F, const FreeModule *G); +M2_bool rawFreeModuleIsEqual(const FreeModule *F, const FreeModule *G); /* drg: connected === */ /* Determines if F and G are the same graded module. If one has a * Schreyer order and one does not, but their ranks and degrees are the * same, then they are considered equal by this routine. */ -const FreeModule /* or null */ *IM2_FreeModule_sum(const FreeModule *F, +const FreeModule /* or null */ *rawFreeModuleSum(const FreeModule *F, const FreeModule *G); /* drg: connected rawDirectSum */ /* The direct sum of two free modules over the same ring, or NULL. * If F or G has a Schreyer order, then so does their direct sum */ -const FreeModule /* or null */ *IM2_FreeModule_tensor(const FreeModule *F, +const FreeModule /* or null */ *rawFreeModuleTensor(const FreeModule *F, const FreeModule *G); /* drg: connected rawTensor*/ /* The tensor product of two free modules over the same ring, or NULL. @@ -108,7 +108,7 @@ const FreeModule /* or null */ *IM2_FreeModule_tensor(const FreeModule *F, * At the moment, the answer is almost yes... */ -const FreeModule /* or null */ *IM2_FreeModule_dual(const FreeModule *F); +const FreeModule /* or null */ *rawFreeModuleDual(const FreeModule *F); /* drg: connected rawDual*/ /* Returns the graded dual F^* of F: if F has basis {f_1,...,f_r}, * with degrees {d_1, ..., d_r}, then F^* has rank r, with @@ -116,7 +116,7 @@ const FreeModule /* or null */ *IM2_FreeModule_dual(const FreeModule *F); * Schreyer order (even if F does). */ -const FreeModule *IM2_FreeModule_symm(int n, const FreeModule *F); +const FreeModule *rawFreeModuleSymm(int n, const FreeModule *F); /* drg: connected rawSymmetricPower*/ /* Returns the n th symmetric power G of F. * If F has basis {f_1,...,f_r}, then G has basis @@ -125,7 +125,7 @@ const FreeModule *IM2_FreeModule_symm(int n, const FreeModule *F); * If F has a Schreyer order, then G is set to have one as well. */ -const FreeModule *IM2_FreeModule_exterior(int n, const FreeModule *F); +const FreeModule *rawFreeModuleExterior(int n, const FreeModule *F); /* drg: connected rawExteriorPower*/ /* Returns the n th exterior power G of F. * If F has basis {f_1,...,f_r}, then G has basis @@ -134,7 +134,7 @@ const FreeModule *IM2_FreeModule_exterior(int n, const FreeModule *F); * If F has a Schreyer order, then G is set to have one as well. */ -const FreeModule /* or null */ *IM2_FreeModule_submodule(const FreeModule *F, +const FreeModule /* or null */ *rawFreeModuleSubmodule(const FreeModule *F, M2_arrayint selection); /* drg: connected rawSubmodule*/ /* Returns a free module obtained by choosing basis elements of F: diff --git a/M2/Macaulay2/e/interface/groebner.cpp b/M2/Macaulay2/e/interface/groebner.cpp index 1ce4d27f610..dbf061d8109 100644 --- a/M2/Macaulay2/e/interface/groebner.cpp +++ b/M2/Macaulay2/e/interface/groebner.cpp @@ -55,7 +55,7 @@ void test_over_RR_or_CC(const Ring *R) } //////////////////////////////////// -const RingElement /* or null */ *IM2_Matrix_Hilbert(const Matrix *M) +const RingElement /* or null */ *rawMatrixHilbert(const Matrix *M) /* This routine computes the numerator of the Hilbert series for coker leadterms(M), using the degrees of the rows of M. nullptr is returned if the ring is not appropriate for @@ -87,7 +87,7 @@ const Matrix *rawKernelOfGB(const Matrix *M) ///////// The following will be removed once the new code is functional ///////////// /////////////////////////////////////////////////////////////////////////////////// -Computation /* or null */ *IM2_GB_make( +Computation /* or null */ *rawGBMake( const Matrix *m, M2_bool collect_syz, int n_rows_to_keep, @@ -121,7 +121,7 @@ Computation /* or null */ *IM2_GB_make( } } -Computation /* or null */ *IM2_res_make(const Matrix *m, +Computation /* or null */ *rawResMake(const Matrix *m, M2_bool resolve_cokernel, int max_level, M2_bool use_max_slanted_degree, @@ -158,7 +158,7 @@ Computation /* or null */ *IM2_res_make(const Matrix *m, } } -Computation /* or null */ *IM2_GB_set_hilbert_function(Computation *C, +Computation /* or null */ *rawGBSetHilbertFunction(Computation *C, const RingElement *h) { try @@ -180,7 +180,7 @@ Computation /* or null */ *IM2_GB_set_hilbert_function(Computation *C, } } -Computation /* or null */ *IM2_GB_force( +Computation /* or null */ *rawGBForce( const Matrix *m, /* trimmed or minimal gens, may be the same as gb */ const Matrix *gb, const Matrix *change, /* same number of columns as 'gb', if not 0 */ @@ -228,7 +228,7 @@ Computation /* or null */ *rawGroebnerWalk(const Matrix *gb, } } -Computation /* or null */ *IM2_Computation_set_stop( +Computation /* or null */ *rawComputationSetStop( Computation *G, M2_bool always_stop, M2_arrayint degree_limit, @@ -447,7 +447,7 @@ const Matrix /* or null */ *rawGBMatrixRemainder(Computation *C, } } -M2_bool IM2_GB_matrix_lift(Computation *C, +M2_bool rawGBMatrixLift(Computation *C, const Matrix *m, const Matrix /* or null */ **result_remainder, const Matrix /* or null */ **result_quotient) @@ -467,7 +467,7 @@ M2_bool IM2_GB_matrix_lift(Computation *C, return false; } -int IM2_GB_contains(Computation *C, const Matrix *m) +int rawGBContains(Computation *C, const Matrix *m) { try { @@ -576,7 +576,7 @@ const FreeModule /* or null */ *rawResolutionGetFree(Computation *C, int level) } } -int IM2_Resolution_status(Computation *C, +int rawResolutionStatus(Computation *C, int *complete_up_through_this_degree, int *complete_up_through_this_level) { @@ -584,13 +584,13 @@ int IM2_Resolution_status(Computation *C, (void) complete_up_through_this_degree; (void) complete_up_through_this_level; #ifdef DEVELOPMENT -#warning "IM2_Resolution_status to be written" +#warning "rawResolutionStatus to be written" #endif ERROR("not re-implemented yet"); return -1; } -enum ComputationStatusCode IM2_Resolution_status_level( +enum ComputationStatusCode rawResolutionStatusLevel( Computation *C, int level, M2_bool minimize, @@ -601,7 +601,7 @@ enum ComputationStatusCode IM2_Resolution_status_level( (void) minimize; (void) complete_up_through_this_degree; #ifdef DEVELOPMENT -#warning "IM2_Resolution_status to be written" +#warning "rawResolutionStatus to be written" #endif ERROR("not re-implemented yet"); return COMP_ERROR; @@ -630,7 +630,7 @@ M2_arrayintOrNull rawResolutionBetti(Computation *C, int type) } } -M2_string IM2_GB_to_string(Computation *C) +M2_string rawGBToString(Computation *C) /* TODO */ { buffer o; diff --git a/M2/Macaulay2/e/interface/groebner.h b/M2/Macaulay2/e/interface/groebner.h index cb5d9e4f1e4..e4649a466f0 100644 --- a/M2/Macaulay2/e/interface/groebner.h +++ b/M2/Macaulay2/e/interface/groebner.h @@ -50,7 +50,7 @@ M2_arrayint rawMinimalBetti(Computation *G, M2_arrayint length_limit); /* connected: rawMinimalBetti */ -const RingElement /* or null */ *IM2_Matrix_Hilbert(const Matrix *M); +const RingElement /* or null */ *rawMatrixHilbert(const Matrix *M); /* This routine computes the numerator of the Hilbert series for coker leadterms(M), using the degrees of the rows of M. NULL is returned if the ring is not appropriate for @@ -64,7 +64,7 @@ const Matrix *rawKernelOfGB(const Matrix *M); /////////////////////////////////////////////////////////////////////////////// /////// The following will be removed once the new code is functional /////// /////////////////////////////////////////////////////////////////////////////// -Computation /* or null */ *IM2_GB_make( +Computation /* or null */ *rawGBMake( const Matrix *m, M2_bool collect_syz, int n_rows_to_keep, @@ -75,7 +75,7 @@ Computation /* or null */ *IM2_GB_make( int strategy, int max_reduction_count); /* drg: connected rawGB */ -Computation /* or null */ *IM2_res_make(const Matrix *m, +Computation /* or null */ *rawResMake(const Matrix *m, M2_bool resolve_cokernel, int max_level, M2_bool use_max_slanted_degree, @@ -85,11 +85,11 @@ Computation /* or null */ *IM2_res_make(const Matrix *m, M2_bool parallelizeByDegree); /* rawGBSetHilbertFunction */ -Computation /* or null */ *IM2_GB_set_hilbert_function(Computation *C, +Computation /* or null */ *rawGBSetHilbertFunction(Computation *C, const RingElement *h); /* rawGBForce */ -Computation /* or null */ *IM2_GB_force( +Computation /* or null */ *rawGBForce( const Matrix *m, /* trimmed or minimal gens, may be the same as gb */ const Matrix *gb, const Matrix *change, /* same number of columns as 'gb', if not 0 */ @@ -114,7 +114,7 @@ Computation /* or null */ *rawMarkedGB( Computation /* or null */ *rawGroebnerWalk(const Matrix *gb, const MonomialOrdering *order1); -Computation /* or null */ *IM2_Computation_set_stop( +Computation /* or null */ *rawComputationSetStop( Computation *G, M2_bool always_stop, M2_arrayint degree_limit, @@ -171,7 +171,7 @@ const Matrix /* or null */ *rawGBMatrixRemainder(Computation *C, /* rawGBMatrixLift: false is returned if there is an error or if the remainder is NON-zero */ -M2_bool IM2_GB_matrix_lift(Computation *C, +M2_bool rawGBMatrixLift(Computation *C, const Matrix *m, const Matrix /* or null */ **result_remainder, const Matrix /* or null */ **result_quotient); @@ -180,7 +180,7 @@ M2_bool IM2_GB_matrix_lift(Computation *C, -2 for errors -1 for containment of span(m) in span(C) i for the smallest index of a column of m not contained in span(C) */ -int IM2_GB_contains(Computation *C, const Matrix *m); +int rawGBContains(Computation *C, const Matrix *m); /******************************************* * Noncommutative Groebner bases *********** @@ -241,12 +241,12 @@ MutableMatrix /* or null */ *rawResolutionGetMutableMatrix2B( const FreeModule /* or null */ *rawResolutionGetFree(Computation *C, int level); /* TODO */ -int IM2_Resolution_status(Computation *C, +int rawResolutionStatus(Computation *C, int *complete_up_through_this_degree, int *complete_up_through_this_level); /* TODO */ -enum ComputationStatusCode IM2_Resolution_status_level( +enum ComputationStatusCode rawResolutionStatusLevel( Computation *C, int level, M2_bool minimize, @@ -255,7 +255,7 @@ enum ComputationStatusCode IM2_Resolution_status_level( M2_arrayintOrNull rawResolutionBetti(Computation *C, int type); /* see engine.h for description of what 'type' should be */ -M2_string IM2_GB_to_string(Computation *C); +M2_string rawGBToString(Computation *C); /* TODO */ unsigned int rawComputationHash(const Computation *C); diff --git a/M2/Macaulay2/e/interface/matrix.cpp b/M2/Macaulay2/e/interface/matrix.cpp index 4617aae7cce..5450717e91e 100644 --- a/M2/Macaulay2/e/interface/matrix.cpp +++ b/M2/Macaulay2/e/interface/matrix.cpp @@ -15,7 +15,6 @@ #include "interface/gmp-util.h" #include "interface/monoid.h" #include "mat.hpp" -#include "matrix-con.hpp" #include "matrix.hpp" #include "mutablemat-defs.hpp" #include "relem.hpp" @@ -27,17 +26,17 @@ namespace M2 { class ARingCC; } -const FreeModule *IM2_Matrix_get_target(const Matrix *M) { return M->rows(); } -const FreeModule *IM2_Matrix_get_source(const Matrix *M) { return M->cols(); } -int IM2_Matrix_n_rows(const Matrix *M) { return M->n_rows(); } -int IM2_Matrix_n_cols(const Matrix *M) { return M->n_cols(); } +const FreeModule *rawMatrixTarget(const Matrix *M) { return M->rows(); } +const FreeModule *rawMatrixSource(const Matrix *M) { return M->cols(); } +int rawMatrixNumRows(const Matrix *M) { return M->n_rows(); } +int rawMatrixNumColumns(const Matrix *M) { return M->n_cols(); } -M2_arrayint IM2_Matrix_get_degree(const Matrix *M) +M2_arrayint rawMatrixDegree(const Matrix *M) { return to_degree_vector(M->get_ring()->degree_monoid(), M->degree_shift()); } -M2_string IM2_Matrix_to_string(const Matrix *M) +M2_string rawMatrixToString(const Matrix *M) { buffer o; try @@ -52,27 +51,13 @@ M2_string IM2_Matrix_to_string(const Matrix *M) } unsigned int rawMatrixHash(const Matrix *M) { return M->hash(); } -const RingElement /* or null */ *IM2_Matrix_get_entry(const Matrix *M, +const RingElement /* or null */ *rawMatrixEntry(const Matrix *M, int r, int c) { try { - if (r < 0 || r >= M->n_rows()) - { - ERROR("matrix row index %d out of range 0 .. %d", r, M->n_rows() - 1); - return nullptr; - } - if (c < 0 || c >= M->n_cols()) - { - ERROR("matrix column index %d out of range 0 .. %d", - c, - M->n_cols() - 1); - return nullptr; - } - ring_elem result; - result = M->elem(r, c); - return RingElement::make_raw(M->get_ring(), result); + return M->entry(r, c); } catch (const exc::engine_error& e) { ERROR(e.what()); @@ -80,54 +65,13 @@ const RingElement /* or null */ *IM2_Matrix_get_entry(const Matrix *M, } } -/* Returns the entries of the matrix in a flat array in row major order +/* Returns the entries of the matrix as an array of rows. */ -engine_RawRingElementArrayArrayOrNull IM2_Matrix_get_entries(const Matrix *M) +engine_RawRingElementArrayArrayOrNull rawMatrixEntries(const Matrix *M) { try { - int ncols = M->n_cols(); - int nrows = M->n_rows(); - if(nrows < 0 || ncols < 0) - { - ERROR("internal error: matrix has a negative size %d by %d", - nrows, - ncols); - return nullptr; - } - engine_RawRingElementArrayArray entries = - getmemarraytype(engine_RawRingElementArrayArray, nrows); - entries->len = nrows; - RingElement *zero = - RingElement::make_raw(M->get_ring(), M->get_ring()->zero()); - for(int r = 0; r < nrows; r++) - { - engine_RawRingElementArray currRow = - getmemarraytype(engine_RawRingElementArray, ncols); - currRow->len = ncols; - std::fill_n(currRow->array, ncols, zero); - entries->array[r] = currRow; - } - //walk through the columns - for(int c = 0; c < ncols; c++) - { - const vec &column = M->elem(c); - for(const vecterm &term : column) - { - if(term.comp < 0 || term.comp >= nrows) - { - ERROR("internal error: matrix contains invalid entries:" - "row index %d out of range 0 .. %d", - term.comp, - nrows - 1); - //Ignoring the entry and continuing - continue; - } - entries->array[term.comp]->array[c] = - RingElement::make_raw(M->get_ring(), term.coeff); - } - } - return entries; + return M->entries(); } catch (const exc::engine_error &e) { ERROR(e.what()); @@ -136,7 +80,7 @@ engine_RawRingElementArrayArrayOrNull IM2_Matrix_get_entries(const Matrix *M) return nullptr; } -const Matrix *IM2_Matrix_identity(const FreeModule *F, int preference) +const Matrix *rawMatrixIdentity(const FreeModule *F, int preference) { (void) preference; #ifdef DEVELOPMENT @@ -145,7 +89,7 @@ const Matrix *IM2_Matrix_identity(const FreeModule *F, int preference) return Matrix::identity(F); } -const Matrix /* or null */ *IM2_Matrix_zero(const FreeModule *F, +const Matrix /* or null */ *rawMatrixZero(const FreeModule *F, const FreeModule *G, int preference) { @@ -156,7 +100,7 @@ const Matrix /* or null */ *IM2_Matrix_zero(const FreeModule *F, return Matrix::zero(F, G); } -const Matrix /* or null */ *IM2_Matrix_make1(const FreeModule *target, +const Matrix /* or null */ *rawMatrix1(const FreeModule *target, int ncols, const engine_RawRingElementArray M, int preference) @@ -168,7 +112,7 @@ const Matrix /* or null */ *IM2_Matrix_make1(const FreeModule *target, return Matrix::make(target, ncols, M); } -const Matrix /* or null */ *IM2_Matrix_make2(const FreeModule *target, +const Matrix /* or null */ *rawMatrix2(const FreeModule *target, const FreeModule *source, M2_arrayint deg, const engine_RawRingElementArray M, @@ -181,7 +125,7 @@ const Matrix /* or null */ *IM2_Matrix_make2(const FreeModule *target, return Matrix::make(target, source, deg, M); } -const Matrix /* or null */ *IM2_Matrix_make_sparse1( +const Matrix /* or null */ *rawSparseMatrix1( const FreeModule *target, int ncols, M2_arrayint rows, @@ -196,7 +140,7 @@ const Matrix /* or null */ *IM2_Matrix_make_sparse1( return Matrix::make_sparse(target, ncols, rows, cols, entries); } -const Matrix /* or null */ *IM2_Matrix_make_sparse2( +const Matrix /* or null */ *rawSparseMatrix2( const FreeModule *target, const FreeModule *source, M2_arrayint deg, @@ -212,7 +156,7 @@ const Matrix /* or null */ *IM2_Matrix_make_sparse2( return Matrix::make_sparse(target, source, deg, rows, cols, entries); } -M2_bool IM2_Matrix_is_implemented_as_dense(const Matrix *M) +M2_bool rawMatrixIsDense(const Matrix *M) /* Is the matrix M implemented as dense? */ { (void) M; @@ -222,7 +166,7 @@ M2_bool IM2_Matrix_is_implemented_as_dense(const Matrix *M) return 0; } -const Matrix /* or null */ *IM2_Matrix_remake2(const FreeModule *target, +const Matrix /* or null */ *rawMatrixRemake2(const FreeModule *target, const FreeModule *source, M2_arrayint deg, const Matrix *M, @@ -239,7 +183,7 @@ const Matrix /* or null */ *IM2_Matrix_remake2(const FreeModule *target, return M->remake(target, source, deg); } -const Matrix /* or null */ *IM2_Matrix_remake1(const FreeModule *target, +const Matrix /* or null */ *rawMatrixRemake1(const FreeModule *target, const Matrix *M, int preference) /* Create a new matrix, from M, with new target, @@ -262,7 +206,7 @@ const Matrix /* or null */ *IM2_Matrix_remake1(const FreeModule *target, } } -const Matrix /* or null */ *IM2_Matrix_random( +const Matrix /* or null */ *rawRandomConstantMatrix( const Ring *R, int r, int c, @@ -306,9 +250,9 @@ const Matrix* /* or null */ rawMatrixReadMsolveFile(const Ring* R, M2_string fil } ///////////////////////////////////////////////////////////////////// -M2_bool IM2_Matrix_is_zero(const Matrix *M) { return M->is_zero(); } +M2_bool rawMatrixIsZero(const Matrix *M) { return M->is_zero(); } int // 1 = true, 0 = false, -1 = error - IM2_Matrix_is_equal(const Matrix *M, const Matrix *N) + rawMatrixIsEqual(const Matrix *M, const Matrix *N) { try { @@ -324,41 +268,12 @@ int // 1 = true, 0 = false, -1 = error } } -M2_bool IM2_Matrix_is_graded(const Matrix *M) { return M->is_homogeneous(); } -const Matrix /* or null */ *IM2_Matrix_concat(const engine_RawMatrixArray Ms) +M2_bool rawMatrixIsHomogeneous(const Matrix *M) { return M->is_homogeneous(); } +const Matrix /* or null */ *rawMatrixConcat(const engine_RawMatrixArray Ms) { try { - unsigned int n = Ms->len; - if (n == 0) - { - ERROR("matrix concat: expects at least one matrix"); - return nullptr; - } - const FreeModule *F = Ms->array[0]->rows(); - const Ring *R = F->get_ring(); - MatrixConstructor mat(Ms->array[0]->rows(), 0); - int next = 0; - for (unsigned int i = 0; i < n; i++) - { - const Matrix *M = Ms->array[i]; - if (R != M->get_ring()) - { - ERROR("matrix concat: different base rings"); - return nullptr; - } - if (F->rank() != M->n_rows()) - { - ERROR("matrix concat: row sizes are not equal"); - return nullptr; - } - for (int j = 0; j < M->n_cols(); j++) - { - mat.append(R->copy_vec(M->elem(j))); - mat.set_column_degree(next++, M->cols()->degree(j)); - } - } - return mat.to_matrix(); + return Matrix::concat(Ms->len, Ms->array); } catch (const exc::engine_error& e) { ERROR(e.what()); @@ -366,31 +281,11 @@ const Matrix /* or null */ *IM2_Matrix_concat(const engine_RawMatrixArray Ms) } } -const Matrix /* or null */ *IM2_Matrix_direct_sum( - const engine_RawMatrixArray Ms) +const Matrix /* or null */ *rawMatrixDirectSum(const engine_RawMatrixArray Ms) { try { - // Check that the matrices all have the same ring, and that there is - // at least one matrix. - unsigned int n = Ms->len; - if (n == 0) - { - ERROR("matrix direct sum: expects at least one matrix"); - return nullptr; - } - const Matrix *result = Ms->array[0]; - const Ring *R = result->get_ring(); - for (unsigned int i = 1; i < n; i++) - if (R != Ms->array[i]->get_ring()) - { - ERROR("matrix direct sum: different base rings"); - return nullptr; - } - for (unsigned int i = 1; i < n; i++) - result = result->direct_sum(Ms->array[i]); - - return result; + return Matrix::direct_sum(Ms->len, Ms->array); } catch (const exc::engine_error& e) { ERROR(e.what()); @@ -398,7 +293,7 @@ const Matrix /* or null */ *IM2_Matrix_direct_sum( } } -const Matrix /* or null */ *IM2_Matrix_tensor(const Matrix *M, const Matrix *N) +const Matrix /* or null */ *rawMatrixTensor(const Matrix *M, const Matrix *N) { try { @@ -422,7 +317,7 @@ const Matrix /* or null */ *rawModuleTensor(const Matrix *M, const Matrix *N) } } -const Matrix /* or null */ *IM2_Matrix_transpose(const Matrix *M) +const Matrix /* or null */ *rawMatrixDual(const Matrix *M) { try { @@ -434,7 +329,7 @@ const Matrix /* or null */ *IM2_Matrix_transpose(const Matrix *M) } } -const Matrix /* or null */ *IM2_Matrix_reshape(const Matrix *M, +const Matrix /* or null */ *rawMatrixReshape(const Matrix *M, const FreeModule *F, const FreeModule *G) { @@ -448,7 +343,7 @@ const Matrix /* or null */ *IM2_Matrix_reshape(const Matrix *M, } } -const Matrix /* or null */ *IM2_Matrix_flip(const FreeModule *F, +const Matrix /* or null */ *rawMatrixFlip(const FreeModule *F, const FreeModule *G) { try @@ -476,7 +371,7 @@ const Matrix /* or null */ *rawWedgeProduct(int p, int q, const FreeModule *F) } } -const Matrix /* or null */ *IM2_Matrix_submatrix(const Matrix *M, +const Matrix /* or null */ *rawMatrixSubmatrix(const Matrix *M, M2_arrayint rows, M2_arrayint cols) { @@ -490,7 +385,7 @@ const Matrix /* or null */ *IM2_Matrix_submatrix(const Matrix *M, } } -const Matrix /* or null */ *IM2_Matrix_submatrix1(const Matrix *M, +const Matrix /* or null */ *rawMatrixSubmatrixColumns(const Matrix *M, M2_arrayint cols) { try @@ -503,7 +398,7 @@ const Matrix /* or null */ *IM2_Matrix_submatrix1(const Matrix *M, } } -const Matrix /* or null */ *IM2_Matrix_koszul(int p, const Matrix *M) +const Matrix /* or null */ *rawMatrixKoszul(int p, const Matrix *M) { try { @@ -537,7 +432,7 @@ const Matrix /* or null */ *rawKoszulMonomials(int nskew, } } -const Matrix /* or null */ *IM2_Matrix_symm(int p, const Matrix *M) +const Matrix /* or null */ *rawMatrixSymmetricPower(int p, const Matrix *M) { try { @@ -549,7 +444,7 @@ const Matrix /* or null */ *IM2_Matrix_symm(int p, const Matrix *M) } } -const Matrix /* or null */ *IM2_Matrix_exterior(int p, +const Matrix /* or null */ *rawMatrixExteriorPower(int p, const Matrix *M, int strategy) { @@ -564,7 +459,7 @@ const Matrix /* or null */ *IM2_Matrix_exterior(int p, } -M2_arrayintOrNull IM2_Matrix_sort_columns(const Matrix *M, +M2_arrayintOrNull rawMatrixSortColumns(const Matrix *M, int deg_order, int mon_order) { @@ -578,7 +473,7 @@ M2_arrayintOrNull IM2_Matrix_sort_columns(const Matrix *M, } } -const Matrix /* or null */ *IM2_Matrix_minors(int p, +const Matrix /* or null */ *rawMatrixMinors(int p, const Matrix *M, int strategy) { @@ -617,28 +512,28 @@ const Matrix /* or null */ *rawMinors( } } -const Matrix /* or null */ *IM2_Matrix_pfaffians(int p, const Matrix *M) +const Matrix /* or null */ *rawMatrixPfaffians(int p, const Matrix *M) { return M->pfaffians(p); } -const RingElement /* or null */ *IM2_Matrix_pfaffian(const Matrix *M) +const RingElement /* or null */ *rawMatrixPfaffian(const Matrix *M) { return RingElement::make_raw(M->get_ring(), M->pfaffian()); } -const Matrix /* or null */ *IM2_Matrix_diff(const Matrix *M, const Matrix *N) +const Matrix /* or null */ *rawMatrixDiff(const Matrix *M, const Matrix *N) { return M->diff(N, 1); } -const Matrix /* or null */ *IM2_Matrix_contract(const Matrix *M, +const Matrix /* or null */ *rawMatrixContract(const Matrix *M, const Matrix *N) { return M->diff(N, 0); } -const Matrix /* or null */ *IM2_Matrix_homogenize(const Matrix *M, +const Matrix /* or null */ *rawMatrixHomogenize(const Matrix *M, int var, M2_arrayint wts) { @@ -667,23 +562,23 @@ const Matrix /* or null */ *rawBasis( M2_arrayintOrNull rawMatrixIndices(const Matrix *f) /* The list of indices of variables which occur in f is returned. */ /* currently requires a polynomial ring */ { return f->support(); } -const Matrix /* or null */ *IM2_Matrix_monomials(M2_arrayint vars, +const Matrix /* or null */ *rawMatrixMonomials(M2_arrayint vars, const Matrix *M) { return M->monomials(vars); } -const Matrix *IM2_Matrix_initial(int nparts, const Matrix *M) +const Matrix *rawMatrixInitial(int nparts, const Matrix *M) { return M->lead_term(nparts); } -M2_arrayint IM2_Matrix_elim_vars(int nparts, const Matrix *M) +M2_arrayint rawMatrixEliminateVariables(int nparts, const Matrix *M) { return M->elim_vars(nparts); } -M2_arrayint IM2_Matrix_keep_vars(int nparts, const Matrix *M) +M2_arrayint rawMatrixKeepVariables(int nparts, const Matrix *M) { return M->elim_keep(nparts); } @@ -700,7 +595,7 @@ engine_RawMatrixPairOrNull rawTopCoefficients(const Matrix *M) return result; } -engine_RawMatrixAndInt IM2_Matrix_divide_by_var(const Matrix *M, +engine_RawMatrixAndInt rawMatrixDivideByVariable(const Matrix *M, int var, int maxdegree) /* If M = [v1, ..., vn], and x = 'var'th variable in the ring, @@ -761,27 +656,12 @@ const Matrix /* or null */ *IM2_Matrix_remove_content(const Matrix *M) return nullptr; } -const Matrix /* or null */ *IM2_Matrix_promote(const FreeModule *newTarget, +const Matrix /* or null */ *rawMatrixPromote(const FreeModule *newTarget, const Matrix *f) { try { - ring_elem a; - const Ring *R = f->get_ring(); - const Ring *S = newTarget->get_ring(); - MatrixConstructor mat(newTarget, f->n_cols()); - Matrix::iterator i(f); - for (int c = 0; c < f->n_cols(); c++) - for (i.set(c); i.valid(); i.next()) - if (S->promote(R, i.entry(), a)) - mat.set_entry(i.row(), c, a); - else - { - ERROR("cannot promote given matrix"); - return nullptr; - } - mat.compute_column_degrees(); - return mat.to_matrix(); + return f->promote(newTarget); } catch (const exc::engine_error& e) { ERROR(e.what()); @@ -789,29 +669,17 @@ const Matrix /* or null */ *IM2_Matrix_promote(const FreeModule *newTarget, } } -const Matrix /* or null */ *IM2_Matrix_lift(int *success_return, +const Matrix /* or null */ *rawMatrixLift(int *success_return, const FreeModule *newTarget, const Matrix *f) { + *success_return = 0; try { - ring_elem a; - const Ring *R = f->get_ring(); - const Ring *S = newTarget->get_ring(); - MatrixConstructor mat(newTarget, f->n_cols()); - Matrix::iterator i(f); - for (int c = 0; c < f->n_cols(); c++) - for (i.set(c); i.valid(); i.next()) - if (R->lift(S, i.entry(), a)) - mat.set_entry(i.row(), c, a); - else - { - // ERROR("cannot lift given matrix"); - return nullptr; - } - mat.compute_column_degrees(); + const Matrix *result = f->lift(newTarget); + if (result == nullptr) return nullptr; *success_return = 1; - return mat.to_matrix(); + return result; } catch (const exc::engine_error& e) { ERROR(e.what()); diff --git a/M2/Macaulay2/e/interface/matrix.h b/M2/Macaulay2/e/interface/matrix.h index 61efe9d27de..29be2eadbe1 100644 --- a/M2/Macaulay2/e/interface/matrix.h +++ b/M2/Macaulay2/e/interface/matrix.h @@ -24,11 +24,11 @@ typedef struct RingElement RingElement; extern "C" { # endif -const Matrix /* or null */ *IM2_Matrix_promote(const FreeModule *newTarget, +const Matrix /* or null */ *rawMatrixPromote(const FreeModule *newTarget, const Matrix *f); /* connected to rawPromote*/ -const Matrix /* or null */ *IM2_Matrix_lift(int *success_return, +const Matrix /* or null */ *rawMatrixLift(int *success_return, const FreeModule *newTarget, const Matrix *f); /* connected to rawLift */ @@ -38,54 +38,54 @@ const Matrix /* or null */ *IM2_Matrix_lift(int *success_return, /**** Matrix routines *****************************/ /**************************************************/ -const FreeModule *IM2_Matrix_get_target( +const FreeModule *rawMatrixTarget( const Matrix *M); /* drg: connected rawTarget*/ -const FreeModule *IM2_Matrix_get_source( +const FreeModule *rawMatrixSource( const Matrix *M); /* drg: connected rawSource, used in rawMatrixColumns*/ -int IM2_Matrix_n_rows(const Matrix *M); /* drg: connected rawNumberOfRows*/ +int rawMatrixNumRows(const Matrix *M); /* drg: connected rawNumberOfRows*/ -int IM2_Matrix_n_cols(const Matrix *M); /* drg: connected rawNumberOfColumns*/ +int rawMatrixNumColumns(const Matrix *M); /* drg: connected rawNumberOfColumns*/ -M2_arrayint IM2_Matrix_get_degree( +M2_arrayint rawMatrixDegree( const Matrix *M); /* drg: connected rawMultiDegree*/ -M2_string IM2_Matrix_to_string(const Matrix *M); /* drg: connected */ +M2_string rawMatrixToString(const Matrix *M); /* drg: connected */ unsigned int rawMatrixHash(const Matrix *M); /* drg: connected to "hash" */ -const RingElement /* or null */ *IM2_Matrix_get_entry( +const RingElement /* or null */ *rawMatrixEntry( const Matrix *M, int r, int c); /* drg: connected rawMatrixEntry, OK*/ -engine_RawRingElementArrayArrayOrNull IM2_Matrix_get_entries(const Matrix *M); +engine_RawRingElementArrayArrayOrNull rawMatrixEntries(const Matrix *M); /*******************************************************************************/ -const Matrix *IM2_Matrix_identity( +const Matrix *rawMatrixIdentity( const FreeModule *F, int preference); /* drg: connected rawIdentity, OK*/ -const Matrix /* or null */ *IM2_Matrix_zero( +const Matrix /* or null */ *rawMatrixZero( const FreeModule *F, const FreeModule *G, int preference); /* drg: connected rawZero, OK */ -const Matrix /* or null */ *IM2_Matrix_make1( +const Matrix /* or null */ *rawMatrix1( const FreeModule *target, int ncols, const engine_RawRingElementArray M, int preference); /* drg: connected rawMatrix1, OK */ -const Matrix /* or null */ *IM2_Matrix_make2( +const Matrix /* or null */ *rawMatrix2( const FreeModule *target, const FreeModule *source, M2_arrayint deg, const engine_RawRingElementArray M, int preference); /* drg: connected rawMatrix2, OK */ -const Matrix /* or null */ *IM2_Matrix_make_sparse1( +const Matrix /* or null */ *rawSparseMatrix1( const FreeModule *target, int ncols, M2_arrayint rows, @@ -93,7 +93,7 @@ const Matrix /* or null */ *IM2_Matrix_make_sparse1( const engine_RawRingElementArray entries, int preference); /* drg: connected rawSparseMatrix1, OK */ -const Matrix /* or null */ *IM2_Matrix_make_sparse2( +const Matrix /* or null */ *rawSparseMatrix2( const FreeModule *target, const FreeModule *source, M2_arrayint deg, @@ -102,11 +102,11 @@ const Matrix /* or null */ *IM2_Matrix_make_sparse2( const engine_RawRingElementArray entries, int preference); /* drg: connected rawSparseMatrix2, OK */ -M2_bool IM2_Matrix_is_implemented_as_dense( +M2_bool rawMatrixIsDense( const Matrix *M); /* connected to rawIsDense */ /* Is the matrix M implemented in the engine as a dense matrix? */ -const Matrix /* or null */ *IM2_Matrix_remake1( +const Matrix /* or null */ *rawMatrixRemake1( const FreeModule *target, const Matrix *M, int preference); /* drg: connected rawMatrixRemake1, OK */ @@ -116,7 +116,7 @@ const Matrix /* or null */ *IM2_Matrix_remake1( columns of the matrix. */ -const Matrix /* or null */ *IM2_Matrix_remake2( +const Matrix /* or null */ *rawMatrixRemake2( const FreeModule *target, const FreeModule *source, M2_arrayint deg, @@ -127,7 +127,7 @@ const Matrix /* or null */ *IM2_Matrix_remake2( the expected rank. */ -const Matrix /* or null */ *IM2_Matrix_random( +const Matrix /* or null */ *rawRandomConstantMatrix( const Ring *R, int r, int c, @@ -141,9 +141,9 @@ const Matrix* /* or null */ rawMatrixReadMsolveFile(const Ring* R, M2_string fil /**********************************************************************************/ -M2_bool IM2_Matrix_is_zero(const Matrix *M); /* drg: connected rawIsZero*/ +M2_bool rawMatrixIsZero(const Matrix *M); /* drg: connected rawIsZero*/ -int IM2_Matrix_is_equal( +int rawMatrixIsEqual( const Matrix *M, const Matrix *N); /* drg: connected === and to rawIsEqual for use with == */ // 1 = true, 0 = false, -1 = error @@ -152,28 +152,28 @@ int IM2_Matrix_is_equal( Therefore, it can happen that M-N == 0, but M != N. */ -M2_bool IM2_Matrix_is_graded( +M2_bool rawMatrixIsHomogeneous( const Matrix *M); /* drg: connected rawIsHomogeneous*/ -const Matrix /* or null */ *IM2_Matrix_concat( +const Matrix /* or null */ *rawMatrixConcat( const engine_RawMatrixArray Ms); /* drg: connected rawConcat*/ -const Matrix /* or null */ *IM2_Matrix_direct_sum( +const Matrix /* or null */ *rawMatrixDirectSum( const engine_RawMatrixArray Ms); /* drg: connected rawDirectSum*/ -const Matrix /* or null */ *IM2_Matrix_tensor( +const Matrix /* or null */ *rawMatrixTensor( const Matrix *M, const Matrix *N); /* drg: connected rawTensor*/ -const Matrix /* or null */ *IM2_Matrix_transpose( +const Matrix /* or null */ *rawMatrixDual( const Matrix *M); /* drg: connected rawDual*/ -const Matrix /* or null */ *IM2_Matrix_reshape( +const Matrix /* or null */ *rawMatrixReshape( const Matrix *M, const FreeModule *F, const FreeModule *G); /* drg: connected rawReshape*/ -const Matrix /* or null */ *IM2_Matrix_flip( +const Matrix /* or null */ *rawMatrixFlip( const FreeModule *F, const FreeModule *G); /* drg: connected rawFlip*/ @@ -185,16 +185,16 @@ const Matrix /* or null */ *rawWedgeProduct( exterior(p,F) ** exterior(q,F) --> exterior(p+q,F) */ -const Matrix /* or null */ *IM2_Matrix_submatrix( +const Matrix /* or null */ *rawMatrixSubmatrix( const Matrix *M, M2_arrayint rows, M2_arrayint cols); /* drg: connected rawSubmatrix*/ -const Matrix /* or null */ *IM2_Matrix_submatrix1( +const Matrix /* or null */ *rawMatrixSubmatrixColumns( const Matrix *M, M2_arrayint cols); /* drg: connected rawSubmatrix*/ -const Matrix /* or null */ *IM2_Matrix_koszul( +const Matrix /* or null */ *rawMatrixKoszul( int p, const Matrix *M); /* drg: connected rawKoszul*/ @@ -209,22 +209,22 @@ const Matrix /* or null */ *rawKoszulMonomials( exterior algebra (on this set of variables). The actual commutativity of the common ring of M and N is ignored. */ -const Matrix /* or null */ *IM2_Matrix_symm( +const Matrix /* or null */ *rawMatrixSymmetricPower( int p, const Matrix *M); /* drg: connected rawSymmetricPower*/ -const Matrix /* or null */ *IM2_Matrix_exterior( +const Matrix /* or null */ *rawMatrixExteriorPower( int p, const Matrix *M, int strategy); /* drg: connected rawExteriorPower*/ -M2_arrayint IM2_Matrix_sort_columns( +M2_arrayint rawMatrixSortColumns( const Matrix *M, int deg_order, int mon_order); /* drg: connected rawSortColumns*/ const Matrix /* or null */ * -IM2_Matrix_minors(int p, const Matrix *M, int strategy); /* drg: unconnected*/ +rawMatrixMinors(int p, const Matrix *M, int strategy); /* drg: unconnected*/ const Matrix /* or null */ *rawMinors( int p, @@ -239,11 +239,11 @@ const Matrix /* or null */ *rawMinors( if given, otherwise starting at the first (0..p-1,0..p-1). */ -const Matrix /* or null */ *IM2_Matrix_pfaffians( +const Matrix /* or null */ *rawMatrixPfaffians( int p, const Matrix *M); /* drg: connected rawPfaffians*/ -const RingElement /* or null */ *IM2_Matrix_pfaffian( +const RingElement /* or null */ *rawMatrixPfaffian( const Matrix *M); const Matrix *rawMatrixCompress( @@ -297,15 +297,15 @@ const Matrix /* or null */ *IM2_Matrix_remove_content( /* Routines for use when the base ring is a polynomial ring of some sort */ -const Matrix /* or null */ *IM2_Matrix_diff( +const Matrix /* or null */ *rawMatrixDiff( const Matrix *M, const Matrix *N); /* drg: connected rawMatrixDiff*/ -const Matrix /* or null */ *IM2_Matrix_contract( +const Matrix /* or null */ *rawMatrixContract( const Matrix *M, const Matrix *N); /* drg: connected rawMatrixContract*/ -const Matrix /* or null */ *IM2_Matrix_homogenize( +const Matrix /* or null */ *rawMatrixHomogenize( const Matrix *M, int var, M2_arrayint wts); /* drg: connected rawHomogenize*/ @@ -334,23 +334,23 @@ const Matrix /* or null */ *rawCoefficients( * will be used (which one is left undefined) */ -const Matrix /* or null */ *IM2_Matrix_monomials( +const Matrix /* or null */ *rawMatrixMonomials( M2_arrayint vars, const Matrix *M); /* drg: connected rawMonomials*/ -const Matrix *IM2_Matrix_initial( +const Matrix *rawMatrixInitial( int nparts, const Matrix *M); /* drg: connected rawInitial*/ -M2_arrayint IM2_Matrix_elim_vars( +M2_arrayint rawMatrixEliminateVariables( int nparts, const Matrix *M); /* drg: connected rawEliminateVariables*/ -M2_arrayint IM2_Matrix_keep_vars( +M2_arrayint rawMatrixKeepVariables( int nparts, const Matrix *M); /* drg: connected rawKeepVariables*/ -engine_RawMatrixAndInt IM2_Matrix_divide_by_var( +engine_RawMatrixAndInt rawMatrixDivideByVariable( const Matrix *M, int var, int maxdegree); /* drg: connected rawDivideByVariable*/ diff --git a/M2/Macaulay2/e/matrix.cpp b/M2/Macaulay2/e/matrix.cpp index 2832d6ec33a..396b75d2c3b 100644 --- a/M2/Macaulay2/e/matrix.cpp +++ b/M2/Macaulay2/e/matrix.cpp @@ -13,6 +13,7 @@ #include "style.hpp" #include "text-io.hpp" #include "ring.hpp" +#include "ringelem.hpp" #include "comb.hpp" #include "polyring.hpp" #include "assprime.hpp" @@ -57,6 +58,70 @@ unsigned int Matrix::computeHashValue() const return hashval; } +const RingElement /* or null */ *Matrix::entry(int r, int c) const +{ + if (r < 0 || r >= n_rows()) + { + ERROR("matrix row index %d out of range 0 .. %d", r, n_rows() - 1); + return nullptr; + } + if (c < 0 || c >= n_cols()) + { + ERROR("matrix column index %d out of range 0 .. %d", + c, + n_cols() - 1); + return nullptr; + } + return RingElement::make_raw(get_ring(), elem(r, c)); +} + +engine_RawRingElementArrayArrayOrNull Matrix::entries() const +{ + int ncols = n_cols(); + int nrows = n_rows(); + if (nrows < 0 || ncols < 0) + { + ERROR("internal error: matrix has a negative size %d by %d", + nrows, + ncols); + return nullptr; + } + + engine_RawRingElementArrayArray entries = + getmemarraytype(engine_RawRingElementArrayArray, nrows); + entries->len = nrows; + + const Ring *R = get_ring(); + RingElement *zero = RingElement::make_raw(R, R->zero()); + for (int r = 0; r < nrows; r++) + { + engine_RawRingElementArray currRow = + getmemarraytype(engine_RawRingElementArray, ncols); + currRow->len = ncols; + std::fill_n(currRow->array, ncols, zero); + entries->array[r] = currRow; + } + + for (int c = 0; c < ncols; c++) + { + const vec &column = elem(c); + for (const vecterm &term : column) + { + if (term.comp < 0 || term.comp >= nrows) + { + ERROR("internal error: matrix contains invalid entries:" + "row index %d out of range 0 .. %d", + term.comp, + nrows - 1); + continue; + } + entries->array[term.comp]->array[c] = + RingElement::make_raw(R, term.coeff); + } + } + return entries; +} + const Matrix /* or null */ *Matrix::make(const FreeModule *target, int ncols, const engine_RawRingElementArray M) @@ -278,6 +343,50 @@ const Matrix /* or null */ *Matrix::remake(const FreeModule *target) const return mat.to_matrix(); } +const Matrix /* or null */ *Matrix::promote(const FreeModule *target) const +{ + ring_elem a; + const Ring *R = get_ring(); + const Ring *S = target->get_ring(); + MatrixConstructor mat(target, n_cols()); + Matrix::iterator i(this); + for (int c = 0; c < n_cols(); c++) + for (i.set(c); i.valid(); i.next()) + if (S->promote(R, i.entry(), a)) + mat.set_entry(i.row(), c, a); + else + { + ERROR("first error occurred while promoting matrix entry at row %d, column %d", + i.row(), + c); + return nullptr; + } + mat.compute_column_degrees(); + return mat.to_matrix(); +} + +const Matrix /* or null */ *Matrix::lift(const FreeModule *target) const +{ + ring_elem a; + const Ring *R = get_ring(); + const Ring *S = target->get_ring(); + MatrixConstructor mat(target, n_cols()); + Matrix::iterator i(this); + for (int c = 0; c < n_cols(); c++) + for (i.set(c); i.valid(); i.next()) + if (R->lift(S, i.entry(), a)) + mat.set_entry(i.row(), c, a); + else + { + ERROR("first error occurred while lifting matrix entry at row %d, column %d", + i.row(), + c); + return nullptr; + } + mat.compute_column_degrees(); + return mat.to_matrix(); +} + const Matrix /* or null */ *Matrix::make(const MonomialIdeal *mi) { const PolynomialRing *P = mi->get_ring()->cast_to_PolynomialRing(); @@ -626,6 +735,81 @@ Matrix *Matrix::concat(const Matrix &m) const return mat.to_matrix(); } +Matrix *Matrix::concat(unsigned int n, const Matrix *const matrices[]) +{ + if (n == 0) + { + ERROR("matrix concat: expects at least one matrix"); + return nullptr; + } + + const FreeModule *F = matrices[0]->rows(); + const Ring *R = F->get_ring(); + MatrixConstructor mat(F, 0); + int next = 0; + for (unsigned int i = 0; i < n; i++) + { + const Matrix *M = matrices[i]; + if (R != M->get_ring()) + { + ERROR("matrix concat: different base rings"); + return nullptr; + } + if (F->rank() != M->n_rows()) + { + ERROR("matrix concat: row sizes are not equal"); + return nullptr; + } + for (int j = 0; j < M->n_cols(); j++) + { + mat.append(R->copy_vec(M->elem(j))); + mat.set_column_degree(next++, M->cols()->degree(j)); + } + } + return mat.to_matrix(); +} + +const Matrix *Matrix::direct_sum(unsigned int n, const Matrix *const matrices[]) +{ + if (n == 0) + { + ERROR("matrix direct sum: expects at least one matrix"); + return nullptr; + } + + const Ring *R = matrices[0]->get_ring(); + for (unsigned int i = 1; i < n; i++) + if (R != matrices[i]->get_ring()) + { + ERROR("matrix direct sum: different base rings"); + return nullptr; + } + + if (n == 1) return matrices[0]; + + FreeModule *F = R->make_FreeModule(); + FreeModule *G = R->make_FreeModule(); + for (unsigned int i = 0; i < n; i++) + { + F->direct_sum_to(matrices[i]->rows()); + G->direct_sum_to(matrices[i]->cols()); + } + + MatrixConstructor mat(F, G, nullptr); + int row_offset = 0; + int col_offset = 0; + for (unsigned int i = 0; i < n; i++) + { + const Matrix *M = matrices[i]; + for (int j = 0; j < M->n_cols(); j++) + mat.set_column(col_offset + j, + R->component_shift(row_offset, M->elem(j))); + row_offset += M->n_rows(); + col_offset += M->n_cols(); + } + return mat.to_matrix(); +} + Matrix *Matrix::direct_sum(const Matrix *m) const { auto R = get_ring(); diff --git a/M2/Macaulay2/e/matrix.hpp b/M2/Macaulay2/e/matrix.hpp index b7d2fce93f8..7e545dc2fc3 100644 --- a/M2/Macaulay2/e/matrix.hpp +++ b/M2/Macaulay2/e/matrix.hpp @@ -10,6 +10,7 @@ #include class MatrixConstructor; +class RingElement; /** * \ingroup matrices @@ -80,6 +81,9 @@ class Matrix : public EngineObject const Matrix /* or null */ *remake(const FreeModule *target) const; + const Matrix /* or null */ *promote(const FreeModule *target) const; + const Matrix /* or null */ *lift(const FreeModule *target) const; + static const Matrix *make(const MonomialIdeal *mi); const Ring *get_ring() const { return rows()->get_ring(); } @@ -89,6 +93,8 @@ class Matrix : public EngineObject ring_elem elem(int i, int j) const; vec &elem(int i) { return mEntries[i]; } const vec &elem(int i) const { return mEntries[i]; } + const RingElement /* or null */ *entry(int r, int c) const; + engine_RawRingElementArrayArrayOrNull entries() const; /*****************************************/ /* The non-const versions of these will go away */ @@ -117,7 +123,8 @@ class Matrix : public EngineObject Matrix *scalar_mult(const ring_elem r, bool opposite_mult) const; Matrix *mult(const Matrix *m, bool opposite_mult) const; Matrix *concat(const Matrix &m) const; - + static Matrix /* or null */ *concat(unsigned int n, + const Matrix *const matrices[]); static Matrix *identity(const FreeModule *F); static Matrix /* or null */ *zero(const FreeModule *F, const FreeModule *G); @@ -131,6 +138,8 @@ class Matrix : public EngineObject static Matrix /* or null */ *flip(const FreeModule *G, const FreeModule *H); Matrix /* or null */ *direct_sum(const Matrix *m) const; + static const Matrix /* or null */ *direct_sum(unsigned int n, + const Matrix *const matrices[]); Matrix /* or null */ *module_tensor(const Matrix *m) const; Matrix /* or null */ *tensor(const Matrix *m) const; Matrix /* or null */ *diff(const Matrix *m, int use_coef) const; diff --git a/M2/Macaulay2/e/monoid.hpp b/M2/Macaulay2/e/monoid.hpp index 418adcf136d..e47e87a3863 100644 --- a/M2/Macaulay2/e/monoid.hpp +++ b/M2/Macaulay2/e/monoid.hpp @@ -122,6 +122,7 @@ class Monoid : public MutableEngineObject const MonomialOrdering *getMonomialOrdering() const { return mo_; } const PolynomialRing *get_degree_ring() const { return mDegreeRing; } const Monoid *degree_monoid() const { return mDegreeMonoid; } + const std::vector &variableNames() const { return mVariableNames; } const_monomial degree_of_var(int v) const { return mDegreeOfVar[v]; } int primary_degree_of_var(int v) const { return mHeftDegrees[v]; } const std::vector &primary_degree_of_vars() const { return mHeftDegrees; } diff --git a/M2/Macaulay2/e/ring.cpp b/M2/Macaulay2/e/ring.cpp index a5d7d28306c..c3fcc6253e3 100644 --- a/M2/Macaulay2/e/ring.cpp +++ b/M2/Macaulay2/e/ring.cpp @@ -65,6 +65,36 @@ FreeModule *Ring::make_FreeModule(int n) const return new FreeModule(this, n, false); } +FreeModule *Ring::make_FreeModule(int ndegrees, int *degrees) const +{ + auto D = this->degree_monoid(); + unsigned int eachdeg = D->n_vars(); + + if (eachdeg == 0) + { + //ERROR("rawFreeModule: degree rank 0, but sequence of degrees given"); + throw exc::engine_error("rawFreeModule: degree rank 0, but sequence of degrees given"); + return nullptr; + } + + unsigned int rank = ndegrees / eachdeg; + if (rank * eachdeg != ndegrees) + { + //ERROR("inappropriate number of degrees"); + throw exc::engine_error("inappropriate number of degrees"); + return nullptr; + } + + monomial deg = D->make_one(); + FreeModule *F = new FreeModule(this, 0, false); + for (unsigned int i = 0; i < rank; i++) + { + D->from_expvector(degrees + i * eachdeg, deg); + F->append(deg); + } + return F; +} + bool Ring::is_field() const { return _isfield == 1; } bool Ring::declare_field() { diff --git a/M2/Macaulay2/e/ring.hpp b/M2/Macaulay2/e/ring.hpp index e3c209db27d..024e65c25a5 100644 --- a/M2/Macaulay2/e/ring.hpp +++ b/M2/Macaulay2/e/ring.hpp @@ -306,6 +306,7 @@ class Ring : public MutableEngineObject virtual FreeModule *make_FreeModule() const; virtual FreeModule *make_Schreyer_FreeModule() const; virtual FreeModule *make_FreeModule(int n) const; + virtual FreeModule *make_FreeModule(int ndegrees, int *degrees) const; virtual SumCollector *make_SumCollector() const; diff --git a/M2/Macaulay2/e/sagbi.hpp b/M2/Macaulay2/e/sagbi.hpp index 17deeaef39e..f6f3282c8ab 100644 --- a/M2/Macaulay2/e/sagbi.hpp +++ b/M2/Macaulay2/e/sagbi.hpp @@ -7,31 +7,41 @@ #include "comp-gb.hpp" /** - @ingroup comp - - @brief Helper routines for computing Sagbi bases. Not currently functional? +* A basic class that implements SAGBI bases (aka canonical subalgebra bases, or Khovanskii bases.) */ - class sagbi { public: - static ring_elem subduct(int numslots, - const PolyRing *R, - ring_elem f, - const RingMap *phi, - GBComputation *J); - static Matrix *subduct(int numparts, const Matrix *m, const RingMap *phi, GBComputation *J); + /** + * A subduction routine for a single ring element: unclear if/how this is ever used. + */ + static ring_elem subduct(int numslots, /**< the number of blocks in a monomial order */ + const PolyRing *R, /**< a polynomial ring */ + ring_elem f, /**< the ring element to be subducted */ + const RingMap *phi, /**< a ring map presenting an algebra */ + GBComputation *J /**< a Groebner basis computation object */ + ); + + /** + * A subduction routine for multiple ring elements in a matrix: unclear if/how this is ever used. + */ + static Matrix *subduct(int numparts, /**< the number of blocks in a monomial order */ + const Matrix *m, /**< a matrix of ring elements to be subducted */ + const RingMap *phi,/**< a ring map presenting an algebra */ + GBComputation *J /**< a Groebner basis computation object */ + ); static ring_elem subduct1(int numslots, - const PolyRing *T, // this is the tensor ring - const PolyRing *S, // this is the poly ring - ring_elem a, - const RingMap *inclusionAmbient, - const RingMap *fullSubstitution, - const RingMap *substitutionInclusion, - GBComputation *gbI, - GBComputation *gbReductionIdeal); + const PolyRing *T, /**< this is the tensor ring (containing original variables for subalgebra and extra variables tagging algebra generators) */ + const PolyRing *S, /**< polynomial ring containing original variables*/ + ring_elem a, /**< this is the ring element to be subducteda*/ + const RingMap *inclusionAmbient, /**< the inclusion map reprensenting a subalgebra */ + const RingMap *fullSubstitution, /**< substitution map sending tag variables to generators */ + const RingMap *substitutionInclusion, /**< combined */ + GBComputation *gbI, /**< Groebner basis of ideal in S encoding quotient ideal */ + GBComputation *gbReductionIdeal /**< Groebner basis of ideal in T encoding generator reductions and lead term of I */ + ); static Matrix *subduct1(int numparts, const Ring *rawT, diff --git a/M2/Macaulay2/e/unit-tests/ARingGFTest.cpp b/M2/Macaulay2/e/unit-tests/ARingGFTest.cpp index 1917a3ee34d..59f4fce54e7 100644 --- a/M2/Macaulay2/e/unit-tests/ARingGFTest.cpp +++ b/M2/Macaulay2/e/unit-tests/ARingGFTest.cpp @@ -47,13 +47,17 @@ static int randomVals[nelements] = { template <> void getElement(const M2::ARingGFFlint& R, int index, - M2::ARingGFGivaro::ElementType& result) + M2::ARingGFFlint::ElementType& result) { M2::ARingGFFlint::ElementType gen; R.init(gen); R.getGenerator(gen); if (index >= nelements) - R.power(result, gen, rawRandomInt(static_cast(R.cardinality()))); + { + long card = 1; + for (int i = 0; i < R.dimension(); i++) card *= R.characteristic(); + R.power(result, gen, rawRandomInt(static_cast(card))); + } else R.power(result, gen, randomVals[index]); R.clear(gen); @@ -64,7 +68,8 @@ TEST(ARingGFFlint, create) M2::ARingGFFlint R(5, 3); EXPECT_EQ(ringName(R), "GF(5,3,Flint)"); - EXPECT_EQ(R.cardinality(), 125); + // TODO: implement cardinality() in ARingGFFlint + // EXPECT_EQ(R.cardinality(), 125); EXPECT_EQ(R.characteristic(), 5); M2_arrayint gen_modpoly = R.getModPolynomialCoeffs(); diff --git a/M2/Macaulay2/e/unit-tests/GivaroTest.cpp b/M2/Macaulay2/e/unit-tests/GivaroTest.cpp deleted file mode 100644 index 60ac45f908e..00000000000 --- a/M2/Macaulay2/e/unit-tests/GivaroTest.cpp +++ /dev/null @@ -1,197 +0,0 @@ -// Copyright (c) 994-2009 by The Givaro group -// This file is part of Givaro. -// Givaro is governed by the CeCILL-B license under French law -// and abiding by the rules of distribution of free software. -// see the COPYRIGHT file for more details. -#include -#include -#include -#include -#include -using namespace Givaro; -template -void TestField(const Field& F) -{ - std::cerr << "Within "; - F.write(std::cerr); - std::cerr << " : " << std::flush; - typename Field::Element a, b, c, d; - F.init(a, 7U); - F.init(b, -29.3); - F.init(c); // empty constructor - F.init(d); // empty constructor - F.add(c, a, b); // c = a+b - // Separate output writing - F.write(std::cout, a) << " + " << std::flush; - F.write(std::cout, b) << " = " << std::flush; - F.write(std::cerr, c) << std::endl; - F.mul(c, a, b); // c = a*b - F.axpy(d, a, b, c); // d = a*b + c; - // Writing all outputs in a single command line - F.write(std::cerr << "Within ") << " : " << std::flush; - F.write(F.write(F.write(F.write(std::cout, c) << " + ", a) << " * ", b) - << " = ", - d) - << std::endl; - { - typename Field::Element e; - F.init(e); - F.assign(e, d); - F.maxpy(e, a, b, d); // e = d-a*b - // Writing all outputs in a single command line - F.write(std::cerr << "Within ") << " : " << std::flush; - F.write(F.write(F.write(F.write(std::cout, d) << " - ", a) << " * ", b) - << " = ", - e) - << std::endl; - } - { - typename Field::Element e; - F.init(e); - F.assign(e, d); - F.maxpyin(e, a, b); // e = d - a*b; - // Writing all outputs in a single command line - F.write(std::cerr << "Within ") << " : " << std::flush; - F.write(F.write(F.write(F.write(std::cout, d) << " - ", a) << " * ", b) - << " = ", - e) - << std::endl; - } - { - typename Field::Element e; - F.init(e); - F.assign(e, d); - F.axmy(e, a, b, d); // e = a*b -d; - // Writing all outputs in a single command line - F.write(std::cerr << "Within ") << " : " << std::flush; - F.write(F.write(F.write(F.write(std::cout, a) << " * ", b) << " - ", d) - << " = ", - e) - << std::endl; - } - { - typename Field::Element e; - F.init(e); - F.assign(e, d); - F.maxpyin(e, a, b); // e = d - a*b; - // Writing all outputs in a single command line - F.write(std::cerr << "Within ") << " : " << std::flush; - F.write(F.write(F.write(F.write(std::cout, d) << " - ", a) << " * ", b) - << " = ", - e) - << std::endl; - } - // Four operations - F.write(F.write(std::cout, a) << " += ", b) << " is "; - F.write(std::cout, F.addin(a, b)) << " ; "; - F.write(F.write(std::cout, a) << " -= ", b) << " is "; - F.write(std::cout, F.subin(a, b)) << " ; "; - F.write(F.write(std::cout, a) << " *= ", b) << " is "; - F.write(std::cout, F.mulin(a, b)) << " ; "; - F.write(F.write(std::cout, a) << " /= ", b) << " is "; - F.write(std::cout, F.divin(a, b)) << std::endl; - F.init(a, 22996); - F.inv(b, a); - F.write(F.write(std::cout << "1/", a) << " is ", b) << std::endl; - F.mul(c, b, a); - F.write(std::cout << "1 is ", c) << std::endl; - F.init(a, 22996); - F.init(b, 22996); - F.write(std::cout << "1/", a) << " is "; - F.invin(a); - F.write(std::cout, a) << std::endl; - F.mulin(a, b); - F.write(std::cout << "1 is ", a) << std::endl; - F.init(a, 37403); - F.inv(b, a); - F.write(F.write(std::cout << "1/", a) << " is ", b) << std::endl; - F.mul(c, b, a); - F.write(std::cout << "1 is ", c) << std::endl; - F.init(a, 37403); - F.init(b, 37403); - F.write(std::cout << "1/", a) << " is "; - F.invin(a); - F.write(std::cout, a) << std::endl; - F.mulin(a, b); - F.write(std::cout << "1 is ", a) << std::endl; -} -extern "C" { -#include -#include -} -int main(int argc, char** argv) -{ - // modulo 13 over 16 bits - Modular C13(13); - TestField(C13); - // modulo 13 over 32 bits - Modular Z13(13); - TestField(Z13); - // modulo 13 over unsigned 32 bits - Modular U13(13); - TestField(U13); -#ifdef __USE_Givaro_SIXTYFOUR__ - // modulo 13 over 64 bits - Modular LL13(13U); - TestField(LL13); -#endif - // modulo 13 fully tabulated - Modular L13(13); - TestField(L13); - // modulo 13 over 32 bits with Montgomery reduction - Montgomery M13(13); - TestField(M13); - Montgomery M3(39989); - TestField(M3); - // modulo 13 with primitive root representation - GFqDom GF13(13); - TestField(GF13); - // modulo 13 over arbitrary size - Modular IntZ13(13); - TestField(IntZ13); - // Zech log finite field with 5^4 elements - GFqDom GF625(5, 4); - TestField(GF625); -#if 0 // Possibly related: https://github.com/cr-marcstevens/m4gb/issues/8 - // Zech log finite field with 256 elements - // and prescribed irreducible polynomial - std::vector::Residu_t> Irred(9); - Irred[0] = 1; - Irred[1] = 1; - Irred[2] = 0; - Irred[3] = 1; - Irred[4] = 1; - Irred[5] = 0; - Irred[6] = 0; - Irred[7] = 0; - Irred[8] = 1; - GFqDom F256(2, 8, Irred); - TestField(F256); -#endif - // Zech log finite field with 3^4 elements - // Using the Q-adic Transform - GFqExt GF81(3, 4); - TestField(GF81); - // Zech log finite field with 2Mb tables - struct rusage tmp1; - getrusage(RUSAGE_SELF, &tmp1); - // user time - double tim = (double)tmp1.ru_utime.tv_sec + - ((double)tmp1.ru_utime.tv_usec) / (1000000.0); - ; - getrusage(RUSAGE_SELF, &tmp1); - tim = (double)tmp1.ru_utime.tv_sec + - ((double)tmp1.ru_utime.tv_usec) / (1000000.0) - tim; - std::cerr << "Initialization took " << tim - << " cpu seconds and : " << std::endl; - std::cerr << tmp1.ru_maxrss << " maximum resident set size" << std::endl - << tmp1.ru_ixrss << " integral shared memory size" << std::endl - << tmp1.ru_idrss << " integral unshared data size" << std::endl - << tmp1.ru_isrss << " integral unshared stack size" << std::endl - << tmp1.ru_minflt << " page reclaims" << std::endl - << tmp1.ru_majflt << " page faults" << std::endl - << tmp1.ru_nswap << " swaps" << std::endl - << tmp1.ru_inblock << " block input operations" << std::endl - << tmp1.ru_oublock << " block output operations" << std::endl; - return 0; -} diff --git a/M2/Macaulay2/e/unit-tests/Makefile.files b/M2/Macaulay2/e/unit-tests/Makefile.files index cc4a8b5e3f8..4f02840bd4f 100644 --- a/M2/Macaulay2/e/unit-tests/Makefile.files +++ b/M2/Macaulay2/e/unit-tests/Makefile.files @@ -3,7 +3,8 @@ SHARED_UNITTEST_CCFILES := \ M2-cpp-replacement \ basics-test \ fromStream \ - util-polyring-creation + util-polyring-creation \ + RingElem UNITTEST_CCFILES := \ $(SHARED_UNITTEST_CCFILES) \ @@ -29,12 +30,12 @@ UNITTEST_CCFILES := \ NewF4Test \ MonoidTest \ MatrixIOTest \ + WeylAlgebraTest \ + QuotientRingTest \ RingRRRTest \ RingCCCTest -# WeylAlgebraTest \ # TODO: add in this file - -# ARingGFTest +# ARingGFTest \ UNITTEST_TARGET := testMain diff --git a/M2/Macaulay2/e/unit-tests/MatrixIOTest.cpp b/M2/Macaulay2/e/unit-tests/MatrixIOTest.cpp index 0d1b0fb5394..73f5fda8830 100644 --- a/M2/Macaulay2/e/unit-tests/MatrixIOTest.cpp +++ b/M2/Macaulay2/e/unit-tests/MatrixIOTest.cpp @@ -3,9 +3,12 @@ #include #include "util-polyring-creation.hpp" +#include "matrix-con.hpp" #include "matrix.hpp" #include "BasicPolyList.hpp" #include "BasicPolyListParser.hpp" +#include "relem.hpp" +#include "error.h" #include "gb-f4/PolynomialList.hpp" #include "gb-f4/GBF4Interface.hpp" #include "VectorArithmetic.hpp" @@ -95,6 +98,296 @@ TEST(MatrixIO, readMsolve) EXPECT_TRUE(M->n_cols() == 4); } +TEST(Matrix, entriesFromSparseColumns) +{ + const Ring* R = simplePolynomialRing(101, {"x"}); + const FreeModule* target = R->make_FreeModule(5); + MatrixConstructor mat(target, 15); + + mat.set_entry(1, 3, R->from_long(11)); + mat.set_entry(3, 2, R->from_long(22)); + mat.set_entry(4, 1, R->from_long(33)); + mat.compute_column_degrees(); + + const Matrix* M = mat.to_matrix(); + engine_RawRingElementArrayArray entries = M->entries(); + + ASSERT_NE(entries, nullptr); + EXPECT_EQ(entries->len, 5); + + for (int r = 0; r < 5; r++) + { + ASSERT_NE(entries->array[r], nullptr); + EXPECT_EQ(entries->array[r]->len, 15); + + for (int c = 0; c < 15; c++) + { + ASSERT_NE(entries->array[r]->array[c], nullptr); + EXPECT_EQ(entries->array[r]->array[c]->get_ring(), R); + + ring_elem expected = R->zero(); + if (r == 1 && c == 3) + expected = R->from_long(11); + else if (r == 3 && c == 2) + expected = R->from_long(22); + else if (r == 4 && c == 1) + expected = R->from_long(33); + + EXPECT_TRUE(R->is_equal(entries->array[r]->array[c]->get_value(), + expected)) + << "entry (" << r << ", " << c << ")"; + } + } +} + +TEST(Matrix, entry) +{ + const Ring* R = simplePolynomialRing(101, {"x"}); + const FreeModule* target = R->make_FreeModule(3); + MatrixConstructor mat(target, 4); + + mat.set_entry(1, 2, R->from_long(17)); + mat.compute_column_degrees(); + + const Matrix* M = mat.to_matrix(); + const RingElement* nonzero = M->entry(1, 2); + ASSERT_NE(nonzero, nullptr); + EXPECT_EQ(nonzero->get_ring(), R); + EXPECT_TRUE(R->is_equal(nonzero->get_value(), R->from_long(17))); + + const RingElement* zero = M->entry(0, 0); + ASSERT_NE(zero, nullptr); + EXPECT_EQ(zero->get_ring(), R); + EXPECT_TRUE(R->is_zero(zero->get_value())); + + EXPECT_EQ(M->entry(3, 0), nullptr); + EXPECT_TRUE(error()); + EXPECT_STREQ(error_message(), "matrix row index 3 out of range 0 .. 2"); + + EXPECT_EQ(M->entry(0, 4), nullptr); + EXPECT_TRUE(error()); + EXPECT_STREQ(error_message(), "matrix column index 4 out of range 0 .. 3"); +} + +TEST(Matrix, concatArray) +{ + const PolynomialRing* R = simplePolynomialRing(101, {"x", "y"}); + const FreeModule* target = R->make_FreeModule(2); + + MatrixConstructor left(target, 2); + left.set_entry(0, 0, R->from_long(2)); + left.set_entry(1, 1, R->var(0)); + left.compute_column_degrees(); + const Matrix* A = left.to_matrix(); + + MatrixConstructor empty(target, 0); + const Matrix* E = empty.to_matrix(); + + MatrixConstructor right(target, 2); + right.set_entry(1, 0, R->from_long(5)); + right.set_entry(0, 1, R->var(1)); + right.compute_column_degrees(); + const Matrix* B = right.to_matrix(); + + const Matrix* const matrices[] = {A, E, B}; + const Matrix* C = Matrix::concat(3, matrices); + + ASSERT_NE(C, nullptr); + EXPECT_EQ(C->rows(), target); + EXPECT_EQ(C->n_rows(), 2); + EXPECT_EQ(C->n_cols(), 4); + + EXPECT_TRUE(R->is_equal(C->elem(0, 0), R->from_long(2))); + EXPECT_TRUE(R->is_zero(C->elem(1, 0))); + EXPECT_TRUE(R->is_zero(C->elem(0, 1))); + EXPECT_TRUE(R->is_equal(C->elem(1, 1), R->var(0))); + EXPECT_TRUE(R->is_zero(C->elem(0, 2))); + EXPECT_TRUE(R->is_equal(C->elem(1, 2), R->from_long(5))); + EXPECT_TRUE(R->is_equal(C->elem(0, 3), R->var(1))); + EXPECT_TRUE(R->is_zero(C->elem(1, 3))); + + const Monoid* D = R->degree_monoid(); + EXPECT_EQ(D->compare(C->cols()->degree(0), A->cols()->degree(0)), 0); + EXPECT_EQ(D->compare(C->cols()->degree(1), A->cols()->degree(1)), 0); + EXPECT_EQ(D->compare(C->cols()->degree(2), B->cols()->degree(0)), 0); + EXPECT_EQ(D->compare(C->cols()->degree(3), B->cols()->degree(1)), 0); +} + +TEST(Matrix, directSum) +{ + const PolynomialRing* R = simplePolynomialRing(101, {"x", "y"}); + const FreeModule* target = R->make_FreeModule(2); + + MatrixConstructor left(target, 2); + left.set_entry(0, 0, R->from_long(2)); + left.set_entry(1, 1, R->var(0)); + left.compute_column_degrees(); + const Matrix* A = left.to_matrix(); + + MatrixConstructor mid(target, 1); + mid.set_entry(0, 0, R->from_long(7)); + mid.compute_column_degrees(); + const Matrix* M = mid.to_matrix(); + + MatrixConstructor right(target, 2); + right.set_entry(1, 0, R->from_long(5)); + right.set_entry(0, 1, R->var(1)); + right.compute_column_degrees(); + const Matrix* B = right.to_matrix(); + + const Matrix* const singleton[] = {A}; + EXPECT_EQ(Matrix::direct_sum(1, singleton), A); + + const Matrix* const matrices[] = {A, M, B}; + const Matrix* C = Matrix::direct_sum(3, matrices); + + ASSERT_NE(C, nullptr); + EXPECT_EQ(C->n_rows(), 6); + EXPECT_EQ(C->n_cols(), 5); + + EXPECT_TRUE(R->is_equal(C->elem(0, 0), R->from_long(2))); // from A + EXPECT_TRUE(R->is_zero(C->elem(1, 0))); + EXPECT_TRUE(R->is_zero(C->elem(2, 0))); + EXPECT_TRUE(R->is_zero(C->elem(3, 0))); + EXPECT_TRUE(R->is_zero(C->elem(4, 0))); + EXPECT_TRUE(R->is_zero(C->elem(5, 0))); + + EXPECT_TRUE(R->is_zero(C->elem(0, 1))); + EXPECT_TRUE(R->is_equal(C->elem(1, 1), R->var(0))); // from A + EXPECT_TRUE(R->is_zero(C->elem(2, 1))); + EXPECT_TRUE(R->is_zero(C->elem(3, 1))); + EXPECT_TRUE(R->is_zero(C->elem(4, 1))); + EXPECT_TRUE(R->is_zero(C->elem(5, 1))); + + EXPECT_TRUE(R->is_zero(C->elem(0, 2))); + EXPECT_TRUE(R->is_zero(C->elem(1, 2))); + EXPECT_TRUE(R->is_equal(C->elem(2, 2), R->from_long(7))); // from M + EXPECT_TRUE(R->is_zero(C->elem(3, 2))); + EXPECT_TRUE(R->is_zero(C->elem(4, 2))); + EXPECT_TRUE(R->is_zero(C->elem(5, 2))); + + EXPECT_TRUE(R->is_zero(C->elem(0, 3))); + EXPECT_TRUE(R->is_zero(C->elem(1, 3))); + EXPECT_TRUE(R->is_zero(C->elem(2, 3))); + EXPECT_TRUE(R->is_zero(C->elem(3, 3))); + EXPECT_TRUE(R->is_zero(C->elem(4, 3))); + EXPECT_TRUE(R->is_equal(C->elem(5, 3), R->from_long(5))); // from B + + EXPECT_TRUE(R->is_zero(C->elem(0, 4))); + EXPECT_TRUE(R->is_zero(C->elem(1, 4))); + EXPECT_TRUE(R->is_zero(C->elem(2, 4))); + EXPECT_TRUE(R->is_zero(C->elem(3, 4))); + EXPECT_TRUE(R->is_equal(C->elem(4, 4), R->var(1))); // from B + EXPECT_TRUE(R->is_zero(C->elem(5, 4))); + + const Monoid* D = R->degree_monoid(); + EXPECT_EQ(D->compare(C->rows()->degree(0), A->rows()->degree(0)), 0); + EXPECT_EQ(D->compare(C->rows()->degree(1), A->rows()->degree(1)), 0); + EXPECT_EQ(D->compare(C->rows()->degree(2), M->rows()->degree(0)), 0); + EXPECT_EQ(D->compare(C->rows()->degree(3), M->rows()->degree(1)), 0); + EXPECT_EQ(D->compare(C->rows()->degree(4), B->rows()->degree(0)), 0); + EXPECT_EQ(D->compare(C->rows()->degree(5), B->rows()->degree(1)), 0); + EXPECT_EQ(D->compare(C->cols()->degree(0), A->cols()->degree(0)), 0); + EXPECT_EQ(D->compare(C->cols()->degree(1), A->cols()->degree(1)), 0); + EXPECT_EQ(D->compare(C->cols()->degree(2), M->cols()->degree(0)), 0); + EXPECT_EQ(D->compare(C->cols()->degree(3), B->cols()->degree(0)), 0); + EXPECT_EQ(D->compare(C->cols()->degree(4), B->cols()->degree(1)), 0); +} + +TEST(Matrix, promote) +{ + const Ring* ZZ = globalZZ; + const PolynomialRing* P = degreeRing({"x"}); + + const FreeModule* zzTarget = ZZ->make_FreeModule(3); + MatrixConstructor mat(zzTarget, 4); + + mat.set_entry(0, 0, ZZ->from_long(7)); + mat.set_entry(2, 1, ZZ->from_long(-3)); + mat.set_entry(1, 3, ZZ->from_long(11)); + mat.compute_column_degrees(); + + const Matrix* M = mat.to_matrix(); + const FreeModule* polynomialTarget = P->make_FreeModule(3); + const Matrix* promoted = M->promote(polynomialTarget); + + ASSERT_NE(promoted, nullptr); + EXPECT_EQ(promoted->rows(), polynomialTarget); + EXPECT_EQ(promoted->get_ring(), P); + EXPECT_EQ(promoted->n_rows(), 3); + EXPECT_EQ(promoted->n_cols(), 4); + + EXPECT_TRUE(P->is_equal(promoted->elem(0, 0), P->from_long(7))); + EXPECT_TRUE(P->is_equal(promoted->elem(2, 1), P->from_long(-3))); + EXPECT_TRUE(P->is_equal(promoted->elem(1, 3), P->from_long(11))); + EXPECT_TRUE(P->is_zero(promoted->elem(1, 1))); +} + +TEST(Matrix, promoteFailureReportsFirstEntry) +{ + const Ring* ZZ = globalZZ; + const PolynomialRing* P = degreeRing({"x"}); + + const FreeModule* polynomialTarget = P->make_FreeModule(3); + MatrixConstructor mat(polynomialTarget, 2); + + mat.set_entry(2, 0, P->from_long(7)); + mat.set_entry(0, 1, P->from_long(11)); + mat.compute_column_degrees(); + + const Matrix* M = mat.to_matrix(); + const Matrix* promoted = M->promote(ZZ->make_FreeModule(3)); + + EXPECT_EQ(promoted, nullptr); + EXPECT_TRUE(error()); + EXPECT_STREQ(error_message(), + "first error occurred while promoting matrix entry at row 2, column 0"); +} + +TEST(Matrix, liftFromPolynomialRingToCoefficientRing) +{ + const Ring* ZZ = globalZZ; + const PolynomialRing* P = degreeRing({"x"}); + + const FreeModule* polynomialTarget = P->make_FreeModule(3); + MatrixConstructor mat(polynomialTarget, 4); + + mat.set_entry(0, 0, P->from_long(7)); + mat.set_entry(2, 1, P->from_long(-3)); + mat.set_entry(1, 3, P->from_long(11)); + mat.compute_column_degrees(); + + const Matrix* M = mat.to_matrix(); + const Matrix* lifted = M->lift(ZZ->make_FreeModule(3)); + + ASSERT_NE(lifted, nullptr); + EXPECT_EQ(lifted->get_ring(), ZZ); + EXPECT_EQ(lifted->n_rows(), 3); + EXPECT_EQ(lifted->n_cols(), 4); + + EXPECT_TRUE(ZZ->is_equal(lifted->elem(0, 0), ZZ->from_long(7))); + EXPECT_TRUE(ZZ->is_equal(lifted->elem(2, 1), ZZ->from_long(-3))); + EXPECT_TRUE(ZZ->is_equal(lifted->elem(1, 3), ZZ->from_long(11))); + EXPECT_TRUE(ZZ->is_zero(lifted->elem(1, 1))); +} + +TEST(Matrix, liftFailureReturnsNull) +{ + const Ring* ZZ = globalZZ; + const PolynomialRing* P = degreeRing({"x"}); + + const FreeModule* polynomialTarget = P->make_FreeModule(2); + MatrixConstructor mat(polynomialTarget, 1); + mat.set_entry(0, 0, P->var(0)); + mat.compute_column_degrees(); + + const Matrix* M = mat.to_matrix(); + EXPECT_EQ(M->lift(ZZ->make_FreeModule(2)), nullptr); + EXPECT_TRUE(error()); + EXPECT_STREQ(error_message(), + "first error occurred while lifting matrix entry at row 0, column 0"); +} + #if 0 TEST(MatrixIO, readMsolveBig1) { diff --git a/M2/Macaulay2/e/unit-tests/QuotientRingTest.cpp b/M2/Macaulay2/e/unit-tests/QuotientRingTest.cpp new file mode 100644 index 00000000000..6363146159e --- /dev/null +++ b/M2/Macaulay2/e/unit-tests/QuotientRingTest.cpp @@ -0,0 +1,65 @@ +// Copyright 2026, The Macaulay2 Authors. +// Tests for ideal creation, Groebner bases, and quotient rings. + +#include +#include "util-polyring-creation.hpp" +#include "RingElem.hpp" +#include "matrix.hpp" +#include "debug.hpp" +// Step 1: idealFromStrings +TEST(IdealCreation, fromStrings) +{ + const PolynomialRing* R = simplePolynomialRing(101, {"x", "y", "z"}); + const Matrix* I = idealFromStrings(R, {"x^2+y", "y^2-z"}); + dmatrix(I); + ASSERT_NE(I, nullptr); + EXPECT_EQ(I->n_rows(), 1); + EXPECT_EQ(I->n_cols(), 2); +} + +// Step 2: computeGB +TEST(GroebnerBasis, simple) +{ + const PolynomialRing* R = simplePolynomialRing(101, {"x", "y"}); + const Matrix* I = idealFromStrings(R, {"x^2-y", "x*y-1"}); + const Matrix* gb = computeGB(I); + ASSERT_NE(gb, nullptr); + dmatrix(gb); + EXPECT_GE(gb->n_cols(), 2); +} + +// Step 3: simpleQuotientRing +TEST(QuotientRing, arithmetic) +{ + const PolynomialRing* R = simplePolynomialRing(101, {"x", "y"}); + const Ring* Q = simpleQuotientRing(R, {"x^2-1"}); + ASSERT_NE(Q, nullptr); + auto x = RingElem::var(Q, 0); + auto one = RingElem::fromInt(Q, 1); + EXPECT_EQ(x * x, one); // x^2 == 1 in Q +} + +TEST(QuotientRing, sphere) +{ + // QQ[x,y,z] / (x^2+y^2+z^2-1) + const PolynomialRing* R = simplePolynomialRing(0, {"x", "y", "z"}); + const Ring* Q = simpleQuotientRing(R, {"x^2+y^2+z^2-1"}); + ASSERT_NE(Q, nullptr); + + auto x = RingElem::var(Q, 0); + auto y = RingElem::var(Q, 1); + auto z = RingElem::var(Q, 2); + auto one = RingElem::fromInt(Q, 1); + + // The defining relation: x^2 + y^2 + z^2 == 1 + EXPECT_EQ(x*x + y*y + z*z, one); + + // Consequence: (x+y+z)^2 == 1 + 2*(x*y + x*z + y*z) + auto lhs = (x + y + z).power(2); + auto rhs = one + 2 * (x*y + x*z + y*z); + EXPECT_EQ(lhs, rhs); +} + +// Local Variables: +// indent-tabs-mode: nil +// End: diff --git a/M2/Macaulay2/e/unit-tests/RingElem.cpp b/M2/Macaulay2/e/unit-tests/RingElem.cpp new file mode 100644 index 00000000000..ae8c4bb7ed9 --- /dev/null +++ b/M2/Macaulay2/e/unit-tests/RingElem.cpp @@ -0,0 +1,66 @@ +// Copyright 2026, The Macaulay2 Authors. + +#include "RingElem.hpp" +#include "BasicPoly.hpp" +#include "polyring.hpp" +#include "monoid.hpp" + +// Convert a BasicPoly (parsed polynomial with mpz_class coefficients and +// varpower monomials) into a ring_elem in the given PolynomialRing. +static ring_elem basicPolyToRingElem(const PolynomialRing *P, const BasicPoly &bp) +{ + const Ring *K = P->getCoefficients(); + const Monoid *M = P->getMonoid(); + int nvars = M->n_vars(); + + ring_elem result = P->from_long(0); + std::vector exp(nvars, 0); + monomial monom = M->make_one(); + + int monomStart = 0; + for (size_t i = 0; i < bp.mCoefficients.size(); ++i) + { + // Convert coefficient: mpz_class -> ring_elem in coefficient ring + ring_elem coeff = K->from_int(bp.mCoefficients[i].get_mpz_t()); + + // Convert monomial: varpower format -> exponent vector -> encoded monomial + int monomLen = bp.mMonomials[monomStart]; + std::fill(exp.begin(), exp.end(), 0); + for (int j = monomStart + 1; j < monomStart + monomLen; j += 2) + { + int var = bp.mMonomials[j]; + int e = bp.mMonomials[j + 1]; + if (var >= 0 && var < nvars) + exp[var] = e; + } + M->from_expvector(exp.data(), monom); + + // Create the term and add it to the result + ring_elem term = P->make_flat_term(coeff, monom); + result = P->add(result, term); + + monomStart += monomLen; + } + + M->remove(monom); + return result; +} + +RingElem RingElem::fromString(const Ring *R, const std::string &s) +{ + // Polynomial rings (including WeylAlgebra, skew, quotient rings, etc.) + const PolynomialRing *P = R->cast_to_PolynomialRing(); + if (P != nullptr) + { + const Monoid *M = P->getMonoid(); + std::vector varnames = M->variableNames(); + BasicPoly bp = parseBasicPoly(s, varnames); + ring_elem result = basicPolyToRingElem(P, bp); + return RingElem(R, result); + } + + // Base rings: parse as an integer + // Works for ZZ, ZZ/p, QQ (integers only), and other rings with from_int + mpz_class val(s); + return RingElem(R, R->from_int(val.get_mpz_t())); +} diff --git a/M2/Macaulay2/e/unit-tests/RingElem.hpp b/M2/Macaulay2/e/unit-tests/RingElem.hpp new file mode 100644 index 00000000000..872f0b7ae09 --- /dev/null +++ b/M2/Macaulay2/e/unit-tests/RingElem.hpp @@ -0,0 +1,119 @@ +// Copyright 2026, The Macaulay2 Authors. +// +// RingElem: A lightweight value-semantics wrapper around (const Ring*, ring_elem). +// Unlike RingElement (which is heap-allocated and returns pointers from operators), +// RingElem lives on the stack and operators return values, making test code concise: +// +// auto x = RingElem::var(W, 0); +// auto Dx = RingElem::var(W, 2); +// auto one = RingElem::fromInt(W, 1); +// EXPECT_EQ(Dx * x - x * Dx, one); +// +// auto f = RingElem::fromString(R, "x^2+3*x*y-1"); + +#pragma once + +#include +#include +#include +#include + +#include "buffer.hpp" +#include "ring.hpp" +#include "polyring.hpp" +#include "monoid.hpp" +#include "BasicPoly.hpp" + +class RingElem +{ + const Ring *mRing; + ring_elem mValue; + + public: + // Primary constructor: wraps a pre-computed ring_elem. + // Use the static factories below for common cases. + RingElem(const Ring *R, ring_elem f) : mRing(R), mValue(f) {} + + // Factory: create the i-th variable of the ring + static RingElem var(const Ring *R, int i) { return RingElem(R, R->var(i)); } + + // Factory: create a ring element from an integer + static RingElem fromInt(const Ring *R, long n) + { + return RingElem(R, R->from_long(n)); + } + + // Factory: create a ring element from a string. + // For polynomial rings: parses "x^2+3*x*y-1" using variable names from the ring. + // For base rings (ZZ, ZZ/p): parses an integer. + // Throws parsing_error on failure. + // NOTE: fromString and toString are not yet inverses of each other. + // toString outputs e.g. "x3+2xyz" while fromString expects "x^3+2*x*y*z". + // TODO: make these round-trip compatible. + static RingElem fromString(const Ring *R, const std::string &s); + + // Accessors + const Ring *ring() const { return mRing; } + ring_elem value() const { return mValue; } + + // Predicates + bool isZero() const { return mRing->is_zero(mValue); } + bool isUnit() const { return mRing->is_unit(mValue); } + + // Comparison + bool operator==(const RingElem &b) const + { + assert(mRing == b.mRing && "RingElem comparison requires elements from the same ring"); + return mRing->is_equal(mValue, b.mValue); + } + bool operator!=(const RingElem &b) const { return !(*this == b); } + + // Arithmetic (returns values, not pointers) + RingElem operator-() const { return RingElem(mRing, mRing->negate(mValue)); } + + RingElem operator+(const RingElem &b) const + { + assert(mRing == b.mRing && "RingElem addition requires elements from the same ring"); + return RingElem(mRing, mRing->add(mValue, b.mValue)); + } + + RingElem operator-(const RingElem &b) const + { + assert(mRing == b.mRing && "RingElem subtraction requires elements from the same ring"); + return RingElem(mRing, mRing->subtract(mValue, b.mValue)); + } + + RingElem operator*(const RingElem &b) const + { + assert(mRing == b.mRing && "RingElem multiplication requires elements from the same ring"); + return RingElem(mRing, mRing->mult(mValue, b.mValue)); + } + + RingElem operator/(const RingElem &b) const + { + assert(mRing == b.mRing && "RingElem division requires elements from the same ring"); + return RingElem(mRing, mRing->divide(mValue, b.mValue)); + } + + // Scalar multiplication + RingElem operator*(long n) const + { + return RingElem(mRing, mRing->mult(mRing->from_long(n), mValue)); + } + friend RingElem operator*(long n, const RingElem &f) { return f * n; } + + RingElem power(int n) const { return RingElem(mRing, mRing->power(mValue, n)); } + + // String output (for gtest diagnostics) + std::string toString() const + { + buffer o; + mRing->elem_text_out(o, mValue); + return std::string(o.str()); + } + + friend std::ostream &operator<<(std::ostream &os, const RingElem &f) + { + return os << f.toString(); + } +}; diff --git a/M2/Macaulay2/e/unit-tests/SubsetTest.cpp b/M2/Macaulay2/e/unit-tests/SubsetTest.cpp index 7e031ab7d1a..d9d4f4b9307 100644 --- a/M2/Macaulay2/e/unit-tests/SubsetTest.cpp +++ b/M2/Macaulay2/e/unit-tests/SubsetTest.cpp @@ -7,81 +7,25 @@ #include #include "comb.hpp" - -TEST(Subsets, encode1) -{ - Subsets C(5, 2); - - Subset a(2, 0); - - for (int i = 0; i < 10; i++) - { - C.decode(i, a); - EXPECT_TRUE(C.isValid(a)); - size_t j = C.encode(a); - EXPECT_EQ(i, j); - } +constexpr int binom(int n, int k) { + return (k < 1 || n <= 1 || n <= k) ? 1 : (binom(n-1, k) * n/(n-k)); } -TEST(Subsets, encode2) -{ - Subsets C(12, 6); - - Subset a(6, 0); - - for (int i = 0; i < 924; i++) - { - C.decode(i, a); - EXPECT_TRUE(C.isValid(a)); - size_t j = C.encode(a); - EXPECT_EQ(i, j); - } -} - -TEST(Subsets, encode3) -{ - Subsets C(12, 0); - - Subset a(0, 0); - - for (int i = 0; i < 1; i++) - { - C.decode(i, a); - EXPECT_TRUE(C.isValid(a)); - size_t j = C.encode(a); - EXPECT_EQ(i, j); - } -} - -TEST(Subsets, encode4) -{ - Subsets C(21, 7); - - Subset a(7, 0); - - for (int i = 0; i < 116280; i++) - { - C.decode(i, a); - EXPECT_TRUE(C.isValid(a)); - size_t j = C.encode(a); - EXPECT_EQ(i, j); - } -} - -TEST(Subsets, encode5) -{ - Subsets C(21, 21); - - Subset a(21, 0); - - for (int i = 0; i < 1; i++) - { - C.decode(i, a); - EXPECT_TRUE(C.isValid(a)); - size_t j = C.encode(a); - EXPECT_EQ(i, j); - } -} +#define TESTSubsets(n, k, e) \ +TEST(Subsets, e) { \ + Subsets C(n, k); \ + Subset a(k, 0); \ + for (int i = 0; i < binom(n,k); i++) {\ + C.decode(i, a); \ + EXPECT_TRUE(C.isValid(a)); \ + size_t j = C.encode(a); \ + EXPECT_EQ(i, j); } } + +TESTSubsets(5, 2, encode1) +TESTSubsets(12, 6, encode2) +TESTSubsets(12, 0, encode3) +TESTSubsets(21, 7, encode4) +TESTSubsets(21,21, encode5) bool sameSubset(const Subset &a, const Subset &b) { diff --git a/M2/Macaulay2/e/unit-tests/WeylAlgebraTest.cpp b/M2/Macaulay2/e/unit-tests/WeylAlgebraTest.cpp new file mode 100644 index 00000000000..f365cea6fa2 --- /dev/null +++ b/M2/Macaulay2/e/unit-tests/WeylAlgebraTest.cpp @@ -0,0 +1,208 @@ +// Copyright 2026, The Macaulay2 Authors. + +#include "interface/ring.h" +#include "relem.hpp" +#include "util.hpp" +#include "weylalg.hpp" +#include "RingTest.hpp" +#include "RingElem.hpp" +#include "util-polyring-creation.hpp" + +class WeylAlgebraTestAccessor { + public: + static ring_elem binomial(const WeylAlgebra* W, int top, int bottom) { + return W->binomial(top, bottom); + } + static ring_elem multinomial(const WeylAlgebra* W, + const ring_elem a, + const_exponents exptop, + const_exponents expbottom) { + return W->multinomial(a, exptop, expbottom); + } + static const Ring* coefficientRing(const WeylAlgebra* W) { + return W->getCoefficients(); + } + static int nderivatives(const WeylAlgebra* W) { + return W->_nderivatives; + } +}; + +// Test fixture: creates a WeylAlgebra QQ[x,y,Dx,Dy] +class WeylAlgebraTest : public ::testing::Test { + protected: + const WeylAlgebra* W = nullptr; + const Ring* K = nullptr; // coefficient ring (QQ) + + void SetUp() override { + W = simpleWeylAlgebra(0, {"x", "y", "Dx", "Dy"}, {0, 1}, {2, 3}); + K = WeylAlgebraTestAccessor::coefficientRing(W); + } + + // Helper: check that a ring_elem in K equals a given long value + void expectEqualLong(ring_elem actual, long expected) { + ring_elem exp = K->from_long(expected); + EXPECT_TRUE(K->is_equal(actual, exp)); + } +}; + +TEST_F(WeylAlgebraTest, create) +{ + EXPECT_FALSE(W->is_commutative_ring()); + EXPECT_TRUE(W->is_weyl_algebra()); + EXPECT_EQ(4, W->n_vars()); + EXPECT_EQ(2, WeylAlgebraTestAccessor::nderivatives(W)); + + std::string ans {"WeylAlgebra(QQGMP[x,y,Dx,Dy,\n" + " DegreeLength => 1,\n" + " Degrees => {1, 1, -1, -1},\n" + " Heft => {1},\n" + " MonomialOrder => {\n" + " GRevLex => {1,1,1,1},\n" + " Position => Up\n" + " }\n" + " ])" }; + buffer o; + W->text_out(o); + EXPECT_EQ(ans, std::string(o.str())); +} + +TEST_F(WeylAlgebraTest, commutator) +{ + auto x = RingElem::var(W, 0); + auto y = RingElem::var(W, 1); + auto Dx = RingElem::var(W, 2); + auto Dy = RingElem::var(W, 3); + auto one = RingElem::fromInt(W, 1); + auto zero = RingElem::fromInt(W, 0); + + // Test Dx*x - x*Dx == 1 (the defining Weyl relation) + EXPECT_EQ(Dx * x - x * Dx, one); + + // Test Dy*y - y*Dy == 1 + EXPECT_EQ(Dy * y - y * Dy, one); + + // Test Dx*y - y*Dx == 0 (cross terms commute) + EXPECT_EQ(Dx * y - y * Dx, zero); +} + +TEST_F(WeylAlgebraTest, binomial) +{ + // binomial(n, k) = n! / (k! * (n-k)!) + // Results are ring elements in the coefficient ring K (QQ). + + // Base cases + expectEqualLong(WeylAlgebraTestAccessor::binomial(W, 0, 0), 1); + expectEqualLong(WeylAlgebraTestAccessor::binomial(W, 5, 0), 1); + expectEqualLong(WeylAlgebraTestAccessor::binomial(W, 5, 1), 5); + + // Standard values + expectEqualLong(WeylAlgebraTestAccessor::binomial(W, 5, 2), 10); + expectEqualLong(WeylAlgebraTestAccessor::binomial(W, 5, 3), 10); + expectEqualLong(WeylAlgebraTestAccessor::binomial(W, 5, 5), 1); + expectEqualLong(WeylAlgebraTestAccessor::binomial(W, 6, 3), 20); + expectEqualLong(WeylAlgebraTestAccessor::binomial(W, 10, 5), 252); + + // Larger values (beyond the cached binomtable) + expectEqualLong(WeylAlgebraTestAccessor::binomial(W, 20, 10), 184756); + expectEqualLong(WeylAlgebraTestAccessor::binomial(W, 30, 15), 155117520); +} + +TEST_F(WeylAlgebraTest, multinomial) +{ + // multinomial(c, top, bottom) = c * product_i binomial(top[i], bottom[i]) + // where the product is over the _nderivatives derivative pairs. + // Here _nderivatives = 2 (for Dx, Dy). + + ring_elem one = K->from_long(1); + + // All bottom entries zero: result should be c * 1 = c + { + int top[] = {5, 3}; + int bottom[] = {0, 0}; + expectEqualLong(WeylAlgebraTestAccessor::multinomial(W, one, top, bottom), 1); + } + + // Single nonzero bottom entry: should give binomial(5,2) = 10 + { + int top[] = {5, 3}; + int bottom[] = {2, 0}; + expectEqualLong(WeylAlgebraTestAccessor::multinomial(W, one, top, bottom), 10); + } + + // Both entries nonzero: binomial(5,2) * binomial(3,1) = 10 * 3 = 30 + { + int top[] = {5, 3}; + int bottom[] = {2, 1}; + expectEqualLong(WeylAlgebraTestAccessor::multinomial(W, one, top, bottom), 30); + } + + // With a scalar coefficient: 7 * binomial(4,2) * binomial(6,3) = 7 * 6 * 20 = 840 + { + ring_elem seven = K->from_long(7); + int top[] = {4, 6}; + int bottom[] = {2, 3}; + expectEqualLong(WeylAlgebraTestAccessor::multinomial(W, seven, top, bottom), 840); + } +} + +TEST_F(WeylAlgebraTest, fromString) +{ + // Single variable + auto x = RingElem::var(W, 0); + EXPECT_EQ(RingElem::fromString(W, "x"), x); + + // Monomial with coefficient + auto Dx = RingElem::var(W, 2); + EXPECT_EQ(RingElem::fromString(W, "3*x^2*Dx"), x * x * Dx * 3); + + // Polynomial with multiple terms + auto y = RingElem::var(W, 1); + auto expected = x * x + y * 3 - RingElem::fromInt(W, 1); + EXPECT_EQ(RingElem::fromString(W, "x^2+3*y-1"), expected); + + // Constant + EXPECT_EQ(RingElem::fromString(W, "5"), RingElem::fromInt(W, 5)); + + // Zero + EXPECT_TRUE(RingElem::fromString(W, "0").isZero()); + + // Negative coefficient + EXPECT_EQ(RingElem::fromString(W, "-x"), -x); +} + +TEST(PolyRingFromString, basic) +{ + const PolynomialRing *R = simplePolynomialRing(101, {"x", "y", "z"}); + + auto x = RingElem::var(R, 0); + auto y = RingElem::var(R, 1); + auto z = RingElem::var(R, 2); + auto one = RingElem::fromInt(R, 1); + + // Single variable + EXPECT_EQ(RingElem::fromString(R, "x"), x); + EXPECT_EQ(RingElem::fromString(R, "z"), z); + + // Monomial with coefficient + EXPECT_EQ(RingElem::fromString(R, "3*x^2*y"), x * x * y * 3); + + // Polynomial + auto f = x * x + y * 3 - one; + EXPECT_EQ(RingElem::fromString(R, "x^2+3*y-1"), f); + + // Coefficient reduction mod 101 + EXPECT_EQ(RingElem::fromString(R, "102*x"), x); + + // Constant + EXPECT_EQ(RingElem::fromString(R, "7"), RingElem::fromInt(R, 7)); + + // Zero + EXPECT_TRUE(RingElem::fromString(R, "0").isZero()); + + // Multi-term polynomial + // NOTE: can't round-trip via toString() yet — it outputs "x3+2xyz-y2+z" + // (no ^ or *) which the parser doesn't accept. TODO: make these compatible. + auto g = RingElem::fromString(R, "x^3+2*x*y*z-y^2+z"); + auto g_expected = x.power(3) + x * y * z * 2 - y * y + z; + EXPECT_EQ(g, g_expected); +} diff --git a/M2/Macaulay2/e/unit-tests/util-polyring-creation.cpp b/M2/Macaulay2/e/unit-tests/util-polyring-creation.cpp index e13e9850613..1e42cde8eee 100644 --- a/M2/Macaulay2/e/unit-tests/util-polyring-creation.cpp +++ b/M2/Macaulay2/e/unit-tests/util-polyring-creation.cpp @@ -1,4 +1,27 @@ +// Copyright 2026, The Macaulay2 Authors. + +#include "util.hpp" #include "util-polyring-creation.hpp" +#include "weylalg.hpp" +#include "interface/ring.h" +#include "interface/aring.h" +#include "BasicPoly.hpp" +#include "BasicPolyList.hpp" +#include "monoid.hpp" +#include "freemod.hpp" +#include "interface/groebner.h" +#include "comp-gb.hpp" + +// Creation of ring type objects +// coefficient rings +// degree ring (not degree monoid) +// monoid +// polynomialRing +// weylAlgebra +// exteriorAlgebra +// associativeAlgebra +// groebnerAlgebra + const Monoid* degreeMonoid(const std::vector& names) { @@ -31,6 +54,28 @@ const PolynomialRing* degreeRing(int ndegrees) return degreeRing({"T"}); } +// Monoids and monomial orderings + +const Monoid* simpleMonoid(const std::vector& names, + MonomialOrdering* monorder, + const PolynomialRing* degRing, + const std::vector& degs, + const std::vector& heft) +{ + // a few checks: + // (#vars of degrees ring) * #vars == #degs + // #heft == #gens degreesRing. + // heft of each degree vector for each vector should be > 0, if heft is non-empty. + + return Monoid::create( + monorder, + degRing, + names, + degs, + heft); +} + + const PolynomialRing* simplePolynomialRing(const Ring* kk, const std::vector& names, MonomialOrdering* monorder) @@ -61,7 +106,7 @@ const PolynomialRing* simplePolynomialRing(int p, const std::vector // heft is 1. // monomial order is grevlex. - const Ring *kk = (p > 0 ? rawARingZZpFlint(p) : rawARingQQFlint()); + const Ring *kk = (p > 0 ? rawARingZZpFlint(p) : IM2_Ring_QQ()); if (kk == nullptr) return nullptr; // one of these routines would have made an error. MonomialOrdering* monorder = MonomialOrderings::join @@ -73,6 +118,84 @@ const PolynomialRing* simplePolynomialRing(int p, const std::vector return simplePolynomialRing(kk, names, monorder); } +const WeylAlgebra* simpleWeylAlgebra(long p, + const std::vector names, + const std::vector comms, + const std::vector derivs) +{ + // if p is 0, use QQ. + // degrees are all set to 1. (degree ring has one variable) + // heft is 1. + // monomial order is grevlex. + + const Ring *kk = (p > 0 ? rawARingZZpFlint(p) : IM2_Ring_QQ()); + if (kk == nullptr) return nullptr; // one of these routines would have made an error. + + MonomialOrdering* monorder = MonomialOrderings::join + ({ + MonomialOrderings::GRevLex(names.size()), + MonomialOrderings::PositionUp() + }); + + int n = static_cast(comms.size()); + std::vector degs(2*n, -1); + std::fill_n(degs.begin(), n, 1); + + const Monoid* M = Monoid::create( + monorder, + degreeRing(1), + names, + degs, + {1}); + + M2_arrayint derivs1 = stdvector_to_M2_arrayint(derivs); + M2_arrayint comms1 = stdvector_to_M2_arrayint(comms); + auto W = WeylAlgebra::create(kk, + M, + derivs1, + comms1, + -1); + + return W; +} + + +const Matrix* idealFromStrings(const PolynomialRing* R, + const std::vector& polys) +{ + auto varnames = R->getMonoid()->variableNames(); + BasicPolyList bpList; + for (const auto& s : polys) + bpList.push_back(parseBasicPoly(s, varnames)); + return toMatrix(R->make_FreeModule(1), bpList); +} + +const Matrix* computeGB(const Matrix* M) +{ + M2_arrayint weights = stdvector_to_M2_arrayint(std::vector{}); + Computation* C = rawGBMake(M, + false, // collect_syz + 0, // n_rows_to_keep + weights, + false, // use_max_degree + 0, // max_degree + 0, // algorithm (default) + 0, // strategy (default) + 10); // max_reduction_count (engine default) + if (C == nullptr) return nullptr; + rawStartComputation(C); + return rawGBGetMatrix(C); +} + +const Ring* simpleQuotientRing(const PolynomialRing* R, + const std::vector& generators) +{ + const Matrix* I = idealFromStrings(R, generators); + if (I == nullptr) return nullptr; + const Matrix* gb = computeGB(I); + if (gb == nullptr) return nullptr; + return PolynomialRing::create_quotient(R, gb); +} // Local Variables: // indent-tabs-mode: nil diff --git a/M2/Macaulay2/e/unit-tests/util-polyring-creation.hpp b/M2/Macaulay2/e/unit-tests/util-polyring-creation.hpp index 95f55f69c3b..829bf4cf3a9 100644 --- a/M2/Macaulay2/e/unit-tests/util-polyring-creation.hpp +++ b/M2/Macaulay2/e/unit-tests/util-polyring-creation.hpp @@ -27,6 +27,28 @@ const PolynomialRing* simplePolynomialRing(const Ring* kk, // This create a polynomial ring with all degrees 1, and with GRevLex order const PolynomialRing* simplePolynomialRing(int p, const std::vector& names); +// Creates a Weyl algebra, with degree rank one, GRevLex monomial order. +const WeylAlgebra* simpleWeylAlgebra(long p, + const std::vector varnames, + const std::vector comms, + const std::vector derivs); + +class Matrix; + +// Create a 1-row matrix (ideal generators) from polynomial strings. +// Each string is a polynomial like "x^2+3*x*y-1". +const Matrix* idealFromStrings(const PolynomialRing* R, + const std::vector& polys); + +// Compute a Groebner basis of the ideal given by the 1-row matrix M. +// Returns a 1-row matrix whose columns are the GB elements. +const Matrix* computeGB(const Matrix* M); + +// Create quotient ring R / (generators). +// Computes the GB of the generators, then forms the quotient. +const Ring* simpleQuotientRing(const PolynomialRing* R, + const std::vector& generators); + // Local Variables: // indent-tabs-mode: nil // End: diff --git a/M2/Macaulay2/e/weylalg.hpp b/M2/Macaulay2/e/weylalg.hpp index 4a58d86d635..de65ab1fb19 100644 --- a/M2/Macaulay2/e/weylalg.hpp +++ b/M2/Macaulay2/e/weylalg.hpp @@ -11,6 +11,8 @@ class WeylAlgebra : public PolyRing { + friend class WeylAlgebraTestAccessor; + int _nderivatives; bool _homogeneous_weyl_algebra; int _homog_var; // Only used if 'homogeneous_weyl_algebra' is true. diff --git a/M2/Macaulay2/packages/undistributed-packages/GenerateD.m2 b/M2/Macaulay2/packages/undistributed-packages/GenerateD.m2 index c4521b46598..6808ad86eae 100644 --- a/M2/Macaulay2/packages/undistributed-packages/GenerateD.m2 +++ b/M2/Macaulay2/packages/undistributed-packages/GenerateD.m2 @@ -669,7 +669,7 @@ export rawHomogenizeMatrix(e:Expr):Expr := ( if isSmallInt(s.1) then (varindex := getSmallInt(s.1); if isSequenceOfSmallIntegers(s.2) then (varweights := getSequenceOfSmallIntegers(s.2); toExpr(Ccode(RawMatrixOrNull, - "IM2_Matrix_homogenize(", + "rawMatrixHomogenize(", M, ",", varindex, ",", weights,