(*:Name: NumberTheory`LLLalgorithm` *)

(*:Summary:  The extended Lenstra Lenstra Lovasz algorithm finds     *)
(*           a lattice basis consisting of "short" vectors,          *)
(*           together with a transformation that relates the new     *)
(*           basis with the original set of generators. Instead of   *)
(*           giving a generating set directly, one may also give     *)
(*           its Gram matrix. When the generators are linearly       *)
(*           dependent, a reduced basis of the lattice of relations  *)
(*           ("null space lattice") is also obtained.                *)

(*:Author:        Wilberd van der Kallen, July 95, Nov 99            *)
(*                http://www.math.ruu.nl/people/vdkallen/kallen.html *)
(*                bugs to:   vdkallen@math.ruu.nl                    *)

(*:Mathematica Version: 2.2 *)
(*:Package Version: 1.10.1 *)

(*:Warnings:                                                         *)
(* For us a lattice in euclidean n-space need not be of full rank n. *)
(* The result with Method->RationalArithmetic may be quite different *)
(* from the result with Method->IntegerArithmetic.                   *)

(*:Sources:  For the original LLL paper, see                         *)
(* Mathematische Annalen 261, 515-534 (1982).                        *)
(* For the case of dependent vectors, see                            *)
(* M. Pohst, A modification of the LLL-algorithm,                    *)
(* Journal of Symbolic Computation 4 (1987), 123--128.               *)
(* Much of the code below relies on:                                 *)
(* Henri Cohen, A course in computational Algebraic Number Theory,   *)
(* Graduate Texts in Mathematics 138, Springer 1993.                 *)

(* :Discussion:  This package adds ExtendedLatticeReduce,            *)
(* GramLatticeReduce and LatticeReducedQ to the builtin function     *)
(* LatticeReduce. There are matrices for which the builtin function  *)
(* gives up, while the package gets an answer (taking its time).     *)
(* The main point however is to give the user more options and to    *)
(* return more information.                                          *)
(* By default                                                        *)
(*           ExtendedLatticeReduce[matrix]                           *)
(* computes a reduced basis of two lattices. One is the row space    *)
(* lattice spanned by the rows of the matrix, the other is the null  *)
(* space lattice, consisting of vectors v that have integer entries  *)
(* and satisfy v.matrix=0. (Contrary to the function NullSpace, we   *)
(* are concerned with multiplication from the left here.)            *)
(* By means of ReduceRowSpace->False or ReduceNullSpace->False one   *)
(* may turn off the part of the algorithm that makes the respective  *)
(* basis reduced. (They will still be bases.) Thus both              *)
(*      Take[Last[#],Length[Last[#]]-Length[First[#]]]&[             *)
(*      ExtendedLatticeReduce[matrix]]                               *)
(* and                                                               *)
(*      Take[Last[#],Length[Last[#]]-Length[First[#]]]&[             *)
(*      ExtendedLatticeReduce[matrix,ReduceRowSpace->False]]         *)
(* will yield a reduced basis of the null space lattice, but the     *)
(* latter should be faster.                                          *)
(* Using ReduceNullSpace->False one may similarly try to speed       *)
(* things up, at the expense of getting much larger entries in the   *)
(* transformation matrix                                             *)
(*      Last[ExtendedLatticeReduce[matrix,ReduceNullSpace->False]]   *)
(* A safer method is provided with Method->Automatic. It does not    *)
(* reduce the null space lattice either, but it prevents that the    *)
(* entries of the transformation matrix really explode.              *)

(* :Limitation: This implementation is not fast enough for really    *)
(*  big problems. For those one may try other implementations of the *)
(*  LLL algorithm.                                                   *)

BeginPackage["LLLalgorithm`"]

ExtendedLatticeReduce::usage:=
"ExtendedLatticeReduce[{v1, v2, ...}] gives
{reduced basis,transformation}
where the transformation is an integral matrix T of determinant plus or
minus one such that the rows of T.{v1, v2, ...} form a reduced basis,
possibly preceded by zero vectors. These zero vectors occur if the
integral vectors vi are linearly dependent.  In that case the
corresponding top part of T is a reduced basis of the `null space
lattice'.  See ReduceRowSpace for details.\n
See also: LatticeReduce, GramLatticeReduce, LatticeReducedQ";

GramLatticeReduce::usage:=
"GramLatticeReduce[{{vi.vj}}] gives {rank of lattice,transformation}
where the transformation is an integral matrix T of determinant plus or
minus one such that the rows of T.{v1, v2, ...} form a reduced basis,
possibly preceded by zero vectors. These zero vectors occur if the
vectors vi are linearly dependent.\n
The Gram matrix {{vi.vj}} should have rational entries.\n
See also: ExtendedLatticeReduce";

LatticeReducedQ::usage:="LatticeReducedQ[{v1, v2, ...}] is True if
{v1, v2, ...} is LLL reduced, after removing zero vectors. It uses
rational arithmetic.";

RationalArithmetic::usage:="RationalArithmetic is a value which can be
assigned to the option Method in ExtendedLatticeReduce and in
GramLatticeReduce. It disables the other options.
The answer obtained with Method->RationalArithmetic is similar to
the answer with ReduceNullSpace->False."

IntegerArithmetic::usage:="IntegerArithmetic is a value which can be
assigned to the option Method in ExtendedLatticeReduce and in
GramLatticeReduce.\n
See also: RationalArithmetic"

ReduceRowSpace::usage:="ReduceRowSpace is an option in
ExtendedLatticeReduce and in GramLatticeReduce.
Let the null space lattice of the matrix m consist of the
vectors v which satisfy v.m=0 and have integer entries.
To compute a reduced basis of the null space lattice,
without computing a reduced basis
of the row space lattice, one may use
GramLatticeReduce[m.Transpose[m],ReduceRowSpace->False]//
Take[Last[#],Length[Last[#]]-First[#]]&\n
See also: ReduceNullSpace"

ReduceNullSpace::usage:="ReduceNullSpace is an option in
ExtendedLatticeReduce and in GramLatticeReduce. If the rows of the
matrix are linearly dependent, then the entries of the transformation
Last[ExtendedLatticeReduce[matrix,ReduceNullSpace->False]]
may be much larger than those of
Last[ExtendedLatticeReduce[matrix]].\n
To keep them of intermediate size, use Method->Automatic,
which also does not reduce the null space lattice. "

RatioLLLCondition::usage:="RatioLLLCondition is an option in
ExtendedLatticeReduce and in GramLatticeReduce.
It tells which ratio to use in the criterion for an LLL
reduced basis. The default is 3/4, but for some applications one should
take it closer to one.\n RatioLLLCondition->{3/4,19/20,99/100} causes
three passes, with ratios 3/4, 19/20, 99/100 respectively."

(*:Examples:

In[1]:= <<LLLalgorithm.m

In[2]:= generatingset={{1,2,3},{4,5,6},{7,8,9},{10,11,12}};

In[3]:= {reducedbasis,transformation}=ExtendedLatticeReduce[generatingset]

Out[3]= {{{2, 1, 0}, {-1, 1, 3}}, 
 
>    {{1, -1, -1, 1}, {1, -2, 1, 0}, {-1, -1, 1, 0}, {2, 1, -1, 0}}}

In[4]:= transformation.generatingset

Out[4]= {{0, 0, 0}, {0, 0, 0}, {2, 1, 0}, {-1, 1, 3}}

In[5]:= Take[%,-Length[reducedbasis]]

Out[5]= {{2, 1, 0}, {-1, 1, 3}}

In[6]:= %==reducedbasis

Out[6]= True

In[7]:= LatticeReducedQ[reducedbasis]

Out[7]= True

In[8]:= GramLatticeReduce[generatingset.Transpose[generatingset]]==
	{Length[reducedbasis],transformation}

Out[8]= True

In[9]:= reducedbasisnullspace=Take[transformation,Length[transformation]-
         Length[reducedbasis]]

Out[9]= {{1, -1, -1, 1}, {1, -2, 1, 0}}

In[10]:= reducedbasisnullspace.generatingset

Out[10]= {{0, 0, 0}, {0, 0, 0}}

In[11]:= LatticeReducedQ[reducedbasisnullspace]

Out[11]= True

In[12]:= ExtendedLatticeReduce[{{},{}}]

Out[12]= {{}, {{1, 0}, {0, 1}}}

In[13]:= ExtendedLatticeReduce[{}]

ExtendedLatticeReduce::empt:
   Are you sure you want to reduce an empty generating set? We use Null for
     the matrix with no rows or columns.

Out[13]= {{}, Null}

*)

Unprotect[ExtendedLatticeReduce,GramLatticeReduce,LatticeReducedQ];

Options[ExtendedLatticeReduce]={Method->IntegerArithmetic,
ReduceRowSpace->True,ReduceNullSpace->Automatic,RatioLLLCondition->{1/3,3/4}};

Options[GramLatticeReduce]={Method->IntegerArithmetic,
ReduceRowSpace->True,ReduceNullSpace->Automatic,RatioLLLCondition->{1/3,3/4}};

Begin["LLLalgorithm`Private`"]


(*      vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv cut here, if you wish      *)
(* x--------------------------------------------------------------x *)

(* First the rational algorithm, which is simpler.                  *)
(* This algorithm is partly copied from the original LLL paper,     *)
(* Mathematische Annalen 261, 515-534 (1982).                       *)
(* To deal with possibly dependent vectors, some modifications      *)
(* are made, after Pohst, as in algorithm 2.6.8 of:                 *)
(* Henri Cohen, A course in computational Algebraic Number Theory,  *)
(* Graduate Texts in Mathematics 138, Springer 1993.                *)
(* As in Cohen, we also make a version that accepts the Gram matrix *)
(* of a generating set as input, but we do not do any of the        *)
(* further optimizing described by Cohen. In particular, we use     *)
(* rational arithmetic even if the input entries are integers.      *)
(* Therefore we call this the rational algorithm, as opposed to the *)
(* other one below, which uses integer arithmetic.                  *)
(* Information for debugging/modifying of the rational algorithm:
During the algorithm the following situation is maintained/restored.
<
(1) hh.binitial==b (Here binitial is the input matrix in the case of
		    ExtendedLatticeReduce and an identity matrix in the
		    case of GramLatticeReduce. Yes, GramLatticeReduce
		    does some superfluous things.)
(2) With respect to a suitable orthonormal basis (in an n dimensional
    vector space and w.r.t. the inner product given by "gram" in the case
    of GramLatticeReduce) the vector b[[3]], say, looks like
    {e[[1]]mu[[3,1]],e[[2]]mu[[3,2]],e[[3]],0,0,..} and one has
(2a) bb[[3]]==e[[3]]^2 (Here bb[[3]] is a rational number, but e[[3]],
			which is never computed, need not be rational).
(2b) mu[[i,j]]==0 if e[[j]]==0.
		       (If e[[j]]==0, the lattice is contained in the
			subspace spanned by the subset of the
			orthonormal basis obtained by deleting the j-th
			element.)
>
After any "doratmethod[next]", "swapratmethod[k]" or "reduceratmethod[k,l]"
executing
   (tmp=next; initmubb; next=tmp;)
should make no difference, as the conditions (1), (2a), (2b) described above
hold.
At the end of ExtendedLatticeReduce/GramLatticeReduce one should have
an LLL reduced basis, i.e.
(3) Abs[mu[[i,j]]]<=1/2 for n >= i > j >= 1
and
(4) bb[[k]]>=(3/4-mu[[k,k-1]]^2)bb[[k-1]] for 1 <= k <= n.
Conditions (3), (4) can also be checked with LatticeReducedQ.
*)
(* x---------------------------------------------------------------------x *)
(*            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^  cut here, if you wish   *)


reduceratmethod[k_,l_]:=
(
If[Abs[mu[[k,l]]]>1/2,
  q=Round[mu[[k,l]]];
  b[[k]]=b[[k]]-q b[[l]];
  hh[[k]]=hh[[k]]-q hh[[l]];
  mu[[k,l]]=mu[[k,l]]-q;
  Do[mu[[k,i]]=mu[[k,i]]-q mu[[l,i]],{i,l-1}];
  ])

swapratmethod[k_]:=
(
{b[[k]],b[[k-1]]}={b[[k-1]],b[[k]]};
{hh[[k]],hh[[k-1]]}={hh[[k-1]],hh[[k]]};
Do[{mu[[k-1,j]],mu[[k,j]]}={mu[[k,j]],mu[[k-1,j]]};,{j,k-2}];
q=mu[[k,k-1]];
bbb=bb[[k]]+q^2 bb[[k-1]];
If[bbb==0,
   {bb[[k]],bb[[k-1]]}={bb[[k-1]],bb[[k]]};
   Do[{mu[[i,k-1]],mu[[i,k]]}={mu[[i,k]],mu[[i,k-1]]};,{i,k+1,n}];
  ,
   If[bb[[k]]==0,
      bb[[k-1]]=bbb;
      mu[[k,k-1]]=1/q;
      Do[mu[[i,k-1]]=mu[[i,k-1]]/q;,{i,k+1,n}];
      ,
      t=bb[[k-1]]/bbb;(* this t should be nonzero *)
      mu[[k,k-1]]=q t;
      bb[[k]]=bb[[k]]t;
      bb[[k-1]]=bbb;
      Do[
	t=mu[[i,k]];
	mu[[i,k]]=mu[[i,k-1]]-q t;
	mu[[i,k-1]]=t+mu[[k,k-1]]mu[[i,k]];
	,{i,k+1,n}
	];
     ];
  ];
)

doratmethod[testLLLcondition]:=
(
reduceratmethod[k,k-1];
If[bb[[k]]<(3/4-mu[[k,k-1]]^2)bb[[k-1]],
   swapratmethod[k];
   k=Max[2,k-1]; (* The  greedy way, referred to below, would be 
                    If[bb[[k]]!=0||bb[[k-1]]==0,k=Max[2,k-1]];
                  *)
   next=testLLLcondition;
   ,
   For[l=k-2,l>0,l--,reduceratmethod[k,l]];
   k++;
   next=finishedq;
   ];
)

doratmethod[finishedq]:=
(
If[k<=n,
  next=testLLLcondition
  ,
  finished=True;
  ]
)

initmubb:=
(
bb=Range[n];
mu=Table[0,{i,n},{j,i-1}];
bst=b;
Do[
  Do[
     mu[[i,j]]=If[bb[[j]]==0,0,(b[[i]].bst[[j]])/bb[[j]]];
     bst[[i]]=bst[[i]]-mu[[i,j]] bst[[j]];
     ,{j,i-1}
    ];
  bb[[i]]=bst[[i]].bst[[i]];
  If[bb[[i]]<0,Message[GramLatticeReduce::neg,gram];Abort[]];
  ,{i,n}
  ];
bst=.;
next=testLLLcondition;
)

empty:=Sequence[];

extendedlatticereduce[m_,Method->RationalArithmetic]:=CheckAbort[
(
b=m;
(* begin check argument *)
error:=(Message[ExtendedLatticeReduce::argrat,b];Abort[]);
If[MatrixQ[b],
   Map[(If[ #[[0]]===Integer || #[[0]]===Rational,Null,error;])&,b,{2}]
  ,
  error;
  ];
(* end check argument *)
n=Length[b];
zero=0 b[[1]];
k=2;
hh=IdentityMatrix[n];
If[b[[1]]=={},n=0]; (* silly exceptional case: m={{},{}} passes the test. *)
If[n>1,
  For[finished=False;initmubb;,!finished,Null,doratmethod[next]] (* main loop *)
  ];
answer={If[#===zero,empty,#]&/@b,hh}; (* delete zero vectors from b *)
(* "answer" is to be passed on to ExtendedLatticeReduce *)
Clear[q,mu,b,hh,bb,bbb,t,k,l,n,zero];
),$Failed                                                       ];

gramlatticereduce[m_,Method->RationalArithmetic]:=CheckAbort[
(
gram=m;
(* begin check argument *)
error:=(Message[GramLatticeReduce::argrat,gram];Abort[]);
If[MatrixQ[gram],
   Map[(If[ #[[0]]===Integer || #[[0]]===Rational,Null,error;])&,gram,{2}]
  ,
  error;
  ];
If[Length[gram]!=Length[gram[[1]]],error];
If[Length[Union@@(gram-Transpose[gram])]>1,error];
(* end check argument *)
n=Length[gram];
k=2;
b=hh=IdentityMatrix[n];
graminitmubb=OwnValues[initmubb]/. v_ . w_ ->v.gram.w;(* just to show off *)
If[n>1,
  (* main loop *)
  For[finished=False;graminitmubb[[1,2]];,!finished,Null,doratmethod[next]];
  ,
  If[gram[[1,1]]<0,Message[GramLatticeReduce::neg,gram];Abort[]];
  bb=gram[[1]];
  ];
answer={n-Length[Position[bb,0]],hh};
(* to be passed on to GramLatticeReduce *)
Clear[q,mu,b,hh,bb,bbb,t,k,l,n,gram,zero,graminitmubb];
),$Failed                                                   ];

latticereducedq[m_]:=CheckAbort[
(
b=m;
(* begin check argument *)
error:=(Message[LatticeReducedQ::argrat,b];Abort[]);
If[MatrixQ[b],
   Map[(If[ #[[0]]===Integer || #[[0]]===Rational,Null,error;])&,b,{2}]
  ,
  error;
  ];
(* end check argument *)
zero=0 b[[1]];
b=If[#===zero,empty,#]&/@b; (* delete zero vectors *)
n=Length[b];
bb=Range[n];
mu=Table[0,{i,n},{j,i-1}];
bst=b;
(* Gram-Schmidt *)
Do[
  Do[
     mu[[i,j]]=If[bb[[j]]==0,0,(b[[i]].bst[[j]])/bb[[j]]];
     bst[[i]]=bst[[i]]-mu[[i,j]] bst[[j]];
     ,{j,i-1}
    ];
  bb[[i]]=bst[[i]].bst[[i]];
  ,{i,n}
  ];
answer=(Max[mu]<=1/2)&&(-1/2<=Min[mu])&&
(And@@(Table[bb[[k]]>=(3/4-mu[[k,k-1]]^2)bb[[k-1]],{k,2,n}]));
(* to be passed on to LatticeReducedQ *)
Clear[b,mu,bb,bst,zero,n];
),$Failed                      ];

extendedlatticereduce[{},Method->RationalArithmetic]:=(
Message[ExtendedLatticeReduce::empt];answer={{},Null});

latticereducedq[{}]:=(answer=True)

latticereducedq[m___]:=(Message[LatticeReducedQ::nargs];$Failed)

(* In case of failure we wish to return calls unevaluated. That is why     *)
(* we pass the results through the variable "answer", and call a           *)
(* procedure which does the actual work through a side effect of rule      *)
(* checking. Note that one could achieve a similar effect using the        *)
(* Module[ {vars}, body /; condition ] construct.                          *)

privatebody:={answer,Clear[answer]}[[1]]; (* The left hand side becomes   *)
               (* visible through Information[ExtendedLatticeReduce] etc. *)

ExtendedLatticeReduce[x___]:=
privatebody/;extendedlatticereduce[x]=!=$Failed

GramLatticeReduce[x___]:=
privatebody/;gramlatticereduce[x]=!=$Failed

LatticeReducedQ[x___]:=
privatebody/;(latticereducedq[x])=!=$Failed

(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(* Error messages, both for rational and integral algorithm.               *)

ExtendedLatticeReduce::argrat:=
"Argument `1` of ExtendedLatticeReduce should be a matrix with rational
entries."

ExtendedLatticeReduce::argint:=
"Argument `1` of ExtendedLatticeReduce should be a matrix with integral
entries."

ExtendedLatticeReduce::empt:=
"Are you sure you want to reduce an empty generating set? We use Null for
the matrix with no rows or columns."

ExtendedLatticeReduce::met:="`1` is not a valid value for the option
Method.";

ExtendedLatticeReduce::argw:="ExtendedLatticeReduce called with wrong
number of arguments. One argument is expected, not counting options.";

GramLatticeReduce::argrat:=
"Argument `1` of GramLatticeReduce should be a symmetric matrix
with rational entries."

GramLatticeReduce::met:="`1` is not a valid value for the option
Method.";

GramLatticeReduce::argw:="GramLatticeReduce called with wrong number of
arguments. One argument is expected, not counting options.";

GramLatticeReduce::neg:="Gram matrix `1` has negative eigenvalue; it should
be positive semi definite."

LatticeReducedQ::argrat:=
"Argument `1` of LatticeReducedQ should be a matrix with rational entries."

LatticeReducedQ::nargs:="LatticeReducedQ called with wrong
number of arguments; one argument is expected";

RatioLLLCondition::argrat:=
"Value `1` of RatioLLLCondition should be a rational number 
between 1/4 and 1, or a list of such numbers."

(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(*            vvvvvvvvvvvvvvvvvvvvvvvvvvvvvv  cut here, if you wish         *)
(* x----------------------------------------------------------------------x *)
(* Next the version with integer arithmetic, which tries to be faster.      *)
(* However, we also make an effort to keep the transformation matrix        *)
(* reasonable. Basically the latter is achieved by using the LLL algorithm  *)
(* (for the  ordinary inner product) on the rows in the "null space part"   *)
(* of the transformation matrix, or for a component of a later row that     *)
(* lies along the span of the "null space part". Thus it is a side effect   *)
(* of our approach that the "null space part" of the transformation matrix  *)
(* is LLL reduced with respect to the ordinary inner product.               *)
(*                                                                          *)
(* When a dependent vector is encountered, we deviate a little from the     *)
(* Pohst algorithm: this dependence is first removed by a change of         *)
(* generators, based on the extended Euclidean algorithm. Just as for the   *)
(* original LLL algorithm, one may give a good complexity estimate for our  *)
(* version (unless you have chosen ReduceNullSpace->False). Indeed it turns *)
(* out that the sizes of numerators la[[i,j]] and denominators d[[j+1]] or  *)
(* diso[[j+1]] are determined mostly by the size of the gram matrix and the *)
(* sizes of its entries, and less by its rank.                              *)
(*                                                                          *)
(* More specifically, one has after each "trickledown" that                 *)
(* 
   diso[[1+isodim]] d[[1+kmax]]==
   Coefficient[Expand[Det[
                          eps IdentityMatrix[kmax]+
                          (Take[#,kmax]&)/@(Take[gram,kmax]) 
                          (*  = upper left kmax by kmax block of gram *)
                         ]
                     ]
              ,
               eps^isodim
              ]
   Now the coefficients of this polynomial in eps are all non-negative, so  *)
(* each coefficient can be estimated by the value at eps=1, i.e. by         *)
(*    Det[ IdentityMatrix[kmax] + (Take[#,kmax]&)/@(Take[gram,kmax]) ].     *)
(*                                                                          *)
(* The idea underlying the proof of the identity is of course to think of   *)
(*    eps IdentityMatrix[kmax] + (Take[#,kmax]&)/@(Take[gram,kmax])         *)
(* as describing a quadratic form on the span of the first kmax vectors.    *)
(* This leads to a picture that is different from the one we employed in    *)
(* the discussion of the rational algorithm. Work with the perturbed form.  *)
(* It is nondegenerate for positive eps. In fact we think of eps as         *)
(* "positive, but infinitely small."  That is, terms with high powers of    *)
(* eps are negligeable in the presence of terms with lower powers. More     *)
(* formally, work with a ring of real formal power series in eps. We now    *)
(* call a vector isotropic if its length is a multiple of eps. One may do   *)
(* the LLL algorithm with the perturbed form, always computing only to a    *)
(* relevant order in eps. That should lead to something very close to our   *)
(* code. In our code eps is not visible. What is there can be thought of as *)
(* the coefficients of appropriate powers of eps. All this is very much in  *)
(* the spirit of Pohst. We push it a little further and do not throw away   *)
(* the isotropic part. We just do the LLL algorithm on everything,          *)
(* including the isotropic part, with the perturbed form. We did not use    *)
(* this alternative picture when actually writing the code. We just treated *)
(* the standard inner product and the inner product given by the gram       *)
(* matrix separately.                                                       *)
(* x----------------------------------------------------------------------x *)
(*            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^  cut here, if you wish    *)

(* For an explanation of notations, see the debugging code below. *)
(* Compare also the discussion of the rational algorithm.         *)


red[l_]:=(* given k *)
(
t=If[l>isodim,d[[l+1]],diso[[l+1]]];
If[Abs[2la[[k,l]]]>t
  ,
  q=Quotient[2la[[k,l]]+t,2t];
  b[[k]]=b[[k]]-q b[[l]];
  la[[k,l]]=la[[k,l]]-q t;
  Do[la[[k,i]]=la[[k,i]]-q la[[l,i]];,{i,l-1}]
  ]
);

SetAttributes[{swaptail,notLLLcond,coefficient},HoldFirst];
(* d and diso are to be called by name, not by value. *)

swaptail[d_]:=(* given k, q *)
(
bbb=(d[[k-1]]d[[k+1]]+q^2)~Quotient~d[[k]];
Do[
   t=la[[i,k]];
   la[[i,k]]=(d[[k+1]]la[[i,k-1]]-q t)~Quotient~d[[k]];
   la[[i,k-1]]=(bbb t+q la[[i,k]])~Quotient~d[[k+1]]
   ,{i,k+1,kmax}
  ];
d[[k]]=bbb
);

swap:=(* given k *)
(
{b[[k]],b[[k-1]]}={b[[k-1]],b[[k]]};
Do[{la[[k-1,j]],la[[k,j]]}={la[[k,j]],la[[k-1,j]]},{j,k-2}];
q=la[[k,k-1]];
If[k>isodim
   ,
   swaptail[d]
   ,
   swaptail[diso]
   ]
);

notLLLcond[d_]:=(* given k *)
four*d[[k+1]]d[[k-1]]<(three*d[[k]]^2-four*la[[k,k-1]]^2);

LLLbranch[cond_]:=(* given k *)
If[cond
  ,
  swap;
  k=Max[2,k-1]
  ,
  For[l=k-2,l>0,l--,red[l]];
  k++;
  If[k<=n,next=incGS,finished=True]
  ];

do[testLLLcondition]:=
(
red[k-1];
If[k>isodim+1
   ,
   LLLbranch[notLLLcond[d]]
   ,
   LLLbranch[notautomatic&&(k!=isodim+1)&&notLLLcond[diso]]
  ]
);

inner[v_,w_]:=v.gram.w;

coefficient[d_,inner_,k_,j_]:=
(
q=inner[b[[k]],b[[j]]];
Do[q=(d[[i+1]]q-la[[k,i]]la[[j,i]])~Quotient~d[[i]],{i,j-1}];
If[j<k,la[[k,j]]=q,d[[k+1]]=q]
);

do[incGS]:=
(
If[k>kmax
  ,
  kmax=k;
  Do[coefficient[d,inner,k,j]
    ,{j,isodim+1,k} (* la[[k,j]] vanishes for j<=isodim *)
    ];
  If[d[[k+1]]<=0,trickledown,rank++]
  ];
next=testLLLcondition
);

matextendedgcd[x_,y_]:=
({gcd,{px,py}}=ExtendedGCD[x,y];{{-y~Quotient~gcd,x~Quotient~gcd},{px,py}});

isoincGS:=
(
Do[coefficient[diso,Dot,isodim,k],{k,isodim}];
Do[coefficient[diso,Dot,k,isodim],{k,isodim+1,kmax}];
Do[red[l],{k,isodim,kmax},{l,k-1,1,-1}] (* Added in Version 1.10. It  *)
(* simplifies the analysis. Tracing of b suggests it keeps b smaller. *)
);


declutch:=
(
Do[la[[isodim,j]]=0,{j,isodim-1}];
Do[la[[k,isodim]]=0,{k,isodim+1,kmax}]
); 

trickledown:=
(
(*      vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv cut here, if you wish      *)
(*                                                                      *)
(*   We take it in our own hands for a while. The zero entry d[[k+1]]   *)
(* should trickle down until b[[k]] joins the isotropic vectors. Only   *)
(* after that do we update "d" and rescale "la" accordingly. This is    *)
(* possible because our k is descending all this time. And we  pay for  *)
(* it by remembering the ratios by which the volumes drop along the     *)
(* way, whence "rat".                                                   *)
(* In the algorithm of Pohst the k is not "descending all this time".   *)
(* A key difference between our integral algorithm and the algorithm    *)
(* of Pohst is that, until the isotropic part is reached, we only       *)
(* decrease k after a swap if bb[[k]]!=0, i.e. only if the k-th vector  *)
(* is independent (modulo isotropic vectors) of the preceding vectors.  *)
(* That is what makes that our k is descending. One may say that we     *)
(* are greedy and do not to lower k because we are sure to gain a large *)
(* ratio on the spot. If you work this out, you find that this greedy   *)
(* strategy amounts to an extended euclidean algorithm, so that is how  *)
(* we actually implement it below.                                      *)
(* Note that Cohen complains about the extensive updating that one runs *)
(* into when converting the algorithm of Pohst more directly to integer *)
(* arithmetic. We have taken this complaint seriously. Of course our    *)
(* complexity analysis becomes more complex, as we have to assess the   *)
(* damage that trickledown may do to the numerators. This damage is     *)
(* limited, in particular with the default options.                     *)
(*                                                                      *)
(*      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^  cut here, if you wish    *)

If[d[[k+1]]<0,Message[GramLatticeReduce::neg,gramorg];Abort[]];
rat=Hold@@(Table[1,{kmax}]); (* Hold@@ is not essential, cf. comment below *)
Do[
   k--;
   mat=matextendedgcd[d[[k+1]],la[[k+1,k]]];
   rat[[k+1]]=mat[[1,2]];
   {b[[k]],b[[k+1]]}=mat.{b[[k]],b[[k+1]]};
   Do[{la[[k,j]],la[[k+1,j]]}=mat.{la[[k,j]],la[[k+1,j]]},{j,k-1}]
   (*      vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv cut here, if you wish      *)
   (*                                                                     *)
   (*  With the same notations as in the debugging code below,            *)
   (*  we should now estimate the new                                     *)
   (*  Max[Abs[Table[mu[[k,i]],{i,k-1}]]]                                 *)
   (*  divided by the old                                                 *)
   (*  Max[Abs[Table[mu[[k+1,i]],{i,k}]]],                                *)
   (*  as this Max is the quantity which at first sight might explode.    *)
   (*  We claim the growth factor is at most twice rat[[k+1]].            *)
   (*  That is an innocent factor, as the product of all these rat[[k+1]] *)
   (*  (the cumulative rat[[kmax]] below) divides the old d[[kmax]].      *)
   (*  Our estimate of the growth factor uses                             *)
   (*  mat[[1,1]] == - mu[[k+1,k]] mat[[2,1]].                            *)
   (*  The other numerators can be estimated in a similar fashion.        *)
   (*  The length of the b[[i]] can also be estimated (in both inner      *)
   (*  products) and apart from a few "doublings of mu", cf. LLL paper,   *)
   (*  these lengths don't get worse before they get better again.        *)
   (*  For more details see my home page.                                 *)
   (*                                                                     *)
   (*      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^  cut here, if you wish     *)
  ,
  {rank}
  ];
isodim++;(* "null space part" = isotropic part of b is increased *)
rat=Rest[FoldList[Times,1,rat]]; (* cumulative volume ratios *)
(* Only now do we update d and la *)
Do[
   d[[k+1]]=(d[[k]])~Quotient~(rat[[k]]^2);
   q=rat[[k]]rat[[k-1]];
   Do[la[[j,k]]=Quotient[la[[j,k-1]],q],{j,k+1,kmax}];
   ,
   {k,kmax,isodim+1,-1}
  ];
d[[isodim+1]]=1;
If[reducenullspace,isoincGS,declutch];
k=Max[isodim,2]
);

ratiolist[r_Rational]:={r};

ratiolist[{r__Rational}]:={r};

ratiolist[r_]:=(Message[RatioLLLCondition::argrat,r];Abort[]);

setratio:=(
LLLratio=First[LLLratiolist];
If[(1/4<LLLratio<1)
  ,
  three=Numerator[LLLratio];
  four=Denominator[LLLratio];
  LLLratiolist=Rest[LLLratiolist]
  ,
  Message[RatioLLLCondition::argrat,LLLratio];Abort[]
  ];
k=2;
finished=False;
If[(reducerowspace||reducenullspace)&&(k<=n)
  ,
   next=incGS
  ,
   finished=True
  ]
);


init:=(
n=Length[gram];
b=Hold@@(IdentityMatrix[n]);
d=Hold@@(Range[n+1]);
la=Hold@@(Table[0,{i,n},{j,i-1}]);
(* The Hold wrapper is used here because Part does not have Attribute  *)
(* HoldFirst. We want to use Part frequently, without re-evaluating    *)
(* all entries each time. This starts to matter when n gets large.     *)
diso=d;
rank=isodim=0;(* isodim is dimension of isotropic subspace *)
LLLratiolist=ratiolist[LLLratio];
If[!reducerowspace 
  ,(* We compute the null space and then start over with it. *)
   kmax=0;
   Do[do[incGS]; 
      (* with possible side effect on k through trickledown *)
      Do[red[l],{k,k,kmax},{l,k-1,1,-1}]
         (*  To see the need for such reducing, try matrices like    *)
         (*        DiagonalMatrix[Table[(2*n + 1 - i)!,{i, 2*n}]].   *)
         (*        Table[Random[Integer,{-1000,1000}],{2*n},{n}]     *)
         (*  for, say, n=8.                                          *)
      ,
      {k,n}
     ];
   n=isodim;(* temporary meaning, for earlier break in main loop *)
   rank=isodim=0;
   la=Hold@@(Table[0,{i,n},{j,i-1}]);(* We re-initialize, partly *)
  ];
kmax=0;
k=1;
do[incGS];
);

post:=
(
b=List@@b;
If[!reducerowspace
   ,rank=Length[gram]-n;n=Length[gram]
    (* restoring n and rank to their default meanings *)
  ]
)

extendedLLL[m_,integerarithmetic]:=CheckAbort[(
basis=m;
(* begin check argument *)
error:=(Message[ExtendedLatticeReduce::argint,basis];Abort[]);
If[MatrixQ[basis],
   Map[(If[ #[[0]]===Integer,Null,error])&,basis,{2}]
  ,
  error
  ];
(* end check argument *)
gram=basis.Transpose[basis];
For[init;
    ,
    LLLratiolist=!={}
    ,
    ,
    For[setratio,!finished,Null,do[next]]; (* main loop *)
    ];
post;
answer={Take[b.basis,-rank],b}; (* to be passed on to 
                                          ExtendedLatticeReduce *)
Clear[d,diso,la,b,gram,basis,q,bbb,t,k,l,n,
       kmax,gcd,px,py,mat,rat,rank,isodim];
),$Failed                                    ];

gramLLL[m_,integerarithmetic]:=CheckAbort[(
gramorg=gram=m;
(* begin check argument *)
error:=(Message[GramLatticeReduce::argrat,gram];Abort[]);
If[MatrixQ[gram],
   Map[(If[ #[[0]]===Integer || #[[0]]===Rational,Null,error])&,gram,{2}]
  ,
  error
  ];
If[Length[gram]!=Length[gram[[1]]],error];
If[Length[Union@@(gram-Transpose[gram])]>1,error];
(* end check argument *)
If[Count[gram,_Rational,{2}]>0,
   gram=LCM@@(Union@@(Map[Denominator,gram,{2}]))gram
  ];
For[init;
    ,
    LLLratiolist=!={}
    ,
    ,
    For[setratio,!finished,Null,do[next]]; (* main loop *)
    ];
post;
answer={rank,b}; (* to be passed on to GramLatticeReduce *)
Clear[d,diso,la,b,gram,gramorg,q,bbb,t,k,l,n,
       kmax,gcd,px,py,mat,rat,rank,isodim]
),$Failed                                ];

(*                MORE OF THE  INTERFACE                         *)

extendedlatticereduce[m_,rules___Rule]:=
extendedLLL[m,Method/.{rules}/.Options[ExtendedLatticeReduce],rules]
(* We will not bother to catch user of unknown options.      *)
(* Wrong values for ReduceRowSpace or ReduceNullSpace will   *)
(* not be caught either.                                     *)

extendedlatticereduce[m__,{rules___},more___]:=
extendedlatticereduce[m,rules,more]

extendedlatticereduce[arg___]:=(
Message[ExtendedLatticeReduce::argw];$Failed);

extendedLLL[{},___]:=(
Message[ExtendedLatticeReduce::empt];answer={{},Null});

extendedLLL[{{},arg___}?MatrixQ,___]:=
answer={{},IdentityMatrix[Length[{{},arg}]]};

extendedLLL[m_,RationalArithmetic,___]:=
extendedlatticereduce[m,Method->RationalArithmetic]

extendedLLL[m_,IntegerArithmetic,rules___]:=
(reducerowspace=
False=!=(ReduceRowSpace/.{rules}/.Options[ExtendedLatticeReduce]);
reducenullspace=
False=!=(ReduceNullSpace/.{rules}/.Options[ExtendedLatticeReduce]);
LLLratio=(RatioLLLCondition/.{rules}/.Options[ExtendedLatticeReduce]);
notautomatic=True;
extendedLLL[m,integerarithmetic])

extendedLLL[m_,Automatic,rules___]:=
(reducerowspace=
False=!=(ReduceRowSpace/.{rules}/.Options[ExtendedLatticeReduce]);
reducenullspace=True; (* even though we do not reduce the null space *)
LLLratio=(RatioLLLCondition/.{rules}/.Options[ExtendedLatticeReduce]);
notautomatic=False;
extendedLLL[m,integerarithmetic])

extendedLLL[_,val_,___]:=(
Message[ExtendedLatticeReduce::met,val];$Failed);

gramlatticereduce[m_,rules___Rule]:=
gramLLL[m,Method/.{rules}/.Options[GramLatticeReduce],rules]
(* We will not bother to catch user of unknown options.      *)
(* Wrong values for ReduceRowSpace or ReduceNullSpace will   *)
(* not be caught either.                                     *)

gramlatticereduce[m__,{rules___},more___]:=
gramlatticereduce[m,rules,more]

gramlatticereduce[arg___]:=(
Message[GramLatticeReduce::argw];$Failed);

gramLLL[m_,RationalArithmetic,___]:=
gramlatticereduce[m,Method->RationalArithmetic]

gramLLL[m_,IntegerArithmetic,rules___]:=
(reducerowspace=
False=!=(ReduceRowSpace/.{rules}/.Options[GramLatticeReduce]);
reducenullspace=
False=!=(ReduceNullSpace/.{rules}/.Options[GramLatticeReduce]);
LLLratio=(RatioLLLCondition/.{rules}/.Options[GramLatticeReduce]);
notautomatic=True;
gramLLL[m,integerarithmetic])

gramLLL[m_,Automatic,rules___]:=
(reducerowspace=
False=!=(ReduceRowSpace/.{rules}/.Options[GramLatticeReduce]);
reducenullspace=True; (* even though we do not reduce the null space *)
LLLratio=(RatioLLLCondition/.{rules}/.Options[GramLatticeReduce]);
notautomatic=False;
gramLLL[m,integerarithmetic])

gramLLL[_,val_,___]:=(Message[GramLatticeReduce::met,val];$Failed);

(*            vvvvvvvvvvvvvvvvvvvvvvvvvvvvvv  cut here, if you wish        *)
(* x---------------------------------------------------------------------x *)
(* To give some idea of the conventions underlying the integral algorithm, *)
(* we now suggest some debugging code (for the default values only).       *)
(*                                                                         *)
(* Replace the main loop by                                                *)
(*
  For[count=0;init;,!finished,Null,do[next];count++;testall];
*)
(* and use definitions like:                                               *)
(*
testmu[k_,j_]:=If[j>isodim
,
If[mu[[k,j]]!=la[[k,j]]/d[[j+1]],
      Print["error: mu not la/d mu=",mu[[k,j]],
            " la/d=",la[[k,j]]/d[[j+1]],
	    "la=",la[[k,j]]," d=",d[[j+1]]," count=",count
	   ];
      Print[" k=",k," j=",j," kmax=",kmax];finished=True;
  ];
,
If[mu[[k,j]]!=la[[k,j]]/diso[[j+1]],
      Print["error: mu not la/diso mu=",mu[[k,j]],
            " la/diso=",la[[k,j]]/diso[[j+1]],
	    "la=",la[[k,j]]," diso=",diso[[j+1]]," count=",count
	   ];
      Print[" k=",k," j=",j," kmax=",kmax];finished=True;
  ];
		 ];

testbb[k_]:=If[j>isodim
,
If[(bb[[k]]!=d[[k+1]]/d[[k]])&&(bb[[k]]!=0),
   Print["error: bb not d/d bb=",bb[[k]]," d=",
         d[[k+1]]," d=",d[[k]],"
         count=",count
        ];finished=True
  ];
,
If[bb[[k]]!=diso[[k+1]]/diso[[k]],
   Print["error: bb not diso/diso bb=",bb[[k]]," diso=",
         diso[[k+1]]," diso=",diso[[k]],"
         count=",count
        ];finished=True
  ];
	      ];

testall:=
(
Print["testing... ",count," k=",k," kmax=",kmax];
bb=Range[n];
mu=Table[0,{i,n},{j,i-1}];
bst=List@@b;
(* Gram-Schmidt for inner product based on gram matrix *)
Do[
  Do[
     mu[[i,j]]=If[bb[[j]]==0,0,(b[[i]]~inner~bst[[j]])/bb[[j]]];
     If[bb[[j]]!=0,testmu[i,j]];
     bst[[i]]=bst[[i]]-mu[[i,j]] bst[[j]];
     ,{j,i-1}
    ];
  bb[[i]]=bst[[i]]~inner~bst[[i]];
  testbb[i];
  ,{i,kmax}
  ];
If[Count[{la,d},_Rational,{2}]>0,
   Print["error: la or d does not consist of integers"," count=",count];
   finished=true;
  ];
If[Position[bb,0]!=Range[kmax-rank],
   Print["error: isotropic vectors did not trickle down correctly"];
   Print["bb=",bb,"kmax=",kmax,"rank=",rank];
  ];
If[kmax!=(isodim+rank),Print["error: dimensions do not add up.",
" kmax=",kmax," isodim=",isodim," rank=",rank," k=",k]];
bst=List@@b;
(* Gram-Schmidt for ordinary inner product  *)
Do[
  Do[
     mu[[i,j]]=(b[[i]].bst[[j]])/bb[[j]];
     testmu[i,j];
     bst[[i]]=bst[[i]]-mu[[i,j]] bst[[j]];
     ,{j,i-1}
    ];
  bb[[i]]=bst[[i]].bst[[i]];
  testbb[i];
  ,{i,isodim}
  ];
Do[
  Do[
     mu[[i,j]]=(b[[i]].bst[[j]])/bb[[j]];
     testmu[i,j];
     ,{j,isodim}
    ];
  ,{i,isodim+1,kmax}
  ];
)
*)
(* x----------------------------------------------------------------------x *)
(*            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^  cut here, if you wish    *)


End[]

Protect[ExtendedLatticeReduce,GramLatticeReduce,LatticeReducedQ,
	RationalArithmetic,IntegerArithmetic,ReduceRowSpace,ReduceNullSpace
       ];

EndPackage[]

Null



