(*: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, 1998                       *)
(*                http://www.math.uu.nl/people/vdkallen/kallen.html *)
(*                bugs to:   vdkallen@math.uu.nl                    *)

(*:Mathematica Version: 3.0 *)
(*:Package Version: 1.111 *)

(*:Warnings:                                                         *)
(* For us a lattice in euclidean n-space need not be of full rank n. *)

(*: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 version is a rather minimal version with few   *)
(* options. It may serve to document the algorithm.                  *)
(*                                                                   *)
(* We add ExtendedLatticeReduce and GramLatticeReduce  to the        *)
(* builtin function LatticeReduce.                                   *)
(*                                                                   *)
(*           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.)            *)

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'.";

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.";

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]:= GramLatticeReduce[generatingset.Transpose[generatingset]]==
	{Length[reducedbasis],transformation}

Out[7]= True

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

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

In[9]:= reducedbasisnullspace.generatingset

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

*)

Unprotect[ExtendedLatticeReduce,GramLatticeReduce];

Options[ExtendedLatticeReduce]={RatioLLLCondition->3/4};

Options[GramLatticeReduce]={RatioLLLCondition->3/4};

Begin["LLLalgorithm`Private`"]

(* The following will be local to this Private block: 
answer           m
b                mat
basis            matextendedgcd
bbb              moreratios
coefficient      n
coefficientiso   
cond             notLLLcond
d                notLLLcondiso
diso             post
                 px
error            py
extendedLLL      q
finished         r
four             rank
gcd              rat
gram             ratiolist
gramLLL          red
gramorg          rules
i                setratio
incGS            swap
init             swaptail
inner            swaptailiso
integerbody      t
isodim           testgramm
                 testLLLcondition
j                testm
k                three
kmax             trickledown
l                v
la               val
LLLbranch        w
LLLratio         x
LLLratiolist     y
*)

(*                     MAIN INTERFACE                               *)

ExtendedLatticeReduce[m_,rules___Rule]:=
(
LLLratio=(RatioLLLCondition/.{rules}/.Options[ExtendedLatticeReduce]);
extendedLLL[m]
)

GramLatticeReduce[m_,rules___Rule]:=
(
LLLratio=(RatioLLLCondition/.{rules}/.Options[GramLatticeReduce]);
gramLLL[m]
)

extendedLLL[m_]:=
CheckAbort[(
           testm[m];
           integerbody;answer={Take[b.basis,-rank],b}
           )
           ,$Failed
          ]

gramLLL[m_]:=
CheckAbort[(
           testgramm[m];
           integerbody;answer={rank,b}
           )
           ,$Failed
          ]

(*                    SUBROUTINES                                  *)

(* ----- Purpose of testm:
Test if main argument of ExtendedLatticeReduce is legal.
Compute gram matrix.
-------  *)

