/* 
     Utrecht, April 1998

     This is program is similar to its GP/PARI counterpart extendedlll.gp,
     but the order of the parameters is different and all five 
     parameters of extendedlll have to be entered.

     The usage of this version of extendedlll is

     extendedlll(m,quality,isgram,reducecolumnspaceoff,reducerelationsoff)

     where
     m is a matrix of integers.
     quality is a rational number, 1/4 < quality < 1.
     (The usual choice is 3/4. If you are in a hurry, try 1/3)
     isgram=1 means that m is a gram matrix.
     isgram=0 means that the columns of m generate the lattice.
     reducecolumnspaceoff=1 means that we do NOT need a reduced lattice basis.
     reducecolumnspaceoff=0 means that we do need a reduced lattice basis.
     reducerelationsoff=1 means that we do NOT care if the first isodim 
                   columns of the transformation matrix are LLL reduced.
                   (These first isodim columns span the kernel of m.)
     reducerelationsoff=0 means that we do care.
     Here isodim = (number of columns of m) - rank (m).
     
     The output of extendedlll is [rank,transformation],
     where rank=rank(m) and the transformation is the one transforming
     the old generators to the new.

     One may work with the standalone program extendedlll-dyn, but then 
     I do not know how to input complicated matrices. 
     One can also install extendedlll inside GP with
     a command like
     install("extendedlll","GGLLL","extendedlll","mypath/libextendedlll.so")
     where mypath tells where libextendedlll.so is to be found.

     Extendedlll is special in that it is designed so that there is a 
     complexity analysis which also estimates the size of the entries of
     the transformation matrix. No similar estimates are known to me for
     the usual implementation.

     Wilberd van der Kallen
     http://www.math.uu.nl/people/vdkallen/kallen.html 

*/

#include "genpari.h"

GEN extendedlll(GEN,GEN,long,long,long);

main()
{
  GEN x,y,z;
  long prec,d,isgram, reducecolumnspaceoff,  reducerelationsoff;
  char s[512];

  init(8000000,2); /* take eight million bytes of memory for the stack */
  prec=3;
  printf("input your matrix in GP format:\n");
  s[0]=0;while(!s[0]) gets(s); x=lisseq(s);
  sor(x,'f',-1,0);
  printf("\ninput quality ratio (a rational number <1 and >1/4 ):\n");
  s[0]=0;while(!s[0]) gets(s); z=lisseq(s);
  sor(z,'f',-1,0);
  printf("\ninput isgram (0 or 1):\n");
  scanf("%d",&d); isgram=d;
  printf("input reducecolumnspaceoff (0 or 1):\n");
  scanf("%d",&d); reducecolumnspaceoff=d;
  printf("input reducerelationsoff (0 or 1):\n");
  scanf("%d",&d); reducerelationsoff=d;
  y=extendedlll(x,z,isgram,reducecolumnspaceoff,reducerelationsoff);
  printf("\n[rank,transformation]=\n");
  sor(y,'f',-1,0);
  if(!isgram) {printf("\nproduct= ");y=gmul(x,(GEN)y[2]);sor(y,'f',-1,0);}
  return 0;
}

/* produit scalaire de x par y sans verification de compatibilite */
static GEN
gscal(GEN x,GEN y)
{
  GEN z,p;
  long i,av,tetpil,lx=lg(x);

  z=gzero; if (lx==1) return z;
  av=avma;
  for (i=1; i<lx; i++)
  {
    p=gmul((GEN)x[i],(GEN)y[i]);
    tetpil=avma; z=gadd(z,p);
  }
  return gerepile(av,tetpil,z);
}





