From d0efdf4d362a68c485d76198d53caeee8f303fa6 Mon Sep 17 00:00:00 2001 From: Max Horn Date: Sat, 21 Dec 2024 02:22:41 +0100 Subject: [PATCH 01/17] WIP: rewrite MatRing to wrap a native matrix --- src/MatRing.jl | 22 -------- src/Matrix.jl | 8 +-- src/generic/GenericTypes.jl | 15 ++---- src/generic/MatRing.jl | 101 ++++++++++++++---------------------- src/generic/Matrix.jl | 10 ++-- 5 files changed, 54 insertions(+), 102 deletions(-) diff --git a/src/MatRing.jl b/src/MatRing.jl index faf2feacb6..ed35143c52 100644 --- a/src/MatRing.jl +++ b/src/MatRing.jl @@ -180,28 +180,6 @@ function show(io::IO, a::MatRing) end end -############################################################################### -# -# Binary operations -# -############################################################################### - -function *(x::MatRingElem{T}, y::MatRingElem{T}) where {T <: NCRingElement} - degree(x) != degree(y) && error("Incompatible matrix degrees") - A = similar(x) - C = base_ring(x)() - for i = 1:nrows(x) - for j = 1:ncols(y) - A[i, j] = base_ring(x)() - for k = 1:ncols(x) - C = mul!(C, x[i, k], y[k, j]) - A[i, j] = add!(A[i, j], C) - end - end - end - return A -end - ############################################################################### # # Ad hoc comparisons diff --git a/src/Matrix.jl b/src/Matrix.jl index 4dc8bf9139..41c03803bc 100644 --- a/src/Matrix.jl +++ b/src/Matrix.jl @@ -823,7 +823,7 @@ end # ############################################################################### -function +(x::MatrixElem{T}, y::MatrixElem{T}) where {T <: NCRingElement} +function +(x::T, y::T) where {T <: MatrixElem} check_parent(x, y) r = similar(x) for i = 1:nrows(x) @@ -834,7 +834,7 @@ function +(x::MatrixElem{T}, y::MatrixElem{T}) where {T <: NCRingElement} return r end -function -(x::MatrixElem{T}, y::MatrixElem{T}) where {T <: NCRingElement} +function -(x::T, y::T) where {T <: MatrixElem} check_parent(x, y) r = similar(x) for i = 1:nrows(x) @@ -845,7 +845,7 @@ function -(x::MatrixElem{T}, y::MatrixElem{T}) where {T <: NCRingElement} return r end -function *(x::MatElem{T}, y::MatElem{T}) where {T <: NCRingElement} +function *(x::T, y::T) where {T <: MatrixElem} ncols(x) != nrows(y) && error("Incompatible matrix dimensions") base_ring(x) !== base_ring(y) && error("Base rings do not match") A = similar(x, nrows(x), ncols(y)) @@ -2628,6 +2628,8 @@ function minors_iterator(M::MatElem, k::Int) return (det(M[rows, cols]) for rows in row_indices for cols in col_indices) end +# TODO: how to get the 'positions' of the minors? + @doc raw""" minors_iterator_with_position(A::MatElem, k::Int) diff --git a/src/generic/GenericTypes.jl b/src/generic/GenericTypes.jl index bdb61fde40..02d2eeaa0a 100644 --- a/src/generic/GenericTypes.jl +++ b/src/generic/GenericTypes.jl @@ -1117,22 +1117,17 @@ struct MatRing{T <: NCRingElement} <: AbstractAlgebra.MatRing{T} end struct MatRingElem{T <: NCRingElement} <: AbstractAlgebra.MatRingElem{T} - base_ring::NCRing - entries::Matrix{T} + data::MatElem{T} - function MatRingElem{T}(R::NCRing, A::Matrix{T}) where T <: NCRingElement - @assert elem_type(R) === T - return new{T}(R, A) + function MatRingElem{T}(A::MatElem{T}) where T <: NCRingElement + return new{T}(A) end end function MatRingElem{T}(R::NCRing, n::Int, A::Vector{T}) where T <: NCRingElement @assert elem_type(R) === T - t = Matrix{T}(undef, n, n) - for i = 1:n, j = 1:n - t[i, j] = A[(i - 1) * n + j] - end - return MatRingElem{T}(R, t) + t = matrix(R, n, n, A) + return MatRingElem{T}(t) end ############################################################################### diff --git a/src/generic/MatRing.jl b/src/generic/MatRing.jl index 9883e1ad61..09886388b0 100644 --- a/src/generic/MatRing.jl +++ b/src/generic/MatRing.jl @@ -16,17 +16,38 @@ elem_type(::Type{MatRing{T}}) where {T <: NCRingElement} = MatRingElem{T} base_ring(a::MatRing{T}) where {T <: NCRingElement} = a.base_ring::parent_type(T) +base_ring(a::MatRingElem{T}) where {T <: NCRingElement} = base_ring(a.data) + @doc raw""" parent(a::MatRingElem{T}) where T <: NCRingElement Return the parent object of the given matrix. """ -parent(a::MatRingElem{T}) where T <: NCRingElement = MatRing{T}(base_ring(a), size(a.entries)[1]) +parent(a::MatRingElem{T}) where T <: NCRingElement = MatRing{T}(base_ring(a), nrows(a.data)) is_exact_type(::Type{MatRingElem{T}}) where T <: NCRingElement = is_exact_type(T) is_domain_type(::Type{MatRingElem{T}}) where T <: NCRingElement = false +############################################################################### +# +# Basic manipulation +# +############################################################################### + +number_of_rows(a::MatRingElem) = nrows(a.data) + +number_of_columns(a::MatRingElem) = ncols(a.data) + +Base.@propagate_inbounds getindex(a::MatRingElem, r::Int, c::Int) = a.data[r, c] + +Base.@propagate_inbounds function setindex!(a::MatRingElem, d::NCRingElement, + r::Int, c::Int) + a.data[r, c] = base_ring(a)(d) +end + +Base.isassigned(a::MatRingElem, i, j) = isassigned(a.data, i, j) + ############################################################################### # # Transpose @@ -34,9 +55,7 @@ is_domain_type(::Type{MatRingElem{T}}) where T <: NCRingElement = false ############################################################################### function transpose(x::MatRingElem{T}) where T <: NCRingElement - arr = permutedims(x.entries) - z = MatRingElem{T}(base_ring(x), arr) - return z + return MatRingElem{T}(transpose(x.data)) end ############################################################################### @@ -48,18 +67,18 @@ end function _can_solve_with_solution_lu(M::MatRingElem{T}, B::MatRingElem{T}) where {T <: RingElement} check_parent(M, B) R = base_ring(M) - MS = MatSpaceElem{T}(R, M.entries) # convert to ordinary matrix - BS = MatSpaceElem{T}(R, B.entries) + MS = M.data # TODO: is this right? + BS = B.data # TODO: is this right? flag, S = _can_solve_with_solution_lu(MS, BS) - SA = MatRingElem{T}(R, S.entries) + SA = MatRingElem{T}(S) return flag, SA end function AbstractAlgebra.can_solve_with_solution(M::MatRingElem{T}, B::MatRingElem{T}) where {T <: RingElement} check_parent(M, B) R = base_ring(M) - MS = MatSpaceElem{T}(R, M.entries) # convert to ordinary matrix - BS = MatSpaceElem{T}(R, B.entries) + MS = M.data # TODO: is this right? + BS = B.data # TODO: is this right? flag, S = can_solve_with_solution(MS, BS) SA = MatRingElem{T}(R, S.entries) return flag, SA @@ -68,10 +87,10 @@ end function _can_solve_with_solution_fflu(M::MatRingElem{T}, B::MatRingElem{T}) where {T <: RingElement} check_parent(M, B) R = base_ring(M) - MS = MatSpaceElem{T}(R, M.entries) # convert to ordinary matrix - BS = MatSpaceElem{T}(R, B.entries) + MS = M.data # convert to ordinary matrix + BS = B.data flag, S, d = _can_solve_with_solution_fflu(MS, BS) - SA = MatRingElem{T}(R, S.entries) + SA = MatRingElem{T}(S) return flag, SA, d end @@ -88,8 +107,7 @@ Return the minimal polynomial $p$ of the matrix $M$. The polynomial ring $S$ of the resulting polynomial must be supplied and the matrix must be square. """ function minpoly(S::Ring, M::MatRingElem{T}, charpoly_only::Bool = false) where {T <: RingElement} - MS = MatSpaceElem{T}(base_ring(M), M.entries) # convert to ordinary matrix - return minpoly(S, MS, charpoly_only) + return minpoly(S, M.data, charpoly_only) end function minpoly(M::MatRingElem{T}, charpoly_only::Bool = false) where {T <: RingElement} @@ -105,7 +123,7 @@ end ############################################################################### function add!(A::MatRingElem{T}, B::MatRingElem{T}) where T <: NCRingElement - A.entries .+= B.entries + A.data = add!(A.data, B.data) return A end @@ -129,48 +147,20 @@ end function (a::MatRing{T})() where {T <: NCRingElement} R = base_ring(a) - entries = Matrix{T}(undef, a.n, a.n) - for i = 1:a.n - for j = 1:a.n - entries[i, j] = zero(R) - end - end - z = MatRingElem{T}(R, entries) + z = MatRingElem{T}(zero_matrix(R, a.n, a.n)) return z end function (a::MatRing{T})(b::S) where {S <: NCRingElement, T <: NCRingElement} R = base_ring(a) - entries = Matrix{T}(undef, a.n, a.n) - rb = R(b) - for i = 1:a.n - for j = 1:a.n - if i != j - entries[i, j] = zero(R) - else - entries[i, j] = rb - end - end - end - z = MatRingElem{T}(R, entries) + z = MatRingElem{T}(scalar_matrix(R, a.n, R(b))) return z end # to resolve ambiguity for MatRing{MatRing{...}} function (a::MatRing{T})(b::T) where {S <: NCRingElement, T <: MatRingElem{S}} R = base_ring(a) - entries = Matrix{T}(undef, a.n, a.n) - rb = R(b) - for i = 1:a.n - for j = 1:a.n - if i != j - entries[i, j] = zero(R) - else - entries[i, j] = rb - end - end - end - z = MatRingElem{T}(R, entries) + z = MatRingElem{T}(scalar_matrix(R, a.n, R(b))) return z end @@ -182,33 +172,20 @@ end function (a::MatRing{T})(b::MatrixElem{S}) where {S <: NCRingElement, T <: NCRingElement} R = base_ring(a) _check_dim(nrows(a), ncols(a), b) - entries = Matrix{T}(undef, nrows(a), ncols(a)) - for i = 1:nrows(a) - for j = 1:ncols(a) - entries[i, j] = R(b[i, j]) - end - end - z = MatRingElem{T}(R, entries) + z = MatRingElem{T}(matrix(R, b)) return z end function (a::MatRing{T})(b::Matrix{S}) where {S <: NCRingElement, T <: NCRingElement} R = base_ring(a) _check_dim(a.n, a.n, b) - entries = Matrix{T}(undef, a.n, a.n) - for i = 1:a.n - for j = 1:a.n - entries[i, j] = R(b[i, j]) - end - end - z = MatRingElem{T}(R, entries) + z = MatRingElem{T}(matrix(R, b)) return z end function (a::MatRing{T})(b::Vector{S}) where {S <: NCRingElement, T <: NCRingElement} _check_dim(a.n, a.n, b) - b = Matrix{S}(transpose(reshape(b, a.n, a.n))) - z = a(b) + z = MatRingElem{T}(matrix(R, a.n, a.n, b)) return z end diff --git a/src/generic/Matrix.jl b/src/generic/Matrix.jl index db4089a490..0ca22035f2 100644 --- a/src/generic/Matrix.jl +++ b/src/generic/Matrix.jl @@ -43,18 +43,18 @@ dense_matrix_type(::Type{T}) where T <: NCRingElement = MatSpaceElem{T} # ############################################################################### -number_of_rows(a::Union{Mat, MatRingElem}) = size(a.entries, 1) +number_of_rows(a::Mat) = size(a.entries, 1) -number_of_columns(a::Union{Mat,MatRingElem}) = size(a.entries, 2) +number_of_columns(a::Mat) = size(a.entries, 2) -Base.@propagate_inbounds getindex(a::Union{Mat, MatRingElem}, r::Int, c::Int) = a.entries[r, c] +Base.@propagate_inbounds getindex(a::Mat, r::Int, c::Int) = a.entries[r, c] -Base.@propagate_inbounds function setindex!(a::Union{Mat, MatRingElem}, d::NCRingElement, +Base.@propagate_inbounds function setindex!(a::Mat, d::NCRingElement, r::Int, c::Int) a.entries[r, c] = base_ring(a)(d) end -Base.isassigned(a::Union{Mat,MatRingElem}, i, j) = isassigned(a.entries, i, j) +Base.isassigned(a::Mat, i, j) = isassigned(a.entries, i, j) ################################################################################ # From 06c5bab01bd4a4004f9c04411e0ace4b838891a3 Mon Sep 17 00:00:00 2001 From: John Abbott Date: Mon, 8 Dec 2025 13:16:28 +0100 Subject: [PATCH 02/17] fixed code so that tests pass --- src/MatRing.jl | 47 ++++++++++++++++++++++++++-------- src/Matrix.jl | 10 ++++---- src/NCPoly.jl | 3 +++ src/generic/GenericTypes.jl | 11 +++++--- src/generic/MatRing.jl | 49 ++++++++++++++++++------------------ test/generic/MatRing-test.jl | 2 +- 6 files changed, 77 insertions(+), 45 deletions(-) diff --git a/src/MatRing.jl b/src/MatRing.jl index ed35143c52..eadcdbe33a 100644 --- a/src/MatRing.jl +++ b/src/MatRing.jl @@ -89,9 +89,10 @@ Create an uninitialized matrix ring element over the given ring and dimension, with defaults based upon the given source matrix ring element `x`. """ function similar(x::MatRingElem, R::NCRing=base_ring(x), n::Int=degree(x)) + @assert (n >= 0) && (n < 2^30) # so that n^2 cannot overflow TT = elem_type(R) - M = Matrix{TT}(undef, (n, n)) - return Generic.MatRingElem{TT}(R, M) + M = fill(zero(R), n^2) + return Generic.MatRingElem(R, n, M) end similar(x::MatRingElem, n::Int) = similar(x, base_ring(x), n) @@ -224,12 +225,33 @@ end ==(x::T, y::MatRingElem{T}) where T <: NCRingElem = y == x +############################################################################### +# +# Basic arithmetic -- delegate to MatElem +# +############################################################################### + +*(x::MatRingElem{T}, y::MatRingElem{T}) where {T <: NCRingElement} = Generic.MatRingElem(matrix(x) * matrix(y)) + ++(x::MatRingElem{T}, y::MatRingElem{T}) where {T <: NCRingElement} = Generic.MatRingElem(matrix(x) + matrix(y)) + +-(x::MatRingElem{T}, y::MatRingElem{T}) where {T <: NCRingElement} = Generic.MatRingElem(matrix(x) - matrix(y)) + +==(x::MatRingElem{T}, y::MatRingElem{T}) where {T <: NCRingElement} = matrix(x) == matrix(y) + + ############################################################################### # # Exact division # ############################################################################### +# TO DO: using pseudo_inv is not ideal +# consider case M = matrix(ZZmod4, 2,2, [2,1,0,1]) +# pseudo_inv(M) gives error, but copuld give matrix(ZZ4, 2,2, [1,-1,0,2]) with denom=2 +# HINT: Consider using solve(f,g;side=:right) or side=:left +# The unused kwargs in the field cases are necessary for generic code to compile. + function divexact_left(f::MatRingElem{T}, g::MatRingElem{T}; check::Bool=true) where T <: RingElement ginv, d = pseudo_inv(g) @@ -237,7 +259,7 @@ function divexact_left(f::MatRingElem{T}, end function divexact_right(f::MatRingElem{T}, - g::MatRingElem{T}; check::Bool=true) where T <: RingElement + g::MatRingElem{T}; check::Bool=true) where T <: RingElement ginv, d = pseudo_inv(g) return divexact(f*ginv, d; check=check) end @@ -297,6 +319,9 @@ function RandomExtensions.make(S::MatRing, vs...) end end +# Sampler for a MatRing not needing arguments (passed via make) +# this allows to obtain the Sampler in simple cases without having to know about make +# (when one can do `rand(M)`, one can expect to be able to do `rand(Sampler(rng, M))`) Random.Sampler(::Type{RNG}, S::MatRing, n::Random.Repetition ) where {RNG<:AbstractRNG} = Random.Sampler(RNG, make(S), n) @@ -402,14 +427,14 @@ end ############################################################################### function identity_matrix(M::MatRingElem{T}, n::Int) where T <: NCRingElement - R = base_ring(M) - arr = Matrix{T}(undef, n, n) - for i in 1:n - for j in 1:n - arr[i, j] = i == j ? one(R) : zero(R) - end - end - z = Generic.MatRingElem{T}(R, arr) + @assert (n >= 0) && (n < 2^30) # so that n^2 cannot overflow + R = base_ring(M) + # ??Simpler?? Generic.MatRingElem(matrix_space(R,n,n)(1)) + arr = fill(zero(R), n^2) + for i in 1:n + arr[i+(i-1)*n] = one(R) + end + z = Generic.MatRingElem(R, n, arr)::MatRingElem{T} return z end diff --git a/src/Matrix.jl b/src/Matrix.jl index 41c03803bc..6d0b12638f 100644 --- a/src/Matrix.jl +++ b/src/Matrix.jl @@ -27,7 +27,7 @@ base_ring(a::MatrixElem{T}) where {T <: NCRingElement} = a.base_ring::parent_typ is_zero_initialized(mat::T) where {T<:MatrixElem} Specify whether the default-constructed matrices of type `T`, via the -`T(R::Ring, ::UndefInitializerm r::Int, c::Int)` constructor, are +`T(R::Ring, ::UndefInitializer, r::Int, c::Int)` constructor, are zero-initialized. The default is `false`, and new matrix types should specialize this method to return `true` if suitable, to enable optimizations. """ @@ -823,7 +823,7 @@ end # ############################################################################### -function +(x::T, y::T) where {T <: MatrixElem} +function +(x::MatElem{T}, y::MatElem{T}) where {T} check_parent(x, y) r = similar(x) for i = 1:nrows(x) @@ -834,7 +834,7 @@ function +(x::T, y::T) where {T <: MatrixElem} return r end -function -(x::T, y::T) where {T <: MatrixElem} +function -(x::MatElem{T}, y::MatElem{T}) where {T} check_parent(x, y) r = similar(x) for i = 1:nrows(x) @@ -845,7 +845,7 @@ function -(x::T, y::T) where {T <: MatrixElem} return r end -function *(x::T, y::T) where {T <: MatrixElem} +function *(x::MatElem{T}, y::MatElem{T}) where {T} ncols(x) != nrows(y) && error("Incompatible matrix dimensions") base_ring(x) !== base_ring(y) && error("Base rings do not match") A = similar(x, nrows(x), ncols(y)) @@ -1197,7 +1197,7 @@ function Base.promote(x::MatrixElem{S}, T <: NCRingElement} U = promote_rule_sym(S, T) if U === S - return x, map(base_ring(x), y) + return x, map(base_ring(x), y)::Vector{S} # Julia needs help here elseif U === T && length(y) != 0 return change_base_ring(parent(y[1]), x), y else diff --git a/src/NCPoly.jl b/src/NCPoly.jl index a5a91b9d0e..4cfa447ca1 100644 --- a/src/NCPoly.jl +++ b/src/NCPoly.jl @@ -697,6 +697,8 @@ end # ############################################################################### +##?? Base.eltype(::Type{M}) where {M<:NCPolyRing} = elem_type(M) + RandomExtensions.maketype(S::NCPolyRing, dr::AbstractUnitRange{Int}, _) = elem_type(S) function RandomExtensions.make(S::NCPolyRing, deg_range::AbstractUnitRange{Int}, vs...) @@ -714,6 +716,7 @@ function rand(rng::AbstractRNG, <:AbstractUnitRange{Int}}}) S, deg_range, v = sp[][1:end] R = base_ring(S) + local f::elem_type(S) # Julia needs some guidance f = S() x = gen(S) for i = 0:rand(rng, deg_range) diff --git a/src/generic/GenericTypes.jl b/src/generic/GenericTypes.jl index 02d2eeaa0a..467da804e9 100644 --- a/src/generic/GenericTypes.jl +++ b/src/generic/GenericTypes.jl @@ -1119,17 +1119,20 @@ end struct MatRingElem{T <: NCRingElement} <: AbstractAlgebra.MatRingElem{T} data::MatElem{T} - function MatRingElem{T}(A::MatElem{T}) where T <: NCRingElement - return new{T}(A) + function MatRingElem(A::MatElem{TT}) where TT <: NCRingElement + nrows(A) == ncols(A) || error("Matrix must be square") + return new{TT}(A) end end -function MatRingElem{T}(R::NCRing, n::Int, A::Vector{T}) where T <: NCRingElement +function MatRingElem(R::NCRing, n::Int, A::Vector{T}) where T <: NCRingElement @assert elem_type(R) === T t = matrix(R, n, n, A) - return MatRingElem{T}(t) + return MatRingElem(t) end +matrix(A::MatRingElem) = A.data + ############################################################################### # # FreeAssociativeAlgebra / FreeAssociativeAlgebraElem diff --git a/src/generic/MatRing.jl b/src/generic/MatRing.jl index 09886388b0..2c017ac8d3 100644 --- a/src/generic/MatRing.jl +++ b/src/generic/MatRing.jl @@ -35,9 +35,9 @@ is_domain_type(::Type{MatRingElem{T}}) where T <: NCRingElement = false # ############################################################################### -number_of_rows(a::MatRingElem) = nrows(a.data) +number_of_rows(a::MatRingElem) = nrows(a.data)::Int # Julia needs guidance here -number_of_columns(a::MatRingElem) = ncols(a.data) +number_of_columns(a::MatRingElem) = ncols(a.data)::Int # Julia needs guidance here Base.@propagate_inbounds getindex(a::MatRingElem, r::Int, c::Int) = a.data[r, c] @@ -54,8 +54,8 @@ Base.isassigned(a::MatRingElem, i, j) = isassigned(a.data, i, j) # ############################################################################### -function transpose(x::MatRingElem{T}) where T <: NCRingElement - return MatRingElem{T}(transpose(x.data)) +function transpose(x::MatRingElem) + return MatRingElem(transpose(matrix(x))) end ############################################################################### @@ -67,30 +67,30 @@ end function _can_solve_with_solution_lu(M::MatRingElem{T}, B::MatRingElem{T}) where {T <: RingElement} check_parent(M, B) R = base_ring(M) - MS = M.data # TODO: is this right? - BS = B.data # TODO: is this right? + MS = matrix(M) + BS = matrix(B) flag, S = _can_solve_with_solution_lu(MS, BS) - SA = MatRingElem{T}(S) + SA = MatRingElem(S) return flag, SA end function AbstractAlgebra.can_solve_with_solution(M::MatRingElem{T}, B::MatRingElem{T}) where {T <: RingElement} check_parent(M, B) R = base_ring(M) - MS = M.data # TODO: is this right? - BS = B.data # TODO: is this right? + MS = matrix(M) + BS = matrix(B) flag, S = can_solve_with_solution(MS, BS) - SA = MatRingElem{T}(R, S.entries) + SA = MatRingElem(S) return flag, SA end function _can_solve_with_solution_fflu(M::MatRingElem{T}, B::MatRingElem{T}) where {T <: RingElement} check_parent(M, B) R = base_ring(M) - MS = M.data # convert to ordinary matrix - BS = B.data + MS = matrix(M) + BS = matrix(B) flag, S, d = _can_solve_with_solution_fflu(MS, BS) - SA = MatRingElem{T}(S) + SA = MatRingElem(S) return flag, SA, d end @@ -107,13 +107,13 @@ Return the minimal polynomial $p$ of the matrix $M$. The polynomial ring $S$ of the resulting polynomial must be supplied and the matrix must be square. """ function minpoly(S::Ring, M::MatRingElem{T}, charpoly_only::Bool = false) where {T <: RingElement} - return minpoly(S, M.data, charpoly_only) + return minpoly(S, matrix(M), charpoly_only) end function minpoly(M::MatRingElem{T}, charpoly_only::Bool = false) where {T <: RingElement} R = base_ring(M) Rx, x = polynomial_ring(R; cached=false) - return minpoly(Rx, M, charpoly_only) + return minpoly(Rx, matrix(M), charpoly_only) end ############################################################################### @@ -123,7 +123,7 @@ end ############################################################################### function add!(A::MatRingElem{T}, B::MatRingElem{T}) where T <: NCRingElement - A.data = add!(A.data, B.data) + #=A.data ==# add!(A.data, B.data) ### !!struct is NOT mutable!! return A end @@ -147,20 +147,20 @@ end function (a::MatRing{T})() where {T <: NCRingElement} R = base_ring(a) - z = MatRingElem{T}(zero_matrix(R, a.n, a.n)) + z = MatRingElem(zero_matrix(R, a.n, a.n)) return z end function (a::MatRing{T})(b::S) where {S <: NCRingElement, T <: NCRingElement} R = base_ring(a) - z = MatRingElem{T}(scalar_matrix(R, a.n, R(b))) + z = MatRingElem(scalar_matrix(R, a.n, R(b))) return z end # to resolve ambiguity for MatRing{MatRing{...}} function (a::MatRing{T})(b::T) where {S <: NCRingElement, T <: MatRingElem{S}} R = base_ring(a) - z = MatRingElem{T}(scalar_matrix(R, a.n, R(b))) + z = MatRingElem(scalar_matrix(R, a.n, R(b))) return z end @@ -172,20 +172,21 @@ end function (a::MatRing{T})(b::MatrixElem{S}) where {S <: NCRingElement, T <: NCRingElement} R = base_ring(a) _check_dim(nrows(a), ncols(a), b) - z = MatRingElem{T}(matrix(R, b)) + z = MatRingElem(matrix(R, b)) return z end function (a::MatRing{T})(b::Matrix{S}) where {S <: NCRingElement, T <: NCRingElement} - R = base_ring(a) _check_dim(a.n, a.n, b) - z = MatRingElem{T}(matrix(R, b)) + R = base_ring(a) + z = MatRingElem(matrix(R, b)) return z end function (a::MatRing{T})(b::Vector{S}) where {S <: NCRingElement, T <: NCRingElement} - _check_dim(a.n, a.n, b) - z = MatRingElem{T}(matrix(R, a.n, a.n, b)) + _check_dim(a.n, a.n, b) + R = base_ring(a) + z = MatRingElem(matrix(R, a.n, a.n, b)) return z end diff --git a/test/generic/MatRing-test.jl b/test/generic/MatRing-test.jl index 83ab00fd93..a047bb77c4 100644 --- a/test/generic/MatRing-test.jl +++ b/test/generic/MatRing-test.jl @@ -72,7 +72,7 @@ @test parent(S()) == S end -@testset "Generic.MatAlg.finitiess" begin +@testset "Generic.MatAlg.finiteness" begin S = matrix_ring(QQ, 3) @test !is_finite(S) @test !is_trivial(S) From 3bbe379812289788e345a6441c4c72ad6bc9059f Mon Sep 17 00:00:00 2001 From: JohnAAbbott <124266874+JohnAAbbott@users.noreply.github.com> Date: Mon, 8 Dec 2025 15:32:33 +0100 Subject: [PATCH 03/17] Update src/NCPoly.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Lars Göttgens --- src/NCPoly.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/NCPoly.jl b/src/NCPoly.jl index 4cfa447ca1..ce59dc9c2d 100644 --- a/src/NCPoly.jl +++ b/src/NCPoly.jl @@ -697,8 +697,6 @@ end # ############################################################################### -##?? Base.eltype(::Type{M}) where {M<:NCPolyRing} = elem_type(M) - RandomExtensions.maketype(S::NCPolyRing, dr::AbstractUnitRange{Int}, _) = elem_type(S) function RandomExtensions.make(S::NCPolyRing, deg_range::AbstractUnitRange{Int}, vs...) From 1f978478ab8f68ad92b8989985d4f7bfb43aea7c Mon Sep 17 00:00:00 2001 From: JohnAAbbott <124266874+JohnAAbbott@users.noreply.github.com> Date: Mon, 8 Dec 2025 15:33:32 +0100 Subject: [PATCH 04/17] Update src/generic/MatRing.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Lars Göttgens --- src/generic/MatRing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/generic/MatRing.jl b/src/generic/MatRing.jl index 2c017ac8d3..0f20aec806 100644 --- a/src/generic/MatRing.jl +++ b/src/generic/MatRing.jl @@ -16,7 +16,7 @@ elem_type(::Type{MatRing{T}}) where {T <: NCRingElement} = MatRingElem{T} base_ring(a::MatRing{T}) where {T <: NCRingElement} = a.base_ring::parent_type(T) -base_ring(a::MatRingElem{T}) where {T <: NCRingElement} = base_ring(a.data) +base_ring(a::MatRingElem{T}) where {T <: NCRingElement} = base_ring(matrix(a)) @doc raw""" parent(a::MatRingElem{T}) where T <: NCRingElement From 9369964e5933aa829c215c81ea68d4c4141a5b25 Mon Sep 17 00:00:00 2001 From: JohnAAbbott <124266874+JohnAAbbott@users.noreply.github.com> Date: Mon, 8 Dec 2025 15:34:22 +0100 Subject: [PATCH 05/17] Update src/generic/MatRing.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Lars Göttgens --- src/generic/MatRing.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/generic/MatRing.jl b/src/generic/MatRing.jl index 0f20aec806..91abeb1f9a 100644 --- a/src/generic/MatRing.jl +++ b/src/generic/MatRing.jl @@ -35,9 +35,9 @@ is_domain_type(::Type{MatRingElem{T}}) where T <: NCRingElement = false # ############################################################################### -number_of_rows(a::MatRingElem) = nrows(a.data)::Int # Julia needs guidance here +number_of_rows(a::MatRingElem) = nrows(matrix(a)) -number_of_columns(a::MatRingElem) = ncols(a.data)::Int # Julia needs guidance here +number_of_columns(a::MatRingElem) = ncols(matrix(a)) Base.@propagate_inbounds getindex(a::MatRingElem, r::Int, c::Int) = a.data[r, c] From da2d2ea2260b49ec3536906f960a17789dd0a120 Mon Sep 17 00:00:00 2001 From: JohnAAbbott <124266874+JohnAAbbott@users.noreply.github.com> Date: Mon, 8 Dec 2025 15:35:09 +0100 Subject: [PATCH 06/17] Update src/generic/MatRing.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Lars Göttgens --- src/generic/MatRing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/generic/MatRing.jl b/src/generic/MatRing.jl index 91abeb1f9a..b73a4fc471 100644 --- a/src/generic/MatRing.jl +++ b/src/generic/MatRing.jl @@ -113,7 +113,7 @@ end function minpoly(M::MatRingElem{T}, charpoly_only::Bool = false) where {T <: RingElement} R = base_ring(M) Rx, x = polynomial_ring(R; cached=false) - return minpoly(Rx, matrix(M), charpoly_only) + return minpoly(Rx, M, charpoly_only) end ############################################################################### From e5a24d4033e1fa6c94a4f03d0380cde1a9c8c38d Mon Sep 17 00:00:00 2001 From: JohnAAbbott <124266874+JohnAAbbott@users.noreply.github.com> Date: Mon, 8 Dec 2025 15:38:47 +0100 Subject: [PATCH 07/17] Update src/generic/MatRing.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Lars Göttgens --- src/generic/MatRing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/generic/MatRing.jl b/src/generic/MatRing.jl index b73a4fc471..7bad3baca6 100644 --- a/src/generic/MatRing.jl +++ b/src/generic/MatRing.jl @@ -186,7 +186,7 @@ end function (a::MatRing{T})(b::Vector{S}) where {S <: NCRingElement, T <: NCRingElement} _check_dim(a.n, a.n, b) R = base_ring(a) - z = MatRingElem(matrix(R, a.n, a.n, b)) + z = MatRingElem(R, a.n, b) return z end From 7437ac038d6bd894ce573c30c51e32ac3d02f52d Mon Sep 17 00:00:00 2001 From: JohnAAbbott <124266874+JohnAAbbott@users.noreply.github.com> Date: Mon, 8 Dec 2025 15:41:20 +0100 Subject: [PATCH 08/17] Update src/generic/GenericTypes.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Lars Göttgens --- src/generic/GenericTypes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/generic/GenericTypes.jl b/src/generic/GenericTypes.jl index 16fd17586a..a59f74d375 100644 --- a/src/generic/GenericTypes.jl +++ b/src/generic/GenericTypes.jl @@ -1178,7 +1178,7 @@ function MatRingElem(R::NCRing, n::Int, A::Vector{T}) where T <: NCRingElement return MatRingElem(t) end -matrix(A::MatRingElem) = A.data +matrix(A::MatRingElem{T}) where {T} = A.data::dense_matrix_type(T) ############################################################################### # From 74045d890fcb50c0143ae4f636ad09b076cac9bc Mon Sep 17 00:00:00 2001 From: John Abbott Date: Mon, 8 Dec 2025 16:46:04 +0100 Subject: [PATCH 09/17] Revision following lgoettgens suggestions --- src/MatRing.jl | 15 ++++----------- src/generic/GenericTypes.jl | 4 ++-- 2 files changed, 6 insertions(+), 13 deletions(-) diff --git a/src/MatRing.jl b/src/MatRing.jl index 24452ffedd..0cb0765f69 100644 --- a/src/MatRing.jl +++ b/src/MatRing.jl @@ -89,10 +89,9 @@ Create an uninitialized matrix ring element over the given ring and dimension, with defaults based upon the given source matrix ring element `x`. """ function similar(x::MatRingElem, R::NCRing=base_ring(x), n::Int=degree(x)) - @assert (n >= 0) && (n < 2^30) # so that n^2 cannot overflow - TT = elem_type(R) - M = fill(zero(R), n^2) - return Generic.MatRingElem(R, n, M) + (n >= 0) || error("Matrix dimension must be non-negative") + (n < 2^30) || error("Matrix dimension is excessively large") + return Generic.MatRingElem(R, n, fill(0,n^2)) # n^2 cannot overflow given check in line above end similar(x::MatRingElem, n::Int) = similar(x, base_ring(x), n) @@ -429,13 +428,7 @@ end function identity_matrix(M::MatRingElem{T}, n::Int) where T <: NCRingElement @assert (n >= 0) && (n < 2^30) # so that n^2 cannot overflow R = base_ring(M) - # ??Simpler?? Generic.MatRingElem(matrix_space(R,n,n)(1)) - arr = fill(zero(R), n^2) - for i in 1:n - arr[i+(i-1)*n] = one(R) - end - z = Generic.MatRingElem(R, n, arr)::MatRingElem{T} - return z + return Generic.MatRingElem(identity_matrix(R,n)) end @doc raw""" diff --git a/src/generic/GenericTypes.jl b/src/generic/GenericTypes.jl index a59f74d375..1de27e4638 100644 --- a/src/generic/GenericTypes.jl +++ b/src/generic/GenericTypes.jl @@ -1173,12 +1173,12 @@ struct MatRingElem{T <: NCRingElement} <: AbstractAlgebra.MatRingElem{T} end function MatRingElem(R::NCRing, n::Int, A::Vector{T}) where T <: NCRingElement - @assert elem_type(R) === T +# @assert elem_type(R) === T t = matrix(R, n, n, A) return MatRingElem(t) end -matrix(A::MatRingElem{T}) where {T} = A.data::dense_matrix_type(T) +matrix(A::MatRingElem{T}) where {T <: NCRingElement} = A.data::dense_matrix_type(T) ############################################################################### # From 15c8c9ab5204e79f0637ae4b02f9b1de1d87af08 Mon Sep 17 00:00:00 2001 From: John Abbott Date: Tue, 9 Dec 2025 12:06:28 +0100 Subject: [PATCH 10/17] Replaced field accessors by accessor-function calls --- src/generic/MatRing.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/generic/MatRing.jl b/src/generic/MatRing.jl index 7bad3baca6..7e512eee34 100644 --- a/src/generic/MatRing.jl +++ b/src/generic/MatRing.jl @@ -23,7 +23,7 @@ base_ring(a::MatRingElem{T}) where {T <: NCRingElement} = base_ring(matrix(a)) Return the parent object of the given matrix. """ -parent(a::MatRingElem{T}) where T <: NCRingElement = MatRing{T}(base_ring(a), nrows(a.data)) +parent(a::MatRingElem{T}) where T <: NCRingElement = MatRing{T}(base_ring(a), nrows(matrix(a))) is_exact_type(::Type{MatRingElem{T}}) where T <: NCRingElement = is_exact_type(T) @@ -39,14 +39,14 @@ number_of_rows(a::MatRingElem) = nrows(matrix(a)) number_of_columns(a::MatRingElem) = ncols(matrix(a)) -Base.@propagate_inbounds getindex(a::MatRingElem, r::Int, c::Int) = a.data[r, c] +Base.@propagate_inbounds getindex(a::MatRingElem, r::Int, c::Int) = matrix(a)[r, c] Base.@propagate_inbounds function setindex!(a::MatRingElem, d::NCRingElement, r::Int, c::Int) - a.data[r, c] = base_ring(a)(d) + matrix(a)[r, c] = base_ring(a)(d) end -Base.isassigned(a::MatRingElem, i, j) = isassigned(a.data, i, j) +Base.isassigned(a::MatRingElem, i, j) = isassigned(matrix(a), i, j) ############################################################################### # @@ -123,7 +123,7 @@ end ############################################################################### function add!(A::MatRingElem{T}, B::MatRingElem{T}) where T <: NCRingElement - #=A.data ==# add!(A.data, B.data) ### !!struct is NOT mutable!! + #=A.data ==# add!(matrix(A), matrix(B)) ### !!struct is NOT mutable!! return A end From b379297d788d3ef59c3399d82167c907d6a09943 Mon Sep 17 00:00:00 2001 From: John Abbott Date: Tue, 9 Dec 2025 13:36:33 +0100 Subject: [PATCH 11/17] indentation --- src/generic/GenericTypes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/generic/GenericTypes.jl b/src/generic/GenericTypes.jl index 1de27e4638..a90c09c0d8 100644 --- a/src/generic/GenericTypes.jl +++ b/src/generic/GenericTypes.jl @@ -1168,7 +1168,7 @@ struct MatRingElem{T <: NCRingElement} <: AbstractAlgebra.MatRingElem{T} function MatRingElem(A::MatElem{TT}) where TT <: NCRingElement nrows(A) == ncols(A) || error("Matrix must be square") - return new{TT}(A) + return new{TT}(A) end end From 12ff76f0b38db9797e6286259566d5807c04d6ad Mon Sep 17 00:00:00 2001 From: John Abbott Date: Tue, 9 Dec 2025 21:32:20 +0100 Subject: [PATCH 12/17] Minor cleaning - thanks to thofma --- src/Matrix.jl | 2 -- src/generic/MatRing.jl | 2 +- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/src/Matrix.jl b/src/Matrix.jl index c0681bfa0c..64be275458 100644 --- a/src/Matrix.jl +++ b/src/Matrix.jl @@ -2619,8 +2619,6 @@ function minors_iterator(M::MatElem, k::Int) return (det(M[rows, cols]) for rows in row_indices for cols in col_indices) end -# TODO: how to get the 'positions' of the minors? - @doc raw""" minors_iterator_with_position(A::MatElem, k::Int) diff --git a/src/generic/MatRing.jl b/src/generic/MatRing.jl index 7e512eee34..911b130fc8 100644 --- a/src/generic/MatRing.jl +++ b/src/generic/MatRing.jl @@ -43,7 +43,7 @@ Base.@propagate_inbounds getindex(a::MatRingElem, r::Int, c::Int) = matrix(a)[r, Base.@propagate_inbounds function setindex!(a::MatRingElem, d::NCRingElement, r::Int, c::Int) - matrix(a)[r, c] = base_ring(a)(d) + matrix(a)[r, c] = d end Base.isassigned(a::MatRingElem, i, j) = isassigned(matrix(a), i, j) From 603a2a689a4222e5de5220ff4237f26e1345c8f9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lars=20G=C3=B6ttgens?= Date: Wed, 10 Dec 2025 11:04:58 +0100 Subject: [PATCH 13/17] Update src/NCPoly.jl --- src/NCPoly.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/NCPoly.jl b/src/NCPoly.jl index ce59dc9c2d..a5a91b9d0e 100644 --- a/src/NCPoly.jl +++ b/src/NCPoly.jl @@ -714,7 +714,6 @@ function rand(rng::AbstractRNG, <:AbstractUnitRange{Int}}}) S, deg_range, v = sp[][1:end] R = base_ring(S) - local f::elem_type(S) # Julia needs some guidance f = S() x = gen(S) for i = 0:rand(rng, deg_range) From 3e114bddedd00e28c66cdb65672c796ad2522c40 Mon Sep 17 00:00:00 2001 From: John Abbott Date: Wed, 10 Dec 2025 14:39:08 +0100 Subject: [PATCH 14/17] Use macro @req --- src/MatRing.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/MatRing.jl b/src/MatRing.jl index 0cb0765f69..5dc418a59d 100644 --- a/src/MatRing.jl +++ b/src/MatRing.jl @@ -89,8 +89,8 @@ Create an uninitialized matrix ring element over the given ring and dimension, with defaults based upon the given source matrix ring element `x`. """ function similar(x::MatRingElem, R::NCRing=base_ring(x), n::Int=degree(x)) - (n >= 0) || error("Matrix dimension must be non-negative") - (n < 2^30) || error("Matrix dimension is excessively large") + @req (n >= 0) "Matrix dimension must be non-negative" + @req (n < 2^30) "Matrix dimension is excessively large" return Generic.MatRingElem(R, n, fill(0,n^2)) # n^2 cannot overflow given check in line above end @@ -426,7 +426,8 @@ end ############################################################################### function identity_matrix(M::MatRingElem{T}, n::Int) where T <: NCRingElement - @assert (n >= 0) && (n < 2^30) # so that n^2 cannot overflow + @req (n >= 0) "Matrix dimension must be non-negative" + @req (n < 2^30) "Matrix dimension is excessively large" R = base_ring(M) return Generic.MatRingElem(identity_matrix(R,n)) end From 3ef99ac1471cb9bdb684f27947f45e641a586d85 Mon Sep 17 00:00:00 2001 From: John Abbott Date: Wed, 10 Dec 2025 14:39:44 +0100 Subject: [PATCH 15/17] Removed cruft --- src/generic/GenericTypes.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/generic/GenericTypes.jl b/src/generic/GenericTypes.jl index a90c09c0d8..dc6d3900e7 100644 --- a/src/generic/GenericTypes.jl +++ b/src/generic/GenericTypes.jl @@ -1173,7 +1173,6 @@ struct MatRingElem{T <: NCRingElement} <: AbstractAlgebra.MatRingElem{T} end function MatRingElem(R::NCRing, n::Int, A::Vector{T}) where T <: NCRingElement -# @assert elem_type(R) === T t = matrix(R, n, n, A) return MatRingElem(t) end From 5d734bd04ef5f1093de2702f1368511b38c3533f Mon Sep 17 00:00:00 2001 From: John Abbott Date: Wed, 10 Dec 2025 14:41:09 +0100 Subject: [PATCH 16/17] Revised add! according to fingolfin's suggestion --- src/generic/MatRing.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/generic/MatRing.jl b/src/generic/MatRing.jl index 911b130fc8..d725a64695 100644 --- a/src/generic/MatRing.jl +++ b/src/generic/MatRing.jl @@ -123,8 +123,7 @@ end ############################################################################### function add!(A::MatRingElem{T}, B::MatRingElem{T}) where T <: NCRingElement - #=A.data ==# add!(matrix(A), matrix(B)) ### !!struct is NOT mutable!! - return A + return MatRingElem(add!(matrix(A), matrix(B))) end ############################################################################### From 677d2d33f97dea18029624ebdbcb71836992e4ba Mon Sep 17 00:00:00 2001 From: John Abbott Date: Wed, 10 Dec 2025 14:42:19 +0100 Subject: [PATCH 17/17] New signatures for + and - following fingolfin suggestion --- src/Matrix.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Matrix.jl b/src/Matrix.jl index 64be275458..0874e46300 100644 --- a/src/Matrix.jl +++ b/src/Matrix.jl @@ -814,7 +814,7 @@ end # ############################################################################### -function +(x::MatElem{T}, y::MatElem{T}) where {T} +function +(x::T, y::T) where {T <: MatElem} check_parent(x, y) r = similar(x) for i = 1:nrows(x) @@ -825,7 +825,7 @@ function +(x::MatElem{T}, y::MatElem{T}) where {T} return r end -function -(x::MatElem{T}, y::MatElem{T}) where {T} +function -(x::T, y::T) where {T <: MatElem} check_parent(x, y) r = similar(x) for i = 1:nrows(x)