testm[m_]:=
(
basis=m;
(* begin check argument *)
error:=(Message[ExtendedLatticeReduce::argint,basis];Abort[]);
If[MatrixQ[basis],
   Map[(If[ #[[0]]===Integer,,error])&,basis,{2}]
  ,
  error
  ];
(* end check argument *)
gram=basis.Transpose[basis];
If[MatrixQ[gram],,error]
(* To catch "matrices" like {{}}, whose transpose is not a matrix. *)
);

(* ----- Purpose of testgramm:
Test if main argument of GramLatticeReduce is legal.
(But positive semi definiteness will have to wait.)
Remove denominators in gram matrix.
-------  *)

testgramm[m_]:=
(
gramorg=gram=m;
(* begin check argument *)
error:=(Message[GramLatticeReduce::argrat,gram];Abort[]);
If[MatrixQ[gram],
   Map[(If[ #[[0]]===Integer || #[[0]]===Rational,,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
  ];
);

(* ----- Purpose of integerbody:
Body of integer arithmetic algorithm.
-------  *)

integerbody:=
(
For[init;
    ,
    moreratios
    ,
    moreratios=(LLLratiolist=!={})
    ,
    For[setratio,!finished,,testLLLcondition]; (* main loop *)
    ];
post;
)

(* ----- Purpose of init:
Initializing for integer arithmetic.
-------  *)

init:=(
n=Length[gram]; (* Number of vectors *)
b=Hold@@(IdentityMatrix[n]); (* Will be the transformation matrix *)
d=Hold@@(Range[n+1]);  (* Will be a list of denominators *)
la=Hold@@(Table[0,{i,n},{j,i-1}]); (* Will be list of numeratos *)
(* 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; (* Will also be a list of denominators *)
rank=isodim=0;(* isodim is dimension of isotropic subspace *)
LLLratiolist=ratiolist[LLLratio];(* The quality ratios. The default is 3/4. *)
moreratios=True;
kmax=0;(* Number of vectors for which tables are filled. *)
k=1; (* The index of the vector on which attention is focussed *)
incGS;
);

(* ----- Purpose of ratiolist:
Interface for the option RatioLLLCondition.
-------  *)

ratiolist[r_Rational]:={r};

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

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

(* ----- Purpose of setratio:
Process the option RatioLLLCondition and help initializing.
-------  *)

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[(k<=n)
  ,
   incGS
  ,
   finished=True
  ]
);

(* ----- Purpose of post:
Removing the Hold wrapper.
-------  *)

post:=
(
b=List@@b;
)

(* ----- Purpose of incGS:
If b[[k]] is new, extend tables.
If existence of other isotropic vector is detected, go find it and
reorganize, changing "everything".  
-------  *)

incGS:=
(
If[k>kmax
  ,
  kmax=k;
  Do[coefficient[k,j]
    ,{j,isodim+1,k} 
    ];
  Do[coefficientiso[k,j]
    ,{j,isodim} 
    ];
  If[d[[k+1]]<=0,trickledown,rank++]
  ];
);

(* ----- Purpose of inner:
The inner product described by the gram matrix.
-------  *)

inner[v_,w_]:=v.gram.w;(* used only in one place *)

(* ----- Purpose of coefficientiso[k,j]:
Computes la[[k,j]] or diso[[k+1]] assuming tables are up to date for smaller
indices. It assumes j<=isodim, as it works with the euclidean inner product
on the rows of the transformation matrix b.
-------  *)

coefficientiso[k_,j_]:=
(
q=Dot[b[[k]],b[[j]]];
Do[q=(diso[[i+1]]q-la[[k,i]]la[[j,i]])~Quotient~diso[[i]],{i,j-1}];
If[j<k,la[[k,j]]=q,diso[[k+1]]=q]
);

(* ----- Purpose of coefficient[k,j]:
Computes la[[k,j]] or d[[k+1]] assuming tables are up to date for smaller
indices. It assumes j>isodim, as it works with the inner product given
by the gram matrix.
-------  *)

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

(* ----- Purpose of testLLLcondition:
After straightening b[[k]] with respect to b[[k-1]], move to testing
the relevant LLL condition, if any.
-------  *)

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

(* ----- Purpose of LLLbranch[cond]:
Depending on cond, take branch in LLL algorithm. 
-------  *)

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,incGS,finished=True]
  ];


(* ----- Purpose of notLLLcond:
Tell if the LLL condition is violated for the inner product given by the
gram matrix. 
   ----- Purpose of notLLLcondiso:
Tell if the LLL condition is violated for the euclidean inner product
on the rows of the transformation matrix b.
-------  *)

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

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

(* ----- Purpose of red[l]:
Subtract appropriate multiple of b[[l]] from b[[k]] and update table la
accordingly.
-------  *) 

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}]
  ]
);

(* ----- Purpose of swap:
Swap b[[k-1]] with b[[k]] and update tables la, d, diso
accordingly.
-------  *) 

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
   ,
   swaptailiso
   ]
);

(* ----- Purpose of swaptail:
Finish updating tables la and d in a swap.
-------  *)

swaptail:=(* 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
);

(* ----- Purpose of swaptailiso:
Finish updating tables la and diso in a swap.
-------  *)

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

(* ----- Purpose of matextendedgcd[x,y]:
Compute the gcd and the matrix implicit in the euclidean algorithm.
-------  *)

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

(* ----- Purpose of trickledown:
To compute another vector which is isotropic with respect to the inner product 
given by the gram matrix and to put it before the nonisotropic vectors.
To adapt b, isodim, k, la, d, diso to the new situation.
-------  *)

trickledown:=
(
(*                                                                      *)
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}];
   For[l=k-1,l>0,l--,red[l]];
  ,
  {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;
Do[coefficientiso[k,isodim],{k,isodim,kmax}];
Do[red[l],{k,isodim,kmax},{l,k-1,1,-1}];
k=Max[isodim,2]
);


(*                MORE OF THE  INTERFACE                         *)


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

ExtendedLatticeReduce[___]:=(
Message[ExtendedLatticeReduce::argw,val];$Failed);

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

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



(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)
(*                          Error messages.                                *)


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

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

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

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

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

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

(* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *)

End[]

Protect[ExtendedLatticeReduce,GramLatticeReduce,RatioLLLCondition];

EndPackage[]

Null