SAGE Worksheet: 2D LLL Toy

#The
LLL_toy2D algorithm is the same as LLL_toy except it has a different
animation which requires that it only take two dimensional vector
spaces as input. The first animation is the same as in the LLL_toy. The
second animation graphs the actual vectors as the reduction process
takes place as opposed to graphing the second LLL condition. The green
vector represents b_0 (or Basis_Matrix[0]) and the blue vector
represents b_1 (or Basis_Matrix[1]).
#

#In this first example, both animations change.  This is not the case in the second and third examples.

a = matrix(2,[-2262,-2,-6,-1]);show(a)

��2262 �6 �2 �1  �

��2262 �6 �2 �1  �

Basis_MatrixA = LLL_toy2D(a)

��6 �60 �1 365  � 
None

��6 �60 �1 365  � 
None

is_LLL_reduced(a)

True

True

#In this example, the meu values always remain zero, but the vectors are still swapped.

b = matrix(2,[-1,3,2,1]);show(b)

��1 2 3 1  �

��1 2 3 1  �

Basis_MatrixB = LLL_toy2D(b)

�2 �1 1 3  � 
None

�2 �1 1 3  � 
None

is_LLL_reduced(b)

True

True

#This matrix is already LLL reduced, so nothing happens.

c = matrix(2,[1,-3,4,1]); show(c)

�1 4 �3 1  �

�1 4 �3 1  �

Basis_MatrixC = LLL_toy2D(c)

�1 4 �3 1  � 
None

�1 4 �3 1  � 
None

is_LLL_reduced(c)

True

True

def RED(k,l,meu,Basis_Matrix,H):
#Input:
#    k -- integer < Basis_Matrix.nrows()= n
#    l -- integer < k
#    mue -- n by n matrix
#    Basis_Matrix -- n by n matrix
#    H -- n by n matrix
#Output:
#    mue -- n by n matrix
#    Basis_Matrix -- n by n matrix
#    H -- n by n matrix
    if abs(meu[k,l]) <= 0.5:
        return(meu,Basis_Matrix,H)
    if abs(meu[k,l]) > 0.5:
        
        q = integer_floor(meu[k][l]+0.5)
        Basis_Matrix[k] = Basis_Matrix[k]-q*Basis_Matrix[l]
        H = H.rows()
        H[k] = H[k]-q*H[l]
        H = matrix(H)

        meu[k,l] = meu[k,l]-q 
        for i in range(l):
            meu[k,i] = meu[k,i]-q*meu[l,i]

    return(meu,Basis_Matrix,H)

def SWAP(k,k_max,meu,B,Basis_Matrix,basis_star,H):
#this is the algorithm that assures the second LLL criteria is reached
#Input:
#    k -- integer < Basis_Matrix.nrows()= n
#    l -- integer < k
#    mue -- n by n matrix
#    B -- list of length n
#    Basis_Matrix -- n by n matrix
#    basis_star -- n by n matrix
#    H -- n by n matrix
#Output:
#    mue -- n by n matrix
#    B -- list of length n
#    basis_star -- n by n matrix
#    Basis_Matrix -- n by n matrix
#    H -- n by n matrix
   
    basis_swap = Basis_Matrix
    Basis_Matrix.swap_rows(k,k-1)

    H1 = H
    H.swap_rows(k,k-1)
    
    meu_swap = copy(meu)
    
    if k > 0:
        for i in range(k-1):
            
            meu[k,i] = meu_swap[k-1,i]
            meu[k-1,i] = meu_swap[k,i]

    meu_entry = meu[k,k-1]
    B_swap = B[k] + meu_entry^2*B[k-1]
    meu[k,k-1] = meu_entry*B[k-1]/B_swap
    

    basis_entry = basis_star[k-1]
    basis_star = basis_star.rows()
    basis_star[k-1] = basis_star[k]+meu_entry*basis_entry
    basis_star[k] = -meu[k,k-1]*basis_star[k] + B[k]/B_swap*basis_entry
    basis_star = matrix(basis_star)
    

    B[k] = B[k-1]*B[k]/B_swap
    B[k-1] = B_swap

    for i in range(k+1,k_max+1):
        t = meu[i,k]
        meu[i,k] = meu[i,k-1] - meu_entry*t
        meu[i,k-1] = t + meu[k,k-1]*meu[i,k]
    
    return(meu,B,basis_star,Basis_Matrix,H)

