Skip to content

Commit

Permalink
Misc cleanup (#1770)
Browse files Browse the repository at this point in the history 
* Fix typo

* Move some ZZMatrix code to fmpz_mat.jl

* Enable 'f(x)' support where f is an arb/acb/real/complex poly

* Remove redundant is_zero_row methods
  • Loading branch information
@fingolfin
fingolfin committed Jun 4, 2024
1 parent 2c69d79 commit c3368ac
Show file tree
Hide file tree
Showing 5 changed files with 208 additions and 228 deletions.
220 changes: 0 additions & 220 deletions src/HeckeMiscMatrix.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,19 +29,6 @@ function Array(a::ZZMatrix; S::Type{T}=ZZRingElem) where {T}
return A
end

function is_zero_row(M::ZZMatrix, i::Int)
GC.@preserve M begin
for j = 1:ncols(M)
m = ccall((:fmpz_mat_entry, libflint), Ptr{ZZRingElem}, (Ref{ZZMatrix}, Int, Int), M, i - 1, j - 1)
fl = ccall((:fmpz_is_zero, libflint), Bool, (Ptr{ZZRingElem},), m)
if !fl
return false
end
end
end
return true
end

function is_zero_row(M::zzModMatrix, i::Int)
zero = UInt(0)
for j in 1:ncols(M)
Expand All @@ -53,26 +40,6 @@ function is_zero_row(M::zzModMatrix, i::Int)
return true
end

function is_zero_row(M::Matrix{ZZRingElem}, i::Int)
for j = 1:Base.size(M, 2)
if !iszero(M[i, j])
return false
end
end
return true
end


function divexact!(a::ZZMatrix, b::ZZMatrix, d::ZZRingElem)
ccall((:fmpz_mat_scalar_divexact_fmpz, libflint), Nothing,
(Ref{ZZMatrix}, Ref{ZZMatrix}, Ref{ZZRingElem}), a, b, d)
end

function mul!(a::ZZMatrix, b::ZZMatrix, c::ZZRingElem)
ccall((:fmpz_mat_scalar_mul_fmpz, libflint), Nothing,
(Ref{ZZMatrix}, Ref{ZZMatrix}, Ref{ZZRingElem}), a, b, c)
end

################################################################################
#
################################################################################
Expand Down Expand Up @@ -178,59 +145,6 @@ end
#
################################################################################

@doc raw"""
mod!(M::ZZMatrix, p::ZZRingElem)
Reduces every entry modulo $p$ in-place, i.e. applies the mod function to every entry.
Positive residue system.
"""
function mod!(M::ZZMatrix, p::ZZRingElem)
GC.@preserve M begin
for i = 1:nrows(M)
for j = 1:ncols(M)
z = ccall((:fmpz_mat_entry, libflint), Ptr{ZZRingElem}, (Ref{ZZMatrix}, Int, Int), M, i - 1, j - 1)
ccall((:fmpz_mod, libflint), Nothing, (Ptr{ZZRingElem}, Ptr{ZZRingElem}, Ref{ZZRingElem}), z, z, p)
end
end
end
return nothing
end

@doc raw"""
mod(M::ZZMatrix, p::ZZRingElem) -> ZZMatrix
Reduces every entry modulo $p$, i.e. applies the mod function to every entry.
"""
function mod(M::ZZMatrix, p::ZZRingElem)
N = deepcopy(M)
mod!(N, p)
return N
end

@doc raw"""
mod_sym!(M::ZZMatrix, p::ZZRingElem)
Reduces every entry modulo $p$ in-place, into the symmetric residue system.
"""
function mod_sym!(M::ZZMatrix, B::ZZRingElem)
@assert !iszero(B)
ccall((:fmpz_mat_scalar_smod, libflint), Nothing, (Ref{ZZMatrix}, Ref{ZZMatrix}, Ref{ZZRingElem}), M, M, B)
end
mod_sym!(M::ZZMatrix, B::Integer) = mod_sym!(M, ZZRingElem(B))

@doc raw"""
mod_sym(M::ZZMatrix, p::ZZRingElem) -> ZZMatrix
Reduces every entry modulo $p$ into the symmetric residue system.
"""
function mod_sym(M::ZZMatrix, B::ZZRingElem)
N = zero_matrix(ZZ, nrows(M), ncols(M))
ccall((:fmpz_mat_scalar_smod, libflint), Nothing, (Ref{ZZMatrix}, Ref{ZZMatrix}, Ref{ZZRingElem}), N, M, B)
return N
end
mod_sym(M::ZZMatrix, B::Integer) = mod_sym(M, ZZRingElem(B))


@doc raw"""
mod_sym!(A::Generic.Mat{AbsSimpleNumFieldElem}, m::ZZRingElem)
Expand All @@ -244,140 +158,6 @@ function mod_sym!(A::Generic.Mat{AbsSimpleNumFieldElem}, m::ZZRingElem)
end
end

################################################################################
#
# In-place HNF
#
################################################################################

function hnf!(x::ZZMatrix)
if nrows(x) * ncols(x) > 100
z = hnf(x)
ccall((:fmpz_mat_set, libflint), Nothing, (Ref{ZZMatrix}, Ref{ZZMatrix}), x, z)

return x
end
ccall((:fmpz_mat_hnf, libflint), Nothing, (Ref{ZZMatrix}, Ref{ZZMatrix}), x, x)
return x
end

################################################################################
#
# Smith normal form with trafo
#
################################################################################

#=
g, e,f = gcdx(a, b)
U = [1 0 ; -divexact(b, g)*f 1]*[1 1; 0 1];
V = [e -divexact(b, g) ; f divexact(a, g)];
then U*[ a 0; 0 b] * V = [g 0 ; 0 l]
=#
@doc raw"""
snf_with_transform(A::ZZMatrix, l::Bool = true, r::Bool = true) -> ZZMatrix, ZZMatrix, ZZMatrix
Given some integer matrix $A$, compute the Smith normal form (elementary
divisor normal form) of $A$. If `l` and/ or `r` are true, then the corresponding
left and/ or right transformation matrices are computed as well.
"""
function snf_with_transform(A::ZZMatrix, l::Bool=true, r::Bool=true)
if r
R = identity_matrix(ZZ, ncols(A))
end

if l
L = identity_matrix(ZZ, nrows(A))
end
# TODO: if only one trafo is required, start with the HNF that does not
# compute the trafo
# Rationale: most of the work is on the 1st HNF..
S = deepcopy(A)
while !is_diagonal(S)
if l
S, T = hnf_with_transform(S)
L = T * L
else
S = hnf!(S)
end

if is_diagonal(S)
break
end
if r
S, T = hnf_with_transform(transpose(S))
R = T * R
else
S = hnf!(transpose(S))
end
S = transpose(S)
end
#this is probably not really optimal...
for i = 1:min(nrows(S), ncols(S))
if S[i, i] == 1
continue
end
for j = i+1:min(nrows(S), ncols(S))
if S[j, j] == 0
continue
end
if S[i, i] != 0 && S[j, j] % S[i, i] == 0
continue
end
g, e, f = gcdx(S[i, i], S[j, j])
a = divexact(S[i, i], g)
S[i, i] = g
b = divexact(S[j, j], g)
S[j, j] *= a
if l
# U = [1 0; -b*f 1] * [ 1 1; 0 1] = [1 1; -b*f -b*f+1]
# so row i and j of L will be transformed. We do it naively
# those 2x2 transformations of 2 rows should be a c-primitive
# or at least a Nemo/Hecke primitive
for k = 1:ncols(L)
x = -b * f
# L[i,k], L[j,k] = L[i,k]+L[j,k], x*L[i,k]+(x+1)*L[j,k]
L[i, k], L[j, k] = L[i, k] + L[j, k], x * (L[i, k] + L[j, k]) + L[j, k]
end
end
if r
# V = [e -b ; f a];
# so col i and j of R will be transformed. We do it naively
# careful: at this point, R is still transposed
for k = 1:nrows(R)
R[i, k], R[j, k] = e * R[i, k] + f * R[j, k], -b * R[i, k] + a * R[j, k]
end
end
end
end

# It might be the case that S was diagonal with negative diagonal entries.
for i in 1:min(nrows(S), ncols(S))
if S[i, i] < 0
if l
multiply_row!(L, ZZRingElem(-1), i)
end
S[i, i] = -S[i, i]
end
end

if l
if r
return S, L, transpose(R)
else
# last is dummy
return S, L, L
end
elseif r
# second is dummy
return S, R, transpose(R)
else
# last two are dummy
return S, S, S
end
end


#Returns a positive integer if A[i, j] > b, negative if A[i, j] < b, 0 otherwise
function compare_index(A::ZZMatrix, i::Int, j::Int, b::ZZRingElem)
a = ccall((:fmpz_mat_entry, libflint), Ptr{ZZRingElem}, (Ref{ZZMatrix}, Int, Int), A, i - 1, j - 1)
Expand Down
4 changes: 0 additions & 4 deletions src/HeckeMiscPoly.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -407,10 +407,6 @@ function mulhigh(a::PolyRingElem{T}, b::PolyRingElem{T}, n::Int) where {T}
return mulhigh_n(a, b, degree(a) + degree(b) - n)
end

function (f::AcbPolyRingElem)(x::AcbFieldElem)
return evaluate(f, x)
end

function mod(f::AbstractAlgebra.PolyRingElem{T}, g::AbstractAlgebra.PolyRingElem{T}) where {T<:RingElem}
check_parent(f, g)
if length(g) == 0
Expand Down
6 changes: 3 additions & 3 deletions src/HeckeMoreStuff.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -698,7 +698,7 @@ function Base.Int128(x::ZZRingElem)
return Base.Int128(BigInt(x))
end

## Make zzModRing iteratible
## Make zzModRing iterable

Base.iterate(R::zzModRing) = (zero(R), zero(UInt))

Expand All @@ -716,7 +716,7 @@ Base.eltype(::Type{zzModRing}) = zzModRingElem
Base.IteratorSize(::Type{zzModRing}) = Base.HasLength()
Base.length(R::zzModRing) = R.n

## Make ZZModRing iteratible
## Make ZZModRing iterable

Base.iterate(R::ZZModRing) = (zero(R), zero(ZZRingElem))

Expand All @@ -734,7 +734,7 @@ Base.eltype(::Type{ZZModRing}) = ZZModRingElem
Base.IteratorSize(::Type{ZZModRing}) = Base.HasLength()
Base.length(R::ZZModRing) = Integer(R.n)

## Make fpField iteratible
## Make fpField iterable

Base.iterate(R::fpField) = (zero(R), zero(UInt))

Expand Down
Loading

0 comments on commit c3368ac

Please sign in to comment.