Numerically stable LLL algorithm using Givens rotations

LLL.jl

#
# Loosely based on the paper at https://core.ac.uk/download/pdf/82729885.pdf but with some significant deviation.
#
#  This implementation should be robust even for bad inputs with large condition numbers. The reduced basis of a lattice L with basis
#  vectors as columns of a matrix M is given by M*lll_unimod(M) .
#  

#
#  Acceptably optimized. Spends just under 50% of its time in the call to LinearAlgebra.qr for 2000 x 2000 matrices.
#  Benchmarks fairly well compared to existing libraries in multiple languages.
#

using LinearAlgebra

function _lll_decrease!(uppertriangular,i::Integer,j::Integer,unimodular)
    γ = Int(round(uppertriangular[i,j]/uppertriangular[i,i]))
    @inbounds @simd for k in 1:size(unimodular,1) # Hot loop.
        unimodular[k,j] -= γ*unimodular[k,i]
        uppertriangular[k,j] -= γ*uppertriangular[k,i]
    end
end

function _swapcols!(X, i::Integer, j::Integer)
    @inbounds @simd for k = 1:size(X,1)
        X[k,i], X[k,j] = X[k,j], X[k,i]
    end
end

function _coeffs(x,y)
    h = hypot(x,y)
    (x/h, y/h)
end

function _lll_swap!(uppertriangular,i::Integer,j::Integer,unimodular) # where j < i
    _swapcols!(unimodular,j,i)
    _swapcols!(uppertriangular,j,i) # Breaks upper triangularity.
    c,s = _coeffs(uppertriangular[j,j],uppertriangular[i,j])
    @inbounds @simd for k = 1:size(uppertriangular,2)  # Fix upper triangularity with a Givens rotation on rows i,j.
        uppertriangular[j,k] , uppertriangular[i,k] = c*uppertriangular[j,k] + s*uppertriangular[i,k] , c*uppertriangular[i,k] - s*uppertriangular[j,k] 
    end
end

#
# Outputs a unimodular integer matrix to change from the bad to the good basis,
# but can also reduce m directly (by mutating it) if you call lll_unimod(m;output=m).
# lll_unimod(m,output=deepcopy(m)) is generally faster than m*lll_unimod(m)
#
function lll_unimod(m::AbstractMatrix, δ=0.75; output= Matrix{Int}(I, size(m)))
    r = LinearAlgebra.qr(m).R # I.e. q' * m = r . The key Invariant of the algorithm as we mutate output and r 
    k = 2 # is that Q*m*output = r where Q (not computed) soaks up all left unitaries. 
    @inbounds while k <= size(r,2)
        if abs(r[k-1,k-1]) < 2*abs(r[k-1,k])    _lll_decrease!(r,k-1,k,output) end
        if r[k,k]^2  + r[k-1,k]^2 < δ*r[k-1,k-1]^2
            _lll_swap!(r,k,k-1,output) 
            k = max(2,k-1)
        else
            for i in (k-2):-1:1
                if abs(r[i,i]) < 2*abs(r[i,k])  _lll_decrease!(r,i,k,output) end
            end
            k += 1
        end
    end
    return output
end