diff --git a/src/Exports.jl b/src/Exports.jl index ee29b244b..0de7b6dfc 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 119829590..ba4586255 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 448700c03..8dbad38a1 100644 --- a/src/matrix.jl +++ b/src/matrix.jl @@ -394,3 +394,296 @@ 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) + throw(ArgumentError("Matrix must be square")) + end + 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 + 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. + 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 + 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 + 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))) + 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 diff --git a/test/matrix-test.jl b/test/matrix-test.jl index 17ce1f7bd..3fe9f771a 100644 --- a/test/matrix-test.jl +++ b/test/matrix-test.jl @@ -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