\\                                                                         \\
\\  Extended Lenstra Lenstra Lovasz algorithm.                             \\
\\                                                                         \\
\\  NAME extendedlll                                                       \\
\\                                                                         \\
\\  DESCRIPTION                                                            \\
\\                                                                         \\
\\  By default extendedlll computes LLL reduced bases of both the lattice  \\
\\  L spanned by the generating set and the lattice R of integer           \\
\\  relations between the generators. It is an extended algorithm in       \\
\\  that it returns the needed transformation as second entry.             \\
\\                                                                         \\
\\  SYNOPSIS                                                               \\
\\                                                                         \\
\\  extendedlll(m[,isgram[,reducecolumnspaceoff[,reducerelationsoff]]])    \\
\\                                                                         \\
\\  Here isgram, reducecolumnspaceoff, reducerelationsoff are 0 or 1.      \\
{\

\\  PROGRAM  \\

extendedlll(X,\
\\                              X is matrix with integer entries           \\

\\         OPTIONS                                                         \\

           isgram,\
\\                              isgram is an option to indicate,           \\
\\                              by means of value 1,                       \\
\\                              that the matrix X should be understood     \\
\\                              as Gram matrix. If isgram is 0 or absent,  \\
\\                              then the columns of the matrix X are       \\
\\                              understood as generators of the lattice L. \\
           reducecolumnspaceoff,\
\\                              option to indicate if nonreduced basis of  \\
\\                              the lattice L is acceptable.               \\
           reducerelationsoff,\
\\                              option to indicate if nonreduced basis of  \\
\\                              the relation lattice R is acceptable.      \\

\\         LOCAL VARIABLES                                                 \\

           B,\
\\            table of denominators                                        \\
           Gram,\
\\               Gram matrix                                               \\
           H,\
\\            transformation matrix                                        \\
           Lam,\
\\              table of numerators                                        \\
           R,
\\           table of volume ratios                                        \\

           A,C,D,Matr11,Matr12,Matr21,Matr22,Q,Q2,T1,T3,Tmp,Tmp2,V1,V3,\
           isodim,k,kk,kmax,l,ll,n,nxt,notfinished,rnk,s,s1,s2)\

\\         LOCAL LOOP VARIABLES ii,j                                       \\

=\
\\  INITIALIZATION  \\

\\ Check matrix argument. \\
if(type(X)!="t_MAT",print("extendedlll: argument is no matrix");1/0,);
Q2=matsize(X);
k=Q2[2];
for(ii=1,Q2[1],\
    for(j=1,k,\
        if(type(X[ii,j])!="t_INT",\
          print("extendedlll: matrix entry is no integer");1/0\
          ,\
          )\
       )\
   );
\\ Get Gram matrix. \\
if(isgram,\
  Gram=X;
  if(Gram==Gram~,\
    ,\
    print("extendedlll: Gram matrix is not symmetric.");1/0\
\\  Actually one may wish to use that the entries below the diagonal do   \\
\\  not affect the computation, and thus forget the symmetry requirement. \\
    );
  ,\
  Gram=(X~)*X\
  );
s1=if((1-reducecolumnspaceoff),1,0);
s2=if((1-reducerelationsoff),1,0);
notfinished=1;
n=matsize(Gram)[1];
H=matid(n);
\\ The array B must be long enough to keep the isotropic part separate. \\
B=vector(n+2,ii,1);
Lam=H;
rnk=0;
isodim=0;
kmax=0;
k=1;
kmax=k;
\\ Gram Schmidt coefficients of first generator. \\
B[k+2]=(H[,k]~)*Gram[,k];
s=sign(B[k+2]);
if(s<=0,\
\\ It is a relation. The rank of the lattice R of relations increases. \\
  isodim=isodim+1;
  if(s,\
    print("extendedlll: Gram matrix is not positive semi-definite.");1/0\
    ,\
    );
  B[isodim+2]=1;
\\ Gram Schmidt coefficients of relation. Here the ordinary inner \\
\\ product is relevant, not the one given by Gram. For the        \\
\\ latter, relations are isotropic. Whence the name isodim.       \\
  B[k+1]=(H[,k]~)*(H[,k]);
  ,\
\\ It is not a relation, and rank of lattice L increases. \\
  rnk=rnk+1\
  );
nxt=0;
k=2;
if((k<=n),nxt=1,notfinished=0);
\\  MAIN LOOP  \\
while(notfinished,\
      if(nxt,\
        if(k>kmax,\
\\        Add generator.     \\
          kmax=k;
\\        Gram Schmidt coefficients of new generator.        \\
          for(j=isodim+1,k,\
              Q=(H[,j]~)*Gram[,k];
              for(ii=isodim+1,j-1,Q=(B[ii+2]*Q-Lam[k,ii]*Lam[j,ii])/B[1+ii]);
              if(j<k,Lam[k,j]=Q,B[k+2]=Q)\
             );
          s=sign(B[k+2]);
          if(s<=0,\
\\          New relation expected.     \\
            if(s,\
              print("extendedlll: Gram matrix has negative eigenvalue.");1/0\
              ,\
              );
            R=vector(kmax,ii,1);
\\          Push generator down into subspaces.          \\
            for(ii=1,rnk,\
                k=k-1;
\\              Extended Euclid.              \\
                A=B[k+2];C=Lam[k+1,k];Matr21=1;D=A;
                if(C,\
                  V1=0;V3=C;
                  while(V3,\
                        Q2=divrem(D,V3);Q=Q2[1];T3=Q2[2];
                        T1=Matr21-Q*V1;Matr21=V1;D=V3;V1=T1;V3=T3\
                       );
                  Matr22=(D-A*Matr21)/C;
                  ,\
                  Matr22=0\
                 );
                Matr11=-C/D;Matr12=A/D;
\\              Matrix which pushes generator one down is            \\
\\              Matr=[Matr11,Matr12;Matr21,Matr22];                  \\
\\              To estimate its 11 entry, note                        \\
\\              Matr11 == - (C/A)*Matr12                             \\

\\              Volume ratio between old and new sublattice.              \\
                R[k+1]=Matr12;
\\              Push one down.              \\
                Tmp=Matr11*H[,k]+Matr12*H[,k+1];
                H[,k+1]=Matr21*H[,k]+Matr22*H[,k+1];H[,k]=Tmp;
\\              Update part of Lam.              \\
                for(j=1,k-1,\
                    Q=Matr11*Lam[k,j]+Matr12*Lam[k+1,j];
                    Lam[k+1,j]=Matr21*Lam[k,j]+Matr22*Lam[k+1,j];Lam[k,j]=Q\
                   )\
               );
\\          Add new relation.          \\
            isodim=isodim+1;
\\          Cumulative volume ratios.          \\
            for(ii=2,kmax,R[ii]=R[ii]*R[ii-1]);
\\          Only now do we update B and Lam.          \\
\\          Thus we avoided a lot of updating.       \\
            forstep(k=kmax,isodim+1,-1,\
                    B[k+2]=(B[1+k])/(R[k]^2);
\\                  These are integers, implying a bound on R[k].          \\
                    Q=R[k]*R[k-1];
                    for(j=k+1,kmax,Lam[j,k]=(Lam[j,k-1])/(Q));
                   );
            B[isodim+2]=1;
\\          Gram Schmidt coefficients of new relation.          \\
            for(k=1,isodim,\
                kk=isodim;ll=k;Q=(H[,kk]~*(H[,ll]));
                for(ii=1,ll-1,Q=(B[ii+1]*Q-Lam[kk,ii]*Lam[ll,ii])/B[ii]);
                if(ll<kk,Lam[kk,ll]=Q,B[kk+1]=Q)\
               );
\\          B[isodim+1] <= (order kmax leading minor of matid(n)+Gram)     \\
            for(k=isodim+1,kmax,\
                kk=k;ll=isodim;Q=(H[,kk]~)*(H[,ll]);
                for(ii=1,ll-1,Q=(B[ii+1]*Q-Lam[kk,ii]*Lam[ll,ii])/B[ii]);
                if(ll<kk,Lam[kk,ll]=Q,B[kk+1]=Q)\
               );
            k=max(isodim,2)\
            ,\
\\          It is not a relation, and rank of lattice L increases.         \\
            rnk=rnk+1\
            )\
          ,\
          );
        nxt=0\
        ,\
\\      Test for swap.      \\
\\      First reduce generator k with generator k-1.      \\
        l=(k-1);A=if(l>isodim,B[l+2],B[l+1]);
        if(abs(2*Lam[k,l])>A,\
          Q=((2*Lam[k,l]+A)\(2*A));
          H[,k]=H[,k]-Q*H[,l];
          Lam[k,l]=Lam[k,l]-Q*A;
          for(ii=1,l-1,Lam[k,ii]=Lam[k,ii]-Q*Lam[l,ii])\
          ,\
          );
\\      Compute if condition for swap is satisfied.      \\
        cond=if(k>isodim+1,\
\\              Test involves no relations.              \\
                s1&&(4*B[k+2]*B[k]<(3*B[1+k]^2-4*Lam[k,k-1]^2))\
                ,\
\\              Test involves relations.              \\
                s2&&((k!=isodim+1)&&\
                4*B[k+1]*B[k-1]<(3*B[k]^2-4*Lam[k,k-1]^2))\
                );
        if(cond,\
\\        Swap.        \\
          Tmp=H[,k-1];H[,k-1]=H[,k];H[,k]=Tmp;
          for(j=1,k-2,Q=Lam[k,j];Lam[k,j]=Lam[k-1,j];Lam[k-1,j]=Q);
          Q=Lam[k,k-1];
          if(k>isodim,\
\\          More Gram Schmidt coefficients to be updated.          \\
            u=(B[k]*B[k+2]+Q^2)/B[1+k];
            for(ii=k+1,kmax,\
                A=Lam[ii,k];
                Lam[ii,k]=(B[k+2]*Lam[ii,k-1]-Q*A)/B[1+k];
                Lam[ii,k-1]=(u*A+Q*Lam[ii,k])/B[k+2]\
               );
            B[1+k]=u\
            ,\
\\          Gram Schmidt coefficients that involve relations.          \\
            u=(B[k-1]*B[k+1]+Q^2)/B[k];
            for(ii=k+1,kmax,\
                A=Lam[ii,k];
                Lam[ii,k]=(B[k+1]*Lam[ii,k-1]-Q*A)/B[k];
                Lam[ii,k-1]=(u*A+Q*Lam[ii,k])/B[k+1]\
               );
            B[k]=u\
            );
          k=max(2,k-1)\
          ,\
\\        No swap required. Reduce generator k with         \\
\\        generators k-2 down to 1.                         \\
          forstep(l=k-2,1,-1,\
                  A=if(l>isodim,B[l+2],B[l+1]);
                  if(abs(2*Lam[k,l])>A,\
                    Q=((2*Lam[k,l]+A)\(2*A));
                    H[,k]=H[,k]-Q*H[,l];
                    Lam[k,l]=Lam[k,l]-Q*A;
                    for(ii=1,l-1,Lam[k,ii]=Lam[k,ii]-Q*Lam[l,ii])\
                    ,\
                    )\
                 );
\\        Go up.        \\
          k=k+1;
          if(k<=n,nxt=1,notfinished=0)\
          )\
        )\
   );
[if(isgram,\
\\  If isgram is 1, then the first entry of the output is the rank of L. \\
    rnk,\
\\  Otherwise the columns of the first entry form a basis of L. \\
\\  This basis is reduced unless reducecolumnspaceoff is 1.     \\
    vecextract(X*H,\
               vector(matsize(X)[1],ii,ii),\
               vector(rnk,ii,n-rnk+ii))\
              )\
,\
\\ The second entry of the output is the required transformation,  \\
\\ as acting on the matrix whose columns span L.                   \\
H
]}