(*: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, January 95 *) (* bugs to: vdkallen@math.ruu.nl *) (*:Mathematica Version: 2.2 *) (*:Package Version: 1.5 *) (*: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). *) (* 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]] *) (* :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,ReduceRowSpace->False]] may be much larger than those of Last[ExtendedLatticeReduce[matrix]]." (*:Examples: In[1]:= < {{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}; Options[GramLatticeReduce]={Method->IntegerArithmetic, ReduceRowSpace->True,ReduceNullSpace->Automatic}; 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, 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]; 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,,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,,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,,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,,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,,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: *) 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"; (* vvvvvvvvvvvvvvvvvvvvvvvvvvvvvv cut here, if you wish *) (* x----------------------------------------------------------------------x *) (* Next the version with integer arithmetic, which tries to be faster. *) (* It does not follow the other algorithm faithfully: When a dependent *) (* vector is encountered, this dependence is first removed by a change of *) (* generators, based on the extended Euclidean algorithm. Hopefully the *) (* algorithm will not blow up at this stage, as some experiments suggest *) (* it may. We make an effort to keep the transformation matrix reasonable. *) (* Basically this 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. *) (* x----------------------------------------------------------------------x *) (* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ cut here, if you wish *) 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]; 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 *) 4d[[k+1]]d[[k-1]]<(3d[[k]]^2-4la[[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[(k!=isodim+1)&¬LLLcond[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[jkmax , 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}] ); declutch:= ( Do[la[[isodim,j]]=0,{j,isodim-1}]; Do[la[[k,isodim]]=0,{k,isodim+1,kmax}] ); trickledown:=(* We take it in our own hands for a while *) ( If[d[[k+1]]<0,Message[GramLatticeReduce::neg,gramorg];Abort[]]; dorg=d; 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-2,l>0,l--,red[l]]; would such reducing help ? *) , {rank} ]; isodim++;(* "null space part" = isotropic part of b is increased *) d[[isodim+1]]=1; Do[ d[[k+1]]=(d[[k]]dorg[[k]])~Quotient~(dorg[[k-1]]rat[[k]]^2) , {k,isodim+1,kmax} ]; Do[ q=Quotient[dorg[[k]],rat[[k]]]; Do[la[[j,k]]=Quotient[la[[j,k-1]]d[[k+1]],q],{j,k+1,kmax}] , {k,n-1,isodim+1,-1} ]; If[reducenullspace,isoincGS,declutch]; k=Max[isodim,2] ); init:=( finished=False; 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. *) rat=diso=d; rank=isodim=0;(* isodim is dimension of isotropic subspace *) If[!reducerowspace ,(* We compute the null space and then start over with it. *) kmax=0; Block[{isoincGS,declutch} (* to make them do nothing *) ,Do[do[incGS],{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]; k=2; If[(reducerowspace||reducenullspace)&&(k<=n) , next=incGS , finished=True ] ); 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,,error])&,basis,{2}] , error ]; (* end check argument *) gram=basis.Transpose[basis]; For[init,!finished,,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,dorg,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,,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,!finished,,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,dorg,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]); 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]); 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,,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