Skip to content
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.

Sign up for GitHub

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

Jump to bottom

Some matrix function cleanup #1434

Merged
merged 6 commits into from
Sep 13, 2023
Merged

Some matrix function cleanup #1434

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
5de5a38
Move `charpoly` from NemoStuff.jl to Matrix.jl; refactor a bit
lgoettgens Sep 13, 2023
3cb8efb
Move `minpoly` from NemoStuff.jl to Matrix.jl; refactor a bit
lgoettgens Sep 13, 2023
c7e8c6d
Remove duplicate `*_kernel` methods
lgoettgens Sep 13, 2023
b6ad6e3
Move `kernel` over different ring from NemoStuff.jl to Matrix.jl
lgoettgens Sep 13, 2023
3561519
Make all `*_kernel` docstrings consistent in variable naming
lgoettgens Sep 13, 2023
9b0a0e6
Fix docs
lgoettgens Sep 13, 2023
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
4 changes: 2 additions & 2 deletions docs/src/matrix.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1208,13 +1208,13 @@ true
### Characteristic polynomial

```@docs
charpoly{T <: RingElem}(::Ring, ::MatElem{T})
charpoly{T <: RingElem}(::PolyRing{T}, ::MatrixElem{T})
```

### Minimal polynomial

```@docs
minpoly{T <: RingElem}(::Ring, ::MatElem{T}, ::Bool)
minpoly{T <: RingElem}(::PolyRing{T}, ::MatElem{T}, ::Bool)
```

### Transforms
Expand Down
113 changes: 77 additions & 36 deletions src/Matrix.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3864,17 +3864,17 @@
###############################################################################

@doc raw"""
left_kernel(a::MatElem{T}) where T <: RingElement
left_kernel(A::MatElem{T}) where T <: RingElement

Return a tuple `n, M` where $M$ is a matrix whose rows generate the kernel
of $M$ and $n$ is the rank of the kernel. The transpose of the output of this
Return a tuple $(n, M)$ where $M$ is a matrix whose rows generate the kernel
of $A$ and $n$ is the rank of the kernel. The transpose of the output of this
function is guaranteed to be in flipped upper triangular format (i.e. upper
triangular format if columns and rows are reversed).
"""
function left_kernel(x::MatElem{T}) where T <: RingElement
!is_domain_type(elem_type(base_ring(x))) && error("Not implemented")
R = base_ring(x)
H, U = hnf_with_transform(x)
function left_kernel(A::MatElem{T}) where T <: RingElement
!is_domain_type(elem_type(base_ring(A))) && error("Not implemented")
R = base_ring(A)
H, U = hnf_with_transform(A)
i = nrows(H)
zero_rows = false
while i > 0 && is_zero_row(H, i)
Expand All @@ -3888,30 +3888,30 @@
end
end

function left_kernel(M::MatElem{T}) where T <: FieldElement
n, N = nullspace(transpose(M))
return n, transpose(N)
function left_kernel(A::MatElem{T}) where T <: FieldElement
n, M = nullspace(transpose(A))
return n, transpose(M)
end

@doc raw"""
right_kernel(a::MatElem{T}) where T <: RingElement
right_kernel(A::MatElem{T}) where T <: RingElement

Return a tuple `n, M` where $M$ is a matrix whose columns generate the
kernel of $a$ and $n$ is the rank of the kernel.
Return a tuple $(n, M)$ where $M$ is a matrix whose columns generate the
kernel of $A$ and $n$ is the rank of the kernel.
"""
function right_kernel(x::MatElem{T}) where T <: RingElement
n, M = left_kernel(transpose(x))
function right_kernel(A::MatElem{T}) where T <: RingElement
n, M = left_kernel(transpose(A))
return n, transpose(M)
end

function right_kernel(M::MatElem{T}) where T <: FieldElement
return nullspace(M)
function right_kernel(A::MatElem{T}) where T <: FieldElement
return nullspace(A)
end

