SAGE Worksheet: General LLL toy with Examples

r"""
Here are some examples of the LLL_toy algorithm implemented on various
matrices given an n by n matrix M of row vectors $b_0,b_1,...,b_(n-1)$,
LLL_toy returns the LLL reduced matrix.
LLL_toy also has two animations associated with it. The first animation
represents the changes in the entries of the matrix mue where the
entries mue[i,j] are as given in Gram-Schmidt, where mue[i,j] =
b_i*b*_j/|b_j| and b*_i = b_i - sum_(j=0)^(i-1)mue[i,j]*b*_j for
(1<= i <=n). Each square of the animation represents the entry of
mue in the corresponding position. Black represents 0 and lighter
colors indicate larger values. As the matrix mue is initialized to the
n x n zero matrix, the animation begins with all the squares black. The
squares appear in lighter shades of gray as the entries' sizes explode,
and finally the matrix returns to its original black state as the first
LLL condition, abs(mue[i,j])<=1/2 for 1<=j<i<=n, is
achieved. The second animation associated with LLL_toy illustrates the
second LLL condition, the Lovasz condition, where we require
|b*_i|>=(3/4-mue[i,i-1]^2)*|b*_(i-1)|. The red lines following the
blue lines represent (3/4-\mue[i,i-1]^2)*|b*_(i-1)| and the blue lines
|b*_i|, where the x-coordinate is i.
Brief warning: You need to run is_LLL_reduced after LLL_toy to check
that the returned matrix, Basis_Matrix is in fact LLL reduced. The
algorithm seems to be much less consistent with larger n, currently
n>5.
"""

a1 = random_matrix(ZZ,5);a1

[ -1  -1  -2  -4  -1]
[ -2  -2   4  -4   2]
[  4   2  -1   1   4]
[  1  -1   1  -1  -1]
[ -1   1  -1   1 -30]

[ -1  -1  -2  -4  -1]
[ -2  -2   4  -4   2]
[  4   2  -1   1   4]
[  1  -1   1  -1  -1]
[ -1   1  -1   1 -30]

a5 = matrix(5,[2,2,1,8,-1,1,9,-2,1,-6,3,-1,2,1,1,1,-1,-2,-3,1,-2,2,5,-1,-1]); show(a5)

0 B B B @  2 1 3 1 �2 2 9 �1 �1 2 1 �2 2 �2 5 8 1 1 �3 �1 �1 �6 1 1 �1  1 C C C A

0 B B B @  2 1 3 1 �2 2 9 �1 �1 2 1 �2 2 �2 5 8 1 1 �3 �1 �1 �6 1 1 �1  1 C C C A

Basis_Matrix = LLL_toy(a5)

0 B B B @  �1 2 1 �1 1 3 2 �1 3 9 0 2 �2 �1 �2 0 1 �3 7 1 �1 0 1 �2 �6  1 C C C A   
None

0 B B B @  �1 2 1 �1 1 3 2 �1 3 9 0 2 �2 �1 �2 0 1 �3 7 1 �1 0 1 �2 �6  1 C C C A   
None

is_LLL_reduced(Basis_Matrix)

True

True

#this is the actual command for the LLL reduction algorithm
#notice that this the reduction is not unique
show(a5.LLL())

0 B B B @  �1 2 1 0 1 3 2 �1 0 1 0 2 �2 1 �5 0 1 �3 �1 1 �1 0 1 �4 0  1 C C C A

0 B B B @  �1 2 1 0 1 3 2 �1 0 1 0 2 �2 1 �5 0 1 �3 �1 1 �1 0 1 �4 0  1 C C C A

a3 = matrix(3,[1,4,-1,4,18,-1,-1,-5,-1]); show(a3)

�1 4 �1 4 18 �5 �1 �1 �1  !

�1 4 �1 4 18 �5 �1 �1 �1  !

Basis_Matrix = LLL_toy(a3)

�0 0 1 �1 0 4 0 1 �1  ! 
None

�0 0 1 �1 0 4 0 1 �1  ! 
None

is_LLL_reduced(a3)

True

True

show(a3.LLL())

�0 0 1 �1 0 0 0 1 0  !

�0 0 1 �1 0 0 0 1 0  !

a6 = random_matrix(ZZ,6); show(a6)

