Rewrite MatRing to wrap a native matrix

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)
@@ -180,28 +181,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 +225,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 +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)
@@ -424,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
 

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::MatElem{T}, y::MatElem{T}) where {T}
    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::MatElem{T}, y::MatElem{T}) where {T}
    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
@@ -2619,6 +2619,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)
 

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)

src/generic/GenericTypes.jl

@@ -1164,24 +1164,22 @@ 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
+function MatRingElem(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)
 end
 
+matrix(A::MatRingElem) = A.data
+
 ###############################################################################
 #
 #   FreeAssociativeAlgebra / FreeAssociativeAlgebraElem