Skip to content

Howell form #1735

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 6 commits into from
Jul 2, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
252 changes: 250 additions & 2 deletions src/MatrixNormalForms.jl
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@

################################################################################
#
# Matrix normal form over (euclidean) rings (hermite_form)
# Matrix normal form over (euclidean) domains (hermite_form)
#
################################################################################

Expand All @@ -87,7 +87,7 @@
# 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
Expand Down Expand Up @@ -143,3 +143,251 @@
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."))

Check warning on line 328 in src/MatrixNormalForms.jl

View check run for this annotation

Codecov / codecov/patch

src/MatrixNormalForms.jl#L328

Added line #L328 was not covered by tests
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."))

Check warning on line 369 in src/MatrixNormalForms.jl

View check run for this annotation

Codecov / codecov/patch

src/MatrixNormalForms.jl#L369

Added line #L369 was not covered by tests
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
28 changes: 28 additions & 0 deletions src/Residue.jl
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,8 @@
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))

Expand Down Expand Up @@ -379,6 +381,32 @@
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)

Check warning on line 408 in src/Residue.jl

View check run for this annotation

Codecov / codecov/patch

src/Residue.jl#L408

Added line #L408 was not covered by tests

###############################################################################
#
# Unsafe functions
Expand Down
2 changes: 2 additions & 0 deletions src/exports.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
50 changes: 50 additions & 0 deletions test/MatrixNormalForms-test.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Loading