From aebe4464d542818b38dac61b624783d6e4a3856d Mon Sep 17 00:00:00 2001 From: John Abbott Date: Tue, 14 Jan 2025 15:57:23 +0100 Subject: [PATCH 1/8] Added new func is_unimodular for square integer matrices --- src/Exports.jl | 1 + src/Nemo.jl | 3 + src/matrix.jl | 298 ++++++++++++++++++++++++++++++++++++++++++++ test/matrix-test.jl | 37 ++++++ 4 files changed, 339 insertions(+) diff --git a/src/Exports.jl b/src/Exports.jl index ee29b244b5..0de7b6dfc7 100644 --- a/src/Exports.jl +++ b/src/Exports.jl @@ -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 diff --git a/src/Nemo.jl b/src/Nemo.jl index 1198295905..ba45862551 100644 --- a/src/Nemo.jl +++ b/src/Nemo.jl @@ -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 diff --git a/src/matrix.jl b/src/matrix.jl index 448700c03c..a740add563 100644 --- a/src/matrix.jl +++ b/src/matrix.jl @@ -394,3 +394,301 @@ function reduce_mod(A::MatElem{T}, B::MatElem{T}) where {T<:FieldElem} 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? +####################################################### + +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) + 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 + +@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) + error("Matrix must be square"); + end + if !(algorithm in [:auto, :CRT, :pauderis_storjohann]) + error("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 + 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 + end + det_mod_m = det_mod_p; + mul!(M, 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) + # 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. + if !is_square(A) + error("Matrix must be square"); + end + if abs(det_mod_m) != 1 + error("det_mod_m must be +1 or -1") + end + if M <= 2^30 + error("modulus must be greater than 2^30") + end + if !(algorithm in [:auto, :CRT, :pauderis_storjohann]) + error("algorithm must be one of [:CRT, :pauderis_storjohann, :auto]") + end + # Deal with two trivial cases -- does this make sense here??? + if nrows(A) == 0 + return true; + end + if nrows(A) == 1 + return (abs(A[1,1]) == 1); + end + @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; + end + while true p = next_prime(p+rand(1:stride)); if !is_divisible_by(M,p) break; end; end + 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 + 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))); + #=Hecke.=#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 + #=Hecke.=#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 + #=Hecke.=#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) + #=Hecke.=#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)); + #=Hecke.=#mod_sym!(R_bar, m) + mul!(M, B0_modm, R_bar) + #=Hecke.=#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 diff --git a/test/matrix-test.jl b/test/matrix-test.jl index 17ce1f7bdc..8eb5f8fb5c 100644 --- a/test/matrix-test.jl +++ b/test/matrix-test.jl @@ -101,3 +101,40 @@ end # eigenvalues_simple not defined for integer matrices, so not tested. end + + +@testset "is_unimodular" begin + 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]) + + + for k in 0:5 + @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 From 19926164fd56806b76e49597d3f41ff0c50d91cf Mon Sep 17 00:00:00 2001 From: John Abbott Date: Tue, 14 Jan 2025 17:14:17 +0100 Subject: [PATCH 2/8] More test coverage --- test/matrix-test.jl | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/test/matrix-test.jl b/test/matrix-test.jl index 8eb5f8fb5c..f12f224472 100644 --- a/test/matrix-test.jl +++ b/test/matrix-test.jl @@ -104,6 +104,16 @@ 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_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]) @@ -117,7 +127,7 @@ end 0, -2, 71, 77, 999925, 1]) - for k in 0:5 + for k in 0:9 @test is_unimodular(U^k) @test is_unimodular(U^k; algorithm=:CRT) @test is_unimodular(U^k; algorithm=:pauderis_storjohann) From 564c2824403c063dd81b2be70736aa2f224d6e3b Mon Sep 17 00:00:00 2001 From: John Abbott Date: Wed, 15 Jan 2025 10:48:48 +0100 Subject: [PATCH 3/8] Fixed syntax error --- test/matrix-test.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/matrix-test.jl b/test/matrix-test.jl index f12f224472..17d4710803 100644 --- a/test/matrix-test.jl +++ b/test/matrix-test.jl @@ -106,7 +106,7 @@ end @testset "is_unimodular" begin # Some trivial inputs - @test is_unimodular(matrix(ZZ,0,0,[]) + @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])) From 2062d74a08d692a2f09a8d2c8e86575d3a85b3f0 Mon Sep 17 00:00:00 2001 From: John Abbott Date: Wed, 15 Jan 2025 11:15:03 +0100 Subject: [PATCH 4/8] Removed lots of semicolons (o/w someone will kick up a fuss) --- src/matrix.jl | 169 ++++++++++++++++++++++++++------------------------ 1 file changed, 88 insertions(+), 81 deletions(-) diff --git a/src/matrix.jl b/src/matrix.jl index a740add563..7ce03f0cc5 100644 --- a/src/matrix.jl +++ b/src/matrix.jl @@ -474,40 +474,41 @@ function is_unimodular(A::ZZMatrix; algorithm=:auto) # 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) - error("Matrix must be square"); + throw(ArgumentError("Matrix must be square")) end if !(algorithm in [:auto, :CRT, :pauderis_storjohann]) - error("algorithm must be one of [:CRT, :pauderis_storjohann, :auto]") + throw(ArgumentError("algorithm must be one of [:CRT, :pauderis_storjohann, :auto]")) end # Deal with two trivial cases if nrows(A) == 0 - return true; + 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! + 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 + 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)); + @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 - det_mod_p -= p; + det_mod_p -= p end if abs(det_mod_p) != 1 - return false; # obviously not unimodular + return false # obviously not unimodular end if det_mod_m != 0 && det_mod_m != det_mod_p - return false; # obviously not unimodular + return false # obviously not unimodular end - det_mod_m = det_mod_p; + det_mod_m = det_mod_p mul!(M, M, p) end return is_unimodular_given_det_mod_m(A, det_mod_m, M; algorithm=algorithm) @@ -523,48 +524,48 @@ function is_unimodular_given_det_mod_m(A::ZZMatrix, det_mod_m::Int, M::ZZRingEle # M is a modulus satisfying M > 2^30; # det_mod_m is det(A) mod M [!not checked!]: it is either +1 or -1. if !is_square(A) - error("Matrix must be square"); + throw(ArgumentError("Matrix must be square")) end if abs(det_mod_m) != 1 - error("det_mod_m must be +1 or -1") + throw(ArgumentError("det_mod_m must be +1 or -1")) end if M <= 2^30 - error("modulus must be greater than 2^30") + throw(ArgumentError("modulus must be greater than 2^30")) end if !(algorithm in [:auto, :CRT, :pauderis_storjohann]) - error("algorithm must be one of [:CRT, :pauderis_storjohann, :auto]") + throw(ArgumentError("algorithm must be one of [:CRT, :pauderis_storjohann, :auto]")) end # Deal with two trivial cases -- does this make sense here??? if nrows(A) == 0 - return true; + return true end if nrows(A) == 1 - return (abs(A[1,1]) == 1); + return (abs(A[1,1]) == 1) end - @vprintln(:UnimodVerif,1,"is_unimodular_given_det_mod_m starting"); + @vprintln(:UnimodVerif,1,"is_unimodular_given_det_mod_m starting") if algorithm == :pauderis_storjohann - @vprintln(:UnimodVerif,1,"User specified Pauderis_Storjohann --> delegating"); + @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)"); + 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]"); + 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 @@ -572,34 +573,40 @@ function is_unimodular_given_det_mod_m(A::ZZMatrix, det_mod_m::Int, M::ZZRingEle end end if algorithm == :CRT - @vprintln(:UnimodVerif,1,"User specified CRT --> delegating"); + @vprintln(:UnimodVerif,1,"User specified CRT --> delegating") end - @vprintln(:UnimodVerif,1,"UnimodularVerification: CRT loop"); + @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,"."); + @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,"!!"); + @vprint(:UnimodVerif,2,"!!") p = 2^56; stride = 2^25; # p = 2^56 is a good choice on my standard 64-bit machine - continue; + continue + 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 - while true p = next_prime(p+rand(1:stride)); if !is_divisible_by(M,p) break; end; end - ZZmodP = Nemo.Native.GF(p; cached = false, check = false); - map_entries!(A_mod_p, ZZmodP, A); - det_mod_p = _det(A_mod_p); + 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; + det_mod_p -= p end if abs(det_mod_p) != 1 || det_mod_m != det_mod_p - return false; # obviously not unimodular + return false # obviously not unimodular end M *= p end - return true; + return true end # THIS FUNCTION NOT OFFICIALLY DOCUMENTED -- not intended for public use. @@ -608,32 +615,32 @@ end # 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); + 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)"); + 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))); - #=Hecke.=#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"); + 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 - #=Hecke.=#mod_sym!(A_mod_m,mm); + 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 @@ -641,26 +648,26 @@ function is_unimodular_Pauderis_Storjohann_Hensel(A::ZZMatrix) mul!(A_mod_m, m) sub!(B0, B0, A_mod_m) m = mm - #=Hecke.=#mod_sym!(B0, m); + 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))))"); + @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 = (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)"); + @vprintln(:UnimodVerif,1,"Stage 2: max num iters is k=$(k)") - ZZmodM,_ = residue_ring(ZZ,m; cached = false); # m is probably NOT prime + 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; + return true end #ZZMatrix and ZZModMatrix are the same type. The modulus is @@ -668,20 +675,20 @@ function is_unimodular_Pauderis_Storjohann_Hensel(A::ZZMatrix) #mul! for ZZModMatrix is mul, followed by reduce. B0_modm = deepcopy(B0) - #=Hecke.=#mod_sym!(B0_modm, m) + 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"); + @vprintln(:UnimodVerif,1,"Starting stage 2 lifting") for i in 0:k-1 - @vprintln(:UnimodVerif,1,"Stage 2 loop: i=$(i)"); + @vprintln(:UnimodVerif,1,"Stage 2 loop: i=$(i)") mul!(R_bar, R, R) #Next 3 lines do: M = lift(B0_modm*MatModM(R_bar)); - #=Hecke.=#mod_sym!(R_bar, m) + mod_sym!(R_bar, m) mul!(M, B0_modm, R_bar) - #=Hecke.=#mod_sym!(M, m) + 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) From f771c9d4040f23b535e7d7110e747c3fe2f086f1 Mon Sep 17 00:00:00 2001 From: John Abbott Date: Wed, 15 Jan 2025 12:37:25 +0100 Subject: [PATCH 5/8] Minor revision to reduce code coverage complaints --- src/matrix.jl | 24 ++++++------------------ 1 file changed, 6 insertions(+), 18 deletions(-) diff --git a/src/matrix.jl b/src/matrix.jl index 7ce03f0cc5..8dbad38a1f 100644 --- a/src/matrix.jl +++ b/src/matrix.jl @@ -523,25 +523,13 @@ function is_unimodular_given_det_mod_m(A::ZZMatrix, det_mod_m::Int, M::ZZRingEle # (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. - if !is_square(A) - throw(ArgumentError("Matrix must be square")) - end - if abs(det_mod_m) != 1 - throw(ArgumentError("det_mod_m must be +1 or -1")) - end - if M <= 2^30 - throw(ArgumentError("modulus must be greater than 2^30")) - end - if !(algorithm in [:auto, :CRT, :pauderis_storjohann]) - throw(ArgumentError("algorithm must be one of [:CRT, :pauderis_storjohann, :auto]")) - end + 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??? - if nrows(A) == 0 - return true - end - if nrows(A) == 1 - return (abs(A[1,1]) == 1) - end + 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") From 9a8c48d9e468bec6c78db3ac46db6869744edc02 Mon Sep 17 00:00:00 2001 From: John Abbott Date: Wed, 15 Jan 2025 12:38:42 +0100 Subject: [PATCH 6/8] Increased test coverage --- test/matrix-test.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/matrix-test.jl b/test/matrix-test.jl index 17d4710803..1100d3293c 100644 --- a/test/matrix-test.jl +++ b/test/matrix-test.jl @@ -111,6 +111,10 @@ end @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)) + @test_throws ArgumentError is_unimodular(matrix(ZZ,0,1,[])) @test_throws ArgumentError is_unimodular(matrix(ZZ,0,0,[]); algorithm=:WRONG) From d62285ca8bcdab24d1ea5a15cfc6adde52d5756c Mon Sep 17 00:00:00 2001 From: John Abbott Date: Wed, 15 Jan 2025 13:32:05 +0100 Subject: [PATCH 7/8] Increase test coverage --- test/matrix-test.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/matrix-test.jl b/test/matrix-test.jl index 1100d3293c..dfad50da82 100644 --- a/test/matrix-test.jl +++ b/test/matrix-test.jl @@ -115,6 +115,11 @@ end @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) From f6c2d6733f3db76f49ddec7818402dc46f4945e6 Mon Sep 17 00:00:00 2001 From: John Abbott Date: Wed, 15 Jan 2025 13:36:06 +0100 Subject: [PATCH 8/8] Increase test coverage --- test/matrix-test.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/matrix-test.jl b/test/matrix-test.jl index dfad50da82..3fe9f771ac 100644 --- a/test/matrix-test.jl +++ b/test/matrix-test.jl @@ -136,7 +136,8 @@ end 0, -2, 71, 77, 999925, 1]) - for k in 0:9 + # 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)