def Initialize(Basis_Matrix):
#Input:
#    Basis_Matrix -- n by n matrix
#Output:
#    basis_star -- n by n matrix
#    mue -- n by n matrix
#    k -- integer
#    k_max -- integer
#    B -- list of length n
#    H -- n by n identity matrix 
    
    B = [0 for i in range(Basis_Matrix.nrows())] #this is the list of dot products of rows
    B[0] = (Basis_Matrix[0])*(Basis_Matrix[0]) #initializing the first entry

    basis_star = matrix(QQ,Basis_Matrix.nrows(),Basis_Matrix.ncols())
    basis_star[0] = Basis_Matrix[0]
    basis_star = matrix(basis_star)

    meu = matrix(QQ,Basis_Matrix.nrows(),Basis_Matrix.ncols())

    k = 1
    k_max = 0

    H = identity_matrix(Basis_Matrix.nrows())
    
    return(basis_star,meu,k,k_max,B,H)

def LLL_toy2D(Basis_Matrix):
#Input:
#    Basis_Matrix -- n by n matrix
#Output:
#    Basis_Matrix -- n by n matrix
#Warning: this only works with Basis_Matrix of size 2x2.

#This algorithm is the basic LLL algorithm found on page 88 of Cohen's Computational
#Number Theory text, specifically algorithm 2.6.3

#Background
# we will be using the typical Gram-Schmidt notation for the Gram-Schmidt process
# with meu for the greek letter mue and basis_star for b*
# a basis is called LLL reduced if 
#     1. abs(meu[i,j])<=1/2 for 0<=j<i<=n 
# and 2. basis_star[i]*basis_star[i]>=(3/4-meu[i,i-1]^2)*basis_star[i-1]*basis_star[i-1]

#Important points about LLL
#  This algorithm only computes Gram-Schmidt coefficients meu[k,j] and B[k]
#  when needed, thus step 2 is just going through one iteration of the Gram-Schmidt process
#  Which have been computed is recorded in the value k_max, which is the larges k-value attained

#note: the following are some assumption which are made when the algorithm is run
# First: that the input matrix does in fact represent a basis, this is first checked
#  in the very first if statement and makes the B[k]=0 redundant
#  Cohen states that the assumption of linearly independent row vectors is can be avoided
#  and is in fact in his algorithm 2.6.8 (Cohen: Computational Number Theory)
# Second: the matrix entries are rational, this is just for ease of my own checking and
#   for greater precision and could be replaced once the algorithm is consistently working    
#preliminary check

#check if this is in fact a basis because the algorithm does not work if it isn't
    if Basis_Matrix.rank() != Basis_Matrix.nrows():
        raise NotImplementedError,"Not a Basis"
    if Basis_Matrix.nrows()>2:
        raise NotImplementedError,"Only implemented for 2x2"
#step 1
# in this step the matrices meu,B,basis_star (b* in Gram-Schmidt notation), and H are initialized
# meu,basis_star, and B are all zeros and H is the identity matrix
# the first Gram-Schmidt values of basis_star and B are then computed (i.e. for the 0th row)
    basis_star,meu,k,k_max,B,H = Initialize(Basis_Matrix)
#Mp collects the matrix objects used to graphically compare the values of meu 
    Mp = [matrix_plot(meu)]  #for animation of meu values as they converge
    Lp0 = [line([(0,0),Basis_Matrix[0]],rgbcolor=(0.5,1,0.5))]
    Lp1 = [line([(0,0),Basis_Matrix[1]])]
#step 2 
# this is just one step of the Gram-Schmidt process, i.e. here we are computing only what we need
# to check whether the basis is or is not reduced
# we use k_max to keep track of where we are in the Gram-Schmidt process, keeping k_max = to the maximum
# value of k which we have currently obtained
    for k in range(1,Basis_Matrix.nrows()):
        
        if k > k_max:
            k_max = k
            basis_star = swap_rows(k,k,basis_star,Basis_Matrix) #sets the k'th row of basis_star = k'th of Basis_Matrix
            
            for j in range(k):
                meu[k,j] = Basis_Matrix[k]*basis_star[j]/B[j]
                basis_star = basis_star.rows();
                basis_star[k] = basis_star[k]-meu[k,j]*basis_star[j]
                basis_star = matrix(basis_star)
                

            B[k] = basis_star[k]*basis_star[k]
            
            if B[k]==0.0:
                return('error Bk = 0, not a basis')
        Mp.append(matrix_plot(meu))
        Lp0.append(line([(0,0),Basis_Matrix[0]],rgbcolor=(0.5,1,0.5)))
        Lp1.append(line([(0,0),Basis_Matrix[1]]))
    
#step 3 [This is where all the work is done]
# in this step we continually run the RED subroutine
# the RED subroutine checks if meu[k,l] is reduced, i.e. less than or equal to 1/2 and if not
# it reduces meu[k,l] by q = floor(meu+1/2) or the closest integer to meu[k,l] ensuring that 
# meu[k,l] is then less than or equal to 1/2
# now notice two things, first, as we are working inductively, we can assume that for l<j<k, meu[k,j] <= 1/2 
# and basis_star[j] is orthogonal to Basis_Matrix[l].  Second, replacing Basis_Matrix[k] by 
# Basis_Matrix[k]-q*Basis_Matrix[l] is equivalent to replacing meu[k,l] by meu[k,l]-q, 
# and for j>l, meu[k,j] is not modified, finally, for 0<=i<l, the algorithm updates meu[k,i] = meu[k,i]-q*meu[l,i]
# so that the matrix meu is consistent with the matrices Basis_Matrix and basis_star as in the Gram-Schmidt process

