diff --git a/src/flint/fmpz_mat.jl b/src/flint/fmpz_mat.jl index 3613f25b4..82c0fc61d 100644 --- a/src/flint/fmpz_mat.jl +++ b/src/flint/fmpz_mat.jl @@ -957,10 +957,7 @@ input. """ function lll_with_transform(x::ZZMatrix, ctx::LLLContext = LLLContext(0.99, 0.51)) z = deepcopy(x) - u = similar(x, nrows(x), nrows(x)) - for i in 1:nrows(u) - u[i, i] = 1 - end + u = identity_matrix(ZZ, nrows(x)) @ccall libflint.fmpz_lll(z::Ref{ZZMatrix}, u::Ref{ZZMatrix}, ctx::Ref{LLLContext})::Nothing return z, u end @@ -1003,6 +1000,34 @@ function lll!(x::ZZMatrix, ctx::LLLContext = LLLContext(0.99, 0.51)) return x end +# for lll_gram, we use the following helpers to deal with the degenerate case. + +function _is_definitely_full_rank(x::ZZMatrix) + if nrows(x) == 0 + return x + end + if Int == Int64 + p = 1099511627791 % Int64 + else + p = 1048583 % Int32 + end + F = Nemo.Native.GF(p) + xmodp = F.(x) + # try to recognize the definite case + return rank(xmodp) == nrows(x) +end + +function _complete_to_basis(x::ZZMatrix) + @assert nrows(x) <= ncols(x) + # compute column HNF and take the inverse + h, t = hnf_with_transform(transpose(x)) + for i in 1:nrows(x) + @assert is_one(h[i, i]) + end + transpose!(t) + return inv!(t) +end + @doc raw""" lll_gram_with_transform(x::ZZMatrix, ctx::LLLContext = LLLContext(0.99, 0.51, :gram)) @@ -1015,13 +1040,45 @@ See [`lll_gram`](@ref) for the used default parameters which can be overridden b supplying an optional context object. """ function lll_gram_with_transform(x::ZZMatrix, ctx::LLLContext = LLLContext(0.99, 0.51, :gram)) - z = deepcopy(x) - u = similar(x, nrows(x), nrows(x)) - for i in 1:nrows(u) - u[i, i] = 1 + @req is_square(x) && is_symmetric(x) "The matrix must be a symmetric square matrix" + # try to recognize the definite case + if _is_definitely_full_rank(x) + if x[1, 1] < 0 + y = neg(x) + u = _lll_gram_with_transform!(y, ctx) + return neg!(y), u + else + y = deepcopy(x) + u = _lll_gram_with_transform!(y, ctx) + return y, u + end end - @ccall libflint.fmpz_lll(z::Ref{ZZMatrix}, u::Ref{ZZMatrix}, ctx::Ref{LLLContext})::Nothing - return z, u + # remove the radical using a unimodular transformation + L = kernel(x; side = :left) + n = nrows(L) + r = nrows(x) - n + C = _complete_to_basis(L) + # move kernel to bottom, so that + # C * x * transpose(C) = [ * | 0 ] + # [-------] + # [ 0 | 0 ] + uu = deepcopy(reverse_rows!(C)) + mul!(C, C, x * transpose(C)) + u = identity_matrix(ZZ, nrows(x)) + if C[1, 1] < 0 + _lll_gram_with_transform!(neg!(@view(C[1:r, 1:r])), @view(u[1:r, 1:r]), ctx) + @show u + neg!(C) + else + _lll_gram_with_transform!(@view(C[1:r, 1:r]), @view(u[1:r, 1:r]), ctx) + end + return C, u * uu +end + +function _lll_gram_with_transform!(x::ZZMatrix, u::ZZMatrix, ctx::LLLContext = LLLContext(0.99, 0.51, :gram)) + # modifies x and u + @ccall libflint.fmpz_lll(x::Ref{ZZMatrix}, u::Ref{ZZMatrix}, ctx::Ref{LLLContext})::Nothing + return x, u end @doc raw""" @@ -1029,29 +1086,60 @@ end Return the Gram matrix $L$ of an LLL-reduced basis of the lattice given by the Gram matrix $x$. -The matrix $x$ must be symmetric and non-singular. +The matrix $x$ must be symmetric and semidefinite. If an indefinite matrix is +provided, the output is undefined. By default, the LLL is performed with reduction parameters $\delta = 0.99$ and $\eta = 0.51$. These defaults can be overridden by specifying an optional context object. """ function lll_gram(x::ZZMatrix, ctx::LLLContext = LLLContext(0.99, 0.51, :gram)) + @req is_square(x) && is_symmetric(x) "The matrix must be a symmetric square matrix" z = deepcopy(x) return lll_gram!(z) end - @doc raw""" lll_gram!(x::ZZMatrix, ctx::LLLContext = LLLContext(0.99, 0.51, :gram)) Compute the Gram matrix of an LLL-reduced basis of the lattice given by the Gram matrix $x$ inplace. -The matrix $x$ must be symmetric and non-singular. +The matrix $x$ must be symmetric and semidefinite. By default, the LLL is performed with reduction parameters $\delta = 0.99$ and $\eta = 0.51$. These defaults can be overridden by specifying an optional context object. """ function lll_gram!(x::ZZMatrix, ctx::LLLContext = LLLContext(0.99, 0.51, :gram)) + # try to recognize the definite case + if _is_definitely_full_rank(x) + if x[1, 1] < 0 + return neg!(_lll_gram!(!neg(x), ctx)) + else + return _lll_gram(x, ctx) + end + end + # remove the radical using a unimodular transformation + L = kernel(x; side = :left) + n = nrows(L) + r = nrows(x) - n + C = _complete_to_basis(L) + # move kernel to bottom, so that + # C * x * transpose(C) = [ * | 0 ] + # [-------] + # [ 0 | 0 ] + reverse_rows!(C) + mul!(x, x, transpose(C)) + C = mul!(x, C, x) + @show C + if C[1, 1] < 0 + neg!(_lll_gram!(neg!(@view(x[1:r, 1:r])), ctx)) + else + _lll_gram!(@view(x[1:r, 1:r]), ctx) + end + return x +end + +function _lll_gram!(x::ZZMatrix, ctx::LLLContext = LLLContext(0.99, 0.51, :gram)) @ccall libflint.fmpz_lll(x::Ref{ZZMatrix}, C_NULL::Ptr{Nothing}, ctx::Ref{LLLContext})::Nothing return x end diff --git a/test/flint/fmpz_mat-test.jl b/test/flint/fmpz_mat-test.jl index f361bda45..763173239 100644 --- a/test/flint/fmpz_mat-test.jl +++ b/test/flint/fmpz_mat-test.jl @@ -612,6 +612,31 @@ end B = gram(A) lll_gram!(B) @test B == lll_gram(gram(A)) + + # Semidefinite cases + A = ZZ[17070 17380 -4930 5840 28160; + 17380 17975 -5135 6015 28835; + -4930 -5135 1485 -1715 -8155; + 5840 6015 -1715 2015 9675; + 28160 28835 -8155 9675 46705] + G = lll_gram(A) + GG = ZZ[5 0 0 0 0; + 0 10 0 0 0; + 0 0 10 0 0; + 0 0 0 0 0; + 0 0 0 0 0] + @test G == GG + @test lll_gram(-A) == -GG + + G, T = lll_gram_with_transform(A) + @test T * A * transpose(T) == G + G, T = lll_gram_with_transform(-A) + @test T * (-A) * transpose(T) == G + + @test_throws ArgumentError lll_gram(ZZ[1 0]) + @test_throws ArgumentError lll_gram_with_transform(ZZ[1 0]) + @test_throws ArgumentError lll_gram(ZZ[1 0; 1 1]) + @test_throws ArgumentError lll_gram_with_transform(ZZ[1 0; 1 1]) end @testset "ZZMatrix.nullspace" begin