Add Singular.jl extension and tests

- Implement Singular-backed methods for HomalgMatrix, syzygies (rows/columns and relative), bases, zero-decision, and left/right division (2- and 3-arg).
- Add singular-ext-test.jl and include in test runner.

Author: kamalsaleh

Project.toml

@@ -6,13 +6,21 @@ version = "0.1.6"
 [deps]
 Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
 
+[weakdeps]
+Singular = "bcd08a7b-43d2-5ff7-b6d4-c458787f915c"
+
+[extensions]
+MatricesForHomalgSingularExt = "Singular"
+
 [compat]
 Nemo = "0.48.0, 0.49, 0.50, 0.51, 0.52, 0.53"
+Singular = "0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28"
 julia = "1.9, 1.10, 1.11"
 
 [extras]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Singular = "bcd08a7b-43d2-5ff7-b6d4-c458787f915c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Documenter"]
+test = ["Test", "Documenter", "Singular"]

ext/MatricesForHomalgSingularExt.jl

@@ -0,0 +1,308 @@
+module MatricesForHomalgSingularExt
+
+using MatricesForHomalg
+import Singular
+
+"""
+    HomalgMatrix(entries, r, c, R)
+Construct a Singular matrix of size r x c over the polynomial ring R,
+with entries given by the nested array `entries` (a vector of vectors).
+The entries are converted to elements of R using the constructor R(...).
+
+# Examples
+```jldoctest
+julia> using Singular, MatricesForHomalg
+julia> R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+(Singular polynomial ring (QQ),(x,y),(dp(2),C), spoly{Singular.n_Q}[x, y])
+julia> entries = [[x, y], [x^2, y^2]]
+2-element Vector{Vector{Singular.n_Q}}:
+    [x, y]
+    [x^2, y^2]
+julia> A = HomalgMatrix(entries, 2, 2, R)
+[ x  y
+x^2  y^2 ]
+```
+"""
+function MatricesForHomalg.HomalgMatrix(entries, nr_rows, nr_cols, R::Singular.PolyRing)
+    m = MatricesForHomalg.zero_matrix(R, nr_rows, nr_cols)
+    flat = vcat(entries...)
+    for i in 1:nr_rows, j in 1:nr_cols
+        m[i, j] = R(flat[(i - 1) * nr_cols + j])
+    end
+    return m
+end
+
+"""
+    SyzygiesOfRows(A::Singular.smatrix)
+
+Compute the matrix of row syzygies of A using Singular.
+Returns a matrix whose rows span the left kernel of A,
+i.e. all rows X satisfying X * A = 0.
+
+# Examples
+```jldoctest
+julia> using Singular, MatricesForHomalg
+
+julia> R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+(Singular polynomial ring (QQ),(x,y),(dp(2),C), spoly{Singular.n_Q}[x, y])
+
+julia> A = Singular.Matrix(R, [x y; x^2 y^2; x*y x*y])
+x*y
+x^2*y^2
+x*y*x*y
+
+julia> S = SyzygiesOfRows(A)
+x*(-y+x)*1
+(-1)*(y+x)*x
+
+julia> iszero(S * A)
+true
+```
+"""
+function MatricesForHomalg.SyzygiesOfRows(A::Singular.smatrix)
+    R = Singular.base_ring(A)
+    # Columns of transpose(A) are the rows of A.
+    # Module from those columns, syz finds relations among them,
+    # i.e. vectors (a1,...,am) with a1*row1 + ... + am*rowm = 0.
+    M = Singular.Module(Singular.transpose(A))
+    S = Singular.syz(M)
+    # Matrix(S) has syzygy vectors as columns; transpose to get rows.
+    return Singular.transpose(Singular.Matrix(S))
+end
+
+"""
+    SyzygiesOfColumns(A::Singular.smatrix)
+
+Compute the matrix of column syzygies of A using Singular.
+Returns a matrix whose columns span the right kernel of A,
+i.e. all columns X satisfying A * X = 0.
+
+# Examples
+```jldoctest
+julia> using Singular, MatricesForHomalg
+
+julia> R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+(Singular polynomial ring (QQ),(x,y),(dp(2),C), spoly{Singular.n_Q}[x, y])
+
+julia> A = Singular.Matrix(R, [x x^2 x*y; y y^2 x*y])
+x, x^2, x*y
+y, y^2, x*y
+
+julia> S = SyzygiesOfColumns(A)
+...
+
+julia> iszero(A * S)
+true
+```
+"""
+function MatricesForHomalg.SyzygiesOfColumns(A::Singular.smatrix)
+    # Columns of A are generators of a submodule of R^m.
+    # syz finds vectors (a1,...,an) with a1*col1 + ... + an*coln = 0,
+    # which is exactly A * [a1,...,an]^T = 0.
+    M = Singular.Module(A)
+    S = Singular.syz(M)
+    if Singular.iszero(S)
+        # syz returns zero module if there are no syzygies, but we want a zero matrix.
+        R = Singular.base_ring(A)
+        return Singular.zero_matrix(R, Singular.ncols(A), 0)
+    else
+        return Singular.Matrix(S)
+    end
+end
+
+"""
+    SyzygiesOfRows(A::Singular.smatrix, N::Singular.smatrix)
+
+Compute the matrix of relative row syzygies of A modulo N using Singular.
+Returns a matrix whose rows K satisfy K * A + L * N = 0 for some L.
+
+# Examples
+```jldoctest
+julia> using Singular, MatricesForHomalg
+
+julia> R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+(Singular polynomial ring (QQ),(x,y),(dp(2),C), spoly{Singular.n_Q}[x, y])
+
+julia> v1 = Singular.vector(R, x, y)
+x*gen(1)+y*gen(2)
+
+julia> v2 = Singular.vector(R, y, x)
+y*gen(1)+x*gen(2)
+
+julia> A = Singular.Matrix(Singular.Module(R, v1, v2))
+x*y
+y*x
+
+julia> N = Singular.Matrix(Singular.Module(R, v1))
+x
+y
+
+julia> K = SyzygiesOfRows(A, N)
+...
+
+julia> # verify K * A is in the row span of N
+```
+"""
+function MatricesForHomalg.SyzygiesOfRows(A::Singular.smatrix, N::Singular.smatrix)
+    At = Singular.transpose(A)
+    Nt = Singular.transpose(N)
+    St = MatricesForHomalg.SyzygiesOfColumns(At, Nt)
+    return Singular.transpose(St)
+end
+
+"""
+    SyzygiesOfColumns(A::Singular.smatrix, N::Singular.smatrix)
+
+Compute the matrix of relative column syzygies of A modulo N using Singular.
+Returns a matrix whose columns K satisfy A * K + N * L = 0 for some L.
+
+# Examples
+```jldoctest
+julia> using Singular, MatricesForHomalg
+
+julia> R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+(Singular polynomial ring (QQ),(x,y),(dp(2),C), spoly{Singular.n_Q}[x, y])
+
+julia> v1 = Singular.vector(R, x, y)
+x*gen(1)+y*gen(2)
+
+julia> v2 = Singular.vector(R, y, x)
+y*gen(1)+x*gen(2)
+
+julia> A = Singular.Matrix(Singular.Module(R, v1, v2))
+x*y
+y*x
+
+julia> N = Singular.Matrix(Singular.Module(R, v1))
+x
+y
+
+julia> K = SyzygiesOfColumns(A, N)
+...
+
+julia> # verify A * K is in the column span of N
+```
+"""
+function MatricesForHomalg.SyzygiesOfColumns(A::Singular.smatrix, N::Singular.smatrix)
+    # modulo(M_A, M_N) computes the kernel of R^n -> M_A / (M_A ∩ M_N),
+    # i.e. vectors (a1,...,an) such that a1*col1(A) + ... + an*coln(A) ∈ Im(N).
+    M_A = Singular.Module(A)
+    M_N = Singular.Module(N)
+    S = Singular.modulo(M_A, M_N)
+    if Singular.iszero(S)
+        R = Singular.base_ring(A)
+        return Singular.zero_matrix(R, Singular.ncols(A), 0)
+    else
+        return Singular.Matrix(S)
+    end
+end
+
+## BasisOfRows / BasisOfColumns
+
+function MatricesForHomalg.BasisOfColumns(A::Singular.smatrix)
+    M = Singular.Module(A)
+    G = Singular.std(M)
+    return Singular.Matrix(G)
+end
+
+function MatricesForHomalg.BasisOfRows(A::Singular.smatrix)
+    return Singular.transpose(MatricesForHomalg.BasisOfColumns(Singular.transpose(A)))
+end
+
+## DecideZeroRows / DecideZeroColumns
+
+function MatricesForHomalg.DecideZeroColumns(B::Singular.smatrix, A::Singular.smatrix)
+    M_A = Singular.Module(A)
+    G = Singular.std(M_A)
+    M_B = Singular.Module(B)
+    R = Singular.reduce(M_B, G)
+    return Singular.Matrix(R)
+end
+
+function MatricesForHomalg.DecideZeroRows(B::Singular.smatrix, A::Singular.smatrix)
+    return Singular.transpose(MatricesForHomalg.DecideZeroColumns(Singular.transpose(B), Singular.transpose(A)))
+end
+
+## LeftDivide / RightDivide (two-argument: solve AX = B)
+
+function MatricesForHomalg.SafeLeftDivide(A::Singular.smatrix, B::Singular.smatrix)
+    M_A = Singular.Module(A)
+    M_B = Singular.Module(B)
+    T, rest = Singular.lift(M_A, M_B)
+    if !Singular.iszero(rest)
+        error("Unable to solve linear system")
+    end
+    return Singular.Matrix(T)
+end
+
+function MatricesForHomalg.LeftDivide(A::Singular.smatrix, B::Singular.smatrix)
+    try
+        return MatricesForHomalg.SafeLeftDivide(A, B)
+    catch
+        return "fail"
+    end
+end
+
+function MatricesForHomalg.SafeRightDivide(B::Singular.smatrix, A::Singular.smatrix)
+    return Singular.transpose(MatricesForHomalg.SafeLeftDivide(Singular.transpose(A), Singular.transpose(B)))
+end
+
+function MatricesForHomalg.RightDivide(B::Singular.smatrix, A::Singular.smatrix)
+    try
+        return MatricesForHomalg.SafeRightDivide(B, A)
+    catch
+        return "fail"
+    end
+end
+
+## LeftDivide / RightDivide (three-argument: solve AX + LY = B)
+
+function MatricesForHomalg.SafeLeftDivide(A::Singular.smatrix, B::Singular.smatrix, L::Singular.smatrix)
+    R = Singular.base_ring(A)
+    nr_rows = Singular.nrows(A)
+    nr_cols_a = Singular.ncols(A)
+    nr_cols_l = Singular.ncols(L)
+    # Combine columns of A and L into one matrix
+    AL = Singular.zero_matrix(R, nr_rows, nr_cols_a + nr_cols_l)
+    for i in 1:nr_rows, j in 1:nr_cols_a
+        AL[i, j] = A[i, j]
+    end
+    for i in 1:nr_rows, j in 1:nr_cols_l
+        AL[i, nr_cols_a + j] = L[i, j]
+    end
+    M_AL = Singular.Module(AL)
+    M_B = Singular.Module(B)
+    T, rest = Singular.lift(M_AL, M_B)
+    if !Singular.iszero(rest)
+        error("Unable to solve linear system")
+    end
+    T_mat = Singular.Matrix(T)
+    # Extract the first nr_cols_a rows (corresponding to A)
+    result = Singular.zero_matrix(R, nr_cols_a, Singular.ncols(T_mat))
+    for i in 1:nr_cols_a, j in 1:Singular.ncols(T_mat)
+        result[i, j] = T_mat[i, j]
+    end
+    return result
+end
+
+function MatricesForHomalg.LeftDivide(A::Singular.smatrix, B::Singular.smatrix, L::Singular.smatrix)
+    try
+        return MatricesForHomalg.SafeLeftDivide(A, B, L)
+    catch
+        return "fail"
+    end
+end
+
+function MatricesForHomalg.SafeRightDivide(B::Singular.smatrix, A::Singular.smatrix, L::Singular.smatrix)
+    return Singular.transpose(MatricesForHomalg.SafeLeftDivide(Singular.transpose(A), Singular.transpose(B), Singular.transpose(L)))
+end
+
+function MatricesForHomalg.RightDivide(B::Singular.smatrix, A::Singular.smatrix, L::Singular.smatrix)
+    try
+        return MatricesForHomalg.SafeRightDivide(B, A, L)
+    catch
+        return "fail"
+    end
+end
+
+end # module

test/runtests.jl

@@ -9,3 +9,4 @@ include("properties.jl")
 include("attributes.jl")
 include("operations.jl")
 include("testmanual.jl")
+include("singular-ext-test.jl")