0 B B B B B B @  �1 2 �3 �1 �1 �40 1 �2 �1 30 1 1 1 �1 1 �2 �1 1 4 �1 �1 �1 �1 �2 �1 101 2 2 1 6 1 11 3 �1 1 �1  1 C C C C C C A

0 B B B B B B @  �1 2 �3 �1 �1 �40 1 �2 �1 30 1 1 1 �1 1 �2 �1 1 4 �1 �1 �1 �1 �2 �1 101 2 2 1 6 1 11 3 �1 1 �1  1 C C C C C C A

Basis_Matrix6 = LLL_toy(a6)

0 B B B B B B @  �1 �1 �3 �4 �40 �3 1 1 �1 29 1 3 �1 1 1 �1 1 4 �1 4 �1 �2 �2 19 1 �1 2 4 6 96 1 1 3 2 �1 16  1 C C C C C C A  
None

0 B B B B B B @  �1 �1 �3 �4 �40 �3 1 1 �1 29 1 3 �1 1 1 �1 1 4 �1 4 �1 �2 �2 19 1 �1 2 4 6 96 1 1 3 2 �1 16  1 C C C C C C A  
None

is_LLL_reduced(Basis_Matrix6)

True

True

a10 = random_matrix(ZZ,10); show(a10)

0 B B B B B B B B B B B B B B @  �1 2 2 1 �5 1 1 2 �1 2 �1 2 1 �1 �1 3 �1 1 1 �1 �3 8 �1 �4 1 8 �5 �5 �1 2 �8 1 5 �2 1 1 �2 �26 2 �1 �2 �1 1 42 1 �1 �7 1 1 �5 1 1 8 2 �1 �41 1 �2 �3 �1 1 1 7 �1 1 �21 �1 23 �1 1 3 1 1 �140 1 1 8 9 �3 �6 1 1 1 �2 1 �2 �5 �1 �1 2 1 �1 �1 �1 4 �2 �3 �16 99 �5  1 C C C C C C C C C C C C C C A

0 B B B B B B B B B B B B B B @  �1 2 2 1 �5 1 1 2 �1 2 �1 2 1 �1 �1 3 �1 1 1 �1 �3 8 �1 �4 1 8 �5 �5 �1 2 �8 1 5 �2 1 1 �2 �26 2 �1 �2 �1 1 42 1 �1 �7 1 1 �5 1 1 8 2 �1 �41 1 �2 �3 �1 1 1 7 �1 1 �21 �1 23 �1 1 3 1 1 �140 1 1 8 9 �3 �6 1 1 1 �2 1 �2 �5 �1 �1 2 1 �1 �1 �1 4 �2 �3 �16 99 �5  1 C C C C C C C C C C C C C C A

r"""
Here are two examples of 10 by 10 matrices. The first, a10, is
LLL_reduced.
For the second, matrix a10b, my toy implementation fails to reduce the
matrix. You can see this in both animations.
In the first animation, just before the loop ends, everything seems to
be converging, the gray squares are fairly minimal, however there is
one white square. Also notice that this white square occurs where the
algorithm fails, which is noted in is_LLL_reduced.
Additionally, in the second animation, at the end of the loop when the
Lovasz condition is being met for the larger i values, the term
(3/4-mue[i,i-1]^2)*|b*_(i-1)| explodes for a small i value.
"""

Basis_Matrix = LLL_toy(a10)

0 B B B B B B B B B B B B B B @  �5 �1 2 2 2 1 9 �7 3 �4 �1 �1 2 �1 1 �1 7 �6 0 1 1 �3 8 2 �1 �5 4 �9 �5 4 1 �8 1 �1 5 �2 21 �47 3 �1 1 �2 �1 �5 1 �7 3 �2 43 �42 �1 1 1 �1 8 1 �9 7 10 �13 1 1 1 1 7 �1 7 16 6 �7 1 3 1 �6 1 8 5 4 �139 136 1 1 1 2 1 �5 2 �3 �1 0 4 1 �1 �5 �1 �3 �6 �10 �2 101  1 C C C C C C C C C C C C C C A  
None

