Skip to content

Commit f143dcb

Browse files
fiekerfieker
authored
Merge pull request #1991 from JohnAAbbott/JAA/is_unimodular
Added new func is_unimodular for square integer matrices
2 parents 6e20099 + f6c2d67 commit f143dcb

File tree

4 files changed

+354
-0
lines changed

4 files changed

+354
-0
lines changed

‎src/Exports.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -373,6 +373,7 @@ export is_uinf
373373
export is_undefined
374374
export is_unit
375375
export is_unknown
376+
export is_unimodular
376377
export is_upper_triangular
377378
export is_zero_entry
378379
export is_zero_row

‎src/Nemo.jl

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -193,6 +193,7 @@ import AbstractAlgebra: set_attribute!
193193
import AbstractAlgebra: Solve
194194
import AbstractAlgebra: terse
195195
import AbstractAlgebra: truncate!
196+
import AbstractAlgebra: add_verbosity_scope
196197

197198
AbstractAlgebra.@include_deprecated_bindings()
198199

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

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

363+
add_verbosity_scope(:UnimodVerif)
364+
362365
if AbstractAlgebra.should_show_banner() && get(ENV, "NEMO_PRINT_BANNER", "true") != "false"
363366
show_banner()
364367
end

‎src/matrix.jl

Lines changed: 293 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -394,3 +394,296 @@ function reduce_mod(A::MatElem{T}, B::MatElem{T}) where {T<:FieldElem}
394394
reduce_mod!(C, B)
395395
return C
396396
end
397+
398+
399+
##################################################################
400+
## Functions related to unimodularity verification
401+
##################################################################
402+
403+
# ***SOURCE ALGORITHM***
404+
# Refined implementation of method from
405+
# Authors: C. Pauderis + A. Storjohann
406+
# Title: Detrministic Unimodularity Certification
407+
# Where: Proc ISSAC 2012, pp. 281--287
408+
# See also: thesis of Colton Pauderis, 2013, Univ of Waterloo
409+
# Deterministic Unimodularity Certification and Applications for Integer Matrices
410+
411+
412+
# add_verbosity_scope(:UnimodVerif) -- see func __init__ in src/Nemo.jl
413+
414+
#######################################################
415+
# Low-level auxiliary functions -- should be elsewhere?
416+
#######################################################
417+
418+
function _det(a::fpMatrix)
419+
a.r < 10 && return lift(det(a))
420+
#_det avoids a copy: det is computed destructively
421+
r = ccall((:_nmod_mat_det, Nemo.libflint), UInt, (Ref{fpMatrix}, ), a)
422+
return r
423+
end
424+
425+
function map_entries!(a::fpMatrix, k::Nemo.fpField, A::ZZMatrix)
426+
ccall((:nmod_mat_set_mod, Nemo.libflint), Cvoid, (Ref{fpMatrix}, UInt), a, k.n)
427+
ccall((:fmpz_mat_get_nmod_mat, Nemo.libflint), Cvoid, (Ref{fpMatrix}, Ref{ZZMatrix}), a, A)
428+
a.base_ring = k # exploiting that the internal repr is the indep of char
429+
return a
430+
end
431+
432+
function Base.copy!(A::ZZMatrix, B::ZZMatrix)
433+
ccall((:fmpz_mat_set, Nemo.libflint), Cvoid, (Ref{ZZMatrix}, Ref{ZZMatrix}), A, B)
434+
end
435+
436+
function Nemo.sub!(A::ZZMatrix, m::Int)
437+
for i=1:nrows(A)
438+
A_p = Nemo.mat_entry_ptr(A, i, i)
439+
sub!(A_p, A_p, m)
440+
end
441+
return A
442+
end
443+
444+
# No allocations; also faster than *=
445+
function Nemo.mul!(A::ZZMatrix, m::ZZRingElem)
446+
for i=1:nrows(A)
447+
A_p = Nemo.mat_entry_ptr(A, i, 1)
448+
for j=1:ncols(A)
449+
mul!(A_p, A_p, m)
450+
A_p += sizeof(Int)
451+
end
452+
end
453+
return A
454+
end
455+
456+
@doc raw"""
457+
is_unimodular(A::ZZMatrix)
458+
is_unimodular(A::ZZMatrix; algorithm=:auto)
459+
460+
Determine whether `A` is unimodular, i.e. whether `abs(det(A)) == 1`.
461+
The method used is either that of Pauderis--Storjohann or using CRT;
462+
the choice is made based on cost estimates for the two approaches.
463+
464+
The kwarg `algorithm` may be specified to be `:CRT` or `:pauderis_storjohann`
465+
to ensure that the corresponding underlying algorithm is used (after a quick
466+
initial unimodularity check). `:CRT` is better if the matrix is unimodular
467+
and has an inverse containing large entries (or if it is not unimodular);
468+
`:pauderis_storjohann` is asymptotically faster when the matrix is unimodular
469+
and its inverse does not contain large entries, but it requires more space
470+
for intermediate expressions.
471+
"""
472+
function is_unimodular(A::ZZMatrix; algorithm=:auto)
473+
# Call this function when no extra info about the matrix is available.
474+
# It does a preliminary check that det(A) is +/-1 modulo roughly 2^100.
475+
# If so, then delegate the complete check to is_unimodular_given_det_mod_m
476+
if !is_square(A)
477+
throw(ArgumentError("Matrix must be square"))
478+
end
479+
if !(algorithm in [:auto, :CRT, :pauderis_storjohann])
480+
throw(ArgumentError("algorithm must be one of [:CRT, :pauderis_storjohann, :auto]"))
481+
end
482+
# Deal with two trivial cases
483+
if nrows(A) == 0
484+
return true
485+
end
486+
if nrows(A) == 1
487+
return (abs(A[1,1]) == 1);
488+
end
489+
# Compute det mod M -- may later include some more primes
490+
target_bitsize_modulus = 100 # magic number -- should not be too small!
491+
@vprintln(:UnimodVerif,1,"Quick first check: compute det by crt up to modulus with about $(target_bitsize_modulus) bits")
492+
p::Int = 2^20 # start point of primes for initial CRT trial -- det with modulus >= 2^22 is much slower
493+
M = ZZ(1)
494+
det_mod_m = 0 # 0 means uninitialized, o/w value is 1 or -1
495+
A_mod_p = identity_matrix(Nemo.Native.GF(2), nrows(A))
496+
while nbits(M) <= target_bitsize_modulus
497+
@vprint(:UnimodVerif,2,".")
498+
p = next_prime(p + rand(1:2^10)) # increment by a random step to a prime
499+
ZZmodP = Nemo.Native.GF(p; cached = false, check = false)
500+
map_entries!(A_mod_p, ZZmodP, A)
501+
@vtime :UnimodVerif 2 det_mod_p = Int(_det(A_mod_p))
502+
if det_mod_p > p/2
503+
det_mod_p -= p
504+
end
505+
if abs(det_mod_p) != 1
506+
return false # obviously not unimodular
507+
end
508+
if det_mod_m != 0 && det_mod_m != det_mod_p
509+
return false # obviously not unimodular
510+
end
511+
det_mod_m = det_mod_p
512+
mul!(M, M, p)
513+
end
514+
return is_unimodular_given_det_mod_m(A, det_mod_m, M; algorithm=algorithm)
515+
end #function
516+
517+
518+
# THIS FUNCTION NOT OFFICIALLY DOCUMENTED -- not intended for public use.
519+
function is_unimodular_given_det_mod_m(A::ZZMatrix, det_mod_m::Int, M::ZZRingElem; algorithm=:auto)
520+
# This UI is convenient for our sophisticated determinant algorithm.
521+
# We estimate cost of computing det by CRT and cost of doing Pauderis-Storjohann
522+
# unimodular verification then call whichever is likely to be faster
523+
# (but we avoid P-S if the space estimate is too large).
524+
# M is a modulus satisfying M > 2^30;
525+
# det_mod_m is det(A) mod M [!not checked!]: it is either +1 or -1.
526+
is_square(A) || throw(ArgumentError("Matrix must be square"))
527+
abs(det_mod_m) == 1 || throw(ArgumentError("det_mod_m must be +1 or -1"))
528+
M > 2^30 || throw(ArgumentError("modulus must be greater than 2^30"))
529+
(algorithm in [:auto, :CRT, :pauderis_storjohann]) || throw(ArgumentError("algorithm must be one of [:CRT, :pauderis_storjohann, :auto]"))
530+
# Deal with two trivial cases -- does this make sense here???
531+
nrows(A) == 0 && return true
532+
nrows(A) == 1 && return (abs(A[1,1]) == 1)
533+
@vprintln(:UnimodVerif,1,"is_unimodular_given_det_mod_m starting")
534+
if algorithm == :pauderis_storjohann
535+
@vprintln(:UnimodVerif,1,"User specified Pauderis_Storjohann --> delegating")
536+
return is_unimodular_Pauderis_Storjohann_Hensel(A)
537+
end
538+
n = ncols(A)
539+
Hrow = hadamard_bound2(A)
540+
Hcol = hadamard_bound2(transpose(A))
541+
DetBoundSq = min(Hrow, Hcol)
542+
Hbits = 1+div(nbits(DetBoundSq),2)
543+
@vprintln(:UnimodVerif,1,"Hadamard bound in bits $(Hbits)")
544+
if algorithm == :auto
545+
# Estimate whether better to "climb out" with CRT or use Pauderis-Storjohann
546+
CRT_speed_factor = 20.0 # scale factor is JUST A GUESS !!! (how much faster CRT is compared to P-S)
547+
total_entry_size = sum(nbits,A)
548+
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
549+
EntrySize = maximum(abs, A)
550+
entry_size_bits = maximum(nbits, A)
551+
@vprintln(:UnimodVerif,1,"entry_size_bits = $(entry_size_bits) log2(EntrySize) = $(ceil(log2(EntrySize)))")
552+
e = max(16, Int(2+ceil(2*log2(n)+log2(EntrySize)))) # see Condition in Fig 1, p.284 of Pauderis-Storjohann
553+
Estimate_PS = ceil(e*log(e)) # assumes soft-linear ZZ arith
554+
@vprintln(:UnimodVerif,1,"Est_CRT = $(Estimate_CRT) and Est_PS = $(Estimate_PS)")
555+
Space_estimate_PS = ZZ(n)^2 * e # nbits of values, ignoring mem mgt overhead
556+
@vprintln(:UnimodVerif,1,"Space est PS: $(Space_estimate_PS) [might be bytes]")
557+
if Estimate_PS < Estimate_CRT #=&& Space_estimate_PS < 2^32=#
558+
# Expected RAM requirement is not excessive, and
559+
# Pauderis-Storjohann is likely faster, so delegate to UniCertZ_hensel
560+
return is_unimodular_Pauderis_Storjohann_Hensel(A)
561+
end
562+
end
563+
if algorithm == :CRT
564+
@vprintln(:UnimodVerif,1,"User specified CRT --> delegating")
565+
end
566+
@vprintln(:UnimodVerif,1,"UnimodularVerification: CRT loop")
567+
# Climb out with CRT:
568+
# in CRT loop start with 22 bit primes; if we reach threshold (empirical), jump to much larger primes.
569+
p = 2^21; stride = 2^10; threshold = 5000000;
570+
A_mod_p = zero_matrix(Nemo.Native.GF(2), nrows(A), ncols(A))
571+
while nbits(M) < Hbits
572+
@vprint(:UnimodVerif,2,".")
573+
# Next lines increment p (by random step up to stride) to a prime not dividing M
574+
if p > threshold && p < threshold+stride
575+
@vprint(:UnimodVerif,2,"!!")
576+
p = 2^56; stride = 2^25; # p = 2^56 is a good choice on my standard 64-bit machine
577+
continue
578+
end
579+
# advance to another prime which does not divide M:
580+
while true
581+
p = next_prime(p+rand(1:stride))
582+
if !is_divisible_by(M,p)
583+
break
584+
end
585+
end
586+
ZZmodP = Nemo.Native.GF(p; cached = false, check = false)
587+
map_entries!(A_mod_p, ZZmodP, A)
588+
det_mod_p = _det(A_mod_p)
589+
if det_mod_p > p/2
590+
det_mod_p -= p
591+
end
592+
if abs(det_mod_p) != 1 || det_mod_m != det_mod_p
593+
return false # obviously not unimodular
594+
end
595+
M *= p
596+
end
597+
return true
598+
end
599+
600+
# THIS FUNCTION NOT OFFICIALLY DOCUMENTED -- not intended for public use.
601+
# Pauderis-Storjohann unimodular verification/certification.
602+
# We use Hensel lifting to obtain the modular inverse B0
603+
# Seems to be faster than the CRT prototype (not incl. here)
604+
# VERBOSITY via :UnimodVerif
605+
function is_unimodular_Pauderis_Storjohann_Hensel(A::ZZMatrix)
606+
n = nrows(A)
607+
# assume ncols == n
608+
EntrySize = maximum(abs, A)
609+
entry_size_bits = maximum(nbits, A)
610+
e = max(16, 2+ceil(Int, 2*log2(n)+entry_size_bits))
611+
ModulusLWB = ZZ(2)^e
612+
@vprintln(:UnimodVerif,1,"ModulusLWB is 2^$(e)")
613+
# Now look for prime: about 25 bits, and such that p^(2^sthg) is
614+
# just slightly bigger than ModulusLWB. Bit-size of prime matters little.
615+
root_order = exp2(ceil(log2(e/25)))
616+
j = ceil(exp2(e/root_order))
617+
p = next_prime(Int(j))
618+
@vprintln(:UnimodVerif,1,"Hensel prime is p = $(p)")
619+
m = ZZ(p)
620+
ZZmodP = Nemo.Native.GF(p)
621+
MatModP = matrix_space(ZZmodP,n,n)
622+
B0 = lift(inv(MatModP(A)))
623+
mod_sym!(B0, m)
624+
delta = identity_matrix(base_ring(A), n)
625+
A_mod_m = identity_matrix(base_ring(A), n)
626+
@vprintln(:UnimodVerif,1,"Starting stage 1 lifting")
627+
while (m < ModulusLWB)
628+
@vprintln(:UnimodVerif,2,"Loop 1: nbits(m) = $(nbits(m))");
629+
copy!(A_mod_m, A)
630+
mm = m^2
631+
mod_sym!(A_mod_m,mm)
632+
@vtime :UnimodVerif 2 mul!(delta, A_mod_m, B0)
633+
sub!(delta, 1)
634+
divexact!(delta, m) # to make matrix product in line below cheaper
635+
@vtime :UnimodVerif 2 mul!(A_mod_m, B0, delta)
636+
mul!(A_mod_m, m)
637+
sub!(B0, B0, A_mod_m)
638+
m = mm
639+
mod_sym!(B0, m)
640+
end
641+
@vprintln(:UnimodVerif,1,"Got stage 1 modular inverse: modulus has $(nbits(m)) bits; value is $(p)^$(Int(round(log2(m)/log2(p))))")
642+
643+
# Hadamard bound brings little benefit; anyway, almost never lift all the way
644+
### H = min(hadamard_bound2(A), hadamard_bound2(transpose(A)))
645+
### println("nbits(H) = $(nbits(H))");
646+
# Compute max num iters in k
647+
k = (n-1)*(entry_size_bits+log2(n)/2) - (2*log2(n) + entry_size_bits -1)
648+
k = 2+nbits(ceil(Int,k/log2(m)))
649+
@vprintln(:UnimodVerif,1,"Stage 2: max num iters is k=$(k)")
650+
651+
ZZmodM,_ = residue_ring(ZZ,m; cached = false) # m is probably NOT prime
652+
## @assert is_one(lift(MatModM(B0*A)))
653+
# Originally: R = (I-A*B0)/m
654+
R = A*B0
655+
sub!(R, 1) #this is -R, but is_zero test is the same, and the loop starts by squaring
656+
divexact!(R, m)
657+
if is_zero(R)
658+
return true
659+
end
660+
661+
#ZZMatrix and ZZModMatrix are the same type. The modulus is
662+
#always passed in as an extra argument when needed...
663+
#mul! for ZZModMatrix is mul, followed by reduce.
664+
665+
B0_modm = deepcopy(B0)
666+
mod_sym!(B0_modm, m)
667+
668+
R_bar = zero_matrix(ZZ, n, n)
669+
M = zero_matrix(ZZ, n, n)
670+
671+
@vprintln(:UnimodVerif,1,"Starting stage 2 lifting")
672+
for i in 0:k-1
673+
@vprintln(:UnimodVerif,1,"Stage 2 loop: i=$(i)")
674+
675+
mul!(R_bar, R, R)
676+
#Next 3 lines do: M = lift(B0_modm*MatModM(R_bar));
677+
mod_sym!(R_bar, m)
678+
mul!(M, B0_modm, R_bar)
679+
mod_sym!(M, m)
680+
# Next 3 lines do: R = (R_bar - A*M)/m; but with less memory allocation
681+
mul!(R, A, M)
682+
sub!(R, R_bar, R)
683+
divexact!(R, m)
684+
if is_zero(R)
685+
return true
686+
end
687+
end
688+
return false
689+
end

