Rewrite MatRing to wrap a native matrix

src/MatRing.jl

@@ -89,9 +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))
-   TT = elem_type(R)
-   M = Matrix{TT}(undef, (n, n))
-   return Generic.MatRingElem{TT}(R, M)
+   @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
 
 similar(x::MatRingElem, n::Int) = similar(x, base_ring(x), n)
@@ -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
@@ -246,20 +224,41 @@ 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)
    return divexact(ginv*f, d; check=check)
 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
@@ -319,6 +318,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)
@@ -424,15 +426,10 @@ 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)
-   return z
+  @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
 
 @doc raw"""

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.
 """
@@ -814,7 +814,7 @@ end
 #
 ###############################################################################
 
-function +(x::MatrixElem{T}, y::MatrixElem{T}) where {T <: NCRingElement}
+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::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 <: MatElem}
    check_parent(x, y)
    r = similar(x)
    for i = 1:nrows(x)
@@ -836,7 +836,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::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))
@@ -1188,7 +1188,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

src/generic/GenericTypes.jl

@@ -1164,24 +1164,21 @@ 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(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
-   @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)
+function MatRingElem(R::NCRing, n::Int, A::Vector{T}) where T <: NCRingElement
+   t = matrix(R, n, n, A)
+   return MatRingElem(t)
 end
 
+matrix(A::MatRingElem{T}) where {T <: NCRingElement} = A.data::dense_matrix_type(T)
+
 ###############################################################################
 #
 #   FreeAssociativeAlgebra / FreeAssociativeAlgebraElem

src/generic/MatRing.jl

@@ -16,27 +16,46 @@ 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(matrix(a))
+
 @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(matrix(a)))
 
 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(matrix(a))
+
+number_of_columns(a::MatRingElem) = ncols(matrix(a))
+
+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)
+    matrix(a)[r, c] = d
+end
+
+Base.isassigned(a::MatRingElem, i, j) = isassigned(matrix(a), i, j)
+
 ###############################################################################
 #
 #   Transpose
 #
 ###############################################################################
 
-function transpose(x::MatRingElem{T}) where T <: NCRingElement
-   arr = permutedims(x.entries)
-   z = MatRingElem{T}(base_ring(x), arr)
-   return z
+function transpose(x::MatRingElem)
+   return MatRingElem(transpose(matrix(x)))
 end
 
 ###############################################################################
@@ -48,32 +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 = MatSpaceElem{T}(R, M.entries) # convert to ordinary matrix
-   BS = MatSpaceElem{T}(R, B.entries)
+   MS = matrix(M)
+   BS = matrix(B)
    flag, S = _can_solve_with_solution_lu(MS, BS)
-   SA = MatRingElem{T}(R, S.entries)
+   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)
-   # TODO: Once #1955 is resolved, the conversion to matrix and back to MatRingElem
-   # should be done better
-   MS = matrix(R, M.entries) # convert to ordinary matrix
-   BS = matrix(R, B.entries)
+   MS = matrix(M)
+   BS = matrix(B)
    flag, S = can_solve_with_solution(MS, BS)
-   SA = MatRingElem{T}(R, Array(S))
+   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 = MatSpaceElem{T}(R, M.entries) # convert to ordinary matrix
-   BS = MatSpaceElem{T}(R, B.entries)
+   MS = matrix(M)
+   BS = matrix(B)
    flag, S, d = _can_solve_with_solution_fflu(MS, BS)
-   SA = MatRingElem{T}(R, S.entries)
+   SA = MatRingElem(S)
    return flag, SA, d
 end
 
@@ -90,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, matrix(M), charpoly_only)
 end
 
 function minpoly(M::MatRingElem{T}, charpoly_only::Bool = false) where {T <: RingElement}
@@ -107,8 +123,7 @@ end
 ###############################################################################
 
 function add!(A::MatRingElem{T}, B::MatRingElem{T}) where T <: NCRingElement
-   A.entries .+= B.entries
-   return A
+   return MatRingElem(add!(matrix(A), matrix(B)))
 end
 
 ###############################################################################
@@ -131,48 +146,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(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(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(scalar_matrix(R, a.n, R(b)))
    return z
 end
 
@@ -184,33 +171,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)
-   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(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)
+   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)
-   b = Matrix{S}(transpose(reshape(b, a.n, a.n)))
-   z = a(b)
+  _check_dim(a.n, a.n, b)
+   R = base_ring(a)
+   z = MatRingElem(R, a.n, b)
    return z
 end