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.

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

Added new function is_unimodular for square integer matrices #1991

Merged
merged 8 commits into from
Jan 16, 2025
Merged
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
1 change: 1 addition & 0 deletions src/Exports.jl
Original file line number Diff line number Diff line change
@@ -373,6 +373,7 @@ export is_uinf
export is_undefined
export is_unit
export is_unknown
export is_unimodular
export is_upper_triangular
export is_zero_entry
export is_zero_row
3 changes: 3 additions & 0 deletions src/Nemo.jl
Original file line number Diff line number Diff line change
@@ -193,6 +193,7 @@ import AbstractAlgebra: set_attribute!
import AbstractAlgebra: Solve
import AbstractAlgebra: terse
import AbstractAlgebra: truncate!
import AbstractAlgebra: add_verbosity_scope

AbstractAlgebra.@include_deprecated_bindings()

@@ -359,6 +360,8 @@ function __init__()

@ccall libflint.flint_set_abort(@cfunction(flint_abort, Nothing, ())::Ptr{Nothing})::Nothing

add_verbosity_scope(:UnimodVerif)

if AbstractAlgebra.should_show_banner() && get(ENV, "NEMO_PRINT_BANNER", "true") != "false"
show_banner()
end
293 changes: 293 additions & 0 deletions src/matrix.jl
Original file line number Diff line number Diff line change
@@ -394,3 +394,296 @@
reduce_mod!(C, B)
return C
end


##################################################################
## Functions related to unimodularity verification
##################################################################

# ***SOURCE ALGORITHM***
# Refined implementation of method from
# Authors: C. Pauderis + A. Storjohann
# Title: Detrministic Unimodularity Certification
# Where: Proc ISSAC 2012, pp. 281--287
# See also: thesis of Colton Pauderis, 2013, Univ of Waterloo
# Deterministic Unimodularity Certification and Applications for Integer Matrices


# add_verbosity_scope(:UnimodVerif) -- see func __init__ in src/Nemo.jl

#######################################################
# Low-level auxiliary functions -- should be elsewhere?
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, I would say most of the coded added here should actually be in src/flint/fmpz_mat.jl

#######################################################

function _det(a::fpMatrix)
a.r < 10 && return lift(det(a))
#_det avoids a copy: det is computed destructively
r = ccall((:_nmod_mat_det, Nemo.libflint), UInt, (Ref{fpMatrix}, ), a)
return r
end

function map_entries!(a::fpMatrix, k::Nemo.fpField, A::ZZMatrix)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The argument order clashes with the general map_entries! signature as defined in AA which is map_entries!(f, dst::MatrixElem{T}, src::MatrixElem{U})

ccall((:nmod_mat_set_mod, Nemo.libflint), Cvoid, (Ref{fpMatrix}, UInt), a, k.n)
ccall((:fmpz_mat_get_nmod_mat, Nemo.libflint), Cvoid, (Ref{fpMatrix}, Ref{ZZMatrix}), a, A)
a.base_ring = k # exploiting that the internal repr is the indep of char
return a
end

function Base.copy!(A::ZZMatrix, B::ZZMatrix)
ccall((:fmpz_mat_set, Nemo.libflint), Cvoid, (Ref{ZZMatrix}, Ref{ZZMatrix}), A, B)
end

function Nemo.sub!(A::ZZMatrix, m::Int)
for i=1:nrows(A)
A_p = Nemo.mat_entry_ptr(A, i, i)
sub!(A_p, A_p, m)
end
return A
end

# No allocations; also faster than *=
function Nemo.mul!(A::ZZMatrix, m::ZZRingElem)
for i=1:nrows(A)
A_p = Nemo.mat_entry_ptr(A, i, 1)
for j=1:ncols(A)
mul!(A_p, A_p, m)
A_p += sizeof(Int)
end
end
return A
end
Comment on lines +445 to +454
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But we already define such a method in src/flint/fmpz_mat.jl:

function mul!(z::ZZMatrixOrPtr, a::ZZMatrixOrPtr, b::ZZRingElemOrPtr)
  @ccall libflint.fmpz_mat_scalar_mul_fmpz(z::Ref{ZZMatrix}, a::Ref{ZZMatrix}, b::Ref{ZZRingElem})::Nothing
  return z
