diff --git a/Project.toml b/Project.toml
index cf8fbcf..c06cfe7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,18 +1,26 @@
 name = "MatricesForHomalg"
 uuid = "29b9b1b6-efa6-450e-8188-a5a2c25df071"
 authors = ["Mohamed Barakat <mohamed.barakat@uni-siegen.de>", "Johanna Knecht <johanna.knecht@student.uni-siegen.de>"]
-version = "0.1.6"
+version = "0.1.7"
 
 [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"
+Nemo = "0.48.0, 0.49, 0.50, 0.51, 0.52, 0.53, 0.54"
+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"]
diff --git a/ext/MatricesForHomalgSingularExt.jl b/ext/MatricesForHomalgSingularExt.jl
new file mode 100644
index 0000000..c960ec9
--- /dev/null
+++ b/ext/MatricesForHomalgSingularExt.jl
@@ -0,0 +1,629 @@
+module MatricesForHomalgSingularExt
+
+using MatricesForHomalg
+import Nemo
+import Singular
+
+MatricesForHomalg.HomalgRingOfIntegersInSingular() = Singular.ZZ
+MatricesForHomalg.HomalgFieldOfRationalsInSingular() = Singular.QQ
+
+MatricesForHomalg.HomalgRingOfIntegersInSingular(n::Int) = Singular.residue_ring(Singular.ZZ, n)[1]
+
+function Base.getindex(R::Union{Singular.Integers, Singular.Rationals}, s::Union{Char, AbstractString, Symbol})
+    return Singular.polynomial_ring(R, [string(s)])[1]
+end
+
+function Base.getindex(R::Union{Singular.Integers, Singular.Rationals}, s::Union{Char, AbstractString, Symbol}, ss::Union{Char, AbstractString, Symbol}...)
+    return Singular.polynomial_ring(R, [string(s); [string(x) for x in ss]])[1]
+end
+
+## assign generating variables
+function MatricesForHomalg.AssignGeneratingVariables(R::Singular.PolyRing)
+    for (sym, gen) in zip(R.S, Singular.gens(R))
+        Core.eval(Main, Expr(:(=), sym, gen))
+    end
+end
+
+function MatricesForHomalg.RingName(Q::Singular.Rationals)
+    return "Q"
+end
+
+function MatricesForHomalg.RingName(Z::Singular.Integers)
+    return "Z"
+end
+
+function MatricesForHomalg.RingName(R::Singular.PolyRing)
+    base = RingName(Singular.base_ring(R))
+    vars = join(map(s -> string(s), R.S), ",")
+    name = base * "[" * vars * "]"
+    if Singular.is_quotient_ring(R)
+        I = Singular.quotient_ideal(R)
+        rels = join([string(Singular.gens(I)[k]) for k in 1:Singular.ngens(I)], ", ")
+        name *= " / ( " * rels * " )"
+    end
+    return name
+end
+
+# Every commutative ring has the invariant basis property
+function MatricesForHomalg.HasHasInvariantBasisProperty(::Singular.PolyRing)
+    return true
+end
+
+function MatricesForHomalg.HasInvariantBasisProperty(::Singular.PolyRing)
+    return true
+end
+
+"""
+    R / relations
+
+Construct the quotient ring of R by the ideal generated by relations.
+Returns only the ring; the indeterminates are assigned as global variables
+with the same names as in R.
+
+# Examples
+```jldoctest
+julia> using Singular, MatricesForHomalg
+julia> R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+julia> S = R / [x^2 + y^2 - 1]
+```
+"""
+function Base.:/(R::Singular.PolyRing, relations::Vector)
+    I = Singular.Ideal(R, relations...)
+    G = Singular.std(I)
+    S, indets = Singular.QuotientRing(R, G)
+    #for (sym, val) in zip(R.S, indets)
+    #    Core.eval(Main, Expr(:(=), sym, val))
+    #end
+    return S
+end
+
+function MatricesForHomalg.Indeterminates(R::Singular.PolyRing)
+    return Singular.gens(R)
+end
+
+## Parse a polynomial string into an element of the ring R
+
+function _subst_vars(expr, vars::Dict{Symbol})
+    expr isa Symbol && haskey(vars, expr) ? vars[expr] :
+    expr isa Expr ? Expr(expr.head, (_subst_vars(a, vars) for a in expr.args)...) :
+    expr
+end
+
+function (R::Singular.PolyRing)(s::String)
+    vars = Dict(sym => g for (sym, g) in zip(R.S, Singular.gens(R)))
+    return Core.eval(@__MODULE__, _subst_vars(Meta.parse(s), vars))
+end
+
+#function (S::Singular.PolyRing)(f::Singular.spoly)
+#    parent(f) == S && return f
+#    B = Singular.MPolyBuildCtx(S)
+#    for (c, e) in zip(Singular.coefficients(f), Singular.exponent_vectors(f))
+#        Singular.push_term!(B, Singular.base_ring(S)(c), e)
+#    end
+#    return Singular.finish(B)
+#end
+
+"""
+    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 = Singular.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
+
+## HomalgZeroMatrix / HomalgIdentityMatrix for Singular rings
+
+function MatricesForHomalg.HomalgZeroMatrix(r, c, R::Singular.PolyRing)
+    return Singular.zero_matrix(R, r, c)
+end
+
+function MatricesForHomalg.HomalgIdentityMatrix(r, R::Singular.PolyRing)
+    return Singular.identity_matrix(R, Int(r))
+end
+
+## RandomMatrix for Singular rings
+function MatricesForHomalg.RandomMatrix(r, c, R::Singular.PolyRing)
+    m = Singular.zero_matrix(R, r, c)
+    for i in 1:r, j in 1:c
+        m[i, j] = R(rand(-5:5)) * prod([rand([zero(R); one(R); Singular.gens(R)...]) for i in 1:rand([1, 1, 1, 2, 2, 3])])
+    end
+    return m
+end
+
+function Base.:-(mat::Singular.smatrix)
+    R = Singular.base_ring(mat)
+    result = Singular.zero_matrix(R, Singular.nrows(mat), Singular.ncols(mat))
+    for i in 1:Singular.nrows(mat), j in 1:Singular.ncols(mat)
+        result[i, j] = -mat[i, j]
+    end
+    return result
+end
+
+## NumberRows / NumberColumns / HomalgRing / TransposedMatrix for smatrix
+
+function MatricesForHomalg.NumberRows(mat::Singular.smatrix)
+    return Singular.nrows(mat)
+end
+
+function MatricesForHomalg.NumberColumns(mat::Singular.smatrix)
+    return Singular.ncols(mat)
+end
+
+function MatricesForHomalg.HomalgRing(mat::Singular.smatrix)
+    return Singular.base_ring(mat)
+end
+
+function MatricesForHomalg.TransposedMatrix(mat::Singular.smatrix)
+    return Singular.transpose(mat)
+end
+
+## IsZero / IsOne / IsEmptyMatrix / IsSymmetricMatrix for smatrix
+
+function MatricesForHomalg.IsZero(mat::Singular.smatrix)
+    return Singular.iszero(mat)
+end
+
+function MatricesForHomalg.IsOne(mat::Singular.smatrix)
+    R = Singular.base_ring(mat)
+    nr = Singular.nrows(mat)
+    nc = Singular.ncols(mat)
+    nr == nc || return false
+    for i in 1:nr, j in 1:nc
+        if i == j
+            mat[i, j] == R(1) || return false
+        else
+            Singular.iszero(mat[i, j]) || return false
+        end
+    end
+    return true
+end
+
+function MatricesForHomalg.IsEmptyMatrix(mat::Singular.smatrix)
+    return Singular.nrows(mat) == 0 || Singular.ncols(mat) == 0
+end
+
+function MatricesForHomalg.IsSymmetricMatrix(mat::Singular.smatrix)
+    nr = Singular.nrows(mat)
+    nc = Singular.ncols(mat)
+    nr == nc || return false
+    for i in 1:nr, j in (i+1):nc
+        mat[i, j] == mat[j, i] || return false
+    end
+    return true
+end
+
+## CertainRows / CertainColumns for smatrix
+
+function MatricesForHomalg.CertainRows(mat::Singular.smatrix, list)
+    R = Singular.base_ring(mat)
+    nc = Singular.ncols(mat)
+    nr_new = length(list)
+    if nr_new == 0
+        return Singular.zero_matrix(R, 0, nc)
+    end
+    result = Singular.zero_matrix(R, nr_new, nc)
+    for (i, row) in enumerate(list)
+        for j in 1:nc
+            result[i, j] = mat[Int(row), j]
+        end
+    end
+    return result
+end
+
+function MatricesForHomalg.CertainColumns(mat::Singular.smatrix, list)
+    R = Singular.base_ring(mat)
+    nr = Singular.nrows(mat)
+    nc_new = length(list)
+    if nc_new == 0
+        return Singular.zero_matrix(R, nr, 0)
+    end
+    result = Singular.zero_matrix(R, nr, nc_new)
+    for (j, col) in enumerate(list)
+        for i in 1:nr
+            result[i, j] = mat[i, Int(col)]
+        end
+    end
+    return result
+end
+
+## UnionOfRows / UnionOfColumns for Singular rings
+
+function MatricesForHomalg.UnionOfRows(R::Singular.PolyRing, nr_cols, list)
+    if length(list) == 0
+        return Singular.zero_matrix(R, 0, nr_cols)
+    end
+    total_rows = sum(Singular.nrows(m) for m in list)
+    result = Singular.zero_matrix(R, total_rows, nr_cols)
+    row_offset = 0
+    for m in list
+        nr = Singular.nrows(m)
+        for i in 1:nr, j in 1:nr_cols
+            result[row_offset + i, j] = m[i, j]
+        end
+        row_offset += nr
+    end
+    return result
+end
+
+function MatricesForHomalg.UnionOfColumns(R::Singular.PolyRing, nr_rows, list)
+    if length(list) == 0
+        return Singular.zero_matrix(R, nr_rows, 0)
+    end
+    total_cols = sum(Singular.ncols(m) for m in list)
+    result = Singular.zero_matrix(R, nr_rows, total_cols)
+    col_offset = 0
+    for m in list
+        nc = Singular.ncols(m)
+        for i in 1:nr_rows, j in 1:nc
+            result[i, col_offset + j] = m[i, j]
+        end
+        col_offset += nc
+    end
+    return result
+end
+
+## KroneckerMat for smatrix
+
+function MatricesForHomalg.KroneckerMat(mat1::Singular.smatrix, mat2::Singular.smatrix)
+    R = Singular.base_ring(mat1)
+    r1, c1 = Singular.nrows(mat1), Singular.ncols(mat1)
+    r2, c2 = Singular.nrows(mat2), Singular.ncols(mat2)
+    result = Singular.zero_matrix(R, r1 * r2, c1 * c2)
+    for i1 in 1:r1, j1 in 1:c1
+        a = mat1[i1, j1]
+        for i2 in 1:r2, j2 in 1:c2
+            result[(i1 - 1) * r2 + i2, (j1 - 1) * c2 + j2] = a * mat2[i2, j2]
+        end
+    end
+    return result
+end
+
+## HomalgDiagonalMatrix for Singular rings
+
+function MatricesForHomalg.HomalgDiagonalMatrix(diagonal_entries, R::Singular.PolyRing)
+    n = length(diagonal_entries)
+    result = Singular.zero_matrix(R, n, n)
+    for i in 1:n
+        result[i, i] = R(diagonal_entries[i])
+    end
+    return result
+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.
+    if Singular.iszero(S)
+        # syz returns zero module if there are no syzygies, but we want a zero matrix.
+        return Singular.zero_matrix(R, 0, Singular.nrows(A))
+    else
+        return Singular.transpose(Singular.Matrix(S))
+    end
+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)
+    if Singular.iszero(G)
+        R = Singular.base_ring(A)
+        return Singular.zero_matrix(R, Singular.nrows(A), 0)
+    else
+        return Singular.Matrix(G)
+    end
+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
+
+## ReducedSyzygiesOfRows / ReducedSyzygiesOfColumns
+
+function MatricesForHomalg.ReducedSyzygiesOfColumns(A::Singular.smatrix)
+    M = Singular.Module(A)
+    S = Singular.syz(M)
+    if Singular.iszero(S)
+        R = Singular.base_ring(A)
+        return Singular.zero_matrix(R, Singular.ncols(A), 0)
+    end
+    G = Singular.std(S; complete_reduction=true)
+    return Singular.Matrix(G)
+end
+
+function MatricesForHomalg.ReducedSyzygiesOfRows(A::Singular.smatrix)
+    return Singular.transpose(MatricesForHomalg.ReducedSyzygiesOfColumns(Singular.transpose(A)))
+end
+
+function MatricesForHomalg.ReducedSyzygiesOfColumns(A::Singular.smatrix, N::Singular.smatrix)
+    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)
+    end
+    G = Singular.std(S; complete_reduction=true)
+    return Singular.Matrix(G)
+end
+
+function MatricesForHomalg.ReducedSyzygiesOfRows(A::Singular.smatrix, N::Singular.smatrix)
+    return Singular.transpose(MatricesForHomalg.ReducedSyzygiesOfColumns(Singular.transpose(A), Singular.transpose(N)))
+end
+
+end # module
diff --git a/makefile b/makefile
new file mode 100644
index 0000000..5bcf586
--- /dev/null
+++ b/makefile
@@ -0,0 +1,10 @@
+.PHONY: test
+
+install:
+	julia -e 'using Pkg; Pkg.develop(path=".");'
+
+uninstall:
+	julia -e 'using Pkg; Pkg.rm("MatricesForHomalg");'
+
+test:
+	julia -e 'using Pkg; Pkg.test("MatricesForHomalg");'
diff --git a/src/Declarations.jl b/src/Declarations.jl
new file mode 100644
index 0000000..e31695d
--- /dev/null
+++ b/src/Declarations.jl
@@ -0,0 +1,21 @@
+
+
+
+function HomalgRingOfIntegersInSingular end; export HomalgRingOfIntegersInSingular
+
+function HomalgFieldOfRationalsInSingular end; export HomalgFieldOfRationalsInSingular
+
+function AssignGeneratingVariables end; export AssignGeneratingVariables
+
+function HasHasInvariantBasisProperty end; export HasHasInvariantBasisProperty
+
+HasHasInvariantBasisProperty(::TypeOfRingForHomalg) = false
+HasHasInvariantBasisProperty(::Union{typeof(QQ), typeof(ZZ)}) = true
+
+function HasInvariantBasisProperty end; export HasInvariantBasisProperty
+
+HasInvariantBasisProperty(::Union{typeof(QQ), typeof(ZZ)}) = true;
+
+function RandomMatrix end; export RandomMatrix
+
+function Indeterminates end; export Indeterminates
diff --git a/src/MatricesForHomalg.jl b/src/MatricesForHomalg.jl
index 904d171..f578b20 100644
--- a/src/MatricesForHomalg.jl
+++ b/src/MatricesForHomalg.jl
@@ -18,6 +18,8 @@ TypeOfFieldForHomalg = Nemo.Field
 
 TypeOfRingElementForHomalg = Nemo.NCRingElement
 
+include("Declarations.jl")
+
 """
     HomalgRingOfIntegers()
 
@@ -46,6 +48,7 @@ function HomalgFieldOfRationals()::Nemo.QQField
     return Nemo.QQ
 end
 
+
 """
     RingName(ring)
 
@@ -77,7 +80,7 @@ export HomalgRingOfIntegers, HomalgFieldOfRationals, RingName
 
 ## Constructors of homalg matrices
 
-TypeOfMatrixForHomalg = Nemo.MatrixElem
+TypeOfMatrixForHomalg::Type = Union{Nemo.MatrixElem, Nemo.SetElem}
 
 """
     HomalgMatrix(L, r, c, R)
@@ -831,7 +834,8 @@ end
 
 export HomalgRing, NumberRows, NumberColumns, TransposedMatrix, ConvertMatrixToRow, ConvertMatrixToColumn,
     RowReducedEchelonForm, BasisOfRows, BasisOfColumns, ZeroRows, ZeroColumns, FirstZeroRow, FirstZeroColumn,
-    SyzygiesOfRows, SyzygiesOfColumns, RowRankOfMatrix, ColumnRankOfMatrix, StringDisplay
+    SyzygiesOfRows, SyzygiesOfColumns, RowRankOfMatrix, ColumnRankOfMatrix, StringDisplay,
+    ReducedSyzygiesOfRows, ReducedSyzygiesOfColumns
 
 ## Operations of homalg matrices
 
@@ -1651,6 +1655,76 @@ function SyzygiesOfColumns(A, N)
     return TransposedMatrix(SyzygiesOfRows(TransposedMatrix(A), TransposedMatrix(N)))
 end
 
+## ReducedSyzygiesOfRows / ReducedSyzygiesOfColumns
+
+"""
+    ReducedSyzygiesOfRows(M)
+
+Return the reduced row syzygies of M, i.e. a matrix S with
+reduced rows satisfying S * M = 0.
+
+```jldoctest
+julia> A = HomalgMatrix(4:9, 3, 2, ZZ)
+[4   5]
+[6   7]
+[8   9]
+
+julia> S = ReducedSyzygiesOfRows(A)
+[1   -2   1]
+
+julia> S * A
+[0   0]
+```
+"""
+function ReducedSyzygiesOfRows(M)
+    return BasisOfRows(SyzygiesOfRows(M))
+end
+
+"""
+    ReducedSyzygiesOfRows(M, M2)
+
+Return the reduced relative row syzygies of M modulo M2,
+i.e. a matrix K with reduced rows satisfying K * M + L * M2 = 0 for some L.
+"""
+function ReducedSyzygiesOfRows(M, M2)
+    return BasisOfRows(SyzygiesOfRows(M, M2))
+end
+
+"""
+    ReducedSyzygiesOfColumns(M)
+
+Return the reduced column syzygies of M, i.e. a matrix S with
+reduced columns satisfying M * S = 0.
+
+```jldoctest
+julia> A = TransposedMatrix(HomalgMatrix(4:9, 3, 2, ZZ))
+[4   6   8]
+[5   7   9]
+
+julia> S = ReducedSyzygiesOfColumns(A)
+[ 1]
+[-2]
+[ 1]
+
+julia> A * S
+[0]
+[0]
+```
+"""
+function ReducedSyzygiesOfColumns(M)
+    return TransposedMatrix(ReducedSyzygiesOfRows(TransposedMatrix(M)))
+end
+
+"""
+    ReducedSyzygiesOfColumns(M, M2)
+
+Return the reduced relative column syzygies of M modulo M2,
+i.e. a matrix K with reduced columns satisfying M * K + M2 * L = 0 for some L.
+"""
+function ReducedSyzygiesOfColumns(M, M2)
+    return TransposedMatrix(ReducedSyzygiesOfRows(TransposedMatrix(M), TransposedMatrix(M2)))
+end
+
 export UnionOfRows, UnionOfColumns, KroneckerMat, CertainColumns, CertainRows,
     SafeRightDivide, UniqueRightDivide, RightDivide, SafeLeftDivide, UniqueLeftDivide, LeftDivide,
     DecideZeroRows, DecideZeroColumns
diff --git a/test/runtests.jl b/test/runtests.jl
index 74e292e..4b91de0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,3 +9,4 @@ include("properties.jl")
 include("attributes.jl")
 include("operations.jl")
 include("testmanual.jl")
+include("singular-ext-test.jl")
diff --git a/test/singular-ext-test.jl b/test/singular-ext-test.jl
new file mode 100644
index 0000000..e32c248
--- /dev/null
+++ b/test/singular-ext-test.jl
@@ -0,0 +1,280 @@
+using Singular
+using MatricesForHomalg
+using Test
+
+@testset "SyzygiesOfRows for Singular matrices" begin
+    R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+
+    # A simple 3x2 matrix
+    A = Singular.zero_matrix(R, 3, 2)
+    A[1, 1] = x
+    A[1, 2] = y
+    A[2, 1] = x^2
+    A[2, 2] = y^2
+    A[3, 1] = x * y
+    A[3, 2] = x * y
+
+    S = SyzygiesOfRows(A)
+
+    # S * A must be zero
+    @test iszero(S * A)
+
+    # The syzygy matrix should have 3 columns (matching rows of A)
+    @test Singular.ncols(S) == 3
+end
+
+@testset "SyzygiesOfColumns for Singular matrices" begin
+    R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+
+    # A 2x3 matrix
+    A = Singular.zero_matrix(R, 2, 3)
+    A[1, 1] = x
+    A[1, 2] = x^2
+    A[1, 3] = x * y
+    A[2, 1] = y
+    A[2, 2] = y^2
+    A[2, 3] = x * y
+
+    S = SyzygiesOfColumns(A)
+
+    # A * S must be zero
+    @test iszero(A * S)
+
+    # The syzygy matrix should have 3 rows (matching columns of A)
+    @test Singular.nrows(S) == 3
+end
+
+@testset "SyzygiesOfRows - identity matrix has no syzygies" begin
+    R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+
+    I = Singular.identity_matrix(R, 3)
+    S = SyzygiesOfRows(I)
+
+    @test iszero(S)
+end
+
+@testset "SyzygiesOfColumns - identity matrix has no syzygies" begin
+    R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+
+    I = Singular.identity_matrix(R, 3)
+    S = SyzygiesOfColumns(I)
+
+    @test iszero(S)
+end
+
+@testset "Relative SyzygiesOfRows for Singular matrices" begin
+    R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+
+    A = Singular.zero_matrix(R, 2, 2)
+    A[1, 1] = x
+    A[1, 2] = y
+    A[2, 1] = x^2
+    A[2, 2] = y^2
+
+    N = Singular.zero_matrix(R, 1, 2)
+    N[1, 1] = x
+    N[1, 2] = y
+
+    K = SyzygiesOfRows(A, N)
+
+    # K * A should be in the row span of N,
+    # i.e., K * A + L * N = 0 for some L,
+    # equivalently: each row of K*A is a multiple of the row of N.
+    @test Singular.ncols(K) == 2  # same as nrows(A)
+end
+
+@testset "Relative SyzygiesOfColumns for Singular matrices" begin
+    R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+
+    A = Singular.zero_matrix(R, 2, 2)
+    A[1, 1] = x
+    A[1, 2] = x^2
+    A[2, 1] = y
+    A[2, 2] = y^2
+
+    N = Singular.zero_matrix(R, 2, 1)
+    N[1, 1] = x
+    N[2, 1] = y
+
+    K = SyzygiesOfColumns(A, N)
+
+    # A * K should be in the column span of N
+    @test Singular.nrows(K) == 2  # same as ncols(A)
+end
+
+@testset "SyzygiesOfRows via ideal syzygies" begin
+    R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+
+    # Row matrix (1x3) - equivalent to ideal syzygy
+    A = Singular.zero_matrix(R, 1, 3)
+    A[1, 1] = x^2 * y
+    A[1, 2] = x * y^2
+    A[1, 3] = x^2 * y^2
+
+    S = SyzygiesOfRows(A)
+
+    # the result should be trivial: a 0x1 matrix since 1 row has no left kernel
+    # (a single row vector's left kernel is trivially 0 or identity)
+    @test iszero(S * A)
+end
+
+@testset "Consistency: SyzygiesOfRows = transpose of SyzygiesOfColumns of transpose" begin
+    R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+
+    A = Singular.zero_matrix(R, 3, 2)
+    A[1, 1] = x
+    A[1, 2] = y
+    A[2, 1] = x^2
+    A[2, 2] = y^2
+    A[3, 1] = x * y
+    A[3, 2] = x * y
+
+    S_rows = SyzygiesOfRows(A)
+    S_cols = SyzygiesOfColumns(Singular.transpose(A))
+
+    # S_rows * A = 0 and transpose(A) * S_cols = 0
+    @test iszero(S_rows * A)
+    @test iszero(Singular.transpose(A) * S_cols)
+end
+
+@testset "BasisOfColumns for Singular matrices" begin
+    R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+
+    A = Singular.zero_matrix(R, 2, 3)
+    A[1, 1] = x;  A[1, 2] = x^2; A[1, 3] = x * y
+    A[2, 1] = y;  A[2, 2] = y^2; A[2, 3] = x * y
+    B = BasisOfColumns(A)
+
+    # B should have same number of rows as A
+    @test Singular.nrows(B) == 2
+    # Column rank should be at most ncols(A)
+    @test Singular.ncols(B) <= Singular.ncols(A)
+    # Each column of A should be reducible to zero modulo B
+    @test iszero(DecideZeroColumns(A, B))
+end
+
+@testset "BasisOfRows for Singular matrices" begin
+    R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+
+    A = Singular.zero_matrix(R, 3, 2)
+    A[1, 1] = x;  A[2, 1] = x^2; A[3, 1] = x * y
+    A[1, 2] = y;  A[2, 2] = y^2; A[3, 2] = x * y
+    B = BasisOfRows(A)
+
+    @test Singular.ncols(B) == 2
+    @test Singular.nrows(B) <= Singular.nrows(A)
+    @test iszero(DecideZeroRows(A, B))
+end
+
+@testset "DecideZeroColumns for Singular matrices" begin
+    R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+
+    A = Singular.zero_matrix(R, 2, 2)
+    A[1, 1] = x; A[1, 2] = y
+    A[2, 1] = y; A[2, 2] = x
+
+    # B is a column in the span of A
+    B = Singular.zero_matrix(R, 2, 1)
+    B[1, 1] = x + y; B[2, 1] = x + y  # = col1 + col2
+    @test iszero(DecideZeroColumns(B, A))
+
+    # Identity should reduce to zero modulo itself
+    I = Singular.identity_matrix(R, 2)
+    @test iszero(DecideZeroColumns(I, I))
+end
+
+@testset "DecideZeroRows for Singular matrices" begin
+    R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+
+    A = Singular.zero_matrix(R, 2, 2)
+    A[1, 1] = x; A[1, 2] = y
+    A[2, 1] = y; A[2, 2] = x
+
+    @test iszero(DecideZeroRows(A, A))
+end
+
+@testset "LeftDivide (2-arg) for Singular matrices" begin
+    R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+
+    A = Singular.zero_matrix(R, 2, 2)
+    A[1, 1] = x; A[1, 2] = y
+    A[2, 1] = y; A[2, 2] = x
+
+    B = Singular.zero_matrix(R, 2, 1)
+    B[1, 1] = x^2 + y^2; B[2, 1] = R(2) * x * y
+
+    X = LeftDivide(A, B)
+    @test X != "fail"
+    @test A * X == B
+end
+
+@testset "LeftDivide (2-arg) fail case" begin
+    R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+
+    A = Singular.zero_matrix(R, 2, 2)
+    A[1, 1] = x; A[1, 2] = y
+    A[2, 1] = y; A[2, 2] = x
+
+    B = Singular.zero_matrix(R, 2, 1)
+    B[1, 1] = x^3 + R(1); B[2, 1] = y^3 + R(1)
+
+    @test LeftDivide(A, B) == "fail"
+end
+
+@testset "RightDivide (2-arg) for Singular matrices" begin
+    R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+
+    A = Singular.zero_matrix(R, 2, 2)
+    A[1, 1] = x; A[1, 2] = y
+    A[2, 1] = y; A[2, 2] = x
+
+    B = Singular.zero_matrix(R, 1, 2)
+    B[1, 1] = x^2 + y^2; B[1, 2] = R(2) * x * y
+
+    X = RightDivide(B, A)
+    @test X != "fail"
+    @test X * A == B
+end
+
+@testset "LeftDivide (3-arg) for Singular matrices" begin
+    R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+
+    # A*X + L*Y = B with known solution
+    A = Singular.zero_matrix(R, 2, 2)
+    A[1, 1] = x; A[1, 2] = y
+    A[2, 1] = y; A[2, 2] = x
+
+    L = Singular.zero_matrix(R, 2, 1)
+    L[1, 1] = R(1); L[2, 1] = R(0)
+
+    # B = A*[x;y] + L*[1] = [x^2+y^2+1; 2xy]
+    B = Singular.zero_matrix(R, 2, 1)
+    B[1, 1] = x^2 + y^2 + R(1); B[2, 1] = R(2) * x * y
+
+    X = SafeLeftDivide(A, B, L)
+    @test Singular.nrows(X) == 2  # ncols(A) rows
+    @test Singular.ncols(X) == 1  # ncols(B) cols
+    # B - A*X should be in column span of L
+    @test iszero(DecideZeroColumns(B - A * X, L))
+end
+
+@testset "RightDivide (3-arg) for Singular matrices" begin
+    R, (x, y) = Singular.polynomial_ring(Singular.QQ, ["x", "y"])
+
+    A = Singular.zero_matrix(R, 2, 2)
+    A[1, 1] = x; A[1, 2] = y
+    A[2, 1] = y; A[2, 2] = x
+
+    L = Singular.zero_matrix(R, 1, 2)
+    L[1, 1] = R(1); L[1, 2] = R(0)
+
+    # B = [x;y]^T * A + [1] * L = [x^2+y^2+1, 2xy]
+    B = Singular.zero_matrix(R, 1, 2)
+    B[1, 1] = x^2 + y^2 + R(1); B[1, 2] = R(2) * x * y
+
+    X = SafeRightDivide(B, A, L)
+    @test Singular.nrows(X) == 1
+    @test Singular.ncols(X) == 2
+    # B - X*A should be in row span of L
+    @test iszero(DecideZeroRows(B - X * A, L))
+end