0 B B B B B B B B B B B B B B @  �5 �1 2 2 2 1 9 �7 3 �4 �1 �1 2 �1 1 �1 7 �6 0 1 1 �3 8 2 �1 �5 4 �9 �5 4 1 �8 1 �1 5 �2 21 �47 3 �1 1 �2 �1 �5 1 �7 3 �2 43 �42 �1 1 1 �1 8 1 �9 7 10 �13 1 1 1 1 7 �1 7 16 6 �7 1 3 1 �6 1 8 5 4 �139 136 1 1 1 2 1 �5 2 �3 �1 0 4 1 �1 �5 �1 �3 �6 �10 �2 101  1 C C C C C C C C C C C C C C A  
None

is_LLL_reduced(Basis_Matrix)

True

True

a10b = random_matrix(ZZ,10); show(a10b)

0 B B B B B B B B B B B B B B @  1 1 �1 1 �4 4 �1 �1 �7 �3 �3 2 86 �1 �1 2 1 1 1 1 �1 6 �54 1 �22 1 �5 �1 1 �5 �1 1 1 �1 1 �1 �1 4 �1 1 1 �7 5 �2 1 16 �1 1 �1 �1 1 �13 1 1 1 �3 3 �1 �2 9 62 �2 1 �1 �6 �1 1 24 15 �5 �3 3 5 1 �1 �36 �3 �1 �1 1 9 1 �2 �1 9 1 �5 �2 1 1 4 �5 3 1 �1 �2 �1 �1 �15 1  1 C C C C C C C C C C C C C C A

0 B B B B B B B B B B B B B B @  1 1 �1 1 �4 4 �1 �1 �7 �3 �3 2 86 �1 �1 2 1 1 1 1 �1 6 �54 1 �22 1 �5 �1 1 �5 �1 1 1 �1 1 �1 �1 4 �1 1 1 �7 5 �2 1 16 �1 1 �1 �1 1 �13 1 1 1 �3 3 �1 �2 9 62 �2 1 �1 �6 �1 1 24 15 �5 �3 3 5 1 �1 �36 �3 �1 �1 1 9 1 �2 �1 9 1 �5 �2 1 1 4 �5 3 1 �1 �2 �1 �1 �15 1  1 C C C C C C C C C C C C C C A

Basis_Matrix10b = LLL_toy(a10b)

0 B B B B B B B B B B B B B B @  1 �1 �2 1 �4 �5 4 2 �9 �1 �1 1 0 2 �1 0 2 �1 2 86 1 �5 0 6 �22 �23 1 5 �4 �54 �1 �1 2 1 1 5 �1 0 �1 1 �2 �1 0 �7 1 2 16 �6 5 5 1 3 6 �13 1 0 �3 �12 10 1 �1 1 �6 �2 �6 18 �1 60 �45 1 1 �3 4 3 �1 �2 �36 0 �1 5 �1 �5 6 1 9 7 1 10 �9 �2 1 �1 2 �5 �1 �2 �2 �1 �14 3  1 C C C C C C C C C C C C C C A  
None

0 B B B B B B B B B B B B B B @  1 �1 �2 1 �4 �5 4 2 �9 �1 �1 1 0 2 �1 0 2 �1 2 86 1 �5 0 6 �22 �23 1 5 �4 �54 �1 �1 2 1 1 5 �1 0 �1 1 �2 �1 0 �7 1 2 16 �6 5 5 1 3 6 �13 1 0 �3 �12 10 1 �1 1 �6 �2 �6 18 �1 60 �45 1 1 �3 4 3 �1 �2 �36 0 �1 5 �1 �5 6 1 9 7 1 10 �9 �2 1 �1 2 �5 �1 �2 �2 �1 �14 3  1 C C C C C C C C C C C C C C A  
None

is_LLL_reduced(a10b)

('i', 'j', 8, 7)
False

('i', 'j', 8, 7)
False

def RED(k,l,mue,Basis_Matrix,H):
#this algorithm assures the first LLL condition is achieved