@doc raw"""
kernel(a::MatElem{T}; side::Symbol = :right) where T <: RingElement
kernel(A::MatElem{T}; side::Symbol = :right) where T <: RingElement

Return a tuple $(n, M)$, where $n$ is the rank of the kernel of $a$ and $M$ is a
Return a tuple $(n, M)$, where $n$ is the rank of the kernel of $A$ and $M$ is a
basis for it. If side is `:right` or not specified, the right kernel is
computed, i.e. the matrix of columns whose span gives the right kernel
space. If side is `:left`, the left kernel is computed, i.e. the matrix
Expand All @@ -3927,6 +3927,26 @@
end
end

################################################################################
#
# Kernel over different rings
#
################################################################################

@doc raw"""
kernel(A::MatElem{T}; R::Ring, side::Symbol = :right) where T <: RingElement

Return a tuple $(n, M)$, where $n$ is the rank of the kernel of $A$ over $R$ and $M$ is a
basis for it. If side is `:right` or not specified, the right kernel is
computed, i.e. the matrix of columns whose span gives the right kernel
space. If side is `:left`, the left kernel is computed, i.e. the matrix
of rows whose span is the left kernel space.
"""
function kernel(A::MatElem{T}, R::Ring; side::Symbol=:right) where T <: RingElement
AR = change_base_ring(R, A)
return kernel(AR; side)

Check warning on line 3947 in src/Matrix.jl

View check run for this annotation

Codecov / codecov/patch

src/Matrix.jl#L3945-L3947 

Added lines #L3945 - L3947 were not covered by tests
end

###############################################################################
#
# Hessenberg form
Expand Down Expand Up @@ -4260,11 +4280,12 @@
end

@doc raw"""
charpoly(V::Ring, Y::MatrixElem{T}) where {T <: RingElement}
charpoly(Y::MatrixElem{T}) where {T <: RingElement}
charpoly(S::PolyRing{T}, Y::MatrixElem{T}) where {T <: RingElement}

Return the characteristic polynomial $p$ of the matrix $M$. The
polynomial ring $R$ of the resulting polynomial must be supplied
and the matrix is assumed to be square.
Return the characteristic polynomial $p$ of the square matrix $Y$.
If a polynomial ring $S$ over the same base ring as $Y$ is supplied,
the resulting polynomial is an element of it.

# Examples

Expand All @@ -4276,8 +4297,8 @@
Matrix space of 4 rows and 4 columns
over residue ring of integers modulo 7

julia> T, x = polynomial_ring(R, "x")
(Univariate polynomial ring in x over residue ring, x)
julia> T, y = polynomial_ring(R, "y")
(Univariate polynomial ring in y over residue ring, y)

julia> M = S([R(1) R(2) R(4) R(3); R(2) R(5) R(1) R(0);
R(6) R(1) R(3) R(2); R(1) R(1) R(3) R(5)])
Expand All @@ -4287,17 +4308,20 @@
[1 1 3 5]

julia> A = charpoly(T, M)
y^4 + 2*y^2 + 6*y + 2

julia> A = charpoly(M)
x^4 + 2*x^2 + 6*x + 2

```
"""
function charpoly(V::Ring, Y::MatrixElem{T}) where {T <: RingElement}
function charpoly(S::PolyRing{T}, Y::MatrixElem{T}) where {T <: RingElement}
!is_square(Y) && error("Dimensions don't match in charpoly")
R = base_ring(Y)
base_ring(V) != base_ring(Y) && error("Cannot coerce into polynomial ring")
base_ring(S) != base_ring(Y) && error("Cannot coerce into polynomial ring")
n = nrows(Y)
if n == 0
return one(V)
return one(S)
end
F = Vector{elem_type(R)}(undef, n)
A = Vector{elem_type(R)}(undef, n)
Expand Down Expand Up @@ -4336,14 +4360,20 @@
F[j] = -s - A[j]
end
end
z = gen(V)
z = gen(S)
f = z^n
for i = 1:n
f = setcoeff!(f, n - i, F[i])
end
return f
end

function charpoly(Y::MatrixElem)
R = base_ring(Y)
Rx, x = polynomial_ring(R; cached=false)
return charpoly(Rx, Y)
end

###############################################################################
#
# Minimal polynomial
Expand All @@ -4366,7 +4396,7 @@
# charpoly iff it has degree n. Otherwise it is meaningless (but it is
# extremely fast to compute over some fields).

function minpoly(S::Ring, M::MatElem{T}, charpoly_only::Bool = false) where {T <: FieldElement}
function minpoly(S::PolyRing{T}, M::MatElem{T}, charpoly_only::Bool = false) where {T <: FieldElement}
!is_square(M) && error("Not a square matrix in minpoly")
base_ring(S) != base_ring(M) && error("Unable to coerce polynomial")
n = nrows(M)
Expand Down Expand Up @@ -4457,18 +4487,20 @@
end

@doc raw"""
minpoly(S::Ring, M::MatElem{T}, charpoly_only::Bool = false) where {T <: RingElement}
minpoly(M::MatElem{T}, charpoly_only::Bool = false) where {T <: RingElement}
minpoly(S::PolyRing{T}, M::MatElem{T}, charpoly_only::Bool = false) where {T <: RingElement}

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.
Return the minimal polynomial $p$ of the square matrix $M$.
If a polynomial ring $S$ over the same base ring as $Y$ is supplied,
the resulting polynomial is an element of it.

# Examples

```jldoctest; setup = :(using AbstractAlgebra)
julia> R = GF(13)
Finite field F_13