# Next we check the other condition of LLL reduced, (the Lovasz condition) B[k]>=(3/4-meu[k,k-1]^2)*B[k-1]
# if this condition is satisfied we increase k by one and continue on, otherwise we need to use the SWAP algorithm
# The SWAP algorithm swaps Basis_Matrix[k] and Basis_Matrix[k-1], updates meu, basis_star, B, and H and then decreases
# k by one (unless k = 1, then k is reset to k) for the following iteration of step 3 where both conditions are again
# checked.
        
        meu,Basis_Matrix,H = RED(k,k-1,meu,Basis_Matrix,H)
        
        while B[k]<(0.75-meu[k,k-1]^2)*B[k-1]:
#        while B[k]<(0.75-meu[k,k-1]^2)*B[k-1]:
            meu,B,basis_star,Basis_Matrix,H = SWAP(k,k_max,meu,B,Basis_Matrix,basis_star,H)
            k = max(1,k-1)
            Mp.append(matrix_plot(meu))
            Lp0.append(line([(0,0),Basis_Matrix[0]],rgbcolor=(0.5,1,0.5)))
            Lp1.append(line([(0,0),Basis_Matrix[1]]))
            meu,Basis_Matrix,H = RED(k,k-1,meu,Basis_Matrix,H)
       
        if B[k] >= (0.75-meu[k,k-1]^2)*B[k-1]:
            for i in range(2,k-1):
                l = k-i
                meu,Basis_Matrix,H = RED(k,k-1,meu,Basis_Matrix,H) 
            Mp.append(matrix_plot(meu))
            Lp0.append(line([(0,0),Basis_Matrix[0]],rgbcolor=(0.5,1,0.5)))
            Lp1.append(line([(0,0),Basis_Matrix[1]]))

    print(show(Basis_Matrix))
    a = animate(Mp)
    a.show() 
#returns a block showing the transformation the meu's go through to become LLL reduced
    c0 = animate(Lp0)
    c1 = animate(Lp1)
    c = c0+c1
    c.show()
#c.show() returns an animation of the transformation of the basis vectors in Basis_Matrix 
#illustrating the steps they go through to become LLL reduced.
#black squares indicate meu is close to zero
    return(Basis_Matrix)

def swap_rows(i,j,M1,M2):
#Input:
#    i -- integer < n = M1.nrows()
#    j -- integer < n = M1.nrows()
#    M1 -- n by n matrix
#    M2 -- n by n matrix
#Output:
#    M1 -- n by n matrix
    M1 = M1.rows()
    M1[i] = M2.rows()[j]
    M1[j] = M2.rows()[i]
    M1 = matrix(M1)

    return(M1)

def is_LLL_reduced(Basis_Matrix):
#Input:
#    Basis_Matrix -- n by n matrix
#Output:
#     boolean -- returns true if Basis_Matrix is LLL reduced
#this nifty algorithm returns true if the input matrix is LLL reduced
    meu,basis_star = Gram_Schmidt(Basis_Matrix)
    for i in range(1,basis_star.nrows()):
        for j in range(i):
#test condition one
            if abs(meu[i,j])>0.5:
                print('i',i,'j',j)
                return(false)
#test condition two
        if abs(basis_star[i]*basis_star[i]) < (3/4 - meu[i,i-i])*basis_star[i-1]*basis_star[i-1]:
            print('i','j',i,j)
            return(false)
    return(true)

def Gram_Schmidt(Basis_Matrix):
#this is the Gram_Schmidt algorithm
#it is called in the is_LLL_reduced()

#step 1
    basis_star,meu,k,k_max,B,H = Initialize(Basis_Matrix)
#step 2
#one note, the only reason k_max is in this code is so that Initialize is still compatible with LLL_toy()
    for k in range(1,Basis_Matrix.nrows()):
        if k > k_max:
            k_max = k
            basis_star = swap_rows(k,k,basis_star,Basis_Matrix) #sets the k'th row of basis_star = k'th of Basis_Matrix
            
            for j in range(k):
                meu[k][j] = Basis_Matrix[k]*basis_star[j]/B[j]
                
                basis_star = basis_star.rows()
                basis_star[k] = basis_star[k]-meu[k][j]*basis_star[j]
                basis_star = matrix(basis_star)
                

            B[k] = basis_star[k]*basis_star[k]
            
            if B[k]==0.0:
                return('error Bk = 0, not a basis')
    return(meu,basis_star)