# Notation: b[1]..b[n] are vectors, bstar[1]..bstar[n] are their Gram-Schmidt vectors.
# L = Z b[1] + .. + Z b[n], and Determinant(L) = product( |bstar[i]|, i=1..n).

# The LLL algorithm computes a new basis for L by doing two types of operations:
#  (I)  Size-reduction: b[k] := b[k] - r*b[j] for some j<k and integer r.
#  (II) Swap: b[k-1],b[k] := b[k],b[k-1].
# Type (I) does not change any bstar[i]. Type (II) changes bstar[k-1]
# and bstar[k]. Only do (II) if it increases
#     ProgressCounter := sum(i * log2( |bstar[i]| ), i=1..n)
# by at least:

MinProgress := 0.1;

# Write b[k] = bstar[k] + sum(mu[k,j]*bstar[j], j=1..k-1).
# Definition: b[1]..b[n] are LLL-reduced when abs(mu[i,j]) <= 0.51 for all j<i and there
# is no k for which swapping b[k-1],b[k] increases ProgressCounter more than MinProgress.

# Then:  Length(bstar[k-1]) <= 1.28 * Length(bstar[k])
# Hence:   Length(bstar[1]) <= 1.28^(n-1) * Length(bstar[i]) for every i.
# That implies Length(b[1]) <= 1.28^(n-1) * Length(shortest non-zero vector in L).

# The number of swaps during the algorithm is bounded by:
#
#     N_swaps <= (ProgressCounter(output) - ProgressCounter(input)) / MinProgress
#
# To bound ProgressCounter(output), note that the maximum Gram-Schmidt length never
# increases, and the minimum Gram-Schmidt length never decreases during the algorithm:

dot := proc(v,w) local i; add(v[i]*w[i], i=1..nops(v)) end:

MyLLL := proc(L) local B,b,bstar,i,j,k,n,newB_k_1,m,mu,r,delta; global MinProgress;
  n := nops(L);
  delta := 2^(-2*MinProgress);
  for i to n do b[i] := L[i] od; # copy the input
  k := 2;
  while k <= n do

     bstar[1] := b[1];  B[1] := dot(bstar[1], bstar[1]);
     # At this point, b[1]..b[k-1] are LLL-reduced, bstar[1]..bstar[k-1] are their
     # Gram-Schmidt vectors, which have SquareLengths: B[1]..B[k-1].

     # Now compute  bstar[k] = b[k] - add(mu[k,j]*bstar[j], j=1..k-1) while
     # adjusting b[k] in such a way that abs(mu[k,j]) <= 0.51 for all j<k.
     for j from k-1 to 1 by -1 do
	mu[k,j] := dot(b[k], bstar[j])/B[j];
	if abs(mu[k,j]) > 0.51 then
		r := round(mu[k,j]);
		b[k] := b[k] - r*b[j]; # This step is called "Size-reduction"
		mu[k,j] := mu[k,j] - r
	fi
     od;
     bstar[k] := b[k] - add(mu[k,j]*bstar[j], j=1..k-1);

     B[k] := dot(bstar[k], bstar[k]);
     lprint("Current log2(G.S. lengths):", seq(evalf(log(B[i])/log(4),4),i=1..k));

     # If we would swap b[k-1],b[k] then the new bstar[k-1] would be:
     #   bstar[k] + mu[k,k-1]*bstar[k-1] and the new B[k-1] would be:
     newB_k_1 := B[k] + mu[k,k-1]^2*B[k-1];

     # If the new B[k-1] is smaller than the old B[k-1] then swapping
     # makes progress. Only swap if it makes at least MinProgress:
     if newB_k_1 < delta * B[k-1] then
	b[k-1], b[k] := b[k], b[k-1];
	k := max(k-1, 2)
     else            # In both cases b[1]..b[k-1] will be LLL-reduced.
	k := k+1
     fi
  od;
  [seq(b[i], i=1..n)] # k>n so b[1]..b[n] is LLL-reduced.
end:


# In a typical input, the first couple of vectors have large Gram-Schmidt
# length while later vectors have small Gram-Schmidt length. With each swap,
# the LLL algorithm moves G.S. length towards later vectors.

L := [
[ 1, 0, 0, 0, 0, 0, 0,  10000000000, 0],
[ 0, 1, 0, 0, 0, 0, 0,  -5119504669, -11699359581], 
[ 0, 0, 1, 0, 0, 0, 0, -11066568655, 11978985199],
[ 0, 0, 0, 1, 0, 0, 0,  19680180515, 6814529537],
[ 0, 0, 0, 0, 1, 0, 0,  -2102714461, -26513252424],
[ 0, 0, 0, 0, 0, 1, 0, -29942321727, 16033513214],
[ 0, 0, 0, 0, 0, 0, 1,  34087169230, 26822234281]
];

# In this example, the initial log2(G.S. lengths) are:
#   33.22   33.45   1.118   0.832   0.831   0.7575  0.7345
# Each swap moves G.S. length forward (from b[k-1] to b[k] for some k).
# The log2(G.S. lengths) in the output are:
#   6.635   10.52  10.74   10.82   10.68   10.74   10.80
# Note: LLL can not always find the shortest vector (that's NP-hard)
# but when the first vector has the smallest G.S. length then we found
# the shortest non-zero vector in L.

MyLLL(L);