GEN extendedlll(GEN x, GEN z, long isgram, long reducecolumnspaceoff, 
                       long reducerelationsoff)
{
  long av=avma,av1,tetpil,lx=lg(x),tx=typ(x),i,j,k,l,p1,n,lim,dec,q2,\
       s1,s2,isodim,kk,kmax,ll,next,notfinished,rnk,s,cond;
  GEN u,B,lam,q,r,h,p2,p3,p4,p5,gram,            
      a,c,d,Matr11,Matr12,Matr21,Matr22,numz,denz,\
      qq,t1,t3,tmp,tmp2,v1,v3,y, *gptr[5];

  if(tx!=19) err(talker,"extendedlll: argument is no matrix");
  l=lg((GEN)x[1]);
  for(i=1;i<l;i++)
    for(j=1;j<lx;j++)
      {if(typ((GEN)coeff(x,i,j))!=1)
         err(talker,"extendedlll: matrix entry is no integer");
      }
  av=avma;lim=(avma+bot)>>1;
  if(isgram) 
    {gram=x;
    if(lg((GEN)x[1])!=lx)  
      err(talker,"extendedlll: Gram matrix is not square.");
    }
  else
    {gram=(gmul(gtrans(x),x));
    }
  numz=numer(z);
  denz=denom(z);
  if(typ(numz)!=1)
         err(talker,"extendedlll: ratio is not a rational number");
  if(typ(denz)!=1)
         err(talker,"extendedlll: ratio is not a rational number");
  if((cmpii(numz,denz)>=0)) 
         err(talker,"extendedlll: ratio is not <1");
  if((cmpii(denz,shifti(numz,2))>=0)) 
         err(talker,"extendedlll: ratio is not >1/4");
  av1=avma;
  B=cgetg(lx+2,18);
  for(j=1;j<=lx+1;j++) B[j]=un;
  lam=cgetg(lx,19);
  for(j=1;j<lx;j++) lam[j]=lgetg(lx,18);
  for(i=1;i<lx;i++)
    for(j=1;j<lx;j++)
       {coeff(lam,i,j)=zero;
       }
  if(1-reducecolumnspaceoff) {s1=1;} else {s1=0;} 
  if(1-reducerelationsoff) {s2=1;} else {s2=0;}
  notfinished=1;
  n=lx-1;
  h=idmat(n);
  rnk=isodim=0;
  kmax=k=1;
  r=cgetg(lx,18);
  for(i=1;i<lx;i++) r[i]=un;
  B[k+2]=(long)gscal((GEN)h[k],(GEN)gram[k]);
  s=signe((GEN)B[k+2]);
  if(s<=0)
    {isodim=isodim+1;
    if(s)
      err(talker,"extendedlll: Gram matrix has negative eigenvalue.");
    B[isodim+2]=un;
    B[k+1]=(long)gscal((GEN)h[k],(GEN)h[k]);
    }
  else
    {rnk++;
    };
  next=0;
  k=2;
  if(k<=n) {next=1;} else {notfinished=0;}
  while(notfinished)
  {
    if (DEBUGLEVEL>5) { fprintferr("k = %ld",k); flusherr(); };
      if(next) 
        {if(k>kmax) 
           {kmax=k;
           if (DEBUGLEVEL>3) { fprintferr("kmax = %ld\n",k); flusherr(); };
           for(j=isodim+1;j<=k;j++)
              {
              q=gscal((GEN)h[j],(GEN)gram[k]);
              for(i=isodim+1;i<j;i++)
                 {q=divii(subii(mulii(q,(GEN)B[i+2]),
                    mulii(gcoeff(lam,k,i),gcoeff(lam,j,i))),(GEN)B[1+i]);
                 };
              if(j<k) {coeff(lam,k,j)=(long)q;} else {B[k+2]=(long)q;}
              if (low_stack(lim, (av1+bot)>>1))
                {
                 if(DEBUGMEM>1) err(warnmem,"[1]: extendedlll");
                 gptr[0]=&B; gptr[1]=&h; gptr[2]=&r; gptr[3]=&lam;
                 gerepilemany(av1,gptr,4);
                }
              };
           s=signe((GEN)B[k+2]);
           if(s<=0)
             {
             if(s) err(talker,
                      "extendedlll: Gram matrix has negative eigenvalue.");
             for(i=1;i<=kmax;i++) r[i]=un;
             for(i=1;i<=rnk;i++)
                {
                k--;
                a=(GEN)B[k+2];c=gcoeff(lam,k+1,k);Matr21=gun;d=a;
                if(signe(c))
                  {
                  v1=gzero;v3=c;
                  while(signe(v3))
                       {
                        qq=gdiventres(d,v3);
                        q=(GEN)qq[1];t3=(GEN)qq[2];
                        t1=subii(Matr21,mulii(q,v1));
                        Matr21=v1;d=v3;v1=t1;v3=t3;
                       };
                  Matr22=divii(subii(d,mulii(a,Matr21)),c);
                  }
                else
                  {
                  Matr22=gzero;
                  };
                Matr11=gneg(divii(c,d));
                Matr12=divii(a,d);
                r[k+1]=(long)Matr12;
                tmp=gadd(gmul(Matr11,(GEN)h[k]),gmul(Matr12,(GEN)h[k+1]));
                h[k+1]=ladd(gmul(Matr21,(GEN)h[k]),gmul(Matr22,(GEN)h[k+1]));
                h[k]=(long)tmp;
                for(j=1;j<k;j++)
                   {
                    q=addii(mulii(Matr11,gcoeff(lam,k,j)),
                            mulii(Matr12,gcoeff(lam,k+1,j)));
                    coeff(lam,k+1,j)=
                      (long)addii(mulii(Matr21,gcoeff(lam,k,j)),
                                  mulii(Matr22,gcoeff(lam,k+1,j)));
                    coeff(lam,k,j)=(long)q;
                   };
                if (low_stack(lim, (av1+bot)>>1))
                  {
                  if(DEBUGMEM>1) err(warnmem,"[2]: extendedlll");
                  gptr[0]=&B; gptr[1]=&h; gptr[2]=&r; gptr[3]=&lam;
                  gerepilemany(av1,gptr,4);
                  }
	        for(l=(k-1);l>=1;l--)
		     {
                      if (low_stack(lim, (av1+bot)>>1))
                        {
                        if(DEBUGMEM>1) err(warnmem,"[3]: extendedlll");
                        gptr[0]=&B; gptr[1]=&h; gptr[2]=&r; gptr[3]=&lam;
                        gerepilemany(av1,gptr,4);
                         }
                      if(l>isodim) {a=(GEN)B[l+2];} else {a=(GEN)B[l+1];}
         	      if((cmpii(absi(u=shifti(gcoeff(lam,k,l),1)),a)>0))
           	        {q=dvmdii(addii(u,a),shifti(a,1),&y);
	    	         if(signe(y)<0) q=addsi(-1,q);
	                 h[k]=lsub((GEN)h[k],gmul(q,(GEN)h[l]));
	                 coeff(lam,k,l)=lsubii(gcoeff(lam,k,l),mulii(q,a));
	                 for(j=1;j<l;j++)   
                            coeff(lam,k,j)=
                            lsubii(gcoeff(lam,k,j),mulii(q,gcoeff(lam,l,j)));
		        }
                      }
                };
              isodim=isodim+1;
              for(i=isodim+1;i<=kmax;i++) {r[i]=lmul((GEN)r[i],(GEN)r[i-1]);}
              for(k=kmax;k>isodim;k--)
                 {B[k+2]=ldivii((GEN)B[1+k],mulii(u=(GEN)r[k],(GEN)r[k]));
                  q=mulii(u,(GEN)r[k-1]);
                  for(j=k+1;j<=kmax;j++) 
                     {coeff(lam,j,k)=ldivii(gcoeff(lam,j,k-1),q);
                     };
                  if (low_stack(lim, (av1+bot)>>1))
                    {
                    if(DEBUGMEM>1) err(warnmem,"[4]: extendedlll");
                    gptr[0]=&B; gptr[1]=&h; gptr[2]=&r; gptr[3]=&lam;
                    gerepilemany(av1,gptr,4);
                    }
                 };
              B[isodim+2]=un;
              for(k=1;k<=isodim;k++)
                 {kk=isodim;ll=k;
                  q=gscal((GEN)h[kk],(GEN)h[ll]);
                  for(i=1;i<ll;i++)
                     {q=divii(subii(mulii(q,(GEN)B[i+1]),
                          mulii(gcoeff(lam,kk,i),
                                gcoeff(lam,ll,i))),(GEN)B[i]);
                     };
                  if(ll<kk) 
                    {coeff(lam,kk,ll)=(long)q;
                    }
                  else 
                    {B[kk+1]=(long)q;
                    };
                  if (low_stack(lim, (av1+bot)>>1))
                    {
                    if(DEBUGMEM>1) err(warnmem,"[5]: extendedlll");
                    gptr[0]=&B; gptr[1]=&h; gptr[2]=&r; gptr[3]=&lam;
                    gerepilemany(av1,gptr,4);
                    }
                 }
              for(k=isodim+1;k<=kmax;k++)
                 {kk=k;ll=isodim;
                  q=gscal((GEN)h[kk],(GEN)h[ll]);
                  for(i=1;i<ll;i++)
                     {q=divii(subii(mulii(q,(GEN)B[i+1]),
                          mulii(gcoeff(lam,kk,i),
                                gcoeff(lam,ll,i))),(GEN)B[i]);
                     };
                  if(ll<kk) 
                    {coeff(lam,kk,ll)=(long)q;
                    }
                  else 
                    {B[kk+1]=(long)q;
                    };
                 }
	      for(k=isodim;k<=kmax;k++)
	         {
	          for(l=(k-1);l>=1;l--)
		     {
                      if (low_stack(lim, (av1+bot)>>1))
                        {
                        if(DEBUGMEM>1) err(warnmem,"[6]: extendedlll");
                        gptr[0]=&B; gptr[1]=&h; gptr[2]=&r; gptr[3]=&lam;
                        gerepilemany(av1,gptr,4);
                         }
                      if(l>isodim) {a=(GEN)B[l+2];} else {a=(GEN)B[l+1];}
         	      if((cmpii(absi(u=shifti(gcoeff(lam,k,l),1)),a)>0))
           	        {q=dvmdii(addii(u,a),shifti(a,1),&y);
	    	         if(signe(y)<0) q=addsi(-1,q);
	                 h[k]=lsub((GEN)h[k],gmul(q,(GEN)h[l]));
	                 coeff(lam,k,l)=lsubii(gcoeff(lam,k,l),mulii(q,a));
	                 for(i=1;i<l;i++)   
                            coeff(lam,k,i)=
                            lsubii(gcoeff(lam,k,i),mulii(q,gcoeff(lam,l,i)));
		        }
                      }
                  }
              if(isodim>2) {k=isodim;} else {k=2;}
             }
           else
             {rnk++;
             };
           }
        next=0;
        }
      else
        {l=(k-1);
         if(l>isodim) {a=(GEN)B[l+2];} else {a=(GEN)B[l+1];}
         if((cmpii(absi(u=shifti(gcoeff(lam,k,l),1)),a)>0))
           {q=dvmdii(addii(u,a),shifti(a,1),&y);
	    if(signe(y)<0) q=addsi(-1,q);
	    h[k]=lsub((GEN)h[k],gmul(q,(GEN)h[l]));
	    coeff(lam,k,l)=lsubii(gcoeff(lam,k,l),mulii(q,a));
	    for(i=1;i<l;i++)   
               coeff(lam,k,i)=
                 lsubii(gcoeff(lam,k,i),mulii(q,gcoeff(lam,l,i)));
           }
         if (low_stack(lim, (av1+bot)>>1))
           {
            if(DEBUGMEM>1) err(warnmem,"[7]: extendedlll");
            gptr[0]=&B; gptr[1]=&h; gptr[2]=&r; gptr[3]=&lam;
            gerepilemany(av1,gptr,4);
           }
         if(k>isodim+1)
           {p2=(GEN)B[k+2];p5=(GEN)B[k+1];cond=s1&&(
             (cmpii(mulii(numz,mulii(p5,p5)),
                mulii(denz,addii(p3=mulii((GEN)B[k],p2),
                  p4=mulii(q=(GEN)coeff(lam,k,l),
                    gcoeff(lam,k,l)))))>0));
           }
         else
           {p2=(GEN)B[k+1];p5=(GEN)B[k];cond=s2&&((k!=isodim+1)&&
             (cmpii(mulii(numz,mulii(p5,p5)),
                mulii(denz,addii(p3=mulii((GEN)B[l],p2),
                  p4=mulii(q=(GEN)coeff(lam,k,l),
                    gcoeff(lam,k,l)))))>0));
           }
         if(cond)
           {p1=h[l];h[l]=h[k];h[k]=p1;
	    for(j=1;j<=k-2;j++) 
	       {
	        p1=coeff(lam,l,j);coeff(lam,l,j)=coeff(lam,k,j);
	        coeff(lam,k,j)=p1;
	       }
            u=divii(addii(p3,p4),p5);
            for(i=k+1;i<=kmax;i++)
               {   a=gcoeff(lam,i,k);
                    coeff(lam,i,k)=ldivii(
                      subii(mulii(p2,gcoeff(lam,i,l)),mulii(q,a)),
                                         (GEN)p5);
                    coeff(lam,i,l)=ldivii(
                      addii(mulii(q,gcoeff(lam,i,k)),mulii(u,a)),
                                         (GEN)p2);
               }
            if(k>isodim) {B[k+1]=(long)u;} else {B[k]=(long)u;}
            if(k>2) k--;
           }
         else
           {
            for(l=k-2;l>0;l--)
               {
                if (low_stack(lim, (av1+bot)>>1))
                  {
                   if(DEBUGMEM>1) err(warnmem,"[8]: extendedlll");
                   gptr[0]=&B; gptr[1]=&h; gptr[2]=&r; gptr[3]=&lam;
                   gerepilemany(av1,gptr,4);
                  }
                if(l>isodim) {a=(GEN)B[l+2];} else {a=(GEN)B[l+1];}
                if((cmpii(absi(u=shifti(gcoeff(lam,k,l),1)),a)>0))
                  {q=dvmdii(addii(u,a),shifti(a,1),&y);
	           if(signe(y)<0) q=addsi(-1,q);
	           h[k]=lsub((GEN)h[k],gmul(q,(GEN)h[l]));
	           coeff(lam,k,l)=lsubii(gcoeff(lam,k,l),mulii(q,a));
	           for(i=1;i<l;i++)   
                      coeff(lam,k,i)=
                        lsubii(gcoeff(lam,k,i),mulii(q,gcoeff(lam,l,i)));
                  }
               };
             k++;
             if(k<lx) {next=1;} else {notfinished=0;}
           }
        }
    if (low_stack(lim, (av1+bot)>>1))
      {
      if(DEBUGMEM>1) err(warnmem,"[9]: extendedlll");
      gptr[0]=&B; gptr[1]=&h; gptr[2]=&r; gptr[3]=&lam;
      gerepilemany(av1,gptr,4);
      }
  }
  tetpil=avma;
  q=cgetg(lx,19);for(i=1;i<lx;i++) q[i]=lcopy((GEN)h[i]);
  r=cgetg(3,18);r[1]=lcopy(stoi(rnk));r[2]=lcopy(q);
  return gerepile(av,tetpil,r);
}