Merge pull request #10 from ChenZhao44/chen/smith-normal-form

Smith normal form and triangularization for Mod matrices

Project.toml

@@ -4,15 +4,17 @@ version = "0.1.10"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-SimplePolynomials = "cc47b68c-3164-5771-a705-2bc0097375a0"
 Mods = "7475f97c-0381-53b1-977b-4c60186c8d62"
 Permutations = "2ae35dd2-176d-5d53-8349-f30d82d94d4f"
+Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
+SimplePolynomials = "cc47b68c-3164-5771-a705-2bc0097375a0"
 
 [compat]
-julia = "1"
-SimplePolynomials = "0.0.1, 0.0.2, 0.0.3, 0.0.4, 0.0.5, 0.0.6, 0.1, 0.2"
 Mods = "1"
 Permutations = "0.4"
+Primes = "0.5"
+SimplePolynomials = "0.0.1, 0.0.2, 0.0.3, 0.0.4, 0.0.5, 0.0.6, 0.1, 0.2"
+julia = "1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/LinearAlgebraX.jl

@@ -27,6 +27,7 @@ include("nullspacex.jl")
 include("char_poly.jl")
 include("homogeneous.jl")
 include("perm.jl")
-
+include("mod_smith_normal_form.jl")
+include("mod_triangular.jl")
 
 end # module

src/detx.jl

@@ -28,7 +28,7 @@ end
 function detx(A::AbstractMatrix{T}) where {T}
     # @info "Using detx! on matrix of type $T"
     r, c = size(A)
-    B = Matrix{Any}(undef, r, c)
+    B = Matrix{T}(undef, r, c)
     for i = 1:r
         for j = 1:c
             B[i, j] = A[i, j]
@@ -89,3 +89,11 @@ function detx!(A::AbstractMatrix{T}) where {T}
     return detx!(A[2:end, 2:end]) * factor
 
 end
+
+function detx!(A::AbstractMatrix{<:Mod})
+    U, row_ops = upper_triangular!(A)
+    det_C = prod(detx(op) for op in row_ops)
+
+    # U = C * A
+    return inv(det_C) * prod(diag(U))
+end

src/mod_smith_normal_form.jl

