mu :=function(v,u)
  return InnerProduct(v,u)/InnerProduct(u,u);
end function;

Proj:=function(v,u)
  return mu(v,u)*u;
end function;

GramSchmidt:=function (M)
  local N;
  M:=Transpose(M);
  N:=M;
  for i:=1 to Ncols(M) do
    N[i]:=M[i];
    for j:=1 to i-1 do
      N[i]:=N[i]-Proj(M[i],N[j]);
    end for;
  end for;
  return Transpose(N);
end function;

N2:=function(v)
  return InnerProduct(v,v);
end function;

LLLreduced:=function(M)
  local N,i,j;
  N:=GramSchmidt(M);
  M:=Transpose(M);
  N:=Transpose(N);
  for i:=1 to Nrows(M) do
    for j:=1 to i-1 do
      if Abs(mu(M[i],N[j])) gt 1/2 then
        return false,mu(M[i],N[j]);
      end if;
    end for;
  end for;
  for i:=2 to Nrows(M) do
    if N2(N[i]) lt 3/4-mu(M[i],N[i-1])^2*N2(N[i-1]) then
      return false,_;
    end if;
  end for;
  return true,_;
end function;

M:=Matrix(Rationals(),2,2,[[1,2],[3,4]]);
LLL(M);
GramSchmidt(M);
LLLreduced(M);
M:=Matrix(Rationals(),3,3,[[1,2,3],[1,4,9],[1,8,27]]);
GramSchmidt(M);
LLLreduced(M);
