Matrix normal forms (Nemocas#1554)

* Skeleton for matrix normal forms
* Matrix normal forms (by calling existing functions)

src/AbstractAlgebra.jl

@@ -636,6 +636,7 @@ include("MPoly.jl")
 include("UnivPoly.jl")
 include("FreeAssAlgebra.jl")
 include("LaurentMPoly.jl")
+include("MatrixNormalForms.jl")
 
 ###############################################################################
 #
@@ -869,10 +870,12 @@ export divhigh
 export divides
 export domain
 export downscale
-export EuclideanRingResidueField
-export EuclideanRingResidueRing
+export echelon_form
+export echelon_form_with_transformation
 export elem_type
 export enable_cache!
+export EuclideanRingResidueField
+export EuclideanRingResidueRing
 export evaluate
 export exp_gcd
 export exponent
@@ -909,6 +912,8 @@ export has_bottom_neighbor
 export has_gens
 export has_left_neighbor
 export hash
+export hermite_form
+export hermite_form_with_transformation
 export hessenberg
 export hessenberg!
 export hnf

src/MatrixNormalForms.jl

@@ -0,0 +1,145 @@
+################################################################################
+#
+#  Matrix normal form over fields (echelon_form)
+#
+################################################################################
+
+@doc raw"""
+    echelon_form(A::MatElem{<:FieldElement}; reduced::Bool = true, shape::Symbol = :upper, trim::Bool = false)
+
+Return a row echelon form $R$ of $A$.
+
+# Keyword arguments
+* `reduced`: Whether the columns of $R$ which contain a pivot are reduced. The
+  default is `true`.
+* `shape`: Whether $R$ is an upper-right (`:upper`, default) or lower-left
+  (`:lower`) echelon form.
+* `trim`: By default, $R$ will have as many rows as $A$ and potentially involve
+  rows only containing zeros. If `trim = true`, these rows are removed, so that
+  the number of rows of $R$ coincides with the rank of $A$.
+
+See also [`echelon_form_with_transformation`](@ref).
+"""
+function echelon_form(A::MatElem{<:FieldElement}; reduced::Bool = true, shape::Symbol = :upper, trim::Bool = false)
+  R = deepcopy(A)
+  r = echelon_form!(R, reduced = reduced, shape = shape)
+  if trim
+    if shape === :upper
+      return sub(R, 1:r, 1:ncols(R))
+    else
+      return sub(R, nrows(R) - r + 1:nrows(R), 1:ncols(R))
+    end
+  end
+  return R
+end
+
+@doc raw"""
+    echelon_form_with_transformation(A::MatElem{<:FieldElement}; reduced::Bool = true, shape::Symbol = :upper)
+
+Return a row echelon form $R$ of $A$ and an invertible matrix $U$ with $R = UA$.
+
+See [`echelon_form`](@ref) for the keyword arguments.
+"""
+function echelon_form_with_transformation(A::MatElem{<:FieldElement}; reduced::Bool = true, shape::Symbol = :upper)
+  if shape === :upper
+    R = hcat(A, identity_matrix(base_ring(A), nrows(A)))
+  else
+    R = hcat(identity_matrix(base_ring(A), nrows(A)), A)
+  end
+  echelon_form!(R, reduced = reduced, shape = shape)
+  if shape === :upper
+    return sub(R, 1:nrows(A), 1:ncols(A)), sub(R, 1:nrows(A), ncols(A) + 1:ncols(R))
+  else
+    return sub(R, 1:nrows(A), nrows(A) + 1:ncols(R)), sub(R, 1:nrows(A), 1:nrows(A))
+  end
+end
+
+# Return the rank of A and put A into echelon form
+# So far, the `reduced` keyword is ignored (that is, the result is always reduced)
+function echelon_form!(A::MatElem{<:FieldElement}; reduced::Bool = true, shape::Symbol = :upper)
+  if shape !== :upper && shape !== :lower
+    throw(ArgumentError("Unsupported argument :$shape for shape: Must be :upper or :lower."))
+  end
+
+  if shape === :lower
+    reverse_cols!(A)
+    r = echelon_form!(A, reduced = reduced, shape = :upper)
+    reverse_cols!(A)
+    reverse_rows!(A)
+    return r
+  end
+
+  return rref!(A)
+end
+
+################################################################################
+#
+#  Matrix normal form over (euclidean) rings (hermite_form)
+#
+################################################################################
+
+@doc raw"""
+    hermite_form(A::MatElem{<:RingElement}; reduced::Bool = true, shape::Symbol = :upper, trim::Bool = false)
+
+Return a Hermite normal form $H$ of $A$.
+It is assumed that `base_ring(A)` is euclidean.
+
+# Keyword arguments
+* `reduced`: Whether the columns of $H$ which contain a pivot are reduced. The
+  default is `true`.
+* `shape`: Whether $R$ is an upper-right (`:upper`, default) or lower-left
+  (`:lower`) echelon form.
+* `trim`: By default, $H$ will have as many rows as $A$ and potentially involve
+  rows only containing zeros. If `trim = true`, these rows are removed, so that
+  the number of rows of $H$ coincides with the rank of $A$.
+
+See also [`hermite_form_with_transformation`](@ref).
+"""
+function hermite_form(A::MatElem{<:RingElement}; reduced::Bool = true, shape::Symbol = :upper, trim::Bool = false)
+  if shape !== :upper && shape !== :lower
+    throw(ArgumentError("Unsupported argument :$shape for shape: Must be :upper or :lower."))
+  end
+
+  if shape === :lower
+    A = reverse_cols(A)
+  end
+  H = hnf_kb(A)
+  r = nrows(H)
+  # Compute the rank (if necessary)
+  if trim
+    while r > 0 && is_zero_row(H, r)
+      r -= 1
+    end
+    H = sub(H, 1:r, 1:ncols(H))
+  end
+  if shape === :lower
+    reverse_cols!(H)
+    reverse_rows!(H)
+  end
+  return H
+end
+
+@doc raw"""
+    hermite_form_with_transformation(A::MatElem{<:RingElement}; reduced::Bool = true, shape::Symbol = :upper)
+
+Return a Hermite normal form $H$ of $A$ and an invertible matrix $U$ with $H = UA$.
+It is assumed that `base_ring(A)` is euclidean.
+
+See [`hermite_form`](@ref) for the keyword arguments.
+"""
+function hermite_form_with_transformation(A::MatElem{<:RingElement}; reduced::Bool = true, shape::Symbol = :upper)
+  if shape !== :upper && shape !== :lower
+    throw(ArgumentError("Unsupported argument :$shape for shape: Must be :upper or :lower."))
+  end
+
+  if shape === :lower
+    A = reverse_cols(A)
+  end
+  H, U = hnf_kb_with_transform(A)
+  if shape === :lower
+    reverse_cols!(H)
+    reverse_rows!(H)
+    reverse_rows!(U)
+  end
+  return H, U
+end

test/AbstractAlgebra-test.jl

@@ -19,3 +19,4 @@ include("PrettyPrinting-test.jl")
 include("PrintHelper-test.jl")
 include("misc-test.jl")
 include("Solve-test.jl")
+include("MatrixNormalForms-test.jl")

test/MatrixNormalForms-test.jl

@@ -0,0 +1,53 @@
+@testset "Echelon form" begin
+  M = matrix(QQ, [ 1 2 3 4 ; 5 6 7 8 ; 9 10 11 12 ])
+  R = echelon_form(M)
+  @test is_rref(R)
+  @test nrows(R) == 3
+  S = echelon_form(M, shape = :lower)
+  R = reverse_rows(reverse_cols(S))
+  @test is_rref(R)
+  @test nrows(S) == 3
+  R = echelon_form(M, trim = true)
+  @test is_rref(R)
+  @test nrows(R) == 2
+  S = echelon_form(M, shape = :lower, trim = true)
+  R = reverse_rows(reverse_cols(S))
+  @test is_rref(R)
+  @test nrows(S) == 2
+
+  R, U = echelon_form_with_transformation(M)
+  @test is_rref(R)
+  @test U*M == R
+  @test is_invertible(U)
+  S, U = echelon_form_with_transformation(M, shape = :lower)
+  R = reverse_rows(reverse_cols(S))
+  @test is_rref(R)
+  @test U*M == S
+end
+
+@testset "Hermite form" begin
+  M = matrix(ZZ, [ 1 2 3 4 ; 5 6 7 8 ; 9 10 11 12 ])
+  H = hermite_form(M)
+  @test is_hnf(H)
+  @test nrows(H) == 3
+  Hl = hermite_form(M, shape = :lower)
+  H = reverse_rows(reverse_cols(Hl))
+  @test is_hnf(H)
+  @test nrows(Hl) == 3
+  H = hermite_form(M, trim = true)
+  @test is_hnf(H)
+  @test nrows(H) == 2
+  Hl = hermite_form(M, shape = :lower, trim = true)
+  H = reverse_rows(reverse_cols(Hl))
+  @test is_hnf(H)
+  @test nrows(Hl) == 2
+
+  H, U = hermite_form_with_transformation(M)
+  @test is_hnf(H)
+  @test U*M == H
+  @test is_invertible(U)
+  Hl, U = hermite_form_with_transformation(M, shape = :lower)
+  H = reverse_rows(reverse_cols(Hl))
+  @test is_hnf(H)
+  @test U*M == Hl
+end