end

(and then generic code in AA also makes mul!(::ZZMatrix, ::ZZRingElem) available.


@doc raw"""
is_unimodular(A::ZZMatrix)
is_unimodular(A::ZZMatrix; algorithm=:auto)

Determine whether `A` is unimodular, i.e. whether `abs(det(A)) == 1`.
The method used is either that of Pauderis--Storjohann or using CRT;
the choice is made based on cost estimates for the two approaches.

The kwarg `algorithm` may be specified to be `:CRT` or `:pauderis_storjohann`
to ensure that the corresponding underlying algorithm is used (after a quick
initial unimodularity check). `:CRT` is better if the matrix is unimodular
and has an inverse containing large entries (or if it is not unimodular);
`:pauderis_storjohann` is asymptotically faster when the matrix is unimodular
and its inverse does not contain large entries, but it requires more space
for intermediate expressions.
"""
function is_unimodular(A::ZZMatrix; algorithm=:auto)
# Call this function when no extra info about the matrix is available.
# It does a preliminary check that det(A) is +/-1 modulo roughly 2^100.
# If so, then delegate the complete check to is_unimodular_given_det_mod_m
if !is_square(A)
throw(ArgumentError("Matrix must be square"))
end
Comment on lines +476 to +478
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While this is fine, I'd recommend to use the @req macro liberally for such tests:

Suggested change
if !is_square(A)
throw(ArgumentError("Matrix must be square"))
end
@req is_square(A) "Matrix must be square"

if !(algorithm in [:auto, :CRT, :pauderis_storjohann])
throw(ArgumentError("algorithm must be one of [:CRT, :pauderis_storjohann, :auto]"))
end
# Deal with two trivial cases
if nrows(A) == 0
return true
end
if nrows(A) == 1
return (abs(A[1,1]) == 1);
end
# Compute det mod M -- may later include some more primes
target_bitsize_modulus = 100 # magic number -- should not be too small!
@vprintln(:UnimodVerif,1,"Quick first check: compute det by crt up to modulus with about $(target_bitsize_modulus) bits")
p::Int = 2^20 # start point of primes for initial CRT trial -- det with modulus >= 2^22 is much slower
M = ZZ(1)
det_mod_m = 0 # 0 means uninitialized, o/w value is 1 or -1
A_mod_p = identity_matrix(Nemo.Native.GF(2), nrows(A))
while nbits(M) <= target_bitsize_modulus
@vprint(:UnimodVerif,2,".")
p = next_prime(p + rand(1:2^10)) # increment by a random step to a prime
ZZmodP = Nemo.Native.GF(p; cached = false, check = false)
map_entries!(A_mod_p, ZZmodP, A)
@vtime :UnimodVerif 2 det_mod_p = Int(_det(A_mod_p))
if det_mod_p > p/2
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While the cost for converting ints to float etc. here, it could indeed still be avoided here, e.g. via use of >>

Suggested change
if det_mod_p > p/2
if det_mod_p > p>>1

det_mod_p -= p
end
if abs(det_mod_p) != 1
return false # obviously not unimodular
end
if det_mod_m != 0 && det_mod_m != det_mod_p
return false # obviously not unimodular

Check warning on line 509 in src/matrix.jl

Codecov / codecov/patch

src/matrix.jl#L509

Added line #L509 was not covered by tests
end
det_mod_m = det_mod_p
mul!(M, M, p)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
mul!(M, M, p)
M = mul!(M, p)

end
return is_unimodular_given_det_mod_m(A, det_mod_m, M; algorithm=algorithm)
end #function


# THIS FUNCTION NOT OFFICIALLY DOCUMENTED -- not intended for public use.
function is_unimodular_given_det_mod_m(A::ZZMatrix, det_mod_m::Int, M::ZZRingElem; algorithm=:auto)
Comment on lines +518 to +519
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Instead of a code comment which nobody will see, I suggest we start using a prefix like _ to indicate internal functions:

Suggested change
# THIS FUNCTION NOT OFFICIALLY DOCUMENTED -- not intended for public use.
function is_unimodular_given_det_mod_m(A::ZZMatrix, det_mod_m::Int, M::ZZRingElem; algorithm=:auto)
function _is_unimodular_given_det_mod_m(A::ZZMatrix, det_mod_m::Int, M::ZZRingElem; algorithm=:auto)

# This UI is convenient for our sophisticated determinant algorithm.
# We estimate cost of computing det by CRT and cost of doing Pauderis-Storjohann
# unimodular verification then call whichever is likely to be faster
# (but we avoid P-S if the space estimate is too large).
# M is a modulus satisfying M > 2^30;
# det_mod_m is det(A) mod M [!not checked!]: it is either +1 or -1.
is_square(A) || throw(ArgumentError("Matrix must be square"))
abs(det_mod_m) == 1 || throw(ArgumentError("det_mod_m must be +1 or -1"))
M > 2^30 || throw(ArgumentError("modulus must be greater than 2^30"))
(algorithm in [:auto, :CRT, :pauderis_storjohann]) || throw(ArgumentError("algorithm must be one of [:CRT, :pauderis_storjohann, :auto]"))
# Deal with two trivial cases -- does this make sense here???
nrows(A) == 0 && return true
nrows(A) == 1 && return (abs(A[1,1]) == 1)
@vprintln(:UnimodVerif,1,"is_unimodular_given_det_mod_m starting")
if algorithm == :pauderis_storjohann
@vprintln(:UnimodVerif,1,"User specified Pauderis_Storjohann --> delegating")
return is_unimodular_Pauderis_Storjohann_Hensel(A)
end
n = ncols(A)
Hrow = hadamard_bound2(A)
Hcol = hadamard_bound2(transpose(A))
DetBoundSq = min(Hrow, Hcol)
Hbits = 1+div(nbits(DetBoundSq),2)
@vprintln(:UnimodVerif,1,"Hadamard bound in bits $(Hbits)")
if algorithm == :auto
# Estimate whether better to "climb out" with CRT or use Pauderis-Storjohann
CRT_speed_factor = 20.0 # scale factor is JUST A GUESS !!! (how much faster CRT is compared to P-S)
total_entry_size = sum(nbits,A)
Estimate_CRT = ceil(((Hbits - nbits(M))*(2.5+total_entry_size/n^3))/CRT_speed_factor) # total_entry_size/n is proportional to cost of reducing entries mod p
EntrySize = maximum(abs, A)
entry_size_bits = maximum(nbits, A)
@vprintln(:UnimodVerif,1,"entry_size_bits = $(entry_size_bits) log2(EntrySize) = $(ceil(log2(EntrySize)))")
e = max(16, Int(2+ceil(2*log2(n)+log2(EntrySize)))) # see Condition in Fig 1, p.284 of Pauderis-Storjohann
Estimate_PS = ceil(e*log(e)) # assumes soft-linear ZZ arith
@vprintln(:UnimodVerif,1,"Est_CRT = $(Estimate_CRT) and Est_PS = $(Estimate_PS)")
Space_estimate_PS = ZZ(n)^2 * e # nbits of values, ignoring mem mgt overhead
@vprintln(:UnimodVerif,1,"Space est PS: $(Space_estimate_PS) [might be bytes]")
if Estimate_PS < Estimate_CRT #=&& Space_estimate_PS < 2^32=#
# Expected RAM requirement is not excessive, and
# Pauderis-Storjohann is likely faster, so delegate to UniCertZ_hensel
return is_unimodular_Pauderis_Storjohann_Hensel(A)
end
end
if algorithm == :CRT
@vprintln(:UnimodVerif,1,"User specified CRT --> delegating")
end
@vprintln(:UnimodVerif,1,"UnimodularVerification: CRT loop")
# Climb out with CRT:
# in CRT loop start with 22 bit primes; if we reach threshold (empirical), jump to much larger primes.
p = 2^21; stride = 2^10; threshold = 5000000;
A_mod_p = zero_matrix(Nemo.Native.GF(2), nrows(A), ncols(A))
while nbits(M) < Hbits
@vprint(:UnimodVerif,2,".")
# Next lines increment p (by random step up to stride) to a prime not dividing M
if p > threshold && p < threshold+stride
@vprint(:UnimodVerif,2,"!!")
p = 2^56; stride = 2^25; # p = 2^56 is a good choice on my standard 64-bit machine
continue

Check warning on line 577 in src/matrix.jl

Codecov / codecov/patch

src/matrix.jl#L575-L577

Added lines #L575 - L577 were not covered by tests
end
# advance to another prime which does not divide M:
while true
p = next_prime(p+rand(1:stride))
if !is_divisible_by(M,p)
break
end
end

Check warning on line 585 in src/matrix.jl

Codecov / codecov/patch

src/matrix.jl#L585

Added line #L585 was not covered by tests
ZZmodP = Nemo.Native.GF(p; cached = false, check = false)
map_entries!(A_mod_p, ZZmodP, A)
det_mod_p = _det(A_mod_p)
if det_mod_p > p/2
det_mod_p -= p
end
if abs(det_mod_p) != 1 || det_mod_m != det_mod_p
return false # obviously not unimodular

Check warning on line 593 in src/matrix.jl

Codecov / codecov/patch

src/matrix.jl#L593

Added line #L593 was not covered by tests
end
M *= p
end
return true
end

# THIS FUNCTION NOT OFFICIALLY DOCUMENTED -- not intended for public use.
# Pauderis-Storjohann unimodular verification/certification.
# We use Hensel lifting to obtain the modular inverse B0
# Seems to be faster than the CRT prototype (not incl. here)
# VERBOSITY via :UnimodVerif
function is_unimodular_Pauderis_Storjohann_Hensel(A::ZZMatrix)
n = nrows(A)
# assume ncols == n
EntrySize = maximum(abs, A)
entry_size_bits = maximum(nbits, A)
e = max(16, 2+ceil(Int, 2*log2(n)+entry_size_bits))
ModulusLWB = ZZ(2)^e
@vprintln(:UnimodVerif,1,"ModulusLWB is 2^$(e)")
# Now look for prime: about 25 bits, and such that p^(2^sthg) is
# just slightly bigger than ModulusLWB. Bit-size of prime matters little.
root_order = exp2(ceil(log2(e/25)))
j = ceil(exp2(e/root_order))
p = next_prime(Int(j))
@vprintln(:UnimodVerif,1,"Hensel prime is p = $(p)")
m = ZZ(p)
ZZmodP = Nemo.Native.GF(p)
MatModP = matrix_space(ZZmodP,n,n)
B0 = lift(inv(MatModP(A)))
mod_sym!(B0, m)
delta = identity_matrix(base_ring(A), n)
A_mod_m = identity_matrix(base_ring(A), n)
@vprintln(:UnimodVerif,1,"Starting stage 1 lifting")
while (m < ModulusLWB)
@vprintln(:UnimodVerif,2,"Loop 1: nbits(m) = $(nbits(m))");
copy!(A_mod_m, A)
mm = m^2
mod_sym!(A_mod_m,mm)
@vtime :UnimodVerif 2 mul!(delta, A_mod_m, B0)
sub!(delta, 1)
divexact!(delta, m) # to make matrix product in line below cheaper
@vtime :UnimodVerif 2 mul!(A_mod_m, B0, delta)
mul!(A_mod_m, m)
sub!(B0, B0, A_mod_m)
m = mm
mod_sym!(B0, m)
end
@vprintln(:UnimodVerif,1,"Got stage 1 modular inverse: modulus has $(nbits(m)) bits; value is $(p)^$(Int(round(log2(m)/log2(p))))")

# Hadamard bound brings little benefit; anyway, almost never lift all the way
### H = min(hadamard_bound2(A), hadamard_bound2(transpose(A)))
### println("nbits(H) = $(nbits(H))");
# Compute max num iters in k
k = (n-1)*(entry_size_bits+log2(n)/2) - (2*log2(n) + entry_size_bits -1)
k = 2+nbits(ceil(Int,k/log2(m)))
@vprintln(:UnimodVerif,1,"Stage 2: max num iters is k=$(k)")

ZZmodM,_ = residue_ring(ZZ,m; cached = false) # m is probably NOT prime
## @assert is_one(lift(MatModM(B0*A)))
# Originally: R = (I-A*B0)/m
R = A*B0
sub!(R, 1) #this is -R, but is_zero test is the same, and the loop starts by squaring
divexact!(R, m)
if is_zero(R)
return true
end

#ZZMatrix and ZZModMatrix are the same type. The modulus is
#always passed in as an extra argument when needed...
#mul! for ZZModMatrix is mul, followed by reduce.

B0_modm = deepcopy(B0)
mod_sym!(B0_modm, m)

R_bar = zero_matrix(ZZ, n, n)
M = zero_matrix(ZZ, n, n)

@vprintln(:UnimodVerif,1,"Starting stage 2 lifting")
for i in 0:k-1
@vprintln(:UnimodVerif,1,"Stage 2 loop: i=$(i)")

mul!(R_bar, R, R)
#Next 3 lines do: M = lift(B0_modm*MatModM(R_bar));
mod_sym!(R_bar, m)
mul!(M, B0_modm, R_bar)
mod_sym!(M, m)
# Next 3 lines do: R = (R_bar - A*M)/m; but with less memory allocation
mul!(R, A, M)
sub!(R, R_bar, R)
divexact!(R, m)
if is_zero(R)
return true
end
end
return false
end
57 changes: 57 additions & 0 deletions test/matrix-test.jl
Original file line number Diff line number Diff line change
@@ -101,3 +101,60 @@ end

# eigenvalues_simple not defined for integer matrices, so not tested.
end


@testset "is_unimodular" begin

# Some trivial inputs
@test is_unimodular(matrix(ZZ,0,0,[]))
@test is_unimodular(matrix(ZZ,1,1,[1]))
@test is_unimodular(matrix(ZZ,1,1,[-1]))
@test !is_unimodular(matrix(ZZ,1,1,[0]))

@test is_unimodular(identity_matrix(ZZ, 11))
@test !is_unimodular(2*identity_matrix(ZZ, 11))
@test !is_unimodular(zero_matrix(ZZ, 11,11))

# This test increases coverage slightly
N = 1+prod(filter(is_prime, ZZ(2)^20:(ZZ(2)^20+ZZ(2)^11)))
@test !is_unimodular(N*identity_matrix(ZZ, 11))


@test_throws ArgumentError is_unimodular(matrix(ZZ,0,1,[]))
@test_throws ArgumentError is_unimodular(matrix(ZZ,0,0,[]); algorithm=:WRONG)

U = matrix(ZZ, 3,3, [ 3, 22, 46, 5, 35, 73, 4, 27, 56])
M = matrix(ZZ, 3,3, [-3, 22, 46, 5, 35, 73, 4, 27, 56])

# U6 is "small" and triggers the 2nd loop in P-S method
U6 = matrix(ZZ, 6,6,
[1, -76, 87, -57, -41, 999950,
0, 1, 0, 0, 0, 0,
0, 999942, 1, 0, 0, 0,
0, -48, 999963, 1, 0, 0,
0, 91, -24, 1000012, 1, 0,
0, -2, 71, 77, 999925, 1])


# Need larger values of k to increase coverage
for k in 0:25
@test is_unimodular(U^k)
@test is_unimodular(U^k; algorithm=:CRT)
@test is_unimodular(U^k; algorithm=:pauderis_storjohann)
@test is_unimodular(U^k; algorithm=:auto)
@test is_unimodular(-U^k)
@test is_unimodular(-U^k; algorithm=:CRT)
@test is_unimodular(-U^k; algorithm=:pauderis_storjohann)
@test is_unimodular(-U^k; algorithm=:auto)
end

@test is_unimodular(U6)
@test is_unimodular(U6; algorithm=:CRT)
@test is_unimodular(U6; algorithm=:pauderis_storjohann)
@test is_unimodular(U6; algorithm=:auto)

@test !is_unimodular(M)
@test !is_unimodular(M; algorithm=:CRT)
@test !is_unimodular(M; algorithm=:pauderis_storjohann)
@test !is_unimodular(M; algorithm=:auto)
end