-
Notifications
You must be signed in to change notification settings - Fork 62
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
Changes from all commits
aebe446
1992616
564c282
2062d74
f771c9d
9a8c48d
d62285c
f6c2d67
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
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? | ||||||||||
####################################################### | ||||||||||
|
||||||||||
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) | ||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. The argument order clashes with the general |
||||||||||
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
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. But we already define such a method in
(and then generic code in AA also makes |
||||||||||
|
||||||||||
@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
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. While this is fine, I'd recommend to use the
Suggested change
|
||||||||||
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 | ||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe 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
|
||||||||||
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) | ||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Suggested change
|
||||||||||
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
There was a problem hiding this comment. Choose a reason for hiding this commentThe 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
Suggested change
|
||||||||||
# 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 |
There was a problem hiding this comment.
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