#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(mue[k,l]) <= 0.5:
        #this checks the first condition
        #if satisfied, exit subroutine
        return(mue,Basis_Matrix,H)
    if abs(mue[k,l]) > 0.5:
        #if the first condition is not satisfied, reduce mue[k,l] by q to satisfy it
        q = integer_floor(mue[k][l]+0.5)
        
        #now correct Basis_Matrix, H, and mue to make sure that the Gram-Schmidt
        #definitions are kept
        Basis_Matrix[k] = Basis_Matrix[k]-q*Basis_Matrix[l]
        H[k] = H[k]-q*H[l]
    
        mue[k,l] = mue[k,l]-q 
        for i in range(l):
            mue[k,i] = mue[k,i]-q*mue[l,i]

    return(mue,Basis_Matrix,H)

def SWAP(k,k_max,mue,B,Basis_Matrix,basis_star,H):
#this is the algorithm that assures the second LLL criteria is reached
#this subroutine is only entered if the Lovasz condition for k has already failed
#it swaps rows k and k-1 of Basis_Matrix and then updates everything in accordance
#to the Gram-Schmidt definitions before we return to step three to try again
#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_Matrix.swap_rows(k,k-1)
    H.swap_rows(k,k-1)
    mue_swap = copy(mue)
    
    for i in range(k-1):
        mue[k,i] = mue_swap[k-1,i]
        mue[k-1,i] = mue_swap[k,i]
   
    mue_entry = mue[k,k-1]
    B_swap = B[k] + mue_entry^2*B[k-1]
    mue[k,k-1] = mue_entry*B[k-1]/B_swap
    
    basis_entry = basis_star[k-1]
    basis_star = basis_star.rows()
    basis_star[k-1] = basis_star[k]+mue_entry*basis_entry
    basis_star[k] = -mue[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 = mue[i,k]
        mue[i,k] = mue[i,k-1] - mue_entry*t
        mue[i,k-1] = t + mue[k,k-1]*mue[i,k]
        
    return(mue,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    
    #this algorithm initializes all the matrices required for both LLL_toy and Gram_Schmidt
    B = [0 for i in range(Basis_Matrix.nrows())] #this is the list of dot products of rows
    B[0] = (Basis_Matrix.rows()[0])*(Basis_Matrix.rows()[0]) #initializing the first entry

    basis_star = matrix(QQ,Basis_Matrix.nrows(),Basis_Matrix.ncols()) #b* as in Gram-Schmidt
    basis_star[0] = Basis_Matrix[0] #initializing the first row

    mue = matrix(QQ,Basis_Matrix.nrows(),Basis_Matrix.ncols()) #creating zeros matrix for later

    k = 1 #initializing the indexing at 1
    k_max = 0 #initializing the counter for the number of iterations of the Gram-Schmidt portion
              #of the algorithm

    H = identity_matrix(Basis_Matrix.nrows()) #initializing H, the transformation matrix
    
    return(basis_star,mue,k,k_max,B,H)

def
LLL_toy(Basis_Matrix):
#Input:
# Basis_Matrix -- n by n matrix
#Output:
# Basis_Matrix -- n by n matrix
#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 mue for the greek letter mue and basis_star for b*
# a basis is called LLL reduced if # 1. abs(mue[i,j])<=1/2 for
0<=j<i<=n # and 2.
basis_star[i]*basis_star[i]>=(3/4-mue[i,i-1]^2)*basis_star[i-1]*basis_star[i-1]
#Important points about LLL
# This algorithm only computes Gram-Schmidt coefficients mue[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 #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"
#step 1
# in this step the matrices mue,B,basis_star (b* in Gram-Schmidt
notation), and H are initialized
# mue,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,mue,k,k_max,B,H =
Initialize(Basis_Matrix) n = Basis_Matrix.nrows() Lovasz = 0
#Mp collects the matrix objects used to graphically compare the values
of mue #Bl is the collection of b[i]*b[i] to compare how this converges
#BMl is the collection of (3/4-mue[i,i-1]^2)*b[i-1]*b[i-1] Bl = [list()
for i in range(n)] BMl =[list() for i in range(n)] Mp =
[matrix_plot(mue)] #for animation of mue values as they converge for i
in range(Basis_Matrix.nrows()): Bl[i] = [line([(i,0),(i,B[i])])] if i
== 0: Lovasz = 0 else: Lovasz = (3/4 - mue[i,i-1]^2)*B[i-1] BMl[i] =
[line([((i+1/2),0),((i+1/2),Lovasz)],rgbcolor=(1,0,0))]
#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,n): 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): mue[k,j] =
Basis_Matrix[k]*basis_star[j]/B[j] basis_star = basis_star.rows();
basis_star[k] = basis_star[k]-mue[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(mue)) for i
in range(n): Bl[i].append(line([(i,0),(i,B[i])])) if i == 0: Lovasz = 0
else: Lovasz = (3/4 - mue[i,i-1]^2)*B[i-1]
BMl[i].append(line([((i+1/2),0),((i+1/2),Lovasz)],rgbcolor=(1,0,0)))
#step 3 [This is where all the work is done]
# in this step we continually run the RED subroutine
# the RED subroutine checks if mue[k,l] is reduced, i.e. less than or
equal to 1/2 and if not
# it reduces mue[k,l] by q = floor(mue+1/2) or the closest integer to
mue[k,l] ensuring that # mue[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, mue[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 mue[k,l]
by mue[k,l]-q, # and for j>l, mue[k,j] is not modified, finally, for
0<=i<l, the algorithm updates mue[k,i] = mue[k,i]-q*mue[l,i]
# so that the matrix mue 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-mue[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 mue, 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. mue,Basis_Matrix,H = RED(k,k-1,mue,Basis_Matrix,H) while
B[k]<(0.75-mue[k,k-1]^2)*B[k-1]: mue,B,basis_star,Basis_Matrix,H =
SWAP(k,k_max,mue,B,Basis_Matrix,basis_star,H) k = max(1,k-1)
Mp.append(matrix_plot(mue)) for i in range(n):
Bl[i].append(line([(i,0),(i,B[i])])) if i == 0: Lovasz = 0 else: Lovasz
= (3/4 - mue[i,i-1]^2)*B[i-1]
BMl[i].append(line([((i+1/2),0),((i+1/2),Lovasz)],rgbcolor=(1,0,0)))
mue,Basis_Matrix,H = RED(k,k-1,mue,Basis_Matrix,H) if B[k] >=
(0.75-mue[k,k-1]^2)*B[k-1]: for i in range(2,k-1): l = k-i
mue,Basis_Matrix,H = RED(k,k-1,mue,Basis_Matrix,H)
Mp.append(matrix_plot(mue)) for i in range(n):
Bl[i].append(line([(i,0),(i,B[i])])) if i == 0: Lovasz = 0 else: Lovasz
= (3/4 - mue[i,i-1]^2)*B[i-1]
BMl[i].append(line([((i+1/2),0),((i+1/2),Lovasz)],rgbcolor=(1,0,0)))
print(show(Basis_Matrix)) #prints the Basis_Matrix a = animate(Mp)
a.show() c = animate(Bl[0]) d = animate(BMl[0]) e = c+d for i in
range(1,n): c = animate(Bl[i]) d = animate(BMl[i]) e = e+c+d e.show()
return(Basis_Matrix)
r""" A quick note of interest: The meat of the proof that this
algorithm returns an LLL reduced basis is not in showing that the
matrix returned is in fact LLL reduced, but that the algorithm does
terminate. This is proven by showing that the number of times SWAP, the
only place in the algorithm where k can decrease, is executed must be
finite. Taking d_i = prod_(1<=j<=i)|b*_j|, d_i is bounded below
and D = prod_(1<=i<=n-1) d_i is decreasing (as it is continually
multiplied by a factor of at most 3/4). So D is bounded below and
decreasing, so this can only continue for a finite number of i's.
However, D changes only when B_i changes, and B_i only changes in SWAP,
therefore SWAP is only iterated a finite number of times. Therefore the
algorithm must terminate. """

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[i] = M2[j]
    M1[j] = M2[i]

    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
    mue,basis_star = Gram_Schmidt(Basis_Matrix)
    for i in range(1,basis_star.nrows()):
        for j in range(i):
#test condition one
            if abs(mue[i,j])>0.5:
                print('i',i,'j',j, mue[i,j])
                return(false)
#test condition two
        if abs(basis_star[i]*basis_star[i]) < (3/4 - mue[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,mue,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[k] = Basis_Matrix[k]
            
            for j in range(k):
                mue[k][j] = Basis_Matrix[k]*basis_star[j]/B[j]
                
                basis_star = basis_star.rows()
                basis_star[k] = basis_star[k]-mue[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(mue,basis_star)