Howell form (#1735)

* Copy Howell form from Hecke

* Normalize the pivots to make it unique

* Tests

* Use `is_zero_entry`

* Test with another ring

* Rename

Author: joschmitt

src/MatrixNormalForms.jl

@@ -74,7 +74,7 @@ end
 
 ################################################################################
 #
-#  Matrix normal form over (euclidean) rings (hermite_form)
+#  Matrix normal form over (euclidean) domains (hermite_form)
 #
 ################################################################################
 
@@ -87,7 +87,7 @@ 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
+* `shape`: Whether $H$ 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
@@ -143,3 +143,251 @@ function hermite_form_with_transformation(A::MatElem{<:RingElement}; reduced::Bo
   end
   return H, U
 end
+
+################################################################################
+#
+#  Matrix normal form over principal ideal rings (howell_form)
+#
+################################################################################
+
+# Works in theory over any principal ideal ring; internally we require functions
+# annihilator, gcdxx and _div_for_howell_form
+
+# Swap rows so that there is a non-zero entry in A[start_row, col].
+# Return 0 if this is not possible, 1 if no swapping was necessary and -1
+# if rows were swapped.
+function _pivot(A::MatElem, start_row::Int, col::Int)
+  if !is_zero_entry(A, start_row, col)
+    return 1
+  end
+
+  for j in start_row + 1:nrows(A)
+    if !is_zero_entry(A, j, col)
+      swap_rows!(A, j, start_row)
+      return -1
+    end
+  end
+
+  return 0
+end
+
+function triangularize!(A::MatElem{<:RingElement})
+  n = nrows(A)
+  m = ncols(A)
+  d = one(base_ring(A))
+  row = 1
+  col = 1
+  while row <= nrows(A) && col <= ncols(A)
+    t = _pivot(A, row, col)
+    if iszero(t)
+      col = col + 1
+      continue
+    end
+    d = d*t
+    for i in (row + 1):nrows(A)
+      if is_zero_entry(A, i, col)
+        continue
+      end
+
+      b, q = divides(A[i, col], A[row, col])
+
+      if b
+        for k in col:m
+          A[i, k] = A[i, k] - q*A[row, k]
+        end
+      else
+        g, s, t, u, v = gcdxx(A[row, col], A[i, col])
+
+        for k in col:m
+          t1 = s*A[row, k] + t*A[i, k]
+          t2 = u*A[row, k] + v*A[i, k]
+          A[row, k] = t1
+          A[i, k] = t2
+        end
+      end
+    end
+    row = row + 1
+    col = col + 1
+  end
+  return d
+end
+
+# Naive version of inplace strong echelon form
+# It is assumed that A has more rows then columns.
+function strong_echelon_form_naive!(A::MatElem{<:RingElement})
+  n = nrows(A)
+  m = ncols(A)
+
+  @assert n >= m
+
+  triangularize!(A)
+
+  T = zero_matrix(base_ring(A), 1, ncols(A))
+  for j in 1:m
+    if !is_zero_entry(A, j, j)
+      # Normalize/canonicalize the pivot
+      u = canonical_unit(A[j, j])
+      if !is_one(u)
+        uinv = inv(u)
+        for i in j:ncols(A)
+          A[j, i] = uinv*A[j, i]
+        end
+      end
+
+      # This is the reduction
+      for i in 1:j - 1
+        if is_zero_entry(A, i, j)
+          continue
+        end
+        q = _div_for_howell_form(A[i, j], A[j, j])
+        for l in i:m
+          A[i, l] = A[i, l] - q*A[j, l]
+        end
+      end
+
+      a = annihilator(A[j, j])
+
+      for k in 1:m
+        T[1, k] = a*A[j, k]
+      end
+    else
+      for k in 1:m
+        T[1, k] = A[j, k]
+      end
+    end
+
+    for i in j+1:m
+
+      if is_zero_entry(T, 1, i)
+        continue
+      end
+
+      if is_zero_entry(A, i, i)
+        for k in i:m
+          T[1, k], A[i, k] = A[i, k], T[1, k]
+        end
+      else
+        b, q = divides(T[1, i], A[i, i])
+        if b
+          for k in i:m
+            T[1, k] = T[1, k] - q*A[i, k]
+          end
+        else
+          g, s, t, u, v = gcdxx(A[i, i], T[1, i])
+          for k in i:m
+            t1 = s*A[i, k] + t*T[1, k]
+            t2 = u*A[i, k] + v*T[1, k]
+            A[i, k] = t1
+            T[1, k] = t2
+          end
+        end
+      end
+    end
+  end
+  return A
+end
+
+function howell_form!(A::MatElem{<:RingElement})
+  @assert nrows(A) >= ncols(A)
+
+  strong_echelon_form_naive!(A)
+
+  # Move the zero rows to the bottom
+  for i in 1:nrows(A)
+    if !is_zero_row(A, i)
+      continue
+    end
+
+    j = findfirst(l -> !is_zero_row(A, l), i + 1:nrows(A))
+    if isnothing(j)
+      break
+    end
+    swap_rows!(A, i, i + j)
+  end
+  return A
+end
+
+@doc raw"""
+    howell_form(A::MatElem{<:RingElement}; reduced::Bool = true, shape::Symbol = :upper, trim::Bool = false)
+
+Return a Howell form $H$ of $A$.
+It is assumed that `base_ring(A)` is a principal ideal ring.
+
+# Keyword arguments
+* `reduced`: Whether the columns of $H$ which contain a pivot are reduced. The
+  default is `true`.
+* `shape`: Whether $H$ is an upper-right (`:upper`, default) or lower-left
+  (`:lower`) echelon form.
+* `trim`: By default, $H$ will have at least as many rows as $A$ and potentially
+  involve rows only containing zeros. If `trim = true`, these rows are removed.
+
+See also [`howell_form_with_transformation`](@ref).
+"""
+function howell_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
+
+  H = deepcopy(A)
+  if shape === :lower
+    H = reverse_cols!(H)
+  end
+
+  if nrows(H) < ncols(H)
+    H = vcat(H, zero_matrix(base_ring(H), ncols(H) - nrows(H), ncols(H)))
+  end
+
+  howell_form!(H)
+
+  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"""
+    howell_form_with_transformation(A::MatElem{<:RingElement}; reduced::Bool = true, shape::Symbol = :upper)
+
+Return a Howell form $H$ of $A$ and a matrix $U$ with $H = UA$.
+Notice that $H$ may have more rows than $A$ and hence $U$ may not be invertible.
+It is assumed that `base_ring(A)` is a principal ideal ring
+
+See [`hermite_form`](@ref) for the keyword arguments.
+"""
+function howell_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
+    B = hcat(reverse_cols(A), identity_matrix(A, nrows(A)))
+  else
+    B = hcat(A, identity_matrix(A, nrows(A)))
+  end
+  if nrows(B) < ncols(B)
+    B = vcat(B, zero(A, ncols(B) - nrows(B), ncols(B)))
+  end
+
+  howell_form!(B)
+
+  m = max(nrows(A), ncols(A))
+  H = sub(B, 1:m, 1:ncols(A))
+  U = sub(B, 1:m, ncols(A) + 1:ncols(B))
+
+  if shape === :lower
+    reverse_cols!(H)
+    reverse_rows!(H)
+    reverse_rows!(U)
+  end
+  return H, U
+end

src/Residue.jl

@@ -96,6 +96,8 @@ function is_zero_divisor_with_annihilator(a::ResElem)
    return !is_unit(g), parent(a)(b)
 end
 
+annihilator(a::ResElem) = is_zero_divisor_with_annihilator(a)[2]
+
 deepcopy_internal(a::ResElem, dict::IdDict) =
    parent(a)(deepcopy_internal(data(a), dict))
 
@@ -379,6 +381,32 @@ function gcd(a::ResElem{T}, b::ResElem{T}) where {T <: RingElement}
    return parent(a)(gcd(gcd(data(a), modulus(a)), data(b)))
 end
 
+# Return g, s, t, u, v in R with sv - tu a unit and
+#   [s t] [a] = [g]
+#   [u v] [b]   [0]
+# Generic implementation which uses HNF over the base ring.
+# g might not coincide with gcd(a, b) because gcd(a, b) is
+# gcd(gcd(data(a), modulus(a)), data(b)) and g is just
+# gcd(data(a), data(b)).
+function gcdxx(a::ResElem{T}, b::ResElem{T}) where {T <: RingElement}
+  check_parent(a, b)
+  R = parent(a)
+  M = matrix(base_ring(R), 2, 1, [data(a), data(b)])
+  H, U = hermite_form_with_transformation(M)
+  @assert is_zero(H[2, 1])
+  return R(H[1, 1]), R(U[1, 1]), R(U[1, 2]), R(U[2, 1]), R(U[2, 2])
+end
+
+# The operation "Quo" on p. 13 of Storjohann "Algorithms for matrix canonical forms"
+function _div_for_howell_form(a::ResElem{T}, b::ResElem{T}) where {T <: RingElement}
+  check_parent(a, b)
+  return parent(a)(div(data(a), data(b)))
+end
+
+# Fallback for euclidean rings (that is, rings implementing the euclidean ring
+# interface)
+_div_for_howell_form(a::T, b::T) where {T <: RingElement} = div(a, b)
+
 ###############################################################################
 #
 #   Unsafe functions

src/exports.jl

@@ -243,6 +243,8 @@ export hnf_with_transform
 export hom
 export hom_direct_sum
 export hooklength
+export howell_form
+export howell_form_with_transformation
 export ideal
 export identity_map
 export identity_matrix

test/MatrixNormalForms-test.jl

@@ -23,6 +23,7 @@
   R = reverse_rows(reverse_cols(S))
   @test is_rref(R)
   @test U*M == S
+  @test is_invertible(U)
 end
 
 @testset "Hermite form" begin
@@ -50,4 +51,53 @@ end
   H = reverse_rows(reverse_cols(Hl))
   @test is_hnf(H)
   @test U*M == Hl
+  @test is_invertible(U)
+end
+
+@testset "Howell form" begin
+  R, _ = residue_ring(ZZ, 30)
+  M = R[2 3 4 5; 6 9 12 15; 10 15 20 30]
+  H = howell_form(M)
+  @test H == R[2 3 4 0; 0 15 0 0; 0 0 0 5; 0 0 0 0]
+  Hl = howell_form(M, shape = :lower)
+  @test Hl == R[0 0 0 0; 0 15 0 0; 16 9 2 0; 0 0 0 5]
+  H = howell_form(M, trim = true)
+  @test H == R[2 3 4 0; 0 15 0 0; 0 0 0 5]
+  Hl = howell_form(M, shape = :lower, trim = true)
+  @test Hl == R[0 15 0 0; 16 9 2 0; 0 0 0 5]
+
+  H, U = howell_form_with_transformation(M)
+  @test H == R[2 3 4 0; 0 15 0 0; 0 0 0 5; 0 0 0 0]
+  @test U*M == H
+  Hl, U = howell_form_with_transformation(M, shape = :lower)
+  @test Hl == R[0 0 0 0; 0 15 0 0; 16 9 2 0; 0 0 0 5]
+  @test U*M == Hl
+
+  R, _ = residue_ring(ZZ, 12)
+  M = R[4 1 0; 3 0 0; 0 0 5]
+  @test howell_form(M) == R[1 1 0; 0 3 0; 0 0 1]
+
+  a = R[4 1 0; 0 0 5; 0 0 0]
+  b = R[8 5 5; 0 9 8; 0 0 10]
+  c = R[4 1 0; 0 3 0; 0 0 1]
+  @test howell_form(a) == c
+  @test howell_form(b) == c
+
+  Qt, t = QQ["t"]
+  R, _ = residue_ring(Qt, t^5)
+  M = R[t^2 t^2 - 1 1; t t - 1 0]
+  H = howell_form(M)
+  @test H == R[t 0 -1; 0 1 -t^3 - t^2 - t - 1; 0 0 t^4]
+  H = howell_form(M, shape = :lower)
+  @test H == R[0 0 0; -t^4 - t^3 - t^2 - t 1 0; -t 0 1]
+  H = howell_form(M, trim = true)
+  @test H == R[t 0 -1; 0 1 -t^3 - t^2 - t - 1; 0 0 t^4]
+  H = howell_form(M, shape = :lower, trim = true)
+  @test H == R[-t^4 - t^3 - t^2 - t 1 0; -t 0 1]
+  H, U = howell_form_with_transformation(M)
+  @test H == R[t 0 -1; 0 1 -t^3 - t^2 - t - 1; 0 0 t^4]
+  @test U*M == H
+  H, U = howell_form_with_transformation(M, shape = :lower)
+  @test H == R[0 0 0; -t^4 - t^3 - t^2 - t 1 0; -t 0 1]
+  @test U*M == H
 end