‎test/matrix-test.jl

Lines changed: 57 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -101,3 +101,60 @@ end
101101

102102
# eigenvalues_simple not defined for integer matrices, so not tested.
103103
end
104+
105+
106+
@testset "is_unimodular" begin
107+
108+
# Some trivial inputs
109+
@test is_unimodular(matrix(ZZ,0,0,[]))
110+
@test is_unimodular(matrix(ZZ,1,1,[1]))
111+
@test is_unimodular(matrix(ZZ,1,1,[-1]))
112+
@test !is_unimodular(matrix(ZZ,1,1,[0]))
113+
114+
@test is_unimodular(identity_matrix(ZZ, 11))
115+
@test !is_unimodular(2*identity_matrix(ZZ, 11))
116+
@test !is_unimodular(zero_matrix(ZZ, 11,11))
117+
118+
# This test increases coverage slightly
119+
N = 1+prod(filter(is_prime, ZZ(2)^20:(ZZ(2)^20+ZZ(2)^11)))
120+
@test !is_unimodular(N*identity_matrix(ZZ, 11))
121+
122+
123+
@test_throws ArgumentError is_unimodular(matrix(ZZ,0,1,[]))
124+
@test_throws ArgumentError is_unimodular(matrix(ZZ,0,0,[]); algorithm=:WRONG)
125+
126+
U = matrix(ZZ, 3,3, [ 3, 22, 46, 5, 35, 73, 4, 27, 56])
127+
M = matrix(ZZ, 3,3, [-3, 22, 46, 5, 35, 73, 4, 27, 56])
128+
129+
# U6 is "small" and triggers the 2nd loop in P-S method
130+
U6 = matrix(ZZ, 6,6,
131+
[1, -76, 87, -57, -41, 999950,
132+
0, 1, 0, 0, 0, 0,
133+
0, 999942, 1, 0, 0, 0,
134+
0, -48, 999963, 1, 0, 0,
135+
0, 91, -24, 1000012, 1, 0,
136+
0, -2, 71, 77, 999925, 1])
137+
138+
139+
# Need larger values of k to increase coverage
140+
for k in 0:25
141+
@test is_unimodular(U^k)
142+
@test is_unimodular(U^k; algorithm=:CRT)
143+
@test is_unimodular(U^k; algorithm=:pauderis_storjohann)
144+
@test is_unimodular(U^k; algorithm=:auto)
145+
@test is_unimodular(-U^k)
146+
@test is_unimodular(-U^k; algorithm=:CRT)
147+
@test is_unimodular(-U^k; algorithm=:pauderis_storjohann)
148+
@test is_unimodular(-U^k; algorithm=:auto)
149+
end
150+
151+
@test is_unimodular(U6)
152+
@test is_unimodular(U6; algorithm=:CRT)
153+
@test is_unimodular(U6; algorithm=:pauderis_storjohann)
154+
@test is_unimodular(U6; algorithm=:auto)
155+
156+
@test !is_unimodular(M)
157+
@test !is_unimodular(M; algorithm=:CRT)
158+
@test !is_unimodular(M; algorithm=:pauderis_storjohann)
159+
@test !is_unimodular(M; algorithm=:auto)
160+
end

0 commit comments

Comments
 (0)