# 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 bstar[i] for any i in {1..n}
# Type (II) only changes bstar[i] for i in {k-1,k}.

# (Lovasz condition). 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).

#  Length(first LLL vector) <= "LLL fudge factor" * Length(shortest vector)
#
#  In many applications, the fact that we are a "fudge factor" away from finding
#  a shortest vector is not a problem. For instance, if the "second shortest vector"
#  (defined as a shortest vector in L - SPAN(shortest vector)) is more than
#  "fudge factor" times longer than Length(shortest vector), then LLL is guaranteed to
#  find the shortest vector. We can often arrange for such a situation.
#  Suppose L is constructed from some equations mod p^a.  Then Determinant(L)
#  equals (p^a)^(the number of equations). By increasing "a", we increase
#  product( |bstar[i]|, i=1..n). That increases lengths of the vectors that
#  we don't want to find. Once they become longer than "fudge factor" times
#  the vectors that we do want to find, then LLL will separate the vectors
#  we do want to find from those that we don't want to find.

# The number of swaps during the algorithm is bounded by:
#
#     N_swaps <= (ProgressCounter(output) - ProgressCounter(input)) / MinProgress
#
# To bound ProgressCounter(output), observe 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,K,n,newB_k_1,m,mu,r,delta;
  global MinProgress;
  n := nops(L);
  delta := 2^(-2*MinProgress); # The factor 2 is because we compute length-square
  for i to n do b[i] := L[i] od; # copy the input
  k := 1;

  while k <= n do
     # Property L: 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].
     if k=1 then
	bstar[1] := b[1];
	B[1] := dot(bstar[1], bstar[1]);
	k := 2  # Property L holds.
     fi;

     # 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[4](B[i]),4),i=1..k));

     # If we would swap b[k-1],b[k] then the new bstar[k-1] would become
     #   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 current B[k-1] then swapping
     # makes progress. Only swap if it makes at least MinProgress:
     if newB_k_1 < delta * B[k-1] then
	# Lovasz condition holds, so swap two vectors:
	b[k-1], b[k] := b[k], b[k-1];
	k := k-1     # Property L holds.
     else
	k := k+1     # Property L holds.
     fi
  od;
  [seq(b[i], i=1..n)] # k>n so b[1]..b[n] are 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: the "fudge factor" means that LLL does not always find the shortest
# vector (that's NP-hard in general) but when the first vector has the
# smallest G.S. length then we know we found the shortest non-zero vector in L.

MyLLL(L);