@@ -0,0 +1,289 @@
+using Primes
+export smith_normal_form
+
+abstract type AbstractMatrixOperation{R} end
+abstract type AbstractRowOperation{R} <: AbstractMatrixOperation{R} end
+abstract type AbstractColumnOperation{R} <: AbstractMatrixOperation{R} end
+struct RowSwap{R} <: AbstractRowOperation{R}
+    i::Int
+    j::Int
+end
+struct RowAddMult{R} <: AbstractRowOperation{R}
+    i::Int
+    j::Int
+    c::R
+end
+struct RowScale{R} <: AbstractRowOperation{R}
+    i::Int
+    u::R
+end
+struct RowSmith{R} <: AbstractRowOperation{R}
+    i::Int
+    j::Int
+    s::R
+    u::R
+    t::R
+    v::R
+end
+struct ColumnSwap{R} <: AbstractColumnOperation{R}
+    i::Int
+    j::Int
+end
+struct ColumnAddMult{R} <: AbstractColumnOperation{R}
+    i::Int
+    j::Int
+    c::R
+end
+struct ColumnScale{R} <: AbstractColumnOperation{R}
+    i::Int
+    u::R
+end
+struct ColumnSmith{R} <: AbstractColumnOperation{R}
+    i::Int
+    j::Int
+    s::R
+    u::R
+    t::R
+    v::R
+end
+
+detx(op::RowSwap{R}) where {R} = (op.i == op.j) ? one(R) : -one(R)
+detx(::RowAddMult{R}) where {R} = one(R)
+detx(op::RowScale{R}) where {R} = op.u
+detx(op::RowSmith{R}) where {R} = op.s * op.v - op.u * op.t
+detx(op::ColumnSwap{R}) where {R} = (op.i == op.j) ? one(R) : -one(R)
+detx(::ColumnAddMult{R}) where {R} = one(R)
+detx(op::ColumnScale{R}) where {R} = op.u
+detx(op::ColumnSmith{R}) where {R} = op.s * op.v - op.u * op.t
+
+(r::AbstractMatrixOperation{R})(M::AbstractMatrix{R}) where R = apply_matrix_operation!(M, r)
+
+function apply_matrix_operation!(M::AbstractMatrix{R}, op::RowSwap{R}) where {R}
+    M[op.i,:], M[op.j,:] = M[op.j,:], M[op.i,:]
+    return M
+end
+function apply_matrix_operation!(M::AbstractMatrix{R}, op::RowAddMult{R}) where {R}
+    M[op.j,:] += op.c * M[op.i,:]
+    return M
+end
+function apply_matrix_operation!(M::AbstractMatrix{R}, op::RowScale{R}) where {R}
+    M[op.i,:] *= op.u
+    return M
+end
+function apply_matrix_operation!(M::AbstractMatrix{R}, op::RowSmith{R}) where {R}
+    M[[op.i, op.j],:] = [op.s op.t; op.u op.v] * M[[op.i, op.j],:]
+    return M
+end
+function apply_matrix_operation!(M::AbstractMatrix{R}, op::ColumnSwap{R}) where {R}
+    M[:,op.i], M[:,op.j] = M[:,op.j], M[:,op.i]
+    return M
+end
+function apply_matrix_operation!(M::AbstractMatrix{R}, op::ColumnAddMult{R}) where {R}
+    M[:,op.j] += op.c * M[:,op.i]
+    return M
+end
+function apply_matrix_operation!(M::AbstractMatrix{R}, op::ColumnScale{R}) where {R}
+    M[:,op.i] *= op.u
+    return M
+end
+function apply_matrix_operation!(M::AbstractMatrix{R}, op::ColumnSmith{R}) where {R}
+    M[:,[op.i, op.j]] = M[:,[op.i, op.j]] * [op.s op.u; op.t op.v]
+    return M
+end
+
+"""
+smith_coeff(a::Mod{N, T}, b::Mod{N, T}[, prime_powers]) where {N, T}
+
+Return a invertable matrix `M` such that `M*[a; b] = [g; 0]`, 
+where `g` is a gcd of `a` and `b`.
+"""
+function smith_coeff(a::Mod{N, T}, b::Mod{N, T}, prime_powers::Union{Nothing, Vector{<:Integer}} = nothing) where {N, T}
+    prime_powers = isnothing(prime_powers) ? [p^d for (p, d) in factor(N)] : prime_powers
+    return smith_coeff(T, value(a), value(b), prime_powers, N)
+end
+
+function smith_coeff(T, a::Integer, b::Integer, prime_powers::Vector{<:Integer}, N::Integer)
+    as = a .% prime_powers
+    bs = b .% prime_powers
+    zas = iszero.(as)
+    zbs = iszero.(bs)
+    za = Mod{N, T}(CRT(BigInt, zas, prime_powers))
+    zb = Mod{N, T}(CRT(BigInt, zbs, prime_powers))
+    zabs = zas .* zbs
+    as = as + zas .* bs
+    bs = bs + zbs .* as
+    pas = gcd.(as, prime_powers)
+    pbs = gcd.(bs, prime_powers)
+    cas = as .÷ pas
+    cbs = bs .÷ pbs
+    pas = pas .* ((!).(zabs))
+    pbs = pbs .* ((!).(zabs))
+    cas = cas + zabs
+    cbs = cbs + zabs
+    @assert cas .* pas == as
+    @assert cbs .* pbs == bs
+    ca = Mod{N, T}(CRT(BigInt, cas, prime_powers))
+    cb = Mod{N, T}(CRT(BigInt, cbs, prime_powers))
+    αs = pbs .÷ (pas + zabs)
+    βs = pas .÷ (pbs + zabs)
+    γs = iszero.(αs)
+    α = Mod{N, T}(CRT(BigInt, αs, prime_powers))
+    β = Mod{N, T}(CRT(BigInt, βs, prime_powers))
+    γ = Mod{N, T}(CRT(BigInt, γs, prime_powers))
+    M_za = Mod{N, T}[1 za; 0 1]
+    M_zb = Mod{N, T}[1 0; zb 1]
+    M_inv_ca_cb = Mod{N, T}[inv(ca) 0; 0 inv(cb)]
+    M_α = Mod{N, T}[1 0; -α 1]
+    M_β = Mod{N, T}[1 -β; 0 1]
+    M_plus_γ = Mod{N, T}[1 γ; 0 1]
+    M_minus_γ = Mod{N, T}[1 -γ; 0 1]
+    M = M_minus_γ * M_plus_γ * M_β * M_α * M_inv_ca_cb * M_zb * M_za
+    res = Mod{N, T}[a; b]
+    res = M * res
+    @assert iszero(res[2])
+    return M
+end
+
+"""
+    divisible(a::Mod{N, T}, b::Mod{N, T}[, prime_powers]) where {N, T}
+
+Check if a | b mod N, where N is the product of `prime_powers`.
+"""
+function divisible(a::Mod{N, T}, b::Mod{N, T}, prime_powers::Union{Nothing, Vector{<:Integer}} = nothing) where {N, T}
+    prime_powers = isnothing(prime_powers) ? [p^d for (p, d) in factor(N)] : prime_powers
+    _, r = mod_divrem(T, value(b), value(a), prime_powers, N)
+    return iszero(r)
+end
+
+function Base.divrem(a::Mod{N, T}, b::Mod{N, T}, prime_powers::Union{Nothing, Vector{<:Integer}} = nothing) where {N, T}
+    prime_powers = isnothing(prime_powers) ? [p^d for (p, d) in factor(N)] : prime_powers
+    return mod_divrem(T, value(a), value(b), prime_powers, N)
+end
+
+function mod_divrem(T, a::Integer, b::Integer, prime_powers::Vector{<:Integer}, N::Integer)
+    as = a .% prime_powers
+    bs = b .% prime_powers
+    pas = gcd.(as, prime_powers)
+    pbs = gcd.(bs, prime_powers)
+
+    # `a` is not divisible by `b`.
+    !all(iszero, rem.(pas, pbs)) && return Mod{N, T}(0), Mod{N, T}(a)
+
+    # `a` is divisible by `b`: `a = b * c`.
+    zbs = iszero.(bs)
+    αs = as .÷ pas 
+    βs = bs .÷ pbs + zbs
+    inv_βs = [gcdx(β, p)[2] for (β, p) in zip(βs, prime_powers)]
+    cs = inv_βs .* αs .* (pas .÷ pbs)
+    cs[zbs] .= 0
+    # @assert all(rem.(ds, prime_powers) .== 0)
+    c = Mod{N, T}(CRT(BigInt, cs, prime_powers))
+    return c, Mod{N, T}(0)
+end
+
+"""
+    smith_normal_form(M::Matrix{Mod{N, T}}[, prime_powers]) where {N, T}
+
+Compute the Smith normal form of `M` over `Z_N`. It returns matrices `S`, `U`, and `V` 
+such that `S = U * M * V` and `S` is in Smith normal form of `M`.
+"""
+function smith_normal_form(M::AbstractMatrix{Mod{N, T}}, prime_powers::Union{Nothing, Vector{<:Integer}} = nothing) where {N, T}
+    m, n = size(M)
+    S, row_ops, col_ops = smith_normal_form!(copy(M), prime_powers)
+    U = eye(Mod{N, T}, m)
+    V = eye(Mod{N, T}, n)
+    for op in row_ops
+        op(U)
+    end
+    for op in col_ops
+        op(V)
+    end
+    return S, U, V
+end
+
+function smith_normal_form!(M::AbstractMatrix{Mod{N, T}}, prime_powers::Union{Nothing, Vector{<:Integer}} = nothing) where {N, T}
+    prime_powers = isnothing(prime_powers) ? [p^d for (p, d) in factor(N)] : prime_powers
+    R = Mod{N, T}
+    m, n = size(M)
+    row_ops = AbstractMatrixOperation{R}[]
+    col_ops = AbstractMatrixOperation{R}[]
+    for i = 1:min(m, n)
+        smith_elimination!(M, row_ops, col_ops, i, prime_powers)
+    end
+    diag_normalized = false
+    while diag_normalized
+        diag_normalized = true
+        for i = 2:min(m, n)
+            if !divisible(M[i-1,i-1], M[i,i])
+                diag_normalized = false
+                op = ColumnAddMult(i, i-1, one(R))
+                push!(col_ops, op)
+                op(M)
+                smith_elimination!(M, row_ops, col_ops, i-1, prime_powers)
+                break
+            end
+        end
+    end
+    return M, row_ops, col_ops
+end
+
+function smith_elimination!(M::AbstractMatrix{Mod{N, T}}, row_ops, col_ops, i, prime_powers::Vector{<:Integer}) where {N, T}
+    R = Mod{N, T}
+    m, n = size(M)
+    @assert 1 <= i <= min(m, n)
+    finished = false
+    while !finished
+        finished = true
+
+        # Eliminate the i-th column.
+        if iszero(M[i,i])
+            for j = i+1:m
+                if !iszero(M[j,i])
+                    op = RowSwap{R}(i, j)
+                    push!(row_ops, op)
+                    op(M)
+                    break
+                end
+            end
+        end
+        for j = i+1:m
+            c, r = divrem(M[j,i], M[i,i])
+            if iszero(r)
+                op = RowAddMult{R}(i, j, -c)
+                push!(row_ops, op)
+                op(M)
+            else
+                op = RowSmith{R}(i, j, smith_coeff(M[i,i], M[j,i], prime_powers)[:]...)
+                push!(row_ops, op)
+                op(M)
+                finished = false
+            end
+        end
+
+        # Eliminate the i-th row.
+        if iszero(M[i,i])
+            for j = i+1:n
+                if !iszero(M[i,j])
+                    op = ColumnSwap{R}(i, j)
+                    push!(col_ops, op)
+                    op(M)
+                    break
+                end
+            end
+        end
+        for j = i+1:n
+            c, r = divrem(M[i,j], M[i,i])
+            if iszero(r)
+                op = ColumnAddMult{R}(i, j, -c)
+                push!(col_ops, op)
+                op(M)
+            else
+                op = ColumnSmith{R}(i, j, smith_coeff(M[i,i], M[i,j], prime_powers)[:]...)
+                push!(col_ops, op)
+                op(M)
+                finished = false
+            end
+        end
+    end
+    return M, row_ops, col_ops
+end