julia> T, y = polynomial_ring(R, "y")
julia> S, y = polynomial_ring(R, "y")
(Univariate polynomial ring in y over finite field F_13, y)

julia> M = R[7 6 1;
Expand All @@ -4478,12 +4510,15 @@
[7 7 5]
[8 12 5]

julia> A = minpoly(T, M)
julia> A = minpoly(S, M)
y^2 + 10*y

julia> A = minpoly(M)
x^2 + 10*x

```
"""
function minpoly(S::Ring, M::MatElem{T}, charpoly_only::Bool = false) where {T <: RingElement}
function minpoly(S::PolyRing{T}, M::MatElem{T}, charpoly_only::Bool = false) where {T <: RingElement}
!is_square(M) && error("Not a square matrix in minpoly")
base_ring(S) != base_ring(M) && error("Unable to coerce polynomial")
n = nrows(M)
Expand Down Expand Up @@ -4583,6 +4618,12 @@
return divexact(p, canonical_unit(p))
end

function minpoly(M::MatElem{T}, charpoly_only::Bool = false) where {T <: RingElement}
R = base_ring(M)
Rx, x = polynomial_ring(R; cached=false)
return minpoly(Rx, M, charpoly_only)
end

###############################################################################
#
# Hermite Normal Form
Expand Down
42 changes: 0 additions & 42 deletions src/NemoStuff.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -208,24 +208,6 @@ Base.copy(f::Generic.Poly) = deepcopy(f)
Base.copy(a::PolyRingElem) = deepcopy(a)
Base.copy(a::SeriesElem) = deepcopy(a)

################################################################################
#
# Minpoly and Charpoly
#
################################################################################

function minpoly(M::MatElem)
k = base_ring(M)
kx, x = polynomial_ring(k, cached=false)
return minpoly(kx, M)
end

function charpoly(M::MatElem)
k = base_ring(M)
kx, x = polynomial_ring(k, cached=false)
return charpoly(kx, M)
end

###############################################################################
#
# Sub
Expand Down Expand Up @@ -263,30 +245,6 @@ function sub(M::Generic.Mat, rows::AbstractUnitRange{Int}, cols::AbstractUnitRan
return z
end

right_kernel(M::MatElem) = nullspace(M)

function left_kernel(M::MatElem)
rk, M1 = nullspace(transpose(M))
return rk, transpose(M1)
end

################################################################################
#
# Kernel over different rings
#
################################################################################

@doc raw"""
kernel(a::MatrixElem{T}, R::Ring; side::Symbol = :right) -> n, MatElem{elem_type(R)}

It returns a tuple $(n, M)$, where $n$ is the rank of the kernel over $R$ and $M$ is a basis for it. If side is $:right$ or not
specified, the right kernel is computed. If side is $:left$, the left kernel is computed.
"""
function kernel(M::MatrixElem, R::Ring; side::Symbol=:right)
MP = change_base_ring(R, M)
return kernel(MP, side=side)
end

################################################################################
#
# Concatenation of matrices